BZOJ3667:Rabin-Miller算法

传送门

题意：
大素数判断并分解。

题解：
Robin-Miller随机算法。
http://blog.csdn.net/thy_asdf/article/details/51347390

—————————————————————————————————————————
发现以前写得有点水啊，现在来重写一发。
Rabin-Miller检测素数：
首先，检测前9个质数prn−1i=1(modn)
其次，发现除了2，n−1为偶数（否则n不为质数）。那么可以顺便做二次检测：
a2=1(modn)，当n为质数时当且仅当a=1或a=n−1。

Pollard-Rho 大数分解：
设n=p⋅q，x,y均为随机数。
根据生日悖论，随机选取n√个x,y,存在|xi−yj|为p,q的概率为50%。
需要两两枚举n次。
考虑gcd(|xi−yi|,n)为p,q的概率，这样的差p+q−2=n√种，那么只需要选取n14个数，两两枚举n√次。

神奇的地方来了：
我们使用一个函数来生成伪随机。
换句话说，我们不断地 使用函数 使用函数 f 来生成（看上去或者说感觉上像的）随机数。
并不是 所有的函都能够这样做，但一个神奇的函数可以。

f(x)=(x2+a)%n

( 我们可以自己指定a，也可以用 rand()随机函数来生成，这不是重点 )
我们从x1=2 或者其他数开始，让x2=f(x1),x3=f(x2)...，生成规则为xn+1=f(xn)

我们来看下面的伪代码:

x:2
while (.. exit criterion .. )
y := f(x);
p := GCD( | y - x | , N);
if ( p > 1)
return “Found factor: p”; x := y;
end
return “Failed. :-(”

假设 N = 55 , f(x)=(x2+2)%55

xn xn+1 gcd(xn,xn+1,N) 2 6 1 6 38 1 38 16 1 16 36 5

你可以发现对于大部分的数据这个算法能够正常运行，但是对于某些数据，它将会进入无限循环。为什么呢？这是因为存在f环的原因。
当它发生的时候，我们会在一个有限数集中进行无限循环
例如，我们可以构造一个伪随机函数并生成如下伪随机数：

2,10,16,23,29,13,16,23,29,13……

在这个例子中，我们最终将会在16,23,29,13这个圈中无限循环，永远找不到因子 。
那么，如何探测环的出现呢？
一种方法是记录当前产生过的所有的数x1,x2,……xn,并检测是否存在xl=xn(l<n)。
在实际过程中，当 n增长到一定大小时，可能会造成的内存不够用的情况。

另一种方法是由Floyd发明的一个算法，我们举例来说明这个聪明而又有趣的算法。
假设我们在一个很长很长的圆形轨道上行走，我们如何知道我们已经走完了一圈呢？
当然，我们可以像第一种方法那样做，但是更聪明的方法是让A 和B 按照B 的速度是
A 的速度的两倍从同一起点开始往前走，当B 第一次敢上A 时(也就是我们常说的套圈)， 我们便知道，B 已经走了至少一圈了。

a := 2;
b := 2;
while ( b != a )
a = f(a); // a runs once
b = f(f(b)); // b runs twice as fast. p = GCD( | b - a | , N);
p = GCD( | b - a | , N);
if ( p > 1)
return “Found factor: p”;
end
return “Failed. :-(“

发现以前学长写的一篇博客挺好的:关于Pollard_rho的期望复杂度证明

n=p⋅q,那么f1(x)=(x2+c)(modn)的rho环比f2(x)=(x1+c)(modp)的环要大，在p环上相遇的步数是期望n14的，且进入p环时有很大可能没有进入n环，那么此时x≠y,设x=k1∗p+b1,y=k2∗p+b2,b1=b2，那么此时可以求出p这一因数。

#include<bits/stdc++.h>typedef unsigned int uint;typedef long long ll;using namespace std;inline uint unit(){    static uint state0=19491001;    state0^=(state0<<13);    state0^=(state0>>17);    state0^=(state0<<5);    return state0;}const int pr[9]={2,3,5,7,11,13,17,19,23};int T;long long n,mx;inline ll gcd(ll a,ll b){    return (b?gcd(b,a%b):a);}inline ll mul(ll a,ll b,ll mod){    return (a*b-(ll)((long double)a/mod*b)*mod+mod)%mod;}inline ll power(ll a,ll b,ll mod){    ll res=1;    for(;b;b>>=1,a=mul(a,a,mod))if(b&1)res=mul(res,a,mod);    return res;}inline bool mr(ll now){    if(now<=2)return mx=max(mx,now),true;    if(!(now&1))return false;    ll d=now-1,cnt=0;    while(!(d&1))cnt++,d>>=1;    for(int i=0;i<=8;i++){        if(pr[i]==now)return true;        if(!now%pr[i])return false;        ll t=power(pr[i],d,now);        if(t==1||t==now-1)continue;        for(int j=1;j<=cnt;j++){            ll tmp=power(t,2,now);            if(tmp==1&&t!=1&&t!=now-1)return false;            t=tmp;        }        if(t!=1)return false;    }    return true;}inline ll phi_rho(ll n,ll c){    ll x=unit()%n+1,y=x;    while(true){        x=(mul(x,x,n)+c)%n;        x=(mul(x,x,n)+c)%n;        y=(mul(y,y,n)+c)%n;        ll p=(y>x?y-x:x-y);        p=gcd(p,n);        if(p!=1)return p;    }}inline bool solve(ll n){    if(mr(n)){        mx=max(mx,n);        return true;    }    else{        ll d=phi_rho(n,unit()%n);        while(d==n)d=phi_rho(n,unit()%n);        solve(n/d);solve(d);        return false;    }}int main(){    scanf("%d",&T);    while(T--){        mx=0;scanf("%lld",&n);        if(solve(n))puts("Prime");        else printf("%lld\n",mx);    }}