bzoj3667: Rabin-Miller演算法
思路:首先我們說說Miller_Rabin演算法
我們發現了費馬小定理
那它倒過來對不對呢
如果a^(p-1)=1(mod p),那麼p一定是素數嗎?
很不幸,是錯的
雖然出錯概率很低,但是可以被卡
於是我們就給它打補丁
我們又找到了一個二次探測的方法
如果p是質數,那麼x^2=1(mod p)只有兩個解1,p-1 (-1)
那麼它倒過來對不對呢
很不幸,又是錯的
但是兩個錯誤演算法加到一起,出錯概率就很低了
那麼我們先隨機出一些數a[i]
每次拿出一個數a
先用費馬小定理去測試
那麼我們就要算a^(n-1)%n
把n-1拆成2^s*d的形式
這樣我們就可以順便進行二次探測了
先算出a^d次方
然後平方s次不就是a^(n-1)嗎
平方的時候順便檢查一下
最後再用費馬小定理檢測即可
可以證明一次檢測出錯的概率是1/4
那麼很多次後就幾乎不出錯了
然後就是pollard_rho了
設要分解的數是n
如果我們有兩個隨機數x,y
如果gcd(x-y,n)!=1&&gcd(x-y,n)!=n
那麼p=gcd(x-y,n)是n的一個約數
隨機根號n次(1,n)的數,就有很大概率有同樣的數
那麼隨機根號p次,就很有可能有兩個數的差是p的倍數了
這樣我們就會走到一個環上,最後就相遇了、
實現時設計一個隨機函式f(x)
設定k為此次暴力跳的路徑長
每次倍長
x暴力迭代
每次做差求gcd
達到k次後把y賦為x
形象一點就是兩個指標在rho型的東西上走
走到環上相同的點,就可以得到一個p的倍數,p是n的一個因子
然後把這個數和n求gcd,就有可能得到一個約數
先特判n是否為質數
然後因為有可能直接走到n的環,所以如果分解不出n之外的因子那就說明這個隨機函式會使你直接走到n的環上,所以再換一個重試即可
拆出一個因數d後遞迴處理d和n/d即可
還有一點就是快速乘法,這題的模數是longlong的,但是又不想寫高精度
一種處理是把乘法看做多次加法,類似快速冪去做
高階的O(1)做法是:
然後就可以解決這道模板題了
#include<cstdio> #include<cstring> #include<cstdlib> #include<iostream> #include<algorithm> #define abs(a) (a>0?a:-(a)) typedef long long ll; const ll a[]={2,3,5,7,11,13,17,19,23,29}; using namespace std; int cas;ll maxs; void read(ll &x){ char ch; for (ch=getchar();!isdigit(ch);ch=getchar()); for (x=0;isdigit(ch);ch=getchar()) x=x*10+ch-'0'; } ll gcd(ll a,ll b){return !b?a:gcd(b,a%b);} ll mul(ll a,ll b,ll p){ ll d=((long double)a/p*b+1e-8); ll res=a*b-d*p; res=res<0?res+p:res; return res; } ll qpow(ll a,ll b,ll c){ ll res=1; for (;b;b>>=1,a=mul(a,a,c)) if (b&1) res=mul(res,a,c); return res; } bool check(ll a,ll n,ll r,ll s){ ll x=qpow(a,r,n),pre=x; for (int i=1;i<=s;i++){ x=mul(x,x,n); if (x==1&&pre!=1&&pre!=n-1) return 0; pre=x; } if (x!=1) return 0; return 1; } bool MR(ll n){ if (n<=1) return 0; ll r=n-1,s=0; while (!(r&1)) r>>=1,s++; for (int i=0;i<9;i++){ if (a[i]==n) return 1; if (!check(a[i],n,r,s)) return 0; } return 1; } ll pol_rho(ll n,ll c){ //printf("%lld %lld\n",n,c); ll k=2,x=rand()%n,y=x,p=1; for (ll i=1;p==1;i++){ x=(mul(x,x,n)+c)%n; p=y>x?y-x:x-y; p=gcd(n,p); if (i==k) y=x,k+=k; //cout<<" "<<x<<' '<<y<<endl; } return p; } void solve(ll n){ //printf("%lld\n",n); if (n==1) return; if (MR(n)){maxs=max(maxs,n);return;} ll t=n; while (t==n) t=pol_rho(n,rand()%(n-1)); //printf("t=%lld\n",t); solve(t),solve(n/t); } int main(){ srand(1564651598); /*ll a,b,c; scanf("%lld %lld %lld",&a,&b,&c); printf("%lld\n",mul(a,b,c));*/ //for (int i=1;i<=1000;i++) if (MR(i)) printf("%d ",i);puts(""); scanf("%d",&cas); while (cas--){ ll x;maxs=0; read(x),solve(x); if (maxs==x) puts("Prime"); else printf("%lld\n",maxs); } return 0; }