Miller_Rabin测试
如果需要快速测试一个数是否是素数,有筛法与试除法
此处介绍的是一种基于费马小定理的不确定性算法,当然,这种算法的出错率是极其微小的,尤其当选择的测试数较多时,因此在OI比赛中成为一种实用的算法
好的进入正题。
我们知道费马小定理,若
p
为素数,那么对于a为任意正整数且a不为p的倍数时,有
然而它的逆定理并不成立,即若 ap−1≡1(modp) ,p不一定是一个素数
我们称满足费马小定理的p为伪素数
它很大可能是一个素数
那么就有费马测试,多弄几个a,如果都符合,它就很有可能是素数
但是存在这样一些数,无论怎么取与这个数互质的a,它都满足费马小定理
定义这样的数叫做卡米切尔数
第一个卡米切尔数是561,=3*11*17
这样的数在
1
~
这样的错误率是不能接受的
考虑强化这个算法。
我们知道若 a2≡1(modp),(a,p)=1
那么 (a+1)(a−1)≡0(modp)
由于
(a,p)=1
那么要么
a+1≡0(modp)
,要么
a−1≡0(modp)
那么同余方程的解为
a1=1,a2=p−1
既然
ap−1≡1(modp)
考虑将 p−1 写成 d∗2q 的形式
那么就有了新的一种判断方法,将p-1不断除以2,判断它的结果是否仍是1或p-1
若是p-1则退出,若是1则继续除下去,直到除不了了或者出现了不是这两种结果的
这样的每除一次都需要跑一次快速幂,所以我们考虑优化。
可以一次直接除到底,然后再乘上来,看是否是出现了p-1或全部都是1
a可以取多几个,这样在 1018 内的错误率都是极小的,并且如果a随机,几乎是不可能卡掉的
附上代码
LL ti(LL x,LL y,LL m)
{
if(y==0||x==0) return 0;
LL s=ti(x,y/2,m);
return (y%2)?((s+s)%m+x)%m:(s+s)%m;
}
bool miller_rabin(LL k)
{
int d[5]={2,3,5,7,10007};
fo(i,0,4)
{
if(k==d[i]) return 1;
if(k%d[i]==0) return 0;
LL t,n=k-1;
while((n&1)==0) n>>=1;
t=ksm(d[i],n,k);
while(n<k-1&&t!=1&&t!=k-1) t=ti(t,t,k),n<<=1;//这里有可能两个小于10^18的数相乘会爆掉long long,所以要把乘改成类似快速幂的加法
if(!((n&1)||t==k-1)) return 0;
}
return 1;
}
Pollard_Rho算法
pollard_rho是一种基于miller_rabin测试的质因数分解算法
算法流程是若分解一个数k,先判断k是否是质数(Miller_rabin)
然后尝试找到k的一个非1非k的因子p,递归分解p和k/p
关键在于找这个因子p
首先介绍一个玄学的东西
生日“悖论”
一个班有k个同学,至少有两个同学的生日相同的概率P是多少?
k=1,P=0
k=2,P=1−1∗364365
k=3,P=1−1∗364365∗363365
。
。
。
令人难以置信的是
k=23,P=0.507
当
k>55
时,概率就很接近1了
当
k=100,P=0.999999692751072
这与人们的感觉严重不符,因此称之为“悖论”。
事实上这个悖论可以扩展到一般情况
给出N个数,至少有一对数之差等于a的概率
这个增长同上面一样,也是非常快的。
回到正题
那么我们此时并不需要暴力随机因子,而是可以通过用两个数相减的方式计算因子。
定义一个随机函数 f(x)=(x∗x+c)modk ,c是这一次找因子随机出来的数
每次随机出新的x,与之前的x做差,与k求gcd,看这个数是否是合法的因子,不是则继续随机
这个算法有两点需要注意的地方。
- 之前有很多个x,取哪一个?
- 事实上,对于同一个c,x是会循环的
形如希腊字母 ρ (rho),又是pollard最先发明的,因此命名为pollard_rho算法
那么有两种做法。
一种是取一开始
x=y
,每次随机
x=f(x),y=f(f(y))
,y以两倍速度随机,如果最后x=y了,那么说明出现了环。
还有一种做法,我并不太知道这个做法的本质,但是据说这个做法非常快。
每次随机
x=f(x)
,在每第2的幂次随机时,更新y,即y=x,每次x与y做差,当某个时刻x=y了说明出现了环
这种做法具体可以看代码
LL ti(LL x,LL y,LL m)//乘法用类似快速幂的做法累加,防止爆long long
{
if(y==0||x==0) return 0;
LL s=ti(x,y/2,m);
return (y%2)?((s+s)%m+x)%m:(s+s)%m;
}
LL get(LL k)
{
LL c=rd(1e9)+1,r,x,y,i=1,n=2;
x=y=rd(k-1)+1;
while(1)
{
i++;
x=(ti(x,x,k)+c)%k;//随机下一个x
if(x==y) return k;//出现环
r=gcd(abs(x-y),k);
if(1<r&&r<k)
{
return r;//找到因子
}
if(i==n) y=x,n<<=1;//n是2的幂,如果随机次数到了2的某次方则更新y
}
}
void pollard_rho(LL k)
{
if(k==1) return;
if(miller_rabin(k))//这是判断是否是素数
{
pr[++pr[0]]=k;
return;
}
LL p=k;
while(p==k&&p!=1) p=get(k);
pollard_rho(k/p),pollard_rho(p);
}
Pollard_rho的复杂度我不太会分析,据大佬说平均是
O(N1/4)
的
这在OI竞赛中都是绝对够用的了