BZOJ_3667_Rabin-Miller算法_Mille

BZOJ_3667_Rabin-Miller算法_Mille_Rabin+Pollard rho

Description

Input

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

Output

第一行CAS(CAS<=350，代表测试数据的组数)
以下CAS行：每行一个数字，保证是在64位长整形范围内的正数。
对于每组测试数据：输出Prime，代表它是质数，或者输出它最大的质因子，代表它是和数

Sample Input

6
2
13
134
8897
1234567654321
1000000000000

Sample Output

Prime
Prime
67
41
4649
5

HINT

数据范围：

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

思路就是分解出一个因数，然后递归下去，直到这个数为质数，取最大的输出即可。

那么首先有判断质数的方法Miller_Rabin:

根据费马小定理,$a^{p-1}mod{p}=1$，$p$为质数。

以及$x^2mod{p}=1$，$p$为质数的解只有$1$和$p-1$。

这两个同时满足且$p$不是质数的概率就很低了，我们还可以多试几次。

具体地，把$n-1$分解成$a^r*2^s$的形式，然后用$a^r$一直平方，同时验证这两个。

然后有基于生日悖论的快速求约数的方法Pollard rho:

选取$0~n-1$的几个数，用连续的两个数$x,y$求$gcd(x-y,n)$，如果这不等于$1$则找到了一个$n$的约数。

其中$x$和$y$可以用某个函数生成，比如说$f(x)=x^2+c$，但这样可能进到一个环里找不到约数。

这时我们跳出来对$c$进行随机赋值就可以。

如何跳出来？可以让$x$一次跑一步$(x=f(x))$，让$y$一次跑两步$(y=f(f(y)))$，这样如果$x=y$就跳出即可。

(有一个bug：当$n=4$时会一直进去，$Abs(x-y)$永远不会等于$2$)

代码：

#include <cstdio>

#include <cstring>

#include <algorithm>

#include <cmath>

#include <cstdlib>

using namespace std;

typedef double f2;

typedef long long ll;

ll ch(ll a,ll b,ll mod) {ll d=(ll)floor(1.0*a/mod*b+0.5),re=a*b-d*mod; return re<0?re+mod:re;}

ll qp(ll x,ll y,ll mod) {ll re=1;for(;y;y>>=1ll,x=ch(x,x,mod)) if(y&1ll) re=ch(re,x,mod); return re;}

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

ll random(ll x,ll y) {

    return ((rand()*(1ll<<45))+(rand()<<30)+(rand()<<15)+(rand()))%(y-x+1)+x;

}

ll Abs(ll x) {return x>0?x:-x;}

ll maxn,a[]={2,3,5,7,11,13,17,19,23,29};

bool check(ll a,ll n,ll r,ll s) {

    ll x=qp(a,r,n),y=x,i;

    for(i=1;i<=s;i++,y=x) {x=ch(x,x,n); if(x==1&&y!=1&&y!=n-1) return 0;}

    return x==1;

}

bool MR(ll n) {

    if(n<=1) return 0; ll r=n-1,s=0,i;

    for(;!(r&1);r>>=1ll,s++);

    for(i=0;i<9;i++) {

        if(a[i]==n) return 1;

        if(!check(a[i],n,r,s)) return 0;

    }

    return 1;

}

ll f(ll x,ll c,ll mod) {return (ch(x,x,mod)+c)%mod;}

ll PR(ll n,ll c) {

    ll x=random(0,n-1),y=x,p;

    for(p=1;p==1;) {

        x=f(x,c,n); y=f(f(y,c,n),c,n); p=gcd(Abs(x-y),n);

    }

    return p;

}

void solve(ll n) {

    if(n<=1) return ; if(MR(n)) {maxn=max(maxn,n); return ;}

    if(n%2==0) {

        while(n%2==0) n/=2;

        solve(n); maxn=max(maxn,2ll); return ;

    }

    ll tmp=n; while(tmp==n) tmp=PR(n,random(0,n-2));

    solve(tmp); solve(n/tmp);

}

int main() {

    srand(19260817); rand();

    int T; scanf("%d",&T);

    while(T--) {

        ll n; scanf("%lld",&n);

        maxn=0; solve(n);

        if(n==maxn) printf("Prime\n");

        else printf("%lld\n",maxn);

    }

}

巴特西