前言
Miller_Rabin是用来快速判断质数,而Pollard_Rho则是用来快速分解质因数
当然两个算法都很玄学,尤其是Pollard_Rho其实是Bogo失散多年的兄弟
Pollard_Rho很不稳定,例如当n=134579128597851时程序可以跑0.1s~1.0s不等
费马小定理
当p为质数且a和p互质时,有
a
p
−
1
≡
1
 
(
m
o
d
 
p
)
a^{p-1}\equiv1\,(mod\,p)
ap−1≡1(modp)一定成立
(当p为合数时可能成立,如果不成立就一定为合数)
现在要证明这个式子
首先由于p是质数,所以有一个优美的性质:
x
y
≡
0
 
(
m
o
d
 
p
)
xy\equiv 0\,(mod\,p)
xy≡0(modp)(x,y>0)成立时,x和y中必须要有至少一个为p的倍数
显然
(
x
 
m
o
d
 
p
)
(
y
 
m
o
d
 
p
)
=
k
p
(x\,mod\,p)(y\,mod\,p)=kp
(xmodp)(ymodp)=kp,则
(
x
 
m
o
d
 
p
)
(
y
 
m
o
d
 
p
)
p
=
k
\frac{(x\,mod\,p)(y\,mod\,p)}{p}=k
p(xmodp)(ymodp)=k
如果x和y mod p都不为0,则找不出能凑成质数p的因子
(当然如果p是合数就不同了,例如8*9 mod 12=0)
考虑a,2a,3a,…,(p-1)a在mod p意义下的值
有两个性质
①这p-1个值两两不同
反证法证明:
假设ma和na(m≠n)在mod p意义下相同,则
m
a
≡
n
a
 
(
m
o
d
 
p
)
ma\equiv na\,(mod\,p)
ma≡na(modp)
则
(
m
−
n
)
a
≡
0
 
(
m
o
d
 
p
)
(m-n)a\equiv 0\,(mod\,p)
(m−n)a≡0(modp)
因为ap互质,所以只可能(m-n)|p
而因为1≤m,n≤p-1且m≠n,所以m-n不可能是p的倍数
②这p-1个值中不会有0
反证法:
假设某个xa mod p=0,则
x
a
≡
0
 
(
m
o
d
 
p
)
xa\equiv 0\,(mod\,p)
xa≡0(modp)
因为ap互质,所以只可能x mod p=0
因为1≤x≤p-1,所以x mod p≠0
因为各不相同且一共只有p-1个数可选,所以a~(p-1)a mod p这p-1个数只可能把1~p-1的数各取一遍
所以
a
∗
2
a
∗
3
a
∗
.
.
.
∗
(
p
−
1
)
a
≡
1
∗
2
∗
3
∗
.
.
.
∗
(
p
−
1
)
 
(
m
o
d
 
p
)
a*2a*3a*...*(p-1)a\equiv 1*2*3*...*(p-1)\,(mod\,p)
a∗2a∗3a∗...∗(p−1)a≡1∗2∗3∗...∗(p−1)(modp)
约掉
(
p
−
1
)
!
(p-1)!
(p−1)!
a
p
−
1
≡
1
 
(
m
o
d
 
p
)
a^{p-1}\equiv 1\,(mod\,p)
ap−1≡1(modp)
得证。
Miller_Rabin
由于费马小定理实在太好♂用了,所以自然会有人想到用其来判断质数
于是xjb找了一堆合数p和一堆底数a来判断后,发现某些合数也可以通过费马小定理的测试
但是人们发现,当用来判断的a越多时,出错的概率就越小
于是有人猜想:如果把所有小于p的a(ap互质)都测试一遍,则判断出来的一定是质数!
然后很快就被打脸了。
事实上,能通过所有与其互质的最小合数很小,仅为561。
(当然如果把所有小于p的a都测一遍且不考虑互质的话绝对不止这个数,而且也不会判错,因为p与所有小于它的正整数都互质,但这样搞还不如√n判断)
所以还是考虑用某个a来判断p
然后有一个著名的二次探测定理,也是Miller_Rabin的精髓所在。
设
x
2
≡
1
 
