Pollard Rho 求大数因子

传送门:证明

大佬水平   

多种方式详解 //强烈建议

#include<iostream>
#include<cstdlib>
#include<algorithm>
using namespace std;
#define ll long long
#define Times 10
ll fac[1000001];
ll num[1000001];
int cnt;
ll gcd(ll a,ll b)//求最大公约数
{
    return b?gcd(b,a%b):a;
}
ll multi(ll a,ll b,ll mod)//快速乘法法求模
{
    ll ret=0;
    while(b)
    {
        if(b&1)
            ret=(ret+a)%mod,b--;
        a=(a+a)%mod;
        b>>=1;
    }
    return ret;
}
ll power(ll a,ll b,ll mod)
{
    ll ret=1;
    a%=mod;
    while(b)
    {
        if(b&1)
        {
            ret=multi(ret,a,mod);
            b--;
        }
        a=multi(a,a,mod);
        b>>=1;
    }
    return ret;
}
bool miller_rabin(ll n)
{
    if(n==2)
        return true;
    if(n<2||!(n&1))
        return false;
    ll m=n-1;
    ll cnt=0;
    while(!(m&1))
    {
        ++cnt;
        m>>=1;
    }
    for(int i=0; i<Times; ++i)//Times==10
    {
        ll a=rand()%(n-2)+2;
        ll x=power(a,m,n);
        ll y=0;
        for(int j=0; j<cnt; ++j)
        {
            y=multi(x,x,n);
            if(y==1&&x!=1&&x!=n-1)
                return false;
            x=y;
        }
        if(y!=1)
            return false;
    }
    return true;
}
ll pollard_rho(ll n,ll c)
{
    ll i=1;
    ll k=2;
    ll x=rand()%(n-2)+2;
    ll y=x;
    while(true)
    {
        ++i;
        x=(multi(x,x,n)+c)%n;//c取 1 运算
        ll d=gcd((y-x+n)%n,n);
        if(1<d&&d<n)
            return d;
        if(x==y)
            return n;
        if(i==k)
        {
            y=x;
            k<<=1;
        }
    }
}
void find(ll n)
{
    if(n==1)
        return ;
    if(miller_rabin(n))
    {
        fac[cnt++]=n;
        return ;
    }
    ll p=n;
    if(p>=n)
        p=pollard_rho(p,rand()%(n-2)+2);
    find(p);
    find(n/p);
}
int main()
{
    ll n;
    while(cin>>n)
    {
        cnt=0;
        find(n);
        sort(fac,fac+cnt);
        num[0]=1;
        int k=1;
        for(int i=1;i<cnt; ++i)
        {
            if(fac[i]==fac[i-1])
                ++num[k-1];//除去重复的数
            else
            {
                num[k]=1;
                fac[k++]=fac[i];
            }
        }
        cout<<n<<"=";
        for(int i=0;i<k;++i)
            cout<<fac[i]<<"^"<<num[i]<<" ";
        cout<<endl;
    }
    return 0;
}

猜你喜欢

转载自blog.csdn.net/sdau_fangshifeng/article/details/81352031