首页 > 编程知识 正文

求0到100的素数的算法,素数怎么算

时间:2023-05-04 14:06:52 阅读:61949 作者:4144

一、引言

平时做主题和运算时,素数的出现次数总是非常频繁。 这里总结了一般的素数判定和计算某个范围内的素数个数的算法。 有些代码来源于匡宾模板。 嗯,毕竟遵循着这个.

二、朴素判断素数算法

判断素数实际上非常简单。 根据定义,要判断某整数n是否为素数,只要在整数区间[2,n-1]内判断是否有某个数m以使n % m==0即可。 代码可以写如下。

intisprime(intn ) {int i; for(I=2; i n; I ) if(n%I==0)返回0; }返回1; }

其实,这个算法是o(n ),感觉很快,但仍然不能满足需求。 因此有O(SQRT(N ) )的算法。 代码可以写如下。

intisprime(longlongn ) {long long i; for(I=2; i * i=n; I ) if(n%I==0)返回0; }返回1; }

原理很巧妙,只是把代码的i n变成了i * i=n,但优化非常高。 你可能注意到了,在前面的代码中,我定义的n是int型,在下面的代码中,我定义了long long。 这是因为在1s内,上一个代码可以判断的位数为1e8,下一个代码可以判断的位数几乎达到1e16。 其理由是,如果a * b=c,如果a是c的约数,则b也是,如果a不是,则b也不是。

这种方法也称为除法。

三、Miller_Rabin素性测试

上面的o(sqrt ) n ) )的算法非常好,但是如果面临更多订单的“数量”,就会变得力不从心。 此时,Miller_Rabin的优势出现了。

Miller_Rabin的理论基础来源于费马小定理。 顺便说一下,费马小定理是数论的四大定理之一。

费马小定理: n是奇数素数,a是任意整数(1 an-1 )时,a^(n-1 ) (1) modn ) ) ) ) ) ) ) ) ) ) ) ) ) )

要测试n是否为素数,首先将n-1分解为(2^s ) d,每次测试开始时,随机选择(1,n-1 )的整数a,然后对于所有r[0,s-1,a ^ d mmid

所以,这个算法只是一个质数测试,他不能完全保证结果是正确的,但当测试次数约为20次时,可以确保正确率几乎达到97%以上。

代码可以写为:

常数int s=20; lmod_mul(lla、LL b、LL n ) {LL res=0; while(b ) {if ) b1 ) RES=) RESa ) % n; a=(a ) % n; b=1; }返回RES; }llmod_exp(lla、LL b、LL n ) {LL res=1; while(b ) {if ) b1 ) RES=mod_mul ) RES,a,n ); a=mod_mul(a,a,n ); b=1; }返回RES; }boolmiller_rabin(lln ) if (n==2||n==3||| n==5||| n==7||| n==11 )返回真; if(n==1|! (n % 2) |! (n % 3) |! (n % 5) |! (n % 7) |! (n % 11 )返回假; LX,pre,u=n - 1,k=0; while (! (u 1) ) k; u=1; }Srand () ll (time ) ) null ); for(intI=0; i S; I ) )//s次测试,s越大,结果越准确x=rand ) (% ) n-2 ) 2; //[2,n]下随机数if(x%n==0) continue; x=mod_exp(x,u,n ); //计算x^u % npre=x; for(intj=0; j k; j ) ) x=mod_mul(x,x,n ); if(x==1pre!=1 pre!=n - 1 )返回假; pre=x; (if ) x!=1)返回假; }返回真; }

四、筛选法

上介

绍的一些素数判断的算法,在某些问题是基本可以适用的了。但是对于另外一类问题却十分尴尬。比如问你整数区间[1, n]内的素数个数是多少。这个时候如果一个个枚举,一个个判断,对于朴素方法来说,最优也是O(nsqrt(n)),即使是Miller_Rabin算法也无法在O(n)的时间内得到结果。于是就有了埃氏筛选法(娇气的摩托筛法)。

对于筛选整数n以内的素数,算法是这么描述的:先把素数2的倍数全部删除,剩下的数第一个为3,再把素数3的倍数全部删除,剩下的第一个数为5,再把素数5的倍数全部删除······直到把n以内最后一个素数的倍数删除完毕,得到n以内的所有素数。代码可以这么写:

const int maxn = 1e7 + 5;int pri[maxn];void getPrime(int n) {for (int i = 0; i <= n; ++i) pri[i] = i;pri[1] = 0;for (int i = 2; i <= n; ++i) {if (!pri[i]) continue;pri[++pri[0]] = i;for (int j = 2; i * j <= n && j < n; ++j) {pri[i * j] = 0;}}}
而事实上,观察可以发现,上面的这种写法有很多次重复计算,这样显然无法做到线性筛选,而另外一种写法却可以得到线性筛选,达到时间复杂度为O(n),代码可以这么写:

