Pollard-rho算法：模板

最新推荐文章于 2022-07-04 13:30:28 发布

weixin_34297300

最新推荐文章于 2022-07-04 13:30:28 发布

阅读量320

点赞数

文章标签： python

原文链接：http://www.cnblogs.com/autsky-jadek/p/7750708.html

版权

#include<algorithm>
#include<cstdio>
#include<cstdlib>
#define N 5500
using namespace std;
typedef long long ll;
ll ct,cnt;  
ll fac[N],num[N];
const int BASE[]={2,3,5,7,11,13,17,19,23};
ll Quick_Mul(ll a,ll p,ll MOD)
{
    if(!p){
        return 0;
    }
    ll ans=Quick_Mul(a,p>>1,MOD);
    ans=(ans+ans)%MOD;
    if((p&1ll)==1ll){
        ans=ans+a%MOD%MOD;
    }
    return ans;
}
ll Quick_Pow(ll a,ll p,ll MOD)
{
    if(!p){
        return 1;
    }
    ll ans=Quick_Pow(a,p>>1,MOD);
    ans=Quick_Mul(ans,ans,MOD);
    if((p&1ll)==1ll){
        ans=a%MOD*ans%MOD;
    }
    return ans;
}
bool test(ll n,ll a,ll d){
    if(n==2){
        return 1;
    }
    if(n==a){
        return 0;
    }
    if(!(n&1)){
        return 0;
    }
    while(!(d&1ll)){
        d>>=1;
    }
    ll t=Quick_Pow(a,d,n);
    if(t==1){
        return 1;
    }
    while(d!=n-1ll && t!=n-1ll && t!=1ll){
        t=Quick_Mul(t,t,n);
        d<<=1;
    }
    return t==n-1ll;
}
bool Miller_Rabin(ll n){
    if(n==1 || n==3825123056546413051ll){
        return 0;
    }
    for(int i=0;i<9;++i){
        if(n==(ll)BASE[i]){
            return 1;
        }
        if(!test(n,(ll)BASE[i],n-1ll)){
            return 0;
        }
    }
    return 1;
}
ll pollard_rho(ll n,ll c){
    ll i=1,k=2;
    ll x=rand()%(n-1)+1;
    ll y=x;
    while(1){
        i++;
        x=(Quick_Mul(x,x,n)+c)%n;
        ll d=__gcd((y-x+n)%n,n);
        if(1ll<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[ct++]=n;
        return;
    }
    ll p=n;
    ll k=c;
    while(p>=n){
    	p=pollard_rho(p,c--);
    }
    find(p,k);
    find(n/p,k);
}
ll n;
int main(){
	srand(233);
    while(scanf("%I64d",&n)!=EOF){
        ct=0;
        find(n,120);
        sort(fac,fac+ct);
        num[0]=1;
        int k=1;
        for(int i=1;i<ct;++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){
        	printf("%I64d^%I64d\n",fac[i],num[i]);
        }
        puts("");
    }  
    return 0;  
}

转载于:https://www.cnblogs.com/autsky-jadek/p/7750708.html