zoukankan      html  css  js  c++  java
  • 『学习笔记』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\_ \]

  • 相关阅读:
    海伦公式
    简单的博弈
    Hello World 代码
    Hello world
    99999999海岛帝国后传:算法大会
    判断质数
    idea plantUML配置
    测试用例评审
    如何编写有效测试用例
    测试用例设计——场景分析法
  • 原文地址:https://www.cnblogs.com/xixike/p/15716686.html
Copyright © 2011-2022 走看看