(
m
o
d
 
p
)
x^2\equiv1\,(mod\,p)
x2≡1(modp),则显然
(
x
−
1
)
(
x
+
1
)
≡
0
 
(
m
o
d
 
p
)
(x-1)(x+1)\equiv0\,(mod\,p)
(x−1)(x+1)≡0(modp)
若x≠1和p-1,则p不为质数
反证法证明:
若x≠1和p-1且p为质数,则(x-1)和(x+1)mod p都不为0
那么就变成了上面的问题,ab mod p=0只有a或b中有一个为p的倍数(mod p=0)
所以(x-1)(x+1) mod p≠0
则x只能为1或p-1,与条件冲突
所以当x≠1和p-1时,p不为质数
注意这并不代表x=1或p-1时p就为质数,因为二次探测定理同样是不确定算法,只能判掉合数
这个定理有什么用呢?
首先根据费马小定理,ap-1 mod p=1
然后对ap-1不断开方,直到有一个结果不为1或者不能继续开方
如果结果为p-1或是当前指数为奇数,则这个数可能是质数
否则一定不是质数
对于每个数x,如果x=1且√x≠1,而且
(√x)2 mod p=1(x),所以
(√x-1)(√x+1) mod p=1
而√x≠1,则不满足二次探测定理,所以p一定不为质数
但如果√x=p-1,则并不能判断x为合数
事实上,还可以继续开方直到某个x为1位置
但显然,如果前面有某个x为1,则x的若干次平方必然为1
而当前的x=p-1,所以前面不可能存在某个x=1,因此继续判断也没有意义了
实现
一个很显然的方法就是把ap-1算出来然后开方
当然这样也很显然meaningless
还有一个很显然的方法就是算出ax(xk=p-1),然后再不断平方直到指数为p-1,把中途的值都记录下来,最后反着判断
这样做没有问题,但是比较麻烦
事实上,可以不用记录所有的值,只需要依次判断当前值就行了
首先如果ax=1,那么显然平方后也是1
如果当前为1,则判断上一个数是否为p-1,不是则为合数
如果直到ap-1都不为1,则a(p-1)/2肯定不为p-1,不必特判
当然如果a是个质数且p=a的话要特判,因为aa-1 mod a=0
Miller_Rabin的正确率很高,如果用2、3、5、7四个数作为a的话,1~25000000000000以内只有3215031751会出错
如果输入的数很大,则可以用慢速乘来搞
code
#include <iostream>
#include <cstdlib>
#include <cstdio>
#define fo(a,b,c) for (a=b; a<=c; a++)
#define fd(a,b,c) for (a=b; a>=c; a--)
using namespace std;
unsigned long long n;
unsigned long long qpower(unsigned long long a,int b)
{
unsigned long long ans=1;
while (b)
{
if (b&1)
ans=ans*a%n;
a=a*a%n;
b>>=1;
}
return ans;
}
bool Miller_Rabin(long long t,int N)
{
if (t==1) return 0;
if (t==2 || t==3 || t==5 || t==7) return 1;
unsigned long long p=t-1;
unsigned long long s,ls=0;
while (!(p&1)) p>>=1;
s=qpower(N,p);
if (s==1) return 1;
while (p<=t-1 && s!=1)
{
ls=s;
s=s*s%t;
p<<=1;
}
return ls==t-1;
}
bool pd(long long n)
{
if ((n/100000==32150) && (n%100000==31751))
return 0;
else
if (Miller_Rabin(n,2) && Miller_Rabin(n,3) && Miller_Rabin(n,5) && Miller_Rabin(n,7))
return 1;
else
return 0;
}
int main()
{
scanf("%lld",&n);
if (pd(n))
printf("Yes\n");
else
printf("No\n");
//in the scope of 1~25000000000000,only 3215031751 is wrong
}
Pollard_Rho
考虑最原始分解质因数的方法:
每次枚举一个≤√n的数,将其全部除去,若最后剩下的数>1的则为其一个质因子
然而这样的复杂度是稳定O(√n)的,考虑更快的方法
先考虑如何求n的一个因子
有一个很显然的方法:
每次随机random一个2~n-1的数,判断是否为n的因子,是则输出,不是则继续
若n=p*q(pq均为质数),则期望时间复杂度为O(n)
而且c++用的是伪随机数,有可能陷入循环无法自拔
生日悖论
有一个问题:
给出n个人,求n个人当中有任意两个人生日相同的期望
证明不会,结论是当n>√365时,相同的概率>50%
可以推广一下,从[1~x]中选n个数,当n>√x时概率>50%
对于生日悖论来说,选择的两个数的关系是a=b
而对于b-a=p这种情况,实际就相当于a加了某个数,结论还是一样
所以可以推广:从[1~x]中选n个数,当n>√x时b-a=p的概率>50%
而n=p*q(最坏情况),则只需要找√n个数然后两两做差,判断差是否为n的因数,这样成功的概率>50%
然而——
这样做的时间复杂度实际上是O(√n2)=O(n),而且正确率堪忧
真·Pollard_Rho
因为找b-a=p的概率实在是太小了,只有
1
n
\frac{1}{n}
n1
但如果我们改一下:
判断gcd(b-a,n)>1呢?(b>a)
显然如果大于1,那么gcd肯定是n的一个因数(最大公因数也是其中一个的因数)
那么可以发现,原来b-a只可能为p或q
而现在b-a可以为p,2p,3p,…(q-1)p,q,2q,3q,…(p-1)q,因为这些数都包含了p或q
所以一共有p+q-2个数,也就是概率为
p
+
q
−
2
n
\frac{p+q-2}{n}
np+q−2
因为n=p*q,所以显然pq越接近它们的和就越小
最坏情况p≈√n,则概率为
√
n
+
√
n
−
2
n
\frac{√n+√n-2}{n}
n√n+√n−2,也就是
√
n
n
\frac{√n}{n}
n√n
一下提升了√n倍,所以显然只需要找 n 4 \sqrt[4]{n} 4n个数(原来是 n \sqrt{n} n个),所以两两匹配的复杂度为 O ( n ∗ l o g   n ) O(\sqrt{n}*log\,n) O(n∗logn),log n太小了所以一般不考虑
具体实现
但是还有一个问题:
如果n特别大,连
n
4
\sqrt[4]{n}
4n都存不下要怎么搞?
显然可以发现,如果把数全部存下来再搞,就相当于每次随机两个数然后做差判断
于是可以每次只存两个数,不断做差判断即可
这样和上面是类似的,最坏时间复杂度为
O
(
n
)
O(\sqrt{n})
O(n)
对于随机数的选择,有一个很好的函数f(x)=(x2+a)%p(a可以为任意数,也可以随机)
显然这样搞实在模p意义下,难免会出现环的情况
于是就有一个神奇的Floyd判环法:
设ab两个初值相同的数,让b按照a的两倍速度来跑,如果在之后的某个时刻a=b则出现了环
然而这样在某些奇♂妙的数(例如4的倍数)下会陷入死循环,所以可以先搞掉质因数2357再说
也可以每次都随机一个a
然而这样搞还是跟Miller_Rabin没有什么关系
假设给出的数是一个大质数,则Pollard_Rho可能会挂掉
所以就需要Miller_Rabin判断当前的数是否是质数
如果要把一个数分解质因数,那么就递归分解到不是质数位置
code
#include <algorithm>
#include <iostream>
#include <cstdlib>
#include <cstdio>
#include <ctime>
#define fo(a,b,c) for (a=b; a<=c; a++)
#define fd(a,b,c) for (a=b; a>=c; a--)
#define abs(x) ((x>0)?(x):-(x))
using namespace std;
long long ans[101];
unsigned long long P,p,A;
int len,i;
unsigned long long chen(unsigned long long a,int b,long long p)
{
unsigned long long ans=0;
while (b)
{
if (b&1)
ans=(ans+a)%p;
a=a*2%p;
b>>=1;
}
return ans;
}
unsigned long long qpower(unsigned long long a,long long b,long long p)
{
unsigned long long ans=1;
while (b)
{
if (b&1)
ans=chen(ans,a,p);
a=chen(a,a,p);
b>>=1;
}
return ans;
}
bool Miller_Rabin(long long t,int N) //N^(t-1)=1
{
if (t==1) return 0;
if (t==2 || t==3 || t==5 || t==7) return 1;
unsigned long long p=t-1;
unsigned long long s,ls=0;
while (!(p&1)) p>>=1;
s=qpower(N,p,t);
if (s==1) return 1;
while (p<=t-1 && s!=1)
{
ls=s;
s=chen(s,s,t);
p<<=1;
}
return ls==t-1;
//in the scope of 1~25000000000000,only 3215031751 is wrong
}
bool pd(long long n)
{
if ((n/100000==32150) && (n%100000==31751))
return 0;
else
if (Miller_Rabin(n,2) && Miller_Rabin(n,3) && Miller_Rabin(n,5) && Miller_Rabin(n,7))
return 1;
else
return 0;
}
long long gcd(long long n,long long m)
{
long long r=n%m;
while (r)
{
n=m;
m=r;
r=n%m;
}
return m;
}
long long f(long long x,long long p) {return (x*x+A)%p;}
long long find(long long a,long long p)
{
long long b=a;
do{
a=f(a,p);
b=f(f(b,p),p);
long long s=gcd(abs(a-b),p);
if (s>1)
return s;
}while (a!=b);
return 0; //failed
}
void Find(long long t)
{
long long s;
while (t>1 && !pd(t))
{
do{
s=find(rand()+2,t);
A=rand()%233333+1;
}while (!s);
Find(s);
t/=s;
}
if (pd(t))
ans[++len]=t;
}
int main()
{
srand(time(NULL)); rand();
A=rand()%233333+1;
scanf("%lld",&p);P=p;
if (p==1)
{
printf("1=1\n");
return 0;
}
if (!p)
{
printf("FaQ\n");
return 0;
}
while (!(p%2)) ans[++len]=2,p/=2;
while (!(p%3)) ans[++len]=3,p/=3;
while (!(p%5)) ans[++len]=5,p/=5;
while (!(p%7)) ans[++len]=7,p/=7;
if (pd(P))
{
printf("%lld is a prime!\n",P);
return 0;
}
Find(p);
sort(ans+1,ans+len+1);
printf("%lld=",P);
printf("%lld",ans[1]);
fo(i,2,len)
printf("*%lld",ans[i]);
printf("\n");
}
后记
Pollard_Rho真好玩
还有一个艾-鲁算法可以判质数但网上好像没有什么资料(http://www.doc88.com/p-885685528836.html)
是什么?♂
参考资料
https://blog.csdn.net/maxichu/article/details/45459533
https://blog.csdn.net/fisher_jiang/article/details/986654
https://www.cnblogs.com/Doggu/p/MillerRabin_PollardRho.html?utm_source=itdadao&utm_medium=referral
https://files-cdn.cnblogs.com/files/Doggu/Pollard-rho算法详解.pdf
http://www.cnblogs.com/Miracevin/p/9697260.html
https://blog.csdn.net/pi9nc/article/details/27209455
http://www.matrix67.com/blog/archives/234#respond
https://blog.csdn.net/rzo_kqp_orz/article/details/68482000