Miller-Rabin 与 Pollard-Rho

ajcxsu
ajcxsu 2018年09月10日
  • 36 次阅读

一个是素数判断算法,一个是分解质因数算法。
后者相当的玄学。

学习链接

https://www.cnblogs.com/Doggu/p/MillerRabin_PollardRho.html?utm_source=itdadao&utm_medium=referral

Miller-Rabin

对于每个数,进行随机+费马小定理+二次探测定理(欧几里得引理)的Miller-Rabin测试。 每次测试错误几率约为 $25\%$ ,则 $k$ 次测试之后成功率约为 $1-25\% ^k$ 。

费马小定理的逆定理:若 $a^{n-1}\equiv 1 \mod n \ \text{(a,b互质)}$ ,则 $n$ 为素数。 这个逆定理在绝大多数情况下是对的,除了“伪素数”。 在这其中还有一类非常顽强的伪素数,对所有的 $a$ 都有该逆定理成立,也被称为Carmichael数。

使用二次探测定理来进行进一步检测: 若 $x^2 \equiv 1 \mod n$ , $n$ 为质数,则 $x$ 有唯二解: $x_1=1,\ x_2=n-1$ 。

对于long long类型,需要使用快速乘来进行Miller-Rabin,否则会溢出。

模板链接:https://loj.ac/problem/143

Code

// Code by ajcxsu
// Problem: miller rabin

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;

ll qass(ll x, ll y, ll mo) {
    return (__int128)x*y%mo;
}

ll qpow(ll x, ll y, ll mo) {
    ll ret=1;
    x%=mo;
    while(y) {
        if(y&1ll) ret=qass(x, ret, mo);
        x=qass(x, x, mo), y>>=1ll;
    }
    return ret;
}

bool miller_rabin(ll n) {
    if(n==1) return false;
    else if(n==2) return true;
    ll u=n-1, k=0;
    ll pre, x;
    while(!(u&1ll)) u>>=1ll, k++;
    for(int i=0;i<20;i++) {
        x=rand()%(n-2)+2;
        x=qpow(x, u, n);
        pre=x;
        for(int j=0;j<k;j++) {
            x=qass(x,x,n);
            if(x==1 && pre!=1 && pre!=n-1) return false;
            pre=x;
        }
        if(x!=1) return false;
    }
    return true;
}

int main() {
    ios::sync_with_stdio(false);
    srand(time(NULL));
    ll na;
    while(cin>>na) {
        if(miller_rabin(na)) cout<<"Y"<<endl;
        else cout<<"N"<<endl;
    }
    return 0;
}

Pollard-Rho

用于求解 $n$ 的任意一个不为 $1$ 或 $n$ 的约数。

由生日悖论,若有随机数列 ${a_i}$ ,则期望数列在 $\sqrt{n}$ 的大小中存在 $i,j$ 有 $[|a_i-a_j|,n]>1$ 。

这样子的话,我们生成随机数列,再判断是否存在 $[|a_i-a_{i-1}|,n]>1$ 。若有,则约数得出。

对于每个约数都进行进一步分解,由于最多存在 $\log n$ 个约数,大体复杂度可以保证。

对于随机数列的生成法,有两种: $$a_i=({a_{i-1}}^2+c) \% n$$ $$a_i=(a_{i-1}+c) \% n$$

第一种函数在处理较小的数字时得到的随机数列不完全,第二种函数处理效率过低(由实践得出,我也不知道为什么)。
因此两种函数都得用上,才得到一个比较理想时间复杂度的程序。

并不认为考场上能打出来这个算法。

Code

// Code by ajcxsu
// Problem: miller rabin

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;

ll qass(ll x, ll y, ll mo) {
    ll ret=0;
    while(y) {
        if(y&1ll) ret=(ret+x)%mo;
        x=(x+x)%mo, y>>=1ll;
    }
    return ret;
}

ll qpow(ll x, ll y, ll mo) {
    ll ret=1;
    x%=mo;
    while(y) {
        if(y&1ll) ret=qass(x, ret, mo);
        x=qass(x, x, mo), y>>=1ll;
    }
    return ret;
}

bool miller_rabin(ll n) {
    if(n==1) return false;
    else if(n==2) return true;
    ll u=n-1, k=0;
    ll pre, x;
    while(!(u&1)) u>>=1, k++;
    for(int i=0;i<5;i++) {
        x=rand()%(n-2)+2;
        x=qpow(x, u, n);
        pre=x;
        for(int j=0;j<k;j++) {
            x=qass(x,x,n);
            if(x==1 && pre!=1 && pre!=n-1) return false;
            pre=x;
        }
        if(x!=1) return false;
    }
    return true;
}

ll gcd(ll a, ll b) { return b?gcd(b, a%b):a; }

const int N=100;
ll fac[N];
ll stk[N], t;
int p;
ll c;
inline ll f(ll x, ll mo) { return (qass(x,x,mo)+c)%mo; }
inline ll f2(ll x, ll mo) { return (x+c)%mo; }
bool pollard_rho(ll x) {
    if(miller_rabin(x)) {
        fac[p++]=x;
        return true;
    }
    c=rand();
    ll a=rand()%x, b=a, p;
    while(1) {
        a=f(a, x);
        b=f(f(b, x), x);
        if(a==b) break;
        p=gcd(abs(a-b), x);
        if(p!=1 && p) {
            stk[++t]=p, stk[++t]=x/p;
            return true;
        }
    }
    while(1) {
        a=f2(a, x);
        b=f2(f2(b, x), x);
        if(a==b) break;
        p=gcd(abs(a-b), x);
        if(p!=1 && p) {
            stk[++t]=p, stk[++t]=x/p;
            return true;
        }
    }
    return false;
}

int main() {
    ios::sync_with_stdio(false);
    srand(time(NULL));
    int n;
    cin>>n;
    while(n--) {
        ll na;
        cin>>na;
        p=0;
        stk[++t]=na;
        while(t) {
            ll deal=stk[t--];
            while(!pollard_rho(deal));
        }
        for(int i=0;i<p;i++) cout<<fac[i]<<' ';
        cout<<endl;
    }
    return 0;
}

本文链接:https://acxblog.site/archives/miller-rabin-and-pollard-rho.html
本文采用 CC BY-NC-SA 3.0 协议进行许可