数论学习_素数测试
质数(prime number)又称素数,有无限个。质数定义为在大于1的自然数中,除了1和它本身以外不再有其他因数,这样的数称为质数。
对于判定质数,有一个很显然的方法就是判断除了1和他本身之外有没有其他的因数了。
1 bool ok(int N) 2 { 3 if(N==1) return 0; 4 for(int i=2;i<N;++i) 5 if(!N%i) return 0; 6 return 1; 7 }
但是有一个简单粗暴的优化方案是判断到sqrt(N)之内的因子就好了,因为因子是成对出现的,假如a是一个因子,那么N/a也是他的因子,显然当a=sqrt(N)时二者相等,如果再继续找到的因子a那么N/a一定出现在sqrt(N)之前,所以我们只判定到sqrt(N)即可。
bool ok(int N) { if(N==1) return 0; int m=sqrt(n+0.5); for(int i=2;i<=m;++i) if(!N%i) return 0; return 1; }
但是对于大量素数的判定,如果还用这个判定法的话也会很耗时,这是筛法就派上用场了。筛法是一种简单检定素数的算法。据说是古希腊的埃拉托斯特尼(Eratosthenes,约公元前274~194年)发明的,又称埃拉托斯特尼筛法(sieve of Eratosthenes)。
筛法的思想也很简单,枚举所有的素数筛去他们所有的倍数的数(在求解范围之内),显然剩下的数都是素数。但其实,我们不必枚举所有的素数,枚举到sqrt(N)之内的素数即可,这样的正确性在于: 当我们枚举到一个素数x的时候,他的倍数有2,3......x-1,这些倍数如果是质数的话,因为小于x所以之前已经被筛过了,如果是合数根据算术分解定理:算术基本定理可表述为:任何一个大于1的自然数 N,如果N不为质数,那么N可以唯一分解成有限个质数的乘积 N=P1a1P2a2P3a3......Pnan,这里P1<P2<P3......<Pn均为质数,其中指数ai是正整数。所以也能分解出更小的质数,所以在之前也已经筛过。所以直接从x的x倍开始筛,
这样的缺点在于不能讲将所有的素数打表,但是能筛出一个判别数组,判定isp[i]是否为真。
1 bool isp[MAXN]; 2 int prime[MAXN], p = 0; 3 void init(int N) 4 { 5 isp[0] = isp[1] = 1; 6 for (int i = 2;;++i) { 7 if (!isp[i]) { 8 //prime[p++] = i; 9 for (int j = i*i;j <= N;++j) 10 isp[j] = 1; 11 } 12 if (i*i > N) break; 13 } 14 }
接下来说的就是一种线性筛法,欧拉筛法,不仅复杂度更低,还能将判别表和素数表都生成出来。
先上代码:
1 bool isp[MAXN]; 2 int prime[MAXN], p = 0; 3 void init(int N) 4 { 5 isp[0] = isp[1] = 1; 6 for (int i = 2;i <= N;++i) 7 { 8 if (!isp[i]) prime[++p] = i; 9 for (int j = 1;j <= p&&i*prime[j] <= N;++j) { 10 isp[i*prime[j]] = 1; 11 if (i%prime[j] == 0) break; 12 } 13 } 14 }
这个算法保证了用每个数的最小素因子来筛掉他一次且仅一次,所以复杂度是线性的。首先这个方法真实运行起来不一定比埃式要快,因为有取模操作,但是他可以求解一些积性函数问题 ,这个以后再说。之所以说他是线性的是因为对于每个合数只会被他最小的质因子筛去一次。实现这一特性的就在于那句break语句。首先我们知道prime[]数组里的素数是递增的,如果i%prime[j]==0,则我们可以认为i=M*prime[j],后面break掉的就是i*prime[k],他可以表示为M*prime[j]*prime[k],也就是说在后面他会被prime[j]的M*prime[k]>i倍筛掉,所以可以直接break。对于每个筛去的数i*prime[j],prime[j]都是这个数的最小质因子。
讨论了素数的筛法时候,再来说一下素数的计数问题,求解小于n且与n互素的整数个数。给出正整数n的唯一分解式: n=p1a1p2a2p3a3......pkak,求1,2,3......n中与n互素的数的个数,不难想到用容斥定理求解,让n减去与n不互素的数,也就是pi的倍数。先减去不是pi的倍数得数,在加上是pipj倍数得数......这些很好计算就是n/pi,n/pipj。
n-n/p1-n/p2....-n/pk+n/p1p2+n/p1p3......换个思路想就是对于每个因子pi,要么贡献是1要么贡献是-1/pi,得 φ(n)=n*(1-1/p1)*(1-1/p2)......(1-1/pk);这就是著名的欧拉函数。
关于欧拉函数的讨论移步: http://www.cnblogs.com/zzqc/p/7436265.html
对于大素数的判定,尽管对单个素数判定有O(sqrt(N))的算法,但是当N超过一定数量级后也很耗时,此时不妨使用Miller Rabin算法,是一种概率性素数判定方法。
基于费马小定理 a^(p-1)%p=1%p | if(isprime(p)); 算法就是随机出部分的a值来判断qpow(a,p-1,p)==1是否全部成立。实验次数越多结果越准确。
1 #include<iostream> 2 using namespace std; 3 #define LL long long 4 LL qpow(LL a, LL b, LL m) 5 { 6 LL r = 1; for (;b;b >>= 1, a = a*a%m) 7 if (b & 1) r = r*a%m; return r; 8 } 9 bool ok(LL N) 10 { 11 for (int i = 1;i <= 50;++i) 12 { 13 LL x = rand()%(N-2)+2; 14 if (qpow(x%N, N - 1, N) != 1) return 0; 15 }return 1; 16 } 17 int main() 18 { 19 LL N; 20 cin >> N; 21 ok(N) ? puts("Yes") : puts("No"); 22 return 0; 23 }
关于素数的一些应用,
1.分解质因数,给出整数n,找到n所有的质因数。
1 for (int i = 2; i<=sqrt(n); i++) 2 { 3 if (n%i == 0) 4 { 5 fac[p++] = i; 6 while (n%i == 0) n /= i; 7 } 8 } 9 if (n > 1) fac[p++] = n;
试除法,找到一个素因子之后就把这个因子除干净,至于为何i只用循环到sqrt(n),这是因为大于sqrt(n)的质因子至多只有一个,下面证明这句话,反证法假设存在
假设pk是第一个大于sqrt(N)的质因子,pl是大于pk的一个质因子,那么显然pl>pk>sqrt(N),pl*pk>sqrt(N)*sqrt(N)=N,显然质因子的乘积不能大于N,所以我们的结论是正确的。 (ps,想这个结论的证明想了一天,真是好傻啊) 其实还有一个性质是,对于任意的合数P,都至少有一个不大于sqrt(N)的质因子。
2.素数分布,求N以内素数个数
用π(N)表示N以内的素数个数,有一个近似公式 π(N)~N/ln(N) ,称为素数定理。下面是问题描述:
小明是一个聪明的孩子,对数论有着很浓烈的兴趣。他发现求1到正整数10n 之间有多少个素数是一个很难的问题,该问题的难以决定于n 值的大小。现在的问题是,告诉你n的值,让你帮助小明计算小于10n的素数的个数值共有多少位?
只要求出位数,所以不必很精确可以采用近似公式,
ans=log10( N/log(N) ) ,其中N=10n ,利用换底公式化简之后就是: N - log10(N) - log10(log(10))
1 #include<iostream> 2 #include<cmath> 3 using namespace std; 4 int main() 5 { 6 int N; 7 double x; 8 while (cin >> N) { 9 x = N - log10(N) - log10(log(10)); 10 cout << (int)x+ 1 << endl; 11 } 12 return 0; 13 }
3.素因子
有概念: n!的素因子分解中的素数p的幂为: n/p+n/p^2+n/p^3+.........
问 题:从输入中读取一个数n,求出n!中末尾0的个数。
零的个数必然和10有关,而10=2*5,我们如果知道2和5因子的个数便能得出答案,根据上面的性质我们可以直接循环计算得到,但其实5的个数远远低于2的个数,所以直接算出因子5的个数就是答案了。
1 #include<iostream> 2 #include<cmath> 3 using namespace std; 4 int main() 5 { 6 int N, M; 7 cin >> N; 8 while (N--) { 9 cin >> M; 10 int s = 0; 11 while (M) { 12 s += M / 5; 13 M /= 5; 14 } 15 cout << s << endl; 16 } 17 return 0; 18 }