bzoj3667 Rabin-Miller算法

Description

---

第一行：CAS,代表数据组数（不大于350），以下CAS行，每行一个数字，保证在64位长整形范围内，并且没有负数。你需要对于每个数字：第一，检验是否是质数，是质数就输出Prime
第二，如果不是质数，输出它最大的质因子是哪个。

保证cas<=350，保证所有数字均在64位长整形范围内。

Solution

---

第一次接触这种题目，感觉这么裸大概不会经常遇到

Miller_Rabin素数测试：

费马小定理：若p是质数且$a\not|p$则$a^{p−1}\pmod p\equiv1$

一个直观的想法是抽一些[2,p-1] 之间的数字求次幂看是否满足费马晓定理
然鹅有一种神奇的卡迈克尔数满足是合数但是能通过测试，因此需要换一种方法

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

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
这个函数最终会走出这样一个图，形如希腊字母ρ

调了很久不对，原来是方法不对-_-||

Code

---

#include <stdio.h>
#include <string.h>
#include <algorithm>
#define rep(i,st,ed) for (int i=st;i<=ed;++i)

typedef long long LL;
typedef long double LD;

int prime[9]={2,3,5,7,11,13,17,19,23};
LL ans;

LL read() {
    LL x=0,v=1; char ch=getchar();
    for (;ch<'0'||ch>'9';v=(ch=='-')?(-1):(v),ch=getchar());
    for (;ch<='9'&&ch>='0';x=x*10+ch-'0',ch=getchar());
    return x*v;
}

LL mul(LL a,LL b,LL MOD) {
    LL tmp=a*b-(LL)((LD)a/MOD*b+1e-8)*MOD;
    if (tmp<0) tmp+=MOD;
    return tmp;
}

LL ksm(LL x,LL dep,LL MOD) {
    if (dep==0) return 1;
    if (dep==1) return x;
    LL tmp=ksm(x,dep/2,MOD);
    tmp=mul(tmp,tmp,MOD);
    if (dep&1) return mul(tmp,x,MOD);
    return tmp;
}

LL gcd(LL x,LL y) {
    return !y?x:gcd(y,x%y);
}

bool test(LL n) {
    rep(i,0,8) if (prime[i]==n) return true;
    if (n<2) return false;
    if (n%2==0) return false;
    LL m=n-1,k=0;
    while (m%2==0) {
        k++; m/=2;
    }
    rep(i,0,8) {
        LL a=prime[i],w=ksm(a,m,n);
        if (w==1||w==n-1||a==n) continue;
        rep(j,1,k) {
            LL u=mul(w,w,n);
            if (w!=1&&w!=n-1&&u==1) return false;
            w=u;
        }
        if (w!=1) return false;
    }
    return true;
}

LL 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+k;
        }
    }
    return p;
}

void solve(LL n) {
    if (n==1) return ;
    if (test(n)) {
        ans=std:: max(ans,n);
        return ;
    }
    LL t=n;
    while (t==n) t=rho(n,rand()%(n-1)+1);
    solve(t); solve(n/t);
}

int main(void) {
    freopen("data.in","r",stdin);
    freopen("myp.out","w",stdout);
    int T; scanf("%d",&T);
    while (T--) {
        LL n=read();
        ans=0;
        solve(n);
        if (ans==n) puts("Prime");
        else printf("%lld\n", ans);
    }
    return 0;
}