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() ) );
        }
    }
    
  • 相关阅读:
    【WCF--初入江湖】04 WCF通信模式
    【WCF--初入江湖】03 配置服务
    c++输出左右对齐设置
    setw()函数
    clion更改大括号的位置
    emacs org-mode 中文手册精简版(纯小白)
    c++ string 类型 大小写转换 
    C++中string类型的find 函数
    string类型 C++
    统计单词数---单词与空格
  • 原文地址:https://www.cnblogs.com/yijan/p/12384527.html
Copyright © 2011-2022 走看看