miller_rabin + pollard

#include<stdio.h>

#include<stdlib.h>

#include<time.h>

#include<math.h>

#include<algorithm>

#define MAXN 100000

using namespace std;

typedef unsigned long long LL;

LL fac[MAXN],cnt,G,L,m,p;

LL min(LL a,LL b)

{

    return a<b?a:b;

}

LL gcd(LL a,LL b)

{

    return b==0?a:gcd(b,a%b);

}

LL mult_mod(LL a,LL b,LL mod)

{

    LL ans=0;

    while(b)

    {

        if(b&1)

            ans=(ans+a)%mod;

        a=(a<<1)%mod;

        b>>=1;

    }

    return ans;

}

LL pow_mod(LL a,LL b,LL mod)

{

    LL d=1;

    a%=mod;

    while(b)

    {

        if(b&1)

            d=mult_mod(d,a,mod);

        a=mult_mod(a,a,mod);

        b>>=1;

    }

    return d%mod;

}

bool witness(LL a,LL n)

{

    LL u=n-1,t=0;

    while((u&1)==0)

    {

        u>>=1;

        t++;

    }

    LL x,x0=pow_mod(a,u,n);

    for(LL i=1; i<=t; i++)

    {

        x=mult_mod(x0,x0,n);

        if(x==1&&x0!=1&&x0!=(n-1))

            return true;

        x0=x;

    }

    if(x!=1)

        return true;

    return false;

}

bool miller_rabin(LL n)

{

    if(n==2) return true;

    if(n<2||!(n&1)) return false;

    for(int j=1; j<=8; j++)

    {

        LL a=rand()%(n-1)+1;

        if(witness(a,n))

            return false;

    }

    return true;

}

LL pollard_rho(LL n,LL c)

{

    LL i=1,k=2,d,x=2,y=2;

    while(true)

    {

        i++;

        x=(mult_mod(x,x,n)+c)%n;

        d=gcd(y-x,n);

        if(d!=1&&d!=n)

            return d;

        if(x==y) return n;

        if(i==k)

        {

            y=x;

            k<<=1;

        }

    }

}

void find_fac(LL n,LL c)

{

    if(miller_rabin(n)||n<=1)

    {

        fac[cnt++]=n;

        return;

    }

    LL p=pollard_rho(n,c);

    while(p>=n) p=pollard_rho(p,c--);

    find_fac(p,c);

    find_fac(n/p,c);

}

void  dfs( LL step,LL num)

{

    if(step==cnt||num>p)

    {

        if(num<=p&&num>m)

            m=num;

        return;

    }

    dfs(step+1,num*fac[step]);

    dfs(step+1,num);

}

int main()

{

    srand(time(NULL));

    while(scanf("%I64u%I64u",&G,&L)!=EOF)

    {

        cnt=0;

        find_fac(L/G,181);

        sort(fac,fac+cnt);

        LL i=0,t=0;

        while(i<cnt)

        {

            LL ans=1,j=i;

            while(j<cnt&&fac[i]==fac[j])

            {

                ans*=fac[i];

                j++;

            }

            fac[t++]=ans;

            i=j;

        }

        cnt=t,m=1,p=sqrt(0.0+(L/G));

        dfs(0,1);

        printf("%I64u %I64u\n",m*G,L/m);

    }

    return 0;

}
巴特西

miller_rabin + pollard_rho模版

最新文章

热门文章
