Pollard_Rho

注：有参考于LinearODE的博客

引入

现有一个问题：如何求一个大数的一个因子（最大、最小、任意一个），而数据范围是64位带符号整形，并且有100次询问。那么枚举显然是不现实的，我们就需要一个更加高效的算法。在学习Pollard_Rho之前，需要先学会Miller_Rabin素数测试。

算法介绍

注：这一部分若未特殊说明，则我们希望分解$n$，并且所有运算都在$\mod n$意义下进行。

随机采样

Pollard_Rho算法思想依旧是基于枚举因数。但是不能正常枚举，于是考虑随机化。但是直接随机一个数的话，概率未免太低。而考虑选取两个数，看它们的差，那概率就会大很多（因为第二个数$i+k$和$i-k$都可以）。那么当我们选取足够多的数时，判断其中存在两数相减是我们需要的一个因数的概率就很大了（这得益于组合随机采样，类似于生日悖论的原理）。

Pollard的伪随机数列

上述的随机采样似乎是一个很好的点子。但是我们需要枚举两两之间的差，并且为了正确性，又需要足够多的随机数（大约可能会接近$O(n^{\frac{1}{4}})$）。这样的话，并没有实质性地对时间复杂度进行优化（依旧接近$O(\sqrt n)$）。

Pollard设计了这样一个数列：随机两个数$a_0$与$c$，$a_i=a_{i-1}^2+c$。用这个数列代替上面说道的两两枚举差有很好的表现。（LinearODE大佬提到这有关于曼德勃罗集，然而本蒟蒻并不知道，于是就不提了。）

而由于是在模意义下运算，所以这个数列一定会循环。有趣的是，因为这个数列轨迹像$\rho$，所以算法名就叫做Pollard_Rho。

用Floyd和Brent’s对算法进行优化

上面提到，Pollard的伪随机数列一定步数后会循环。那么出现环的时候我们就需要退出。标题中的两个都是不错的算法。

可以参考Floyd判圈，和Brent’s判圈。

下面使用Floyd判圈。

算法实现

通过上面的描述，我相信大家已经可以写出这么一段代码了：

LL Pollard_Rho( LL T, LL Seed ) {
    LL x, y;
    x = y = rand() % T;
    y = ( ( __int128 ) y * y % T + Seed ) % T;
    while( x != y ) {
        LL d = gcd( ( x - y > 0 ) ? x - y : y - x, T );
        if( d > 1 ) return d;
        x = ( ( __int128 ) x * x % T + Seed ) % T;
        y = ( ( __int128 ) y * y % T + Seed ) % T;
        y = ( ( __int128 ) y * y % T + Seed ) % T;
    }
    return T;
}

需要注意的是，$x$和$y$的取值范围需要在$0$至$T-1$之间，并且一开始不要$x=f(x), y=f(f(y))$，只要$y=f(y)$就可以了。否则的话这个算法会在$T=4$的时候出现一点问题（$T=4$时最长的循环节长度只有$2​$，所以必须要最前面两项）。

如果一次不成功，那么改变$Seed$再次尝试即可。

94分

到此为止，已经可以在【模板】Pollard-Rho算法中获得$94$分了。完整代码如下：

#include <bits/stdc++.h>
#define LL long long
using namespace std;

LL FastPow( LL x, LL y, LL Mod ) {
    LL Ans = 1;
    for( ; y; y >>= 1, x = ( __int128 ) x * x % Mod )
        if( y & 1 )
            Ans = ( __int128 ) Ans * x % Mod;
    return Ans;
}

bool Miller_Rabin( LL n, LL Times ) {
    if( n == 2 || n == 3 ) return true;
    if( n < 2 || n % 2 == 0 ) return false;
    LL m = n - 1, p = 0;
    while( !( m & 1 ) ) {
        m >>= 1;
        ++p;
    }
    if( n == 2 || n == 3 ) return true;
    if( n < 2 || n % 2 == 0 ) return false;
    for( LL i = 0; i < Times; ++i ) {
        LL a = rand() % ( n - 2 ) + 2;
        a = FastPow( a, m, n );
        for( LL i = 1; i <= p; ++i ) {
            LL b = a;
            a = ( __int128 ) a * a % n;
            if( a == 1 && b != 1 && b != n - 1 ) return false;
        }
        if( a != 1 ) return false;
    }
    return true;
}

LL gcd( LL x, LL y ) {
    LL m = x % y;
    while( m ) {
        x = y; y = m; m = x % y;
    }
    return y;
}

