BZOJ3667: Rabin-Miller算法 (Miller-Rabin&&pol_rho&&特技快速乘学习笔记)

来源:互联网 发布:歌曲伴奏降调软件 编辑:程序博客网 时间:2024/06/05 09:19

学习了两个新算法,判素数和分解大合数。

Miller-Rabin

Miller-Rabin(以下简称MR)算的时候大概就是这样子算。


那么怎么分解质因数呢?

这两篇博文应该已经说得比较清楚了。
http://www.cnblogs.com/thythy/p/5493624.html http://blog.csdn.net/maxichu/article/details/45459533

rho每次只找一个因数,要找n的因数,每次随机一个x出来,首先令y=x,然后令y不断做如下变化y=y*y+c(c是自己定的一个随机数),然后用abs(x-y)去跟n求gcd,如果不为1,则找到一个,有一个步长扩二倍的优化(详见代码),不太明白为啥,估计是为了防止死循环。

找到一个因数之后递归求解(配合MR),找到所有质因数,期望复杂度O(n1/4)(口胡)



另外一个问题就是计算MOD为64位整数的a*b%MOD,代码给出了一种奇技淫巧,在上面链接的第一篇blog里面有讲解

//QWsin#include<cmath>#include<cstdio>#include<cstring>#include<cstdlib>#include<iostream>#include<algorithm>using namespace std;typedef long long ll;const int a[]={2,3,5,7,11,13,17,19,23,29};ll maxs=0;ll mul(ll a,ll b,ll p){    ll d=(long double)a/p*b+0.5;    ll ret=a*b-d*p;    if(ret<0) ret+=p;    return ret;}ll Pow(ll x,ll p,ll MOD){    ll ret=1;    for(;p;p>>=1,x=mul(x,x,MOD)) if(p&1) ret=mul(ret,x,MOD);    return ret; }inline int MR(ll n){    if(n==2) return 1;    ll t=n-1;int cnt=0;    while((t&1)==0) t>>=1,++cnt;    for(int i=0;i<10;++i)    {        if(a[i]==n) return 1;        ll fac=Pow(a[i],t,n),nxt;        for(int i=1;i<=cnt;++i)        {            nxt=mul(fac,fac,n);            if(nxt==1&&fac!=1&&fac!=n-1) return 0;            fac=nxt;        }        if(fac!=1) return 0;    }    return 1;}ll gcd(ll a,ll b){return b?gcd(b,a%b):a;}ll pol_rho(ll n,ll c){    ll k=2,x=rand()%n,y=x,p=1;    for(ll i=1;p==1;++i){        y=(mul(y,y,n)+c)%n;        p=x>y?x-y:y-x;        p=gcd(p,n);        if(i==k) x=y,k<<=1;    }    return p;}inline void solve(ll n){//  printf("%lld\n",n);    if(MR(n)){maxs=max(maxs,n);return ;}    ll t=n;    while(t==n||t==1)         t=pol_rho(n,rand()%(n-1));    solve(n/t);    solve(t);}int main(){    srand(20010827);    int T;cin>>T;    while(T--) {        ll n;cin>>n;        maxs=0;        solve(n);        if(maxs==n) puts("Prime");        else printf("%lld\n",maxs);    }    return 0;}
0 0
原创粉丝点击