BZOJ 4802: 欧拉函数(大数因数分解算法 Pollard_rho 和素数测试算法 Miller_Rabin)
Description
已知N,求phi(N)
Input
正整數N。N<=10^18
Output
輸出phi(N)
Sample Input
8
Sample Output
4
Solution
-
我們知道 φ(n)=n(1?1p1)(1?1p2)???(1?1pk)φ(n)=n(1-\frac{1}{p_1})(1-\frac{1}{p_2})···(1-\frac{1}{p_k})φ(n)=n(1?p1?1?)(1?p2?1?)???(1?pk?1?)
-
其中 n=∏pikin=\prod{p_i^{k_i}}n=∏piki?? ,pip_ipi? 為互不相同的質數。
-
很自然,我們需要將 nnn 分解質因數,但是 n≤1018n\le10^{18}n≤1018 ,很大。
-
于是我們可以使用 Pollard_rho 算法。
-
其主要思想是找到兩個數 p,qp,qp,q ,使得 n=p?qn=p*qn=p?q ,在不斷繼續分解。
-
設遞歸過程 proc(x)proc(x)proc(x) ,若 xxx 是質數(用 Miller_Rabin 算法判斷),則 xxx 是 nnn 的一個質因數,并退出過程。
-
否則找到一個數 yyy ,使 gcd(x,y)>1gcd(x,y)>1gcd(x,y)>1 ,繼續遞歸 proc(y)proc(y)proc(y) 、proc(xy)proc(\frac{x}{y})proc(yx?) 。
-
那么怎樣找到這個數呢?
-
我們用到了隨機化的思想和生日悖論,即找到兩個數 x1,x2x1,x2x1,x2 使得 gcd(∣x1?x2∣,x)>1gcd(|x1-x2|,x)>1gcd(∣x1?x2∣,x)>1 ,這樣概率就會大得多。
-
設變換函數 f(x)=(x2+c)%mf(x)=(x^2+c)\%mf(x)=(x2+c)%m ,其中 ccc 是一個隨機變量,那么令 x1,x2x1,x2x1,x2 沿著該函數一直走。
-
可這樣走是會成環走回原點的,如下圖:
-
它的形狀形似希臘字母 ρρρ ,該算法名即得于此。
-
判環的話可以用 Floyd 判圈算法:x1=f(x1)x1=f(x1)x1=f(x1)x2=f(f(x2))x2=f(f(x2))x2=f(f(x2))
-
即 x2x2x2 以兩倍速度跑,當 x1=x2x1=x2x1=x2 時重構即可。
-
可證單 Pollard_rho 算法的期望時間復雜度為 O(n14)O(n^{\frac{1}{4}})O(n41?) 。
-
然而其中還有一個算法 Miller_Rabin (素數測試算法)還沒介紹。
-
我們當然可以 O(n)O(n)O(n) 線篩出素數,但是在這種時間空間都不允許的情況下,如何快速判定單個素數呢?
-
根據費馬小定理,當 nnn 為素數時,有: an?1≡1(modn)a^{n-1}\equiv1(mod\ n)an?1≡1(mod?n)
-
但是當正整數 nnn 滿足 an?1≡1(modn)a^{n-1}\equiv1(mod\ n)an?1≡1(mod?n) 時不一定 nnn 就是素數(即逆定理不成立)。
-
但是我們想如果對于很多個 aaa 都滿足上式,是否就能說明 nnn 就是素數了呢?
-
我們發現這樣錯誤率依然很高,這時二次探測定理就派上用場了。
-
若一個數 xxx 滿足: x2≡1(modp)x^2\equiv1(mod\ p)x2≡1(mod?p) ,則合法的解就只有 x=1x=1x=1 或 x=p?1x=p-1x=p?1 。
-
證明很簡單,就是因為 x2?1=(x+1)(x?1)x^2-1=(x+1)(x-1)x2?1=(x+1)(x?1) ,合法解顯然。
-
那么 an?1a^{n-1}an?1 就可以拆成 ad?2ta^{d*2^t}ad?2t ,其中 ddd 是奇數。
-
那么我們先用 ada^dad 檢測,之后每次平方,若 (ax)2≡1{(a^x)}^2\equiv1(ax)2≡1 則 ax≡1或p?1(modp)a^x\equiv1或p-1(mod\ p)ax≡1或p?1(mod?p) 。
-
如此一來,檢測的成功率就大了不少,更容易檢測出是否為質數了。
-
我們只需隨機幾個 aaa (或拿幾個質數)來判斷即可,復雜度約為 log23(n)log_2^3(n)log23?(n)。
-
至此,我們就成功地實現了算法 Miller_Rabin 和 Pollard_rho 。
Code
#include<cstdio> #include<algorithm> #include<iostream> #include<ctime> #include<cstdlib> using namespace std; typedef long long LL; typedef long double LD; const int N=65,P=5,prime[P]={2,3,7,61,24251}; LL n,pn; LL f[N]; inline LL mul(LL a,LL b,LL p) {a%=p,b%=p;LL c=(LD)a*b/p;c=a*b-c*p;if(c<0) c+=p; elseif(c>p) c-=p;return c; } inline LL ksm(LL x,LL y,LL p) {LL s=1;while(y){if(y&1) s=mul(s,x,p);x=mul(x,x,p);y>>=1;}return s; } inline LL twice(LL a,LL p) {LL d=p-1;int t=0;while(!(d&1)) d>>=1,t++;//p-1=d*2^tLL x,y;x=y=ksm(a,d,p);while(t--){y=mul(x,x,p);if(y==1 && x^1 && x^p-1) return 0;x=y;}return y; } inline LL random(LL up) {return (LL)rand()*rand()%up; } LL gcd(LL x,LL y) {return !y?x:gcd(y,x%y); } inline LL abs(LL x,LL y) {return x<y?y-x:x-y; } inline bool Miller_Rabin(LL x) {for(int i=0;i<P;i++){if(x==prime[i]) return true;if(twice(prime[i],x)^1) return false;}return true; } inline LL trans(LL x,LL y,LL z) {return (mul(x,x,z)+y)%z; } void Pollard_rho(LL m) {if(Miller_Rabin(m)){f[++f[0]]=m;return;}LL x1=0,x2=0,c=0,p=1;while(p==1 || p==m){x1=trans(x1,c,m);x2=trans(trans(x2,c,m),c,m);while(x1==x2){c=random(m);x1=x2=random(m);x2=trans(x2,c,m);}p=gcd(abs(x1-x2),m);}Pollard_rho(p);Pollard_rho(m/p); } inline LL getphi(LL m) {sort(f+1,f+1+f[0]);f[0]=unique(f+1,f+1+f[0])-(f+1);LL sum=m;for(int i=1;i<=f[0];i++)sum=sum/f[i]*(f[i]-1);return sum; } int main() {scanf("%lld",&n);if(n==1) return 0&puts("1");srand(time(NULL));Pollard_rho(n);pn=getphi(n);printf("%lld",pn);return 0; }總結
以上是生活随笔為你收集整理的BZOJ 4802: 欧拉函数(大数因数分解算法 Pollard_rho 和素数测试算法 Miller_Rabin)的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: JZOJ 5700. 【gdoi2018
- 下一篇: JZOJ 5794. 2018.08.1