Miller-rabin素性測試與pollard-rho快速質因數分解
~~~~
Miller-Rabin
可以在O(log n)的時間內檢測一個數是否為質數。
理論上是一個概率演算法,但在一定範圍以內的數經過驗證全部可以判出來。
所以OI範圍內可以視作是一個確定演算法。
費馬小定理
若p是質數,且,則有
其逆定理不一定成立,但不符合逆定理的一定不是素數。
根據這個檢測n的素性。
即:嘗試若干個a(a與n不能有公約數(尤拉定理條件),即對於質數n而言a不能是n的倍數,對於合數而言沒關係了)
檢測
然而這還不夠強,會存在判不出的情況
對於Carmichael數,第一個是561,任意與他互質的數判出來都是Yes.
二次探測
顯然地,對於質數n與正整數a<=n
,則有
又因為n是質數,,所以a=1或a=n-1.
與費馬小相同,若不符合此條件則n一定不是質數。
both are ok
兩個加起來,設mr(n,a)表示以a為底數,判定n的返回結果。
將n-1中所有約數2提出來。設去除所有2之後的的數是d.
從d開始二次探測即可。最後考慮一下費馬小。
int mr(ll x,ll a) {
a%=x;
if (a==0) return 1; //注意a不能是x的倍數。
ll d=x-1;
while (!(d&1)) d>>=1;
ll t=ksm(a,d,x);
for (; d<x; d<<=1) {
if (t==1) return 1; //後面所有t都是1,沒必要繼續判了。
if (t!=x-1) {
t=t*t%x;
if (t==1) return 0; //若後一個是1,且當前不是x-1與1,那麼合數。
} else t=t*t%x;
}
return t==1; //要是撐到最後就要看費馬小的結論了。
}
以下資料來自matrix67
1. 如果被測數小於4 759 123 141,那麼只需要測試三個底數2, 7和61就足夠了
2. 如果你每次都用前7個素數(2, 3, 5, 7, 11, 13和17)進行測試,所有不超過341 550 071 728 320的數都是正確的。
3. 如果選用2, 3, 7, 61和24251作為底數,那麼10^16內唯一的強偽素數為46 856 248 255 981。
當然最好需要rand若干個數。
由上可知,檢測一個數n的時間是logn的。
pollar rho
據說時間複雜度是,反正比較玄學吧。
利用生日悖論的結論,作差撞因子會比直接列舉概率大(很多)(超級反直覺)。發現可以利用gcd求約數,用log的時間來做一個炒雞優化。
大概做法
假設我要分解n,那麼
首先每一次用miller rabin判一下n是否質數。
然後定義個隨機函式f(x)=(x^2+c)%c
兩個數a,b,a隨機,b在a的基礎上走一步(取一個f),他們的差是d,將d的絕對值與n求gcd,若不為1或n則找到了一個因子,分下去即可。
若為1或n,那麼要對a,b進行更改,可以用一個隨機函式f(x)=(x^2+c)%mo,c是rand()出來的一個常數。
如果能保證每一對a,b都是有效的,一直做下去很快就能找到因子。(生日悖論)
但會發現可能有迴圈節,這時候無論多少次都是是出不瞭解的。
用一下floyed的判圈思想,b每次走2步,a每次走1步。這樣當a,b重合時(套1~2圈)就不需要繼續判斷了,重新找c(重構隨機函式)。
為什麼a,b初值要間隔一步: 若a,b一開始取相同值,則4是無論如何也判不出來的。 大概和他只有長度為2的圈有關吧。
bzoj4802,注意有longlong相乘。
#include <cstdio>
#include <iostream>
#include <cstdlib>
#include <ctime>
#define abs(a) ((a)>0?(a):-(a))
#include <algorithm>
using namespace std;
typedef long long ll;
typedef long double ld;
ll gcd(ll a,ll b) {return b==0?a:gcd(b,a%b);}
ll _(ll a,ll b,ll c) {
if (a>=c)a%=c; if (b>=c)b%=c;
ll d=(ld)a*b/c;
ll r=a*b-d*c;
return (r%c+c)%c;
}
ll n;
ll ksm(ll a,ll b,ll mo) {
ll ret=1;
for (ll i=1; i<=b; i*=2) {
if (b&i) ret=_(ret,a,mo);
a=_(a,a,mo);
}
return ret;
}
int mr(ll x,ll a) {
a%=x;
if (a==0) return 1;
ll d=x-1;
while (!(d&1)) d>>=1;
ll t=ksm(a,d,x);
for (; d<x; d<<=1) {
if (t==1) return 1;
if (t!=x-1) {
t=_(t,t,x);
if (t==1) return 0;
} else t=_(t,t,x);
}
return t==1;
}
int isprime(ll x) {
return mr(x,2) && mr(x,7) && mr(x,61) && mr(x,rand());
}
ll fo[10000];
int cs=0;
void pr(ll x) {
if (isprime(x)) {
fo[++fo[0]]=x;
return;
}
int flag=fo[0];
while (flag==fo[0]) {
ll c=rand()%x,a=rand()%x;
ll b=(_(a,a,x)+c)%x;
while (a!=b) {
cs++;
ll det=gcd(abs(a-b),x);
if (det>1 && det<x) {
pr(x/det);
pr(det);
return;
}
a=(_(a,a,x)+c)%x;
b=(_(b,b,x)+c)%x;
b=(_(b,b,x)+c)%x;
}
}
}
int main() {
freopen("pol.in","r",stdin);
// freopen("pol.out","w",stdout);
srand(time(0));
rand();
ll nn=0;
cin>>nn;
for (ll n=nn; n<=nn; n++) {
if (n==1) {
cout<<"1"<<endl; return 0;
}
fo[0]=0;
pr(n);
sort(fo+1,fo+1+fo[0]);
ll ans=1;
for (int i=1,cs=0; i<=fo[0]+1; i++) {
if ( i==fo[0]+1 || i!=1 && fo[i-1] != fo[i]) {
ll d = ksm(fo[i-1],cs,n); if (d==0) d=n;
d-=ksm(fo[i-1],cs-1,n);
ans*=d;
cs=1;
} else cs++;
}
cout<<ans<<endl;
}
}