Pollard Rho 跟 Miller Rabin 的合併算法
Miller Rabin 算法：用於快速檢查一個數字是不是質數
假設我們現在在檢查p。
若p是質數，根據費馬小定理，對於每個a，a ^ p === a (mod p)。
兩邊同除a (就自行想像一下吧XDD)，可能a ^ (p-1) === 1 (mod p)。
其中，若x ^ y === 1，x合理的解就只有可能是 +1 或 -1。
所以，可以發現，a ^ ( (p-1)/2 ), a ^ ( ( (p-1)/2 )/2 ), .......，只要是1 或 -1，那這個數字可能就是質數！
複雜度是好的 O(lg N)。
接下來介紹 Pollard Rho 算法，這個算法可以求出一個整數x的因數們。
在進入Pollard Rho算法之前，先用Miller Rabin判斷這個x是不是質數。
如果不是質數的話，我們先定義一個亂數生成函數f(x) = a*x^2 + c。
接著，選擇任意一個開始值s，並且令 f1 = f(s), f2 = f( f(s) )
之後，如果在某個時刻，f2和f1的差和x的最大公因數不是1 ， 那我們就找到一個因數了，遞迴下去求解。
如果在某個時候，f1跟f2一樣了，我們就在換一個起始點 & 生成函數的a跟c值。
要不然，執行f1 = f(f1), f2 = f( f(f2) )，相當於是f1往前走一步，f2往前走兩步
數學家告訴我們，這樣做的話，複雜度大約只有O(N ^ (1/4))
#include <iostream> #include <cstdio> #include <vector> #include <algorithm> #include <cstring> #include <utility> #include <cmath> #include <ctime> #include <cstdlib> #include <queue> #include <stack> #include <set> #include <map> #include <cassert> #include <iomanip> #include <bitset> using namespace std; typedef long long LL; typedef long double ld; typedef pair<int,int> pii; typedef pair<LL,LL> pLL; typedef vector<int> vint; typedef vector<LL> vLL; typedef vector<pii> vpii; typedef vector<pLL> vpLL; #define SZ(x) ((int)(x).size()) #define ALL(x) (x).begin(),(x).end() #define F first #define S second #define MP make_pair #define PB push_back #define Si(x) scanf("%d",&(x)); #define Sii(x,y) scanf("%d %d",&(x),&(y)); #define Siii(x,y,z) scanf("%d %d %d",&(x),&(y),&(z)); #define Siiii(x,y,z,w) scanf("%d %d %d %d",&(x),&(y),&(z),&(w)); #define Siiiii(x,y,z,w,a) scanf("%d %d %d %d %d",&(x),&(y),&(z),&(w),&(a)); #define Siiiiii(x,y,z,w,a,b) scanf("%d %d %d %d %d %d",&(x),&(y),&(z),&(w),&(a),&(b)); #define SL(x) scanf("%lld",&(x)); #define SLL(x,y) scanf("%lld %lld",&(x),&(y)); #define SLLL(x,y,z) scanf("%lld %lld %lld",&(x),&(y),&(z)); #define SLLLL(x,y,z,w) scanf("%lld %lld %lld %lld",&(x),&(y),&(z),&(w)); #define SLLLLL(x,y,z,w,a) scanf("%lld %lld %lld %lld %lld",&(x),&(y),&(z),&(w),&(a)); #define SLLLLLL(x,y,z,w,a,b) scanf("%lld %lld %lld %lld %lld %lld",&(x),&(y),&(z),&(w),&(a),&(b)); #define Pi(x) printf("%d\n",(x)); #define Pii(x,y) printf("%d %d\n",(x),(y)); #define Piii(x,y,z) printf("%d %d %d\n",(x),(y),(z)); #define Piiii(x,y,z,w) printf("%d %d %d %d\n",(x),(y),(z),(w)); #define Piiiii(a,b,c,d,e) printf("%d %d %d %d %d\n",(a),(b),(c),(d),(e)); #define Piiiiii(a,b,c,d,e,f) printf("%d %d %d %d %d %d\n",(a),(b),(c),(d),(e),(f)); #define PL(x) printf("%lld\n",(x)*1LL); #define PLL(x,y) printf("%lld %lld\n",(x)*1LL,(y)*1LL); #define PLLL(x,y,z) printf("%lld %lld %lld\n",(x)*1LL,(y)*1LL,(z)*1LL); #define PLLLL(x,y,z,w) printf("%lld %lld %lld %lld\n",(x)*1LL,(y)*1LL,(z)*1LL,(w)*1LL); #define PLLLLL(a,b,c,d,e) printf("%lld %lld %lld %lld %lld\n",(a)*1LL,(b)*1LL,(c)*1LL,(d)*1LL,(e)*1LL); #define PLLLLLL(a,b,c,d,e,f) printf("%lld %lld %lld %lld %lld %lld\n",(a)*1LL,(b)*1LL,(c)*1LL,(d)*1LL,(e)*1LL,(f)*1LL); #define Pi1(x) printf("%d", (x)); #define PL1(x) printf("%lld",(x)); #define Pspace putchar(' '); #define Pendl puts(""); #define MEM0(x) memset( (x), 0, sizeof( (x) ) ) #define MEM1(x) memset( (x),-1, sizeof( (x) ) ) #define REP1(i,n) for (int i = 1; (n) >= i ; ++i) #define REP0(i,n) for (int i = 0; (n) > i ; ++i) #define IOS ios::sync_with_stdio(0); cin.tie(0); void Parr(int *arr,int L,int R) { for (int i=L;R>=i;i++) { printf("%d%c",arr[i]," \n"[i==R]); } } void Pvec(vint v) { for (int i=0;SZ(v)>i;i++) { printf("%d%c",v[i]," \n"[i==SZ(v)-1]); } } void Sarr(int *arr,int L,int R) { for (int i=L;R>=i;i++) { Si(arr[i]); } } //#define LL addd(LL a,LL b,LL mod) { return ((a)+(b)>=(mod))?(a)+(b)-(mod):(a)+(b); } const int G = (1LL<<31)-1; LL mull(LL a,LL b,LL mod) { if (a<G && b<G) return a*b%mod; LL ret = 0; LL now = a; while (b) { if (b&1) ret = addd(ret,now,mod); now = addd(now,now,mod); b>>=1; } return ret; } LL ppow(LL a,LL n,LL mod) { LL ret = 1; LL now = a; while (n) { if (n&1) ret = mull(ret,now,mod); now = mull(now,now,mod); n>>=1; } return ret; } LL gcd(LL a,LL b) { if (b==0) return a; else return gcd(b,a%b); } const bool PRIME = 1, COMPOSITE = 0; bool miller_rabin(LL n,LL a) { if (gcd(n,a) == n) return PRIME; else if (gcd(n,a) != 1) return COMPOSITE; LL d=n-1,r=0; while (d%2==0) { d >>= 1; ++r; } LL ret = ppow(a,d,n); if (ret == 1 || ret == n-1) return PRIME; while (r--) { ret = mull(ret,ret,n); if (ret == n-1) return PRIME; } return COMPOSITE; } bool isPrime(LL n) { LL as[7] = {2,325,9375,28178,450775,9780504,1795265022}; for (int i=0;7>i;++i) { //cout<<"i = "<<i<<endl; if (miller_rabin(n,as[i]) == COMPOSITE) return COMPOSITE; } return PRIME; } const LL C=2934852462451LL; const LL D=126871905557494LL; LL rnd=98134513458734897LL; LL myRnd() { return rnd = (rnd + C)^D; } LL a,c; LL doo(LL x,LL n) { return addd( mull( a, mull(x,x,n),n ),c ,n); } #define aabs(x) (x)>=0?(x):-(x) LL solve(LL n) { if (isPrime(n)) return n; if (!(n&1)) return 2; a=myRnd()%n; if (!a) a=1; c=myRnd()%n; //cout<<"a = "<<a<<" , c = "<<c<<endl; while (c == 0 || c==2) c=myRnd()%n; LL start = myRnd()%n; //cout<<"start = "<<start<<" , "; LL s1=doo(start,n);//cout<<"s1 = "<<s1<<endl; LL s2=doo(s1,n); //cout<<"n = "<<n<<" : "; while (true) { if (s1 == s2) { start = myRnd()%n; //a=myRnd()+1; a=myRnd()%n; if (!a) a=1; c=myRnd()%n; while (c == 0 || c==2) c=myRnd()%n; s1 = doo(start,n); s2 = doo(s1,n); continue; } //cout<<"s1 = "<<s1<<" , s2 = "<<s2<<endl; LL _=gcd(aabs(s1-s2),n); if (_ != 1) { return min(solve(_),solve(n/_)); } s1 = doo(s1,n); s2 = doo(s2,n); s2 = doo(s2,n); } } int main () { //srand(time(NULL)); int T; Si(T); while (T--) { LL n; SL(n); if (isPrime(n)) { puts("Prime"); } else { PL(solve(n)); } } }
沒有留言:
張貼留言