素数一套:Miller-Rabin 素性检验算法&&Pollard-Rho算法&&线性筛
2021/4/28 22:26:40
本文主要是介绍素数一套:Miller-Rabin 素性检验算法&&Pollard-Rho算法&&线性筛,对大家解决编程问题具有一定的参考价值,需要的程序猿们随着小编来一起学习吧!
部分目录
- Solved 94 / 304 K Gym 100753K Upside down primes
- 高效判断素数
- 快速幂取模
- 继续 Miller-Rabin 素性检验算法
- Unsolved 61 / 182 F Gym 100753F Divisions
- Pollard-Rho算法
- 朴素的算法,试除法
- 生日悖论
- 伪随机数序列
- Floyd判环算法
- 改进log n
- 求最大因数
- 关于template < class T> ,map和vector用法
- AC代码
- 不需要以上算法的方法
- 附:判断立方的方法
- 附:线性筛
Solved 94 / 304 K Gym 100753K Upside down primes
先粘一个我们队的ac代码:
(998ms,我直呼大爷好)
这个判断素数,是线性的,判断一个素数可简单用此,不用线性筛
#include<bits/stdc++.h> using namespace std; int b[20]; bool check(long long a){ while(a){ int tmp=a%10; if(tmp==3||tmp==4||tmp==7) return true; a/=10; } return false; } long long change(long long a){ int cnt=0; long long tmp=1; long long ans=0; while(a){ b[++cnt]=a%10; a/=10; } for(int i=cnt;i>=1;i--){ if(b[i]==6) ans+=9*tmp; else if(b[i]==9) ans+=6*tmp; else ans+=b[i]*tmp; tmp*=10; } return ans; } int main(){ long long a; cin>>a; if(a==1){ cout<<"no"; return 0; } int len=sqrt(a);//由于sqrt(2)=1.414<2,故不会进入循环,不用管 // cout<<"len="<<len<<endl; if(check(a)){ cout<<"no"<<endl; return 0; } for(int i=2;i<=len;i++){ if(a%i==0){ cout<<"no"; return 0; } } a=change(a); // cout<<"a="<<a<<endl; for(int i=2;i<=len;i++){ if(a%i==0){ cout<<"no"; return 0; } } cout<<"yes"<<endl; return 0; }
题意:
判断一个数以及它旋转180度后,是否都是素数
范围:n:1~10e16,还好不大,用最初级的方法没有超时是想不到的。所以说,有希望就要勇于尝试。
题目思路:
字符串读入,
几种情况:
- n=1 ——no
- n中有347之一 ——no
- n不是素数 ——no
- n倒过来不是素数 ——no
- 以上均未no ——yes
接下来的关键是:
介绍高效判断素数的算法
高效判断素数
接下来正式介绍高效判断素数的算法
写一个框架,自顶向下码:
#include<stdio.h> #include<string.h> #include<stdlib.h> #include<math.h> //#include<iostream> //#include<algorithm> #ifdef LOCAL FILE*FP=freopen("text.in","r",stdin); FILE*fp=freopen("text.out","w",stdout); #endif using namespace std; #define mem(a,b) memset(a,b,sizeof(a)) #define ll long long #define _forplus(i,a,b) for(register int i=(a); i<=(b); i++) #define _forsub(i,a,b) for(register int i=(a); i>=(b); i--) char ch[20]={0}; //判断是否素数,注意long long bool isp(ll n){ ...判断素数 } //得正数 inline int check(int a){ if(a==0||a==2||a==5||a==8||a==1)return a; if(a==9)return 6; if(a==6)return 9; } ll change(int tag,char*ch,int len){ ll sum=0; if(tag==1){ _forplus(i,1,len){ sum*=10; sum+=ch[i]-'0'; } }else if(tag==2){ _forsub(i,len,1){ sum*=10; sum+=check(ch[i]-'0'); } } return sum; } //得倒数 int main(){ scanf("%s",ch+1); int len=strlen(ch+1); //判断1 if(len==1&&ch[1]=='1'){ printf("no\n");return 0; } //判断有无347 _forplus(i,1,len){ if(ch[i]=='3'||ch[i]=='4'||ch[i]=='7'){ printf("no\n");return 0; } } //判断正数、倒数,是素数返回1 ,不是返回0 if(!isp(change(1,ch,len))||!isp(change(2,ch,len))){ printf("no\n");return 0; } printf("yes\n"); return 0; }
接下来是判断一个素数如何高效
初版,904ms,和队长想法一样
bool isp(ll n){ int len=(int)sqrt(n); _forplus(i,2,len){ if(n%i==0){ return false; } } return true; }
卡常版,不看2的倍数,速度快一倍,451ms
bool isp(ll n){ if(n==2)return true; if(n%2==0)return false; int len=(int)sqrt(n); for(register int i=3;i<=len;i+=2){ if(n%i==0){ return false; } } return true; }
接下来是大招(甚至帮我判断好了1):
Miller-Rabin 素性检验算法
参考知乎题解
int prime[10]={2, 3, 5, 7, 11, 13, 17, 19, 23, 29}; bool isp(ll n) { if (n <= 1) return false; if (n == 2) return true; int s = 0, t = n - 1; while (!(t % 2)) ++s, t >>= 1; // 求解 n-1=2^s*t for (int i = 0; i < 10 && prime[i] < n; ++i) { int a = prime[i]; int b = 1, m = a, p = t; while (p) { //快速幂,求 b=a^t if (p % 2) b = ((long long) b * m) % n; m = ((long long)m * m) % n; p >>= 1; } if (b == 1) continue; for (int j = 1; j <= s; ++j) { // 进行 s 次二次检验 int k = ((long long)b * b) % n; if(k == 1 && b != n-1) return false; b = k; } if (b != 1) return false; } return true; }
很遗憾的是,也许是我没看懂,也许听说该算法有一定的错误率,wa在了test12…
去看看cf的test12:
Test: #12, time: 0 ms., memory: 0 KB, exit code: 0, checker exit code: 1, verdict: WRONG_ANSWER Input 19 Output no Answer yes Checker Log wrong answer 1st words differ - expected: 'yes', found: 'no'
。。。有点无语,只是被一个小小的19卡到了啊。
看看原理再说
以下来自以上网站,讲的很好:
- 威尔逊定理:对于自然数 p>1,p 是质数当且仅当
- 恶补同余
两个整数a,b,若它们除以整数m所得的余数相等,则称a,b对于模m同余
记作 a ≡ b (mod m) 读作a同余于b模m,或读作a与b关于模m同余。
比如 26 ≡ 14 (mod 12) - 恶补余数
余数不能超过除数,且余数应为正整数
所以-1除以p所得余数为p-1
但是DEVc++上的%结果是-1,可能与简化计算有关。
math.h库里的
double fmod(double x, double y) 返回 x 除以 y 的余数。 fmod(-1,p)
fmod(-1,p)的结果也是-1(注意返回的是double,要么(int)转化)
只对正数有效,判断若为%结果为负,"结果+除数"才是余数
但是这里并不需要算对余数,
只要利用
a%b=c 则相当于a多了个c的零头 (就算计算机算负了,也是少了个一小段) 则(a-c)%b=0 则((p-1)!+1)%p=0 又有步步取余变化,等价于: ((p-1)!%p+1)%p
写出代码:
bool isp(ll n) { ll res=1; for(register ll i=2;i<n;i++){ res*=i; res%=n; } res=(res+1)%n; return res==0; }
虽然效率更慢了(TLE于test25),单独这个并不能减少运算次数。
- 费马小定理:如果 p 是素数,那么有
- 即(a^p-a)%p=0
- 还是((a^(p-1)-1)* a)%p=0
但是不是(a^(p-1)-1)%p=0,如1*2%2=0,1%2=1.
a为任意整数.
虽然费马小定理的逆命题是不成立的,但是不排除它在绝大多数情况下都是成立的。
因此用 p对若干个a符合(a^p-a)%p=0 ,判断p是素数,仍然是不准的。
如 a=2,对于经典的卡迈克数561,它虽然是合数(561=3×11×17),但是会被这个算法判定为质数。
卡迈尔数参考
此时代码
bool 你是质数吗(int n) { if (n <= 1) return false; int t = 1, m = 2, p = n; while(p) { // 快速幂取模 if (p % 2) t = ((long long)t * m) % n; m = (m * m) % n; p >>= 1; } t = (t - 2) % n; return t == 0; }
这里用到了新知识点
快速幂取模
写在:快速幂&&矩阵快速幂&&矩阵乘法&&传矩阵
继续 Miller-Rabin 素性检验算法
二次检验定理:对于质数 p,在0~p-1范围内,满足
的整数只有 1 和 p-1。
简洁证明的博客
- 解释:a|b 数学 a能把b整除
费马二探
如果p是一个质数,而整数a不是p的倍数,则有a^(p-1)≡1(mod p)
- 若a,b,c为任意3个整数,m为正整数,且(m,c)=1,则当a·c≡b·c(mod m)时,有a≡b(mod m)。
证明:a·c≡b·c(mod m)可得ac–bc≡0(mod m)可得(a-b)·c≡0(mod m)。因为(m,c)=1即m,c互质,c可以约去,a– b≡0(mod m)可得a≡b(mod m)。
所以等价的:a ^(p-1)≡1(mod p), a ^(p)≡a(mod p)
参考推导
我不由得想,网上写的东西都没有按人的思维去解释,看了好久,现在在练学习方法。。。也发现,动笔+举例子,提高效率。我觉得现在效率低在一直不了解上,以后多这么学
不是耗时多就好,不是笨,专注更重要
因为我想要去滑冰去,就专心看几分钟,就看完了。。。
专心
其实是这样的
a(p-1)≡1(mod p)——费马定理
条件:整数a不是p的倍数
则,a(p-1)%p应为1
若p为质数,x2≡1(mod p),那么小于p的解只有两个,x=1或x=p−1——二次探测
条件:x2%p为1
则,各x%p应为1或p-1
网上都是代码里,从小到大,一点都不好理解
实际上是从大到小:
测试p,
选定a 为底
以下都要符合,
只要以下有一个不行,就是错:
a(p-1)除以p,余1
p-1化为u*2t(p显然是偶数)
所以上式就是
(a^ (u *2t-1))2余1除以p,
a^ (u *2t-1)应为1或p-1,
若为1或p-1,继续,直到a^u
若为别的,就不是素数
倒过来
x2%p==1&&x!=1&&x!=p-1
则不是素数
若有的
x2%p!=1
则不用看x,进入下一轮,x2当成x看
#include<cstdio> #include<cstring> #include<algorithm> using namespace std; int prime[10]={2,3,5,7,11,13,17,19,23,29}; int Quick_Multiply(int a,int b,int c) //快速积(和快速幂差不多) { long long ans=0,res=a; while(b) { if(b&1) ans=(ans+res)%c; res=(res+res)%c; b>>=1; } return (int)ans; } int Quick_Power(int a,int b,int c) //快速幂,这里就不赘述了 { int ans=1,res=a; while(b) { if(b&1) ans=Quick_Multiply(ans,res,c); res=Quick_Multiply(res,res,c); b>>=1; } return ans; } bool Miller_Rabin(int x) //判断素数 { int i,j,k; int s=0,t=x-1; if(x==2) return true; //2是素数 if(x<2||!(x&1)) return false; //如果x是偶数或者是0,1,那它不是素数 while(!(t&1)) //将x分解成(2^s)*t的样子 { s++; t>>=1; } for(i=0;i<10&&prime[i]<x;++i) //随便选一个素数进行测试 { int a=prime[i]; int b=Quick_Power(a,t,x); //先算出a^t for(j=1;j<=s;++j) //然后进行s次平方 { k=Quick_Multiply(b,b,x); //求b的平方 if(k==1&&b!=1&&b!=x-1) //用二次探测判断 return false; b=k; } if(b!=1) return false; //用费马小定律判断 } return true; //如果进行多次测试都是对的,那么x就很有可能是素数 } int main() { int x; scanf("%d",&x); if(Miller_Rabin(x)) printf("Yes"); else printf("No"); return 0; }
终于看完了
真的挺怀疑自己的,好笨啊
呜呜呜
花了10个小时,搞一个,寒假也是
比常课难多了,为什么这样子。。。
哎,,,
成果:
#include<stdio.h> #include<string.h> #include<stdlib.h> #include<math.h> //#include<iostream> //#include<algorithm> #ifdef LOCAL FILE*FP=freopen("text.in","r",stdin); FILE*fp=freopen("text.out","w",stdout); #endif using namespace std; #define mem(a,b) memset(a,b,sizeof(a)) #define ll long long #define _forplus(i,a,b) for(register int i=(a); i<=(b); i++) #define _forsub(i,a,b) for(register int i=(a); i>=(b); i--) #define N 7 char ch[20]={0}; //判断是否素数,注意long long int prime[N]={2,325,9375,28178,450775,9780504,1795265022}; ll Quick_Multiply(ll a,ll b,ll c) //快速积(和快速幂差不多) { ll ans=0,res=a; while(b) { if(b&1)ans=(ans+res)%c; res=(res+res)%c; b>>=1; } return ans; } ll Quick_Power(ll a,ll b,ll c) //快速幂,这里就不赘述了 { ll ans=1,res=a; while(b) { if(b&1)ans=Quick_Multiply(ans,res,c); res=Quick_Multiply(res,res,c); b>>=1; } return ans; } bool isp(ll x) //Miller_Rabin判断素数 { int i,j; ll k;//注意溢出!(宁愿全ll) ll s=0,t=x-1; if(x==2)return true; if(x<2||!(x&1))return false; while(!(t&1)) { s++; t>>=1; } for(i=0;i<N&&prime[i]<x;++i){ ll a=prime[i]; ll b=Quick_Power(a,t,x); for(j=1;j<=s;++j){ k=Quick_Multiply(b,b,x); if(k==1&&b!=1&&b!=x-1) return false; b=k; } if(b!=1)return false; } return true; } /*bool isp(ll n){ if(n==2)return true; if(n%2==0)return false; int len=(int)sqrt(n); for(register int i=3;i<=len;i+=2){ if(n%i==0){ return false; } } return true; }*/ //得正数 inline int check(int a){ if(a==0||a==2||a==5||a==8||a==1)return a; if(a==9)return 6; if(a==6)return 9; } ll change(int tag,char*ch,int len){ ll sum=0; if(tag==1){ _forplus(i,1,len){ sum*=10; sum+=ch[i]-'0'; } }else if(tag==2){ _forsub(i,len,1){ sum*=10; sum+=check(ch[i]-'0'); } } return sum; } //得倒数 int main(){ scanf("%s",ch+1); int len=strlen(ch+1); //判断1 if(len==1&&ch[1]=='1'){ printf("no\n");return 0; } //判断有无347 _forplus(i,1,len){ if(ch[i]=='3'||ch[i]=='4'||ch[i]=='7'){ printf("no\n");return 0; } } //判断正数、倒数,是素数返回1 ,不是返回0 if(!isp(change(1,ch,len))||!isp(change(2,ch,len))){ printf("no\n");return 0; } printf("yes\n"); return 0; }
Accepted 30 ms
Unsolved 61 / 182 F Gym 100753F Divisions
这是一道分解质因数的题,和上题关系比较密切
给一个数,问它能被多少个数整除?
范围:1~1e18
和上题不太一样,不能快速判断素数了,因为要求有多少个。
那我们有一个定理,只要知道可以分解出的各质数为多少个
如
12=2231
则有
1=2030 3=2031
2=2130 6=2131
4=2230 12=2231
即
n=2q13q2…
答案为
累乘(qi+1) , i from 1 to n.
所以这就是算质因数个数方法
但是打素数表,一个一个判断就超时了
先看看例子:
INPUT | OUTPUT |
---|---|
Sample Input 1 | Sample Output 1 |
12 | 6 |
Sample Input 2 | Sample Output 2 |
999999999999999989 | 2 |
Sample Input 3 | Sample Output 3 |
100000007700000049 | 4 |
学会了快速判断素数的算法后,我似乎有点思路了
有思路就多想一下
对于输入n
如999999999999999989,可以快速判断为素数
则输出2
如100000007700000049,不是素数
则寻找子素数1000000072
先除一个100000007,再看是不是素数?
但是,还是慢了
如果枚举试除的话,到100000007是超时的
翻题解,不出所料,思路是对的
而需要学一个新算法:
Pollard-Rho算法
——John Pollard发明的, 能快速找到大整数的一个非1、非自身的因子的算法。
Pollard-Rho算法
Pollard_rho算法的大致流程是 :
先判断当前数是否是素数(Miller_rabin)了,如果是则直接返回。如果不是素数的话,试图找到当前数的一个因子(可以不是质因子)。然后递归对该因子和约去这个因子的另一个因子进行分解
全B站唯一讲解视频——看博客总是呆懵,看看视频试试是不是更合适
确实理解的多一点
知乎的,他有很多很多算法知识
这个文章关键是要一步一步看,不要跳太快了,一句一句慢慢看,循序渐进
目的:快速找到大整数的一个非1、非自身的因子的算法。
朴素的算法,试除法
方法1:
从2到n-1一个一个试
改进方法:
如果是素数,退出,之后就是2~sqrt(n)
复杂度:O(n)
方法2:
搞笑的方法,随机数看是不是
x = randint(2, n - 1); // 生成2和n-1之间的随机数
不是就再来一次
然而它正是Pollard-Rho算法的基础
O(n)
方法3:
x = randint(2, n - 1); // 生成2和n-1之间的随机数 d = gcd(n, x); // 求gcd
d不为0,返回d
个人改进猜想:
(其实也可以只产生sqrt(n)内的随机数,以后应该也可以只产生sqrt(n)内的随机数而不用GCD)
(后来发现还是要GCD的,GCD不是只把数降到sqrt(n))
生日悖论
如果我们不断在某个范围内生成随机整数,很快便会生成到重复的数,期望大约在根号级别
应用到原题上,即是对于最坏情形 n=p2,如果我们不断在 [1,n-1]间生成随机数,那么期望在生成大约sqrt( p )=n1/4个数后,可以出现两个在模p下相同的数(注意[0,n-1]间的随机数模p大致是[0,p-1]间的随机数)。那么这两个数的差的绝对值 abs,就一定满足abs整除p,则gcd(abs,n)是p的倍数,是因子之一
但这件事意义并没有那么大,因为这个高概率是在两两比较下才成立的。所以我们需要一些技巧。
伪随机数序列
Pollard使用一种特别的伪随机数生成器来生成 [0,n-1]间的伪随机数序列:设序列第一个数为 x , f(x)=(x2+c) mod n,则 x,f(x),f(f(x)),…为一个伪随机数序列。
(也解释了随机数生成原理)
因为每个数都是由前一个数决定的,可以生成的数又是有限的,那么迟早会进入循环。当然,这个循环很可能是混循环
Floyd判环算法
我们在这里使用Floyd判环算法
(跟最短路的Floyd算法是同一个人,但不是同一个算法,这个也叫龟兔赛跑算法)
附:这个UP主的最短路算法学习,这个UP主讲的很好
设置两个变量 t,r,每次判断是否有gcd ( | t - r | , n ) > 1,如果没有,就令 t=f ( t ) , r = f ( f ( r ) ) 。因为 r 跑得更快,如果没有找到答案,最终会与 t 在环上相遇,,t==r,这时退出,换一个 c 重新生成伪随机数。
那么,这有什么好处呢?
注意到若 | t - r | 整除 p , 则 | f ( t ) - f ( r ) | 整除 p ,
由此可得,只要环上距离为 d 的两个数满足条件,那么所有距离为 d 的数都满足条件。在Floyd判环的过程中,每次移动都相当于在检查一个新的距离 d ,这样就不需要进行两两比较了
(这里的 d 是指环上的距离,以相邻为1个单位 )
这个算法的复杂度依赖于这个伪随机数生成器的随机程度,还没有被严格证明。如果它是足够随机的,那么期望复杂度显然是 O ( n1/4log n)。
(我也不知道它是怎么显然的,,,害,可能是:
期望在生成大约sqrt( p )=n1/4个数后,可以出现两个在模p下相同的数)
当前已经比较实用的代码:
ll Pollard_Rho(ll N) { if (N == 4) // 特判4,4找不到解,因为不够随机 return 2; if (is_prime(N)) // 特判质数 return N; while (1) { ll c = randint(1, N - 1); // 生成随机的c auto f = [=](ll x) { return ((lll)x * x + c) % N; }; // lll表示__int128,防溢出 ll t = f(0), r = f(f(0)); while (t != r) { ll d = gcd(abs(t - r), N); if (d > 1) return d; t = f(t), r = f(f(r)); } } }
特判4,4找不到解,因为不够随机
解释:
auto f = [=](ll x) { return ((lll)x * x + c) % N; };
这相当于一个函数,
auto f(ll x){ return ((lll)x * x + c) % N; }
= [=] 我现在还看不懂,,,还没学c++
区别还有最后的 ;就是了。
改进log n
原理:
如果 gcd ( d , n ) > 1 ,
则 gcd ( kd mod n , n ) > 1,
感觉是对的,不知道怎么证明
所以我们可以减少求公因数的次数,即先把一些待选数乘起来,再统一与 n 求公因数
令固定距离C=128
ll Pollard_Rho(ll N) { if (N == 4) return 2; if (is_prime(N)) return N; while (1) { ll c = randint(1, N - 1); auto f = [=](ll x) { return ((lll)x * x + c) % N; }; ll t = 0, r = 0, p = 1, q; do { for (int i = 0; i < 128; ++i) // 令固定距离C=128 { t = f(t), r = f(f(r)); if (t == r || (q = (lll)p * abs(t - r) % N) == 0) // 如果发现环,或者积即将为0,退出 break; p = q; } ll d = gcd(p, N); if (d > 1) return d; } while (t != r); } }
注意,如果积在乘的过程中等于0,那么后面都是0了,可以直接退出,而且这时大概率已经找到了解(除非 | t - r | 整除 n )
(比如,4*3%6==0,则gcd( 4 ,6 )就是解)
求最大因数
给了:
ll max_prime_factor(ll x) { ll fac = Pollard_Rho(x); if (fac == x) return x; else return max(max_prime_factor(fac), max_prime_factor(x / fac)); }
真的这样吗?
比如12
答案应为6
若为
第一层,fac=4
return max(max_prime_factor(4), max_prime_factor(3));
哦,这因子是指
最大素因子
答案为3
关于template < class T> ,map和vector用法
这里
AC代码
特判一下1,注意
#include<stdio.h> #include<string.h> #include<stdlib.h> #include<math.h> #include<time.h> #include<iostream> #include<algorithm> #include<random> #include<map> #ifdef LOCAL FILE*FP=freopen("text.in","r",stdin); FILE*fp=freopen("text.out","w",stdout); #endif using namespace std; #define mem(a,b) memset(a,b,sizeof(a)) #define ll long long #define lll __int128 #define _forplus(i,a,b) for(register int i=(a); i<=(b); i++) #define _forsub(i,a,b) for(register int i=(a); i>=(b); i--) #define INF (1e9+7) #define N 7 map<int,int>p; //判断是否素数,注意long long ll prime[N]={2,325,9375,28178,450775,9780504,1795265022}; ll Quick_Multiply(ll a,ll b,ll c) //快速积(和快速幂差不多) { ll ans=0,res=a; while(b) { if(b&1)ans=(ans+res)%c; res=(res+res)%c; b>>=1; } return ans; } ll Quick_Power(ll a,ll b,ll c) //快速幂,这里就不赘述了 { ll ans=1,res=a; while(b) { if(b&1)ans=Quick_Multiply(ans,res,c); res=Quick_Multiply(res,res,c); b>>=1; } return ans; } bool is_prime(ll x) //Miller_Rabin判断素数 { int i,j; ll k;//注意溢出!(宁愿全ll) ll s=0,t=x-1; if(x==2)return true; if(x<2||!(x&1))return false; while(!(t&1)) { s++; t>>=1; } for(i=0;i<N&&prime[i]<x;++i){ ll a=prime[i]; ll b=Quick_Power(a,t,x); for(j=1;j<=s;++j){ k=Quick_Multiply(b,b,x); if(k==1&&b!=1&&b!=x-1) return false; b=k; } if(b!=1)return false; } return true; } template <class T> T randint(T l, T r = 0) // 生成随机数建议用<random>里的引擎和分布,而不是rand()模数,那样保证是均匀分布 { static mt19937 eng(time(0)); if (l > r) swap(l, r); uniform_int_distribution<T> dis(l, r); return dis(eng); } ll Pollard_Rho(ll n) { if (n == 4) // 特判4,4找不到解,因为不够随机 return 2; if (n == 1)//特判1 ,但是小心1的因子只有1个,其他质数是2个,它不是质数 return 1; if (is_prime(n)) // 特判质数 return n; while (1) { ll c = randint((ll)1, n - 1); // 生成随机的c auto f = [=](ll x) { return ((lll)x * x + c) % n; }; // lll表示__int128,防溢出 ll t = f(0), r = f(f(0)); while (t != r) { ll d = __gcd(abs(t - r), n); if (d > 1) return d; t = f(t), r = f(f(r)); } } } void findp(ll n){ ll k; k=Pollard_Rho(n); if(k==n){ p[k]++; return; } findp(k); findp(n/k); } int main(){ ll n; scanf("%lld",&n); if(n==1){ printf("1\n"); return 0; } findp(n); ll s=1; map<int,int>::iterator i; map<int,int>::iterator end=p.end(); for(i=p.begin();i!=end;i++){ s*=i->second+1; } printf("%d\n",s); return 0; }
不需要以上算法的方法
道高一尺魔高一丈,我找到了不需要以上算法的方法
不得不说,这个方法,用了数论里的“分段”,
乘法是可以分段的,1e18可以化为1e6内分解+剩余几种情况讨论
根本不用直接找因子
并且LL也是1e18,这个方法降到1e6+几个判定,一般是足够的
妙了妙了
可以用于看各素数因子各有多少个,剩余情况Pollard-Rho找也可以
大意是,先把1e6内的素数找到,并分解之,
然后如果剩下的还有,那必由2个以内素数因子组成
因为如果有第三个数,必大于1e6,则整体超过1e18,不对
只要讨论,剩下的是
1,1
2,一个素数
3,因子是2个相同质数
4,因子是2个不同质数
#include <iostream> #include <algorithm> #include <string> #include <cmath> #include <vector> #include <sstream> #include <stack> #include <queue> #include <map> #include <climits> #include <cstdio> #include <set> #include <bitset> using namespace std; #define db(a) (cout << (#a) << " = " << (a) << endl) typedef unsigned long long ll; /** Deterministic Miller-Rabin **/ #define MAX_P (1000000ull) vector<ll> primes; vector<ll> A({2, 3, 5, 7, 11, 13, 17, 19, 23}); // if n < 3,825,123,056,546,413,051, it is enough to test a = 2, 3, 5, 7, 11, 13, 17, 19, and 23. ll N; ll mulmod(ll a, ll b, ll n) { ll erg = 0; ll r = 0; while(b > 0) { // unsigned long long gives enough room for base 10 operations ll x = ((a%n) * (b%10)) % n; for(ll i=0;i<r;i++) x = (x*10)%n; erg = (erg + x) % n; r++; b /= 10; } return erg; } ll fastexp(ll a, ll b, ll n) { if(b == 0) return 1; if(b == 1) return a%n; ll res = 1; while(b > 0) { if(b % 2 == 1) res = mulmod(a, res, n); a = mulmod(a, a, n); b /= 2; } return res; } bool mrtest(ll n) { if(n == 1) return false; ll d = n-1; ll s = 0; while(d % 2 == 0) { s++; d /= 2; } for(ll j=0;j<(ll)A.size();j++) { if(A[j] > n-1) continue; ll ad = fastexp(A[j], d, n); if(ad % n == 1) continue; bool notcomp = false; for(ll r=0;r<=max(0ull,s-1);r++) { ll rr = fastexp(2,r,n); ll ard = fastexp(ad, rr, n); if(ard % n == n-1) {notcomp = true; break;} } if(!notcomp) { return false; } } return true; } bool sieve[MAX_P+1]; void preprimes(ll l) { ll limit = min(MAX_P, l); for(ll i=0;i<=limit;i++) sieve[i] = false; for(ll p=2;p<=limit;p++) { if(sieve[p]) continue; sieve[p] = true; primes.push_back(p); ll i = 2*p; while(i <= limit) { sieve[i] = true; i += p; } } } // how accurate is this? can we generate cases were this fails? bool sqtest(ll n) { ll ns = (ll) sqrt((long double)n); if(ns*ns==n) return true; cout<<ns*ns; ns++; if(ns*ns==n) return true; cout<<ns*ns; return false; } int main() { cin>>N; ll erg = 1; preprimes(N); ll restN = N; for(ll i=0;i<(ll)primes.size();i++) { if(N % primes[i] == 0) { ll tN = N; ll p = 0; while(tN % primes[i] == 0) {tN /= primes[i]; restN /= primes[i]; p++;} // p occurences of primes[i] erg *= p+1; } } //如果剩下的数存在 (必大于MAX_P) if(restN> 1 && sqtest(restN)) erg *= 3; // 如果这个数能整开平方 ,则该素数项指数为2 else if(restN > 1 && !mrtest(restN)) erg *= 4; // 如果这个数是合数 ,则该素数项指数为1,1 //为什么只分解成了2个素数 , //因为如果有第三个数,必大于1e6,则整体超过1e18,不对 else if(restN > 1) erg *= 2; //如果它就是质数,则该项指数1 cout << erg << "\n"; return 0; }
附:判断立方的方法
可以通过i=(ll)pow(x, 1.0/3)求得x的立方根
如同这样如法炮制
pow(i,3)
i++
pow(i,3)
看看其中一个是不是
因为浮点数误差,可能本来是,
比如n,算成n-eps,ll成n-1了
这样子就可以防止之
附:线性筛
for(i=2;i<=sqrt(n);i++){ if(p[i]==0){ for(j=i*i;j<=n;j+=i){ //之所以可以i*i开始,是因为i*(i更小的)即(i更小的)*i,即已经算过了! p[j]=1; } } }
这篇关于素数一套:Miller-Rabin 素性检验算法&&Pollard-Rho算法&&线性筛的文章就介绍到这儿,希望我们推荐的文章对大家有所帮助,也希望大家多多支持为之网!
- 2024-11-23Springboot应用的多环境打包入门
- 2024-11-23Springboot应用的生产发布入门教程
- 2024-11-23Python编程入门指南
- 2024-11-23Java创业入门:从零开始的编程之旅
- 2024-11-23Java创业入门:新手必读的Java编程与创业指南
- 2024-11-23Java对接阿里云智能语音服务入门详解
- 2024-11-23Java对接阿里云智能语音服务入门教程
- 2024-11-23JAVA对接阿里云智能语音服务入门教程
- 2024-11-23Java副业入门:初学者的简单教程
- 2024-11-23JAVA副业入门:初学者的实战指南