zoukankan      html  css  js  c++  java
  • Pollard-Rho 算法

    Pollard-Rho

    一种复杂度大概在 $ O(n^{frac 1 4} log n) $ 的分解质因数方法。

    Miller-Rabin

    给定一个 $ 10^{18} $ 范围的数,判断质数

    由于费马小定理,我们知道当 $ amod p eq 0 $ 且 $ p $ 为 质数 时,$ a^{p-1} equiv 1 pmod p $ ,而如果 $ p $ 不为质数,这东西有一定几率成立。

    所以我们有一个想法,随机一些 $ a $ ,对每个 $ a $ 用这个式子判断一下,如果有一个错误了它就直接是合数了。

    但是有些合数更特殊,对于 $ 1leq a < p $ 都成立这个式子,比如 $ 561 $ 。

    于是有了个二次探测定理。如果 $ x^2 equiv 1 pmod p $ 则有 $ xequiv pm 1 pmod p $ ,其中 $ p $ 为质数。

    证明很显然, $ p | x^2 - 1 o p | (x+1)(x-1) $ ,所以 $ p | (x+1) $ 或 $ p | (x-1) $。

    然而仍然,对于合数这个东西有一定几率不成立。

    这个算法的做法就是,检测 $ a^{p-1} $ 模 $ p $ 下的值,如果不是 1 或者 -1 那么必然不是合数。如果算出来是 1 并且 $ 2 | frac{p-1} 2 $ ,那么递归下去算 $ frac{p-1} 2 $ ,如果算出来是 $ p-1 $ 直接返回为质数。

    这样做并不是一定正确的,所以我们得多拿几个 $ a $ 来测。

    复杂度是 $ O(klog^2 n) $ 的。

    具体而言:

    • $ n leq 4 imes 10^9 $ 时, $ 2,7,61 $ 即可保证正确
    • $ n leq 10^{18} $ 时,$ 2,3,5,7,11,13,17,19,23 $ 可以保证正确
    • $ n leq 10^{19} $ 时,$ 2,3,5,7,11,13,17,19,23,29,31,37 $ 可以保证正确。

    虽然多数时候这个复杂度不会成为瓶颈,但是它仍然可以被优化。可以预处理出 $ p - 1 $ 的分解中 2 的个数 $ k $ 然后先算出 $ a{frac{p-1}{2k}} $ 再一步步平方回去,这样复杂度就可以优化掉 ksm 变成 $ O(klog n) $

    Pollard-Rho

    现在已经完全明白了如何判质数,那么怎么质因数分解呢?

    首先用 Miller-Rabin 判一下它是不是质数,是就直接 return 了。

    否则,它一定有一个 $ leq sqrt n $ 的因子 $ p $

    于是有一个乱搞思路,随机选出 $ n^{frac 1 4} $ 个数,两两做差再求 gcd ,期望情况下可以分解出来,但是这样并不优秀。

    所以考虑设函数 $ f(i) = (f(i-1)^2 + a) pmod n $ 。其中 $ a $ 是一个常数。我们发现,这样的话 $ f $ 必然会循环。它最后的形状很类似一个 $ ho $ ,先经历尾巴再进入一个循环。(同时这也是为啥名字叫 Pollard-Rho)。

    显然,即使是在$ mod p $ 的意义下它仍然是循环的,我们可以认为这个 函数 的数字出现是随机的,所以 $ mod p $ 意义下的最终环的长度是 $ O(p^frac 1 2) = O(n^ frac 1 4) $ 的。

    具体应该怎么做呢? Pollard-Rho 有两种实现方法:

    • Floyd 判环

      一种比较好写好理解的方法,但是貌似不如第二种优秀。

      做法是,每次让 $ s $ 走一步,$ t $ 走两步,然后求一下 $ y-x $ 和 $ n $ 的 $ gcd $ 。但是如果 $ x = y $ 了就说明进入了 $ mod n $ 的循环,就挂掉了,分解失败,需要再次随机继续来。

      如果我们可以在某一时刻发现 $ gcd(y-x , n) eq 1 $ 这就意味着我们找到了一个 $mod p $ 意义下的环,并且没有找到 $ mod n $ 意义下的环。成功进行了一次分解!

    好现在可以总结一下,Pollard-Rho 的本质是在环上跑,目的是在进入 $ mod n $ 的环之前进入 $ mod p $ 的环。如果做到了,那就可以分解出最小质因数 $ p $ 的某个倍数(并且也是 $ n $ 的约数),否则(比如同时进入这两个环)就需要调参重来。

    然后是一种更优秀的判环方法:

    • Brent 判环

      每次 $ x $ 走一步,当它走了 $ 2^k $ 时让 $ y = x $,然后 $ x $ 继续走。同时在 $ x $ 走了 $ 128 $ 的倍数步的时候,我们 check 一次。

      我也不知道为啥 (128) 这个数字很优秀,但是这样做确实很优秀。

    #include "iostream"
    #include "algorithm"
    #include "cstring"
    #include "cstdio"
    #include "cmath"
    #include "vector"
    using namespace std;
    typedef long long ll;
    
    ll mul( ll a , ll b , ll md ) {
        return ( a * b - (ll)( (long double)a / md * b + 0.5 ) * md + md ) % md;
    }
    ll Pow( ll x , ll a , ll md ) {
        ll cur = x % md , ans = 1;
        while( a ) {
            if( a & 1 ) ans = mul( ans , cur , md );
            cur = mul( cur , cur , md ) , a >>= 1;
        }
        return ans;
    }
    
    const int ck[] = {2,3,5,7,11,13,17,19,23} , _l = 8;
    bool miller( ll n ) {
        if( n == 1 ) return false;
        ll t = n - 1; int cn = 0;
        while( !( t & 1 ) ) t >>= 1 , ++ cn;
        for( int i = 0 ; i < _l ; ++ i ) {
            if( n == ck[i] ) return true;
            ll a = Pow( ck[i] , t , n ) , nex = a;
            for( int j = 1 ; j <= cn ; ++ j ) {
                nex = mul( a , a , n );
                if( nex == 1 && a != 1 && a != n - 1 ) return false;
                a = nex;
            }
            if( a != 1 ) return false;
        }
        return true;
    }
    
    inline ll f( ll x , ll c , ll md ) {
        return ( mul( x , x , md ) + c ) % md;
    }
    
    inline ll _rand(  ) {
        return (ll) rand() << 32 | rand();
    }
    inline ll _randw() {
        return (ll)rand() << 48 | (ll)rand() << 32 | rand() << 16 | rand();
    }
    inline ll _abs( ll x ) {
        return x > 0 ? x : -x;
    }
    inline ll gcd( ll a , ll b ) {
        return !b ? a : gcd( b , a % b );
    }
    
    inline ll pollard_rho( ll n ) {
        ll s = 0 , t = 0 , c = _rand() % ( n - 1 ) + 1 , val = 1;
        for( int cir = 1 ; ; cir <<= 1 , s = t , val = 1 ) {
            for( int i = 0 ; i < cir ; ++ i ) {
                t = f( t , c , n ) , val = mul( val , _abs( t - s ) , n );
                if( i % 127 == 0 ) {
                    ll g = gcd( val , n );
                    if( g != 1 ) return g;
                }
            }
            ll g = gcd( val , n );
            if( g != 1 ) return g;
        }
    }
    
    vector<ll> divs;
    inline void analyze( ll n ) {
        if( n == 1 ) return;
        if( miller( n ) ) { divs.push_back( n ); return; }
        ll d = n;
        while( d == n ) d = pollard_rho( n );
        n /= d;
        analyze( n ) , analyze( d );
    }
    
    int main() {
        srand( time(0) );
        int T;cin >> T;
        ll n;
        while( T-- ) {
            scanf("%lld",&n);
            if( miller( n ) ) { puts("Prime"); continue; }
            divs.clear();
            analyze( n );
            printf("%lld
    ",*max_element( divs.begin() , divs.end() ) );
        }
    }
    
  • 相关阅读:
    .NetCore Grpc 客服端 工厂模式配置授权
    DOCKER 拉取 dotnet 镜像太慢 docker pull mcr.microsoft.com too slow
    Introducing .NET 5
    VSCode 出现错误 System.IO.IOException: The configured user limit (128) on the number of inotify instances has been reached.
    Omnisharp VsCode Attaching to remote processes
    zookeeper3.5.5 centos7 完全分布式 搭建随记
    Hadoop2.7.7 centos7 完全分布式 配置与问题随记
    MySQL索引 索引分类 最左前缀原则 覆盖索引 索引下推 联合索引顺序
    SQL基础随记3 范式 键
    MySQL调优 优化需要考虑哪些方面
  • 原文地址:https://www.cnblogs.com/yijan/p/12384527.html
Copyright © 2011-2022 走看看