Miller-Rabin与Pollard-Rho备忘

Miller-Rabin素性测试算法：

　　根据费马小定理当p为素数时成立，所以如果存在一个a使x不满足此定理，则x必然不为素数。

　　但这是充分条件而不是必要条件，所以对于每个a，可能存在满足定理的x，这时就要选取多个a同时检测，这种验证素性的方法即为Miller-Rabin算法。

　　当a取2,3,5,7时，可以直接检测1e13内的所有整数。

　　但是存在非素数通过检测，这时需要进行二次检测。

　　可以证明当p为奇素数时，$x^2equiv 1(mod p)$的解有且仅有两个：1和p-1。根据这个定理可以再次检测出一些非素数。

　　但是仍然存在一些数无法辨别，这些数被称为“强伪素数”，多选取一些数为底数检测即可（一般在[1,p)内选3个左右做二次检测就可以保证一定的正确率了）。

　　概括一下算法的具体流程：

　　1.先特判掉1,2和2的倍数。

　　2.选取3个以上[1,p)的数a。

　　3.对每个a先将p-1中2的因子去除，再逐个加上并实时检测是否出现不合法情况。

　　4.同时要注意判定$a^{p-1}equiv 1(mod p)$

Pollard-Rho质因数分解算法：

　　一种复杂度证明较为复杂的算法，主要思想是先求出n的一个因子p，然后对于p和n/p分别递归下去，如果发现p为素数则停止递归。（这里判断素数需要用到Miller-Rabin）。

　　主要思想是，让a和b同时在$f(x)=x^2+c$的轨迹上走（c需要变化），每2的次幂步进行一次a=b。每次判定若gcd(|a-b|,n)在(1,n)中则返回，当a==b时退出。

　　这样的算法，如果将a和b直到a=b的轨迹画出来，会是一条链加一个环，a每次在上面走1步，b走两步，形如$ ho$。

#include<cstdio>
#include<vector>
#include<algorithm>
#define rep(i,l,r) for (int i=(l); i<=(r); i++)
typedef long long ll;
using namespace std;

vector<ll>ls;
const int A[]={0,2,3,5,7,61,24251,131};
ll T,e,n,c,d,y,p,r,w,mx;
ll Rand(ll l,ll r){ return (((ll)rand()<<31)+rand())%(r-l+1)+l; }
ll gcd(ll a,ll b){ return (b) ? gcd(b,a%b) : a; }

ll ksc(ll a,ll b,ll p)
{
    ll t=a*b-(ll)((long double)a*b/p+0.5)*p;
    return (t<0)?t+p:t; 
}

ll ksm(ll a,ll b,ll mod){
    ll res=1;
    for (; b; a=ksc(a,a,mod),b>>=1)
        if (b & 1) res=ksc(res,a,mod);
    return res;
}

ll find(ll n,int c){
    ll i=1,k=2,x=Rand(0,n-1),y=x,d;
    while (1){
        i++; x=(ksc(x,x,n)+c)%n; d=gcd(abs(y-x),n);
        if (d>1 && d<n) return d;
        if (y==x) return -1;
        if (i==k) y=x,k<<=1;
    }
}

ll Rho(ll n,int c){ ll p=-1; while (p==-1) p=find(n,c--); return p; }

bool chk(ll a,ll n){
    ll m=n-1,x,y; int k=0;
    while (!(m&1)) m>>=1,k++;
    x=ksm(a,m,n);
    rep(i,1,k){
        y=ksm(x,2,n);
        if (y==1 && x!=1 && x!=n-1) return 1;
        x=y;
    }
    return y!=1;
}

bool Miller(ll n){
    if (n==2) return 1;
    if (!(n&1)) return 0;
    rep(i,1,3) if (chk(Rand(1,n-1),n)) return 0;
    return 1;
}

void Fac(ll x,int c){
    if (x==1) return;
    if (Miller(x)) { ls.push_back(x); return; }
    ll p=Rho(x,c); Fac(p,c);
    while (x%p==0) x/=p;
    Fac(x,c);
}

int main(){
    for (scanf("%lld",&T); T--; ){
        scanf("%lld",&n); ls.clear();
        if (Miller(n)) { puts("Prime"); continue; }
        Fac(n,19260817); mx=0;
        for (vector<ll>::iterator it=ls.begin(); it!=ls.end(); it++) mx=max(mx,*it);
        printf("%lld
",mx);
    }
    return 0;
}