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

题意

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

分析

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

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

long long乘long long取模黑科技：

代码

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
#define ll long long
using 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;
}