LL Pollard_Rho( LL T, LL Seed ) {
    LL x, y;
    x = y = rand() % T;
    y = ( ( __int128 ) y * y % T + Seed ) % T;
    while( x != y ) {
        LL d = gcd( ( y - x ) > 0 ? y - x : x - y, T );
        if( d > 1 ) return d;
        x = ( ( __int128 ) x * x % T + Seed ) % T;
        y = ( ( __int128 ) y * y % T + Seed ) % T;
        y = ( ( __int128 ) y * y % T + Seed ) % T;
    }
    return T;
}

LL Ans;

void Find( LL x, LL Seed ) {
    if( x <= Ans ) return;
    if( Miller_Rabin( x, 8 ) ) {
        Ans = x;
        return;
    }
    LL y = x;
    while( x == y )    y = Pollard_Rho( x, Seed-- );
    while( !( x % y ) ) x /= y;
    Find( x, Seed ); Find( y, Seed );
    return;
}

int main() {
    srand( ( unsigned LL ) "我不相信会卡19260817这个模数" );
    LL T; scanf( "%lld", &T );
    for( ; T--; ) {
        LL x; scanf( "%lld", &x );
        Ans = 1;
        Find( x, 19260817 );
        if( Ans == x ) printf( "Prime\n" ); else printf( "%lld\n", Ans );
    }
    return 0;
}

然后剩下的$6$分去哪里了？（实测最后一个点需要6s+）

100分

如果你相信自己的卡常技术，那么可以略过这最后一部分不看了。

我们发现程序中多次求了$gcd$。而$gcd$的速度不是那么友好。于是我们需要优化。又发现如果$gcd(a,b)>1$，那么一定有$gcd(ac\%b,b)>0$。其中$c>0$。所以我们就可以将一些数打包，然后一起求$gcd$。需要注意的是，$c=0$的情况需要特判。

所以最后我们就得到了满分代码：

#include <bits/stdc++.h>
#define LL long long
using namespace std;

LL FastPow( LL x, LL y, LL Mod ) {
    LL Ans = 1;
    for( ; y; y >>= 1, x = ( __int128 ) x * x % Mod )
        if( y & 1 )
            Ans = ( __int128 ) Ans * x % Mod;
    return Ans;
}

bool Miller_Rabin( LL n, LL Times ) {
    if( n == 2 || n == 3 ) return true;
    if( n < 2 || n % 2 == 0 ) return false;
    LL m = n - 1, p = 0;
    while( !( m & 1 ) ) {
        m >>= 1;
        ++p;
    }
    if( n == 2 || n == 3 ) return true;
    if( n < 2 || n % 2 == 0 ) return false;
    for( LL i = 0; i < Times; ++i ) {
        LL a = rand() % ( n - 2 ) + 2;
        a = FastPow( a, m, n );
        for( LL i = 1; i <= p; ++i ) {
            LL b = a;
            a = ( __int128 ) a * a % n;
            if( a == 1 && b != 1 && b != n - 1 ) return false;
        }
        if( a != 1 ) return false;
    }
    return true;
}

LL gcd( LL x, LL y ) {
    LL m = x % y;
    while( m ) {
        x = y; y = m; m = x % y;
    }
    return y;
}

LL Pollard_Rho( LL T, LL Seed ) {
    LL x, y;
    x = y = rand() % T;
    y = ( ( __int128 ) y * y % T + Seed ) % T;
    LL Cnt = 0, Mul = 1;
    while( x != y ) {
        LL Last = Mul;
        Mul = ( __int128 ) Mul * ( ( x - y > 0 ) ? x - y : y - x ) % T; ++Cnt;
        if( Mul == 0 ) {
            LL d = gcd( Last, T );
            return d;
        }
        if( Cnt == 127 ) {
            Cnt = 0;
            LL d = gcd( Mul, T );
            Mul = 1;
            if( d > 1 ) return d;
        }
        x = ( ( __int128 ) x * x % T + Seed ) % T;
        y = ( ( __int128 ) y * y % T + Seed ) % T;
        y = ( ( __int128 ) y * y % T + Seed ) % T;
    }
    if( Mul > 1 ) {
        LL d = gcd( Mul, T );
        if( d > 1 ) return d;
    }
    return T;
}

LL Ans;

void Find( LL x, LL Seed ) {
    if( x <= Ans ) return;
    if( Miller_Rabin( x, 8 ) ) {
        Ans = x;
        return;
    }
    LL y = x;
    while( x == y )    y = Pollard_Rho( x, Seed-- );
    while( !( x % y ) ) x /= y;
    Find( x, Seed ); Find( y, Seed );
    return;
}

int main() {
    srand( ( unsigned LL ) "我不相信会卡19260817这个模数" );
    LL T; scanf( "%lld", &T );
    for( ; T--; ) {
        LL x; scanf( "%lld", &x );
        Ans = 1;
        Find( x, 19260817 );
        if( Ans == x ) printf( "Prime\n" ); else printf( "%lld\n", Ans );
    }
    return 0;
}