数学#素数判定Miller_Rabin+大数因数分解Pollard_rho算法 POJ 1811&2429

素数判定Miller_Rabin算法详解：

http://blog.csdn.net/maxichu/article/details/45458569

大数因数分解Pollard_rho算法详解：

http://blog.csdn.net/maxichu/article/details/45459533

然后是参考了kuangbin的模板：

http://www.cnblogs.com/kuangbin/archive/2012/08/19/2646396.html

模板如下：

//快速乘 (a*b)%mod
//二进制竖式乘法：
//10101*1011=
//10101*1+10101*2^1*1+10101*2^2*0+10101*2^3*1
//位上是1的就加，是0的就不加
ll mult_mod(ll a,ll b,ll mod)
{
    a%=mod;
    b%=mod;
    ll res=0;
    while(b)
    {
        if(b&1)
        {
            res+=a;
            res%=mod;
        }
        a<<=1;
        if(a>=mod) a%=mod;
        b>>=1;
    }
    return res;
}

//快速幂 (x^n)%mod
ll pow_mod(ll x,ll n,ll mod)
{
    if(n==1) return x%mod;
    x%=mod;
    ll t=x;
    ll res=1;
    while(n)
    {
        if(n&1) res=mult_mod(res,t,mod);
        t=mult_mod(t,t,mod);//t平方
        n>>=1;//变成n/2
    }
    return res;
}

//若为合数返回true，不一定返回false
//a^(n-1)=1(mod n) -> a^((2^t)*x)=-1(mod n) == a^x=1(mod n)
bool test(ll a,ll n,ll x,ll t)//miller_rabin算法的核心
{
    ll res=pow_mod(a,x,n);//a^x mod n
    ll last=res;
    for(int i=1;i<=t;i++)
    {
        res=mult_mod(res,res,n);//res=res*res mod n
        if(res==1&&last!=1&&last!=n-1)
            return true;
        last=res;
    }
    if(res!=1) return true;
    return false;
}

//若为素数（或伪素数）返回true，合数返回false
bool miller_rabin(ll n)
{
    if(n<2) return false;
    if(n==2) return true;
    if((n&1)==0) return false; //偶数
    ll x=n-1,t=0;
    while((x&1)==0)//n-1=(2^t)*x;
    {
        x>>=1;
        t++;
    }
    for(int i=0;i<times;i++)//进行随机判定
    {
        ll a=rand()%(n-1)+1;//随机找0~n-1的整数
        if(test(a,n,x,t))//
            return false;
    }
    return true;
}

ll factor[100];//保存质因数分解结果
int tot;//记录质因数个数，下标从0开始

ll gcd(ll a,ll b)
{
    if(a==0) return 1;
    if(a<0) return gcd(-a,b);
    while(b)
    {
        ll c=a%b;
        a=b;
        b=c;
    }
    return a;
}

ll pollard_rho(ll x,ll c)
{
    ll i=1,k=2;
    ll x0=rand()%x;
    ll y=x0;
    while(1)
    {
        i++;
        x0=(mult_mod(x0,x0,x)+c)%x;
        ll d=gcd(y-x0,x);
        if(d!=1&&d!=x) return d;
        if(y==x0) return x;
        if(i==k)
        {
            y=x0;
            k+=k;
        }
    }
}

//对n进行素因子分解
void find_factor(ll n)
{
    if(miller_rabin(n))//若为素数
    {
        factor[tot++]=n;
        return ;
    }
    ll p=n;
    while(p>=n)
        p=pollard_rho(p,rand()%(n-1)+1);
    find_factor(p);
    find_factor(n/p);
}

---

POJ 1811  完完全全的模板题

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<ctime>
#include<algorithm>
#include<cstring>
using namespace std;
typedef long long ll;

const int times=20;//随机算法判定次数

//快速乘 (a*b)%mod
//二进制竖式乘法：
//10101*1011=
//10101*1+10101*2^1*1+10101*2^2*0+10101*2^3*1
//位上是1的就加，是0的就不加
ll mult_mod(ll a,ll b,ll mod)
{
    a%=mod;
    b%=mod;
    ll res=0;
    while(b)
    {
        if(b&1)
        {
            res+=a;
            res%=mod;
        }
        a<<=1;
        if(a>=mod) a%=mod;
        b>>=1;
    }
    return res;
}

//快速幂 (x^n)%mod
ll pow_mod(ll x,ll n,ll mod)
{
    if(n==1) return x%mod;
    x%=mod;
    ll t=x;
    ll res=1;
    while(n)
    {
        if(n&1) res=mult_mod(res,t,mod);
        t=mult_mod(t,t,mod);//t平方
        n>>=1;//变成n/2
    }
    return res;
}

//若为合数返回true，不一定返回false
//a^(n-1)=1(mod n) -> a^((2^t)*x)=-1(mod n) == a^x=1(mod n)
bool test(ll a,ll n,ll x,ll t)//miller_rabin算法的核心
{
    ll res=pow_mod(a,x,n);//a^x mod n
    ll last=res;
    for(int i=1;i<=t;i++)
    {
        res=mult_mod(res,res,n);//res=res*res mod n
        if(res==1&&last!=1&&last!=n-1)
            return true;
        last=res;
    }
    if(res!=1) return true;
    return false;
}

