1. 程式人生 > >Miller-rabin素性測試與pollard-rho快速質因數分解

Miller-rabin素性測試與pollard-rho快速質因數分解

~~~~

Miller-Rabin

可以在O(log n)的時間內檢測一個數是否為質數。
理論上是一個概率演算法,但在一定範圍以內的數經過驗證全部可以判出來。
所以OI範圍內可以視作是一個確定演算法。

費馬小定理

若p是質數,且a<p,則有
ap1=1(modp)

其逆定理不一定成立,但不符合逆定理的一定不是素數。
根據這個檢測n的素性。

即:嘗試若干個a(a與n不能有公約數(尤拉定理條件),即對於質數n而言a不能是n的倍數,對於合數而言沒關係了)
檢測an1modn1

然而這還不夠強,會存在判不出的情況
對於Carmichael數,第一個是561,任意與他互質的數判出來都是Yes.

二次探測

顯然地,對於質數n與正整數a<=n
a2=1(modn),則有n|(a+1)(a1)
又因為n是質數,a<=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

據說時間複雜度是n0.25,反正比較玄學吧。

利用生日悖論的結論,作差撞因子會比直接列舉概率大(很多)(超級反直覺)。發現可以利用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;
    }
}