bzoj 3667: Rabin-Miller算法 素数测试

题意

你需要对于每个数字：第一，检验是否是质数，是质数就输出Prime
第二，如果不是质数，输出它最大的质因子是哪个。 保证是在64位长整形范围内的正数。

分析

Miller_Rabin素数测试：
费马小定理：若p是质数且a!|p则ap−1modp=1
二次探测定理：若a2modp=1且a!=1,a!=p−1则p必为合数
于是我们得出了专业Miller_Rabin算法：
a取 2, 3, 5, 7, 11, 13, 17, 19, 23
设p−1=2sd，我们先计算ad，然后给它平方s次。
若ad=1，则该a无法验证（视为通过）；
若平方过程中=1，则不通过；
若平方过程中=-1，则该a无法验证（视为通过）；
把最后ap−1的判断也视作上面的过程，则程序最后return 不通过。

Pollard_rho大数分解质因数：
设随机函数f(x)=x2+c，我们一开始随机出c、x1和x2，并判断gcd(|x1−x2|,n)是否大于1，若是，则返回该gcd；否则x1=f(x1)，x2=f(f(x2))。若此时x1=x2，则重新随机c、x1和x2。

long long乘long long取模黑科技：

代码

#include<iostream>#include<cstdio>#include<cstdlib>#include<cstring>#include<algorithm>#define ll long longusing namespace std;ll prime[9]={2,3,5,7,11,13,17,19,23},mx;ll read(){    ll x=0,f=1;char ch=getchar();    while (ch<'0'||ch>'9'){if(ch=='0')f=-1;ch=getchar();}    while (ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}    return x*f;}ll mul(ll a,ll b,ll MOD){    ll tmp=a*b-(ll)((long double)a/MOD*b+1e-8)*MOD;    if (tmp<0) tmp+=MOD;    return tmp;}ll ksm(ll x,ll y,ll MOD){    if (!y) return 1;    if (y==1) return x%MOD;    ll w=ksm(x,y/2,MOD);    w=mul(w,w,MOD);    if (y%2==1) w=mul(w,x,MOD);    return w;}ll gcd(ll x,ll y){    if (!y) return x;    if (x%y==0) return y;    else return gcd(y,x%y);}bool MR(ll n){    if (n<2) return 0;    if (n==2) return 1;    if (n%2==0) return 0;    ll m=n-1;int k=0;    while (m%2==0)    {        k++;m/=2;    }    for (int i=0;i<9;i++)    {        ll a=prime[i],w=ksm(a,m,n);        if (w==1||w==n-1||a==n) continue;        for (int j=1;j<=k;j++)        {            ll u=mul(w,w,n);            if (w!=1&&w!=n-1&&u==1) return 0;            w=u;        }        if (w!=1) return 0;    }    return 1;}int rho(ll n,ll c){    ll x1=(ll)(rand()+1)%n,x2=x1,p=1,k=2;    for (ll i=1;p==1;i++)    {        x1=(mul(x1,x1,n)+c)%n;        p=x1>x2?x1-x2:x2-x1;        p=gcd(p,n);        if (i==k)        {            x2=x1;k+=k;        }    }    return p;}void solve(ll n){    if (n==1) return;    if (MR(n))    {        mx=max(mx,n);        return;    }    ll t=n;    while (t==n) t=rho(n,rand()%(n-1)+1);    solve(t);solve(n/t);}int main(){    int t;    scanf("%d",&t);    while (t--)    {        ll x=read();        mx=0;        solve(x);        if (mx==x) printf("Prime\n");        else printf("%lld\n",mx);    }    return 0;}