const int MAXN = 1e7 + 5;void getPrime(){memset(prime, 0, sizeof(prime));for (int i = 2;i <= MAXN;i++) {if (!prime[i])prime[++prime[0]] = i;for (int j = 1;j <= prime[0] && prime[j] <= MAXN / i;j++) {prime[prime[j] * i] = 1;if (i%prime[j] == 0) break;}}}
来自kuangbin的模板。


五、容斥原理

从上面的代码可以发现,显然这种筛法只能应付达到1e7这种数量级的运算,即使是线性的筛选法,也无法满足,因为在ACM竞赛中,1e8的内存是极有可能获得Memery Limit Exceed的。

于是可以考虑容斥原理。

以AHUOJ 557为例,1e8的情况是筛选法完全无法满足的,但是还是考虑a * b = c的情况,1e8只需要考虑10000以内的素数p[10000],然后每次先减去n / p[i],再加上n / (p[i] * p[j])再减去n / (p[i] * p[j] * p[k])以此类推...于是就可以得到正确结果了。

代码如下:

#include <cmath>#include <cstdio>using namespace std; const int maxn = 10005;int sqrn, n, ans = 0;bool vis[maxn];int pri[1500] = {0};void init(){ vis[1] = true; int k = 0; for(int i = 2; i < maxn; i++){ if(!vis[i]) pri[k++] = i; for(int j = 0; j < k && pri[j] * i < maxn; j++){ vis[pri[j] * i] = true; if(i % pri[j] == 0) break; } }}void dfs(int num, int res, int index){ for(int i = index; pri[i] <= sqrn; i++){ if(1LL * res * pri[i] > n){ return; } dfs(num + 1, res * pri[i], i+1); if(num % 2 == 1){ ans -= n / (res * pri[i]); }else{ ans += n / (res * pri[i]); } if(num == 1) ans++; }}int main(){ init(); while(~scanf("%d",&n) && n){ ans = n; sqrn = sqrt((double)n); dfs(1,1,0); printf("%dn",ans-1); } return 0;}
六、Meissel-Lehmer算法

最后介绍的这个算法可以说是黑科技级别的,它能够在3s内求出1e11之内的素数个数。据说还有在300ms内求出1e11的个数的。可以参考wiki里面原理。然后给出来自Codeforces 665F题目里面的代码。

#define MAXN 100 // pre-calc max n for phi(m, n)#define MAXM 10010 // pre-calc max m for phi(m, n)#define MAXP 40000 // max primes counter#define MAX 400010 // max prime#define setbit(ar, i) (((ar[(i) >> 6]) |= (1 << (((i) >> 1) & 31)))) #define chkbit(ar, i) (((ar[(i) >> 6]) & (1 << (((i) >> 1) & 31))))#define isprime(x) (( (x) && ((x)&1) && (!chkbit(ar, (x)))) || ((x) == 2))namespace pcf {long long dp[MAXN][MAXM];unsigned int ar[(MAX >> 6) + 5] = { 0 };int len = 0, primes[MAXP], counter[MAX];void Sieve() {setbit(ar, 0), setbit(ar, 1);for (int i = 3; (i * i) < MAX; i++, i++) {if (!chkbit(ar, i)) {int k = i << 1;for (int j = (i * i); j < MAX; j += k) setbit(ar, j);}}for (int i = 1; i < MAX; i++) {counter[i] = counter[i - 1];if (isprime(i)) primes[len++] = i, counter[i]++;}}void init() {Sieve();for (int n = 0; n < MAXN; n++) {for (int m = 0; m < MAXM; m++) {if (!n) dp[n][m] = m;else dp[n][m] = dp[n - 1][m] - dp[n - 1][m / primes[n - 1]];}}}long long phi(long long m, int n) {if (n == 0) return m;if (primes[n - 1] >= m) return 1;if (m < MAXM && n < MAXN) return dp[n][m];return phi(m, n - 1) - phi(m / primes[n - 1], n - 1);}long long Lehmer(long long m) {if (m < MAX) return counter[m];long long w, res = 0;int i, a, s, c, x, y;s = sqrt(0.9 + m), y = c = cbrt(0.9 + m);a = counter[y], res = phi(m, a) + a - 1;for (i = a; primes[i] <= s; i++) res = res - Lehmer(m / primes[i]) + Lehmer(primes[i]) - 1;return res;}}
七、后记

写这么多...发现我会的也就这么点了....找了不少的模板...嘛...不过把这些知识梳理一遍,思路也变得清晰许多。以后继续加油~~~

版权声明:该文观点仅代表作者本人。处理文章:请发送邮件至 三1五14八八95#扣扣.com 举报,一经查实,本站将立刻删除。