miller_rabin + pollard_rho模版

#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
",m*G,L/m);
    }
    return 0;
}