//若为素数（或伪素数）返回true，合数返回false
bool miller_rabin(ll n)
{
    if(n<2) return false;
    if(n==2) return true;
    if((n&1)==0) return false; //偶数
    ll x=n-1,t=0;
    while((x&1)==0)//n-1=(2^t)*x;
    {
        x>>=1;
        t++;
    }
    for(int i=0;i<times;i++)//进行随机判定
    {
        ll a=rand()%(n-1)+1;//随机找0~n-1的整数
        if(test(a,n,x,t))//
            return false;
    }
    return true;
}

ll factor[100];//保存质因数分解结果
int tot;//记录质因数个数，下标从0开始

ll gcd(ll a,ll b)
{
    if(a==0) return 1;
    if(a<0) return gcd(-a,b);
    while(b)
    {
        ll c=a%b;
        a=b;
        b=c;
    }
    return a;
}

ll pollard_rho(ll x,ll c)
{
    ll i=1,k=2;
    ll x0=rand()%x;
    ll y=x0;
    while(1)
    {
        i++;
        x0=(mult_mod(x0,x0,x)+c)%x;
        ll d=gcd(y-x0,x);
        if(d!=1&&d!=x) return d;
        if(y==x0) return x;
        if(i==k)
        {
            y=x0;
            k+=k;
        }
    }
}

//对n进行素因子分解
void find_factor(ll n)
{
    if(miller_rabin(n))//若为素数
    {
        factor[tot++]=n;
        return ;
    }
    ll p=n;
    while(p>=n)
        p=pollard_rho(p,rand()%(n-1)+1);
    find_factor(p);
    find_factor(n/p);
}

int main()
{
    srand(time(0)); //需要ctime文件，poj的g++无法实现
    int t;
    scanf("%d",&t);
    ll n;
    while(t--)
    {
        scanf("%I64d",&n);
        if(miller_rabin(n))
        {
            printf("Prime
");
            continue;
        }
        tot=0;
        find_factor(n);
        sort(factor,factor+tot);
        printf("%I64d
",factor[0]);
    }
}

---

POJ 2689 

这道题的函数是自己重新打了一遍的，因为理解不透彻，打算边打边想，结果还是没想清楚，尤其是后面的pollard_rho算法，而且还因为把long long型变量声明成了int型，导致一直TLE。。。

AC代码如下：

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<ctime>
#include<algorithm>
#include<cstring>
#include<cmath>
using namespace std;
typedef long long ll;

const int times=20;//随机算法判定次数
ll factor[1000],tot,ans,mina,k;

ll gcd(ll a,ll b)
{
    if(b==0)
        return a;
    return gcd(b,a%b);
}

ll mult_mod(ll a,ll b,ll mod)
{
    a%=mod; b%=mod;
    ll res=0;
    while(b)
    {
        if(b&1)
        {
            res+=a;
            res%=mod;
        }
        a<<=1;
        if(a>=mod) a%=mod;
        b>>=1;
    }
    return res;
}

ll pow_mod(ll x,ll n,ll mod)
{
    if(n==1) return x%mod;
    x%=mod;
    ll t=x,res=1;
    while(n)
    {
        if(n&1) res=mult_mod(res,t,mod);
        t=mult_mod(t,t,mod);
        n>>=1;
    }
    return res;
}

bool test(ll a,ll n)
{
    ll x=n-1,t=0,res,last;
    while((x&1)==0)
    {
        x>>=1;
        t++;
    }
    last=pow_mod(a,x,n);

    for(int i=0;i<t;i++)
    {
        res=pow_mod(last,2,n);
        if(res==1&&last!=1&&last!=n-1)
            return true;
        last=res;
    }
    if(res!=1) return true;
    return false;
}

bool millier_rabin(ll n)
{
    if(n==2) return true;
    if(n==1||(n&1)==0) return false;
    for(int i=0;i<times;i++)
    {
        ll a=rand()%(n-1)+1;
        if(test(a,n))
            return false;
    }
    return true;
}

ll pollard_rho(ll x,ll c)
{
    ll x0,y,i=1,k=2;
    x0=rand()%x;
    y=x0;
    while(1)
    {
        i++;
        x0=(mult_mod(x0,x0,x)+c)%x;
        ll d=gcd(y-x0,x);
        if(d>1&&d<x) return d;
        if(y==x0) break;
        if(i==k)
        {
            y=x0;
            k+=k;
        }
    }
    return x;
}

void find_fac(ll n,int c)
{
    if(n==1) return;
    if(millier_rabin(n))
    {
        factor[tot++]=n;
        return;
    }
    ll p=n;
    while(p>=n)
        p=pollard_rho(p,c--);
    find_fac(p,c);
    find_fac(n/p,c);
}

void dfs(ll i,ll x)
{
    if(i>=tot)
    {
        if(x>mina&&x<=k)
            mina=x;
        return;
    }
    dfs(i+1,x);
    dfs(i+1,x*factor[i]);
}

int main()
{
    ll n,m;
    while(~scanf_s("%I64d%I64d",&n,&m))
    {
        srand(time(0)); //需要ctime文件，poj的g++无法实现
        if(n==m)
        {
            printf("%I64d %I64d
",n,m);
            continue;
        }
        tot=0;
        ans=m/n;
        find_fac(ans,107);
        sort(factor,factor+tot);
        int i,j=0;
        for(i=1;i<tot;i++)
        {
            while(factor[i-1]==factor[i]&&i<tot)
                factor[j]*=factor[i++];
            if(i<tot)
                factor[++j]=factor[i];
        }
        tot=j+1; mina=1;
        k=(ll)sqrt(ans*1.0);
        dfs(0,1);
        printf("%I64d %I64d
",mina*n,ans/mina*n);
    }
    return 0;
}