Pollard-Rho算法求大数质因子

/*
 * 大整数分解到现在都是世界级的难题,但却是一个重要的研究方向,大整数在公共密钥的研究上有着重要的作用
 * Pollard Rho算法的原理就是通过某种方法得到两个整数a和b.而待分解的大整数为n,计算p=gcd(abs(a-b),n),直到p不为1或者a,b出现循环为止.
 * 然后再判断是否为n,如果p==n||p==1,那么返回n时一个质数
 * 否则p就是n的一个因子
 * 那么我们又可以递归地计算Pollard(p)和Pollard(n/p).
 * 最后就可以推出n的所有质因子
 *
 * 具体操作中我们常常使用函数 x[i+1]=(x[i]*x[i]+c)%n 来逐步迭代计算a和b的值,通常c取1,即b=a*a+1,在下一次计算中,将b的值赋给a,再次使用上式来计算新的b的值,当a,b出现循环时即可退出判断.(初值自己确定)
 * 但是这样的话判断循环比较麻烦,这里给出Floyd(没错又是他)发明的一个聪明而又有趣的算法:
 *     假设我们在一个很长很长的圆形轨道上面行走,如何知道自己已经走了一圈了呢?
 *     可以让两个人A和B按照 vb = va<<1 从同一起点开始向前走,当B第一次赶上A时,我们就知道B已经走了两圈
 *     所以我们可以把x当作B,把y当作A,然后进行循环测试
 * 
 * 对于Pollard Rho算法,它可以在O(sqrt(p))的时间复杂度内找到n的一个小因子p,可见效率还是可以的
 * 但是对于一个因子很少或者因子值很大的大整数n来说,这个算法的复杂度依然不是很好
 */
//以下给出Pollard Rho和Miller-Rabin素数测试配合使用的整数分解算法
#include <algorithm>
#include <iostream>
#include <cstdlib>
#include <cstring>
#include <cstdio>
typedef long long ll;
const int counts = 10,N = 5001;

ll tot,cnt,fac[N],num[N];

ll gcd(ll a,ll b)
{return b?gcd(b,a%b):a;}

ll qpow(ll a,ll x,ll p)
{
    ll ret=1;
    for(;x;x>>=1,a=a*a%p)
        if(x&1)
            ret=ret*a%p;
    return ret;
}

ll multi(ll a,ll b,ll p)
{
    ll ans=0;
    a%=p;
    for(;b;b>>=1,a=(a<<1)%p)
        if(b&1)
            ans=(ans+a)%p;
    return ans;
}

bool Miller_Rabin(ll n)
{
    if(n==2)return true;
    if(n<2 || !(n&1))return false;
    ll m=n-1,a,x,y;int k=0;
    while(!(m&1))++k,m>>=1;
    for(int i=0;i<counts;++i)
    {
        a=rand()%(n-1)+1;
        x=qpow(a,m,n);
        y=0;
        for(int j=0;j<k;++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,k=2,x=rand()%(n-1)+1,y=x,d;
    while("fighting")
    {
        ++i;
        x=(multi(x,x,n)+c)%n;
        d=gcd((y-x+n),n);
        if(1<d && d<n)return d;
        if(y == x)return n;
        if(i == k)y=x,k<<=1;
    }
}

void find(ll n,int c)
{
    if(n == 1)return ;
    if(Miller_Rabin(n))
    {
        fac[tot++]=n;
        return ;
    }
    ll p=n,k=c;
    while(p>=n)p=Pollard_Rho(p,c--);
    find(p,k);
    find(n/p,k);
}

int main()
{
    ll n;
    while(std::cin>>n)
    {
        tot=0;
        find(n,120);
        std::sort(fac,fac+tot);
        num[0]=1;
        int k=1;
        for(int i=1;i<tot;++i)
        {
            if(fac[i] == fac[i-1])
                ++num[k-1];
            else
            {
                num[k]=1;
                fac[k++]=fac[i];
            }
        }
        cnt=k;
        for(int i=0;i<cnt;++i)
            std::cout<<fac[i]<<"^"<<num[i]<<" ";
        std::cout<<std::endl;
    }
    return 0;
}