『学习笔记』PollardRho 算法 题解

前言

说是学习笔记,其实窝并没有打算写太多(太麻烦了,而且我的理解还不是特别深,可能也写不清楚),所以打算大概写两句,然后贴个板子。

前置知识

  • \(Miller-Rabin\) 素性测试。
  • 倍增基础应用。

\(Miller-Rabin\) 素性测试

我们知道有费马小定理: \(a^{p - 1} \equiv 1\pmod p\)\(p\) 是质数且 \(a\) 是小于 \(p\) 的正整数。

那么费马小定理的逆定理是否成立呢?

答案是否定的,\(2^{340} \equiv 1 \pmod{341}\),但是 \(341 = 11 \times 13\)

这样的数被称为伪素数。事实上,这样的数有很多。

有人会问:如果我选取小于 \(n\) 的所有质数当作底数都判断一遍,是否能让这个逆定理变得正确呢?

很遗憾,就算取再多的底数总是会有合数通过条件,例如 561,所以我们并不能用费马小定理的逆定理当作判断素数的方法。

这时 \(Miller\)\(Rabbin\) 出现了,他们共同提出了二次探测定理

这个定理是这样的:

假设有 \(x^2 \equiv 1 \pmod p\),那么 \(x = 1\)\(x = p - 1\)

证明很简单,原式即为:

\[x^2 - 1 \equiv 0 \pmod p \\ (x + 1)(x - 1) \equiv 0 \pmod p \]

所以 \(x = 1\)\(x = p - 1\)

注意:\(x\) 等于上述两个数中的任意一个都表示它有可能是质数。

下面我们来演示一下上面的定理如何应用在费马素性测试上:

\[2^{340} \equiv 1 \pmod {341} \]

\(2^{340} = (2^{170})^2\),把上面的定理套进去,并计算 \(2^{170}\ \ mod\ \ 341\) 的值,我们发现是 1。于是继续套用公式,计算 \(2^{85}\ \ mod\ \ 341\),发现其等于 32。这时,341 的真正面目就暴露了出来,它是一个合数!

所以我们的计算过程就是不断地令指数除以 2,然后对模数取模进行判断。

当出现取模之后的结果不为 1 或 \(p - 1\) 时,证明它不是质数。如果指数已经是奇数了,那么它就是一个质数。

需要一提的是,单独使用 2 这个底数进行判断还是会有误差,例如 2047 就能通过以 2 为底数的 \(MR\) 素性测试。

上述就是 \(Miller-Rabin\) 素性测试的原理了。

\(Code:\)

namespace Miller_Rabin{
    const int p[9] = {0, 2, 3, 5, 7, 11, 13, 17, 19}, totp = 8;

    inline ll mul(ll a, ll b, ll mod) {return (a * b - (ll)((long double)a * b / mod) * mod + mod) % mod;}

    inline ll qpow(ll a, ll b, ll mod){
        ll res = 1; a  %= mod;
        while(b){
            if(b & 1) res = mul(res, a, mod);
            a = mul(a, a, mod), b >>= 1;
        }
        return res;
    }

    inline bool check(ll x, ll p){
        if((x % p == 0) || qpow(p, x - 1, x) != 1) return 0;
        ll k = x - 1;
        while(!(k & 1)){
            ll t = qpow(p, k >>= 1, x);
            if(t != 1 && t != (x - 1)) return 0;
            if(t == x - 1) return 1;
        }
        return 1;
    }

    inline bool MR(ll x){
        if(x <= 1) return 0;
        if(x <= 3) return 1;
        if(!(x % 2) || !(x % 3)) return 0;
        for(int i = 1; i <= totp; ++i){
            if(x == p[i]) return 1;
            if(!check(x, p[i])) return 0;
        }
        return 1;
    }
}
using namespace Miller_Rabin;

倍增就不多说了。

\(Pollard-Rho\)

核心思想:每次猜一个数,判断是否是 \(n\) 的约数,并递归猜下去,如果是约数且是质数的话就更新答案。

考虑构造一个随机数列 \(x_i = x_{i - 1}^2 + c\)通过列举不难发现,这个东西是一定会形成环的。

我们可以在上面倍增找点。当找到当前的环上的点,判断其是否是 \(n\) 的约数。

(不想写太多了,各位巨佬不喜勿喷)

代码

#include <bits/stdc++.h>
#define ll long long

using namespace std;

namespace IO{
    inline ll read(){
        ll x = 0;
        char ch = getchar();
        while(!isdigit(ch)) ch = getchar();
        while(isdigit(ch)) x = (x << 3) + (x << 1) + ch - '0', ch = getchar();
        return x;
    }

    template <typename T> inline void write(T x){
        if(x > 9) write(x / 10);
        putchar(x % 10 + '0');
    }
}
using namespace IO;

namespace Miller_Rabin{
    const int p[9] = {0, 2, 3, 5, 7, 11, 13, 17, 19}, totp = 8;

    inline ll mul(ll a, ll b, ll mod) {return (a * b - (ll)((long double)a * b / mod) * mod + mod) % mod;}

    inline ll qpow(ll a, ll b, ll mod){
        ll res = 1; a  %= mod;
        while(b){
            if(b & 1) res = mul(res, a, mod);
            a = mul(a, a, mod), b >>= 1;
        }
        return res;
    }

    inline bool check(ll x, ll p){
        if((x % p == 0) || qpow(p, x - 1, x) != 1) return 0;
        ll k = x - 1;
        while(!(k & 1)){
            ll t = qpow(p, k >>= 1, x);
            if(t != 1 && t != (x - 1)) return 0;
            if(t == x - 1) return 1;
        }
        return 1;
    }

    inline bool MR(ll x){
        if(x <= 1) return 0;
        if(x <= 3) return 1;
        if(!(x % 2) || !(x % 3)) return 0;
        for(int i = 1; i <= totp; ++i){
            if(x == p[i]) return 1;
            if(!check(x, p[i])) return 0;
        }
        return 1;
    }
}
using namespace Miller_Rabin;

namespace Pollard_Rho{
    #define Rand(x) (1ll * rand() * rand() % (x) + 1)

    inline ll add(ll x, ll mod) {return x >= mod ? x - mod : x;}
    inline ll gcd(ll a, ll b) {return !b ? a : gcd(b, a % b);}

    inline ll PR(ll n){
        ll x0 = Rand(n - 1), x = x0, c = Rand(n), s = 1;
        ll k = 1, step = 0, d;
        while(true){
            x = add(mul(x, x, n) + c, n);
            s = mul(s, abs(x - x0), n);
            if(!s || (x == x0)) return n;
            if((++step) == k){
                if((d = gcd(s, n)) != 1) return d;
                x0 = x, k <<= 1;
            }
        }
    }

    inline void solve(ll x, ll &ans){
        if((x == 1) || x < ans) return;
        if(MR(x)) return ans = max(ans, x), void();
        ll y = PR(x);
        while(x == y) y = PR(x);
        while(x % y == 0) x /= y;
        solve(x, ans), solve(y, ans);
    }

    inline ll calc(ll n){
        ll ans = 0;
        solve(n, ans);
        return ans;
    }
}
using namespace Pollard_Rho;

ll T, n;

signed main(){
    srand(20060707);
    T = read();
    while(T--){
        ll n = read();
        if(MR(n)) puts("Prime");
        else write(calc(n)), puts("");
    }
    return 0;
}

\[\_EOF\_ \]

原文地址:https://www.cnblogs.com/xixike/p/15716686.html