学习链接

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://pst.iorinn.moe/archives/miller-rabin-and-pollard-rho.html
许可： https://pst.iorinn.moe/license.html
若无特别说明，博客内的文章默认将采用 CC BY 4.0 许可协议进行许可☆

网站

加入讨论吧...