#include <iostream>
#include<stdlib.h>
#include<cmath>
#include<ctime>
#include<algorithm>
#include<cstring>
#include<stdio.h>
using namespace std;
typedef long long ll;
const ll MAX_INDEEP = 10000;
const ll TABLE_SIZE = 131071;
ll sqrt_table[TABLE_SIZE] = {0};
const ll maxn=100;
ll factor[maxn];
ll cnt;
ll gcd(ll a,ll b){
    while(b){
        ll tem=b;
        b=a%b;
        a=tem;
    }
    return a;
}
ll try_ans(ll n){
    ll sqrt_n=sqrt(n*1.0);
    ll P0=sqrt_n,Q0=1,Q1=n-P0*P0;
    ll B1,P1,Q2;
    if(Q1==0)return P0;
    while(sqrt_table[Q1%TABLE_SIZE]!=Q1){
        B1=(sqrt_n+P0)/Q1;
        P1=B1*Q1-P0;
        Q2=Q0+B1*(P0-P1);
        P0=P1;
        Q0=Q1;
        Q1=Q2;
    }
    B1=(sqrt_n-P0)/sqrt(Q1*1.0);
    P0=B1*sqrt(Q1*1.0)+P0;
    Q0=sqrt(Q1*1.0);
    Q1=(n-P0*P0)/Q0;
    P1=P0;
    P0=0;
    ll step=MAX_INDEEP;
    while(P1!=P0&&step--){
        P0=P1;
        B1=(sqrt_n+P0)/Q1;
        P1=B1*Q1-P0;
        Q2=Q0+B1*(P0-P1);
        Q0=Q1;
        Q1=Q2;
    }
    return P1;
}
ll SQUFOF(ll n){
    ll t=0;
    for(ll k=1;t==0||t==1;k++){
        ll tem=try_ans(k*n);
        if(tem==0)continue;
        t=gcd(tem,n);
    }
    return t;
}
ll  mul_mod(ll  a,ll  b,ll  p){
    ll  res=a%p,ans=0;
    while(b){
        if(b&1)ans+=res;
        if(ans>p)ans-=p;
        res<<=1;
        if(res>p)res-=p;
        b>>=1;
    }
    return ans;
}
ll  mod_exp(ll  a,ll  n,ll  p){
    ll  d=a%p,ans=1;
    while(n>1){
        if(n&1)ans=mul_mod(ans,d,p);
        d=mul_mod(d,d,p);
        n>>=1;
    }
    return mul_mod(ans,d,p);
}
bool miller_rabin(ll n, ll times)
{
    if(n==2)return 1;
    if(n<2||!(n&1))return 0;
    ll a, u=n-1, x, y;
    int t=0;
    while(u%2==0){
        t++;
        u/=2;
    }
    srand(time(0));
    for(int i=0;i<times;i++)
	{
        a = rand() % (n-1) + 1;
        x = mod_exp(a, u, n);
        for(int j=0;j<t;j++)
		{
            y = mul_mod(x, x, n);
            if ( y == 1 && x != 1 && x != n-1 )
                return false; //must not
            x = y;
        }
        if( y!=1) return false;
    }
    return true;
}
void findfactor(ll  n){
    if(n==1)return ;
    if(miller_rabin(n,6)){
        factor[cnt++]=n;
        return ;
    }
    ll  p=SQUFOF(n);
    findfactor(p);
    findfactor(n/p);
}
ll Lcm[maxn];
ll num[maxn][2];
int len2,t2;
int binary(int x){
    int L=0,R=len2;
    while(L<R){
        int mid=(L+R)>>1;
        if(Lcm[mid]<x)L=mid+1;
        else if(Lcm[mid]>x)R=mid;
        else return mid;
    }
    return -1;
}
int A[maxn];
ll ansa,ansb;
ll maxs;
void dfs(int k){
    if(k==t2+1){
        ll aa=1,bb=1;
        for(int i=0;i<=t2;i++){
            aa*=(ll)pow(Lcm[i],num[i][A[i]]*1.0);
            bb*=(ll)pow(Lcm[i],num[i][!A[i]]*1.0);
        }
        if(aa+bb<maxs){
            maxs=aa+bb;
            ansa=aa;
            ansb=bb;
        }
    }
    else{
        for(A[k]=0;A[k]<2;A[k]++){
            dfs(k+1);
        }
    }
}
int main()
{
    //freopen("in.txt","r",stdin);
    for (ll i = 0; i < (1 << 16); i++)
              sqrt_table[i * i % TABLE_SIZE] = i * i;
    ll a, b;
    while(scanf("%I64d%I64d",&a,&b)!=EOF){
        len2=1;
        maxs=ll((ll)1<<63)-1;
        memset(A,0,sizeof(A));
        memset(num,0,sizeof(num));
        t2=0;
        if(miller_rabin(b,6))Lcm[0]=b,num[0][1]=1;
        else{
            cnt=0;
            findfactor(b);
            sort(factor,factor+cnt);
            Lcm[0]=factor[0];
            for(int i=1;i<cnt;i++){
                if(factor[i]!=factor[i-1])Lcm[len2++]=factor[i];
            }
            for(int i=1;i<cnt;i++){
                if(factor[i]==factor[i-1])num[t2][1]++;
                else t2++;
            }
            for(int i=0;i<=t2;i++)num[i][1]++;
        }
        if(miller_rabin(a,6)){
            int index=binary(a);
            num[index][0]=1;
        }
        else{
            cnt=0;
            findfactor(a);
            sort(factor,factor+cnt);
            for(int i=0;i<cnt;i++){
                int index=binary(factor[i]);
                num[index][0]++;
            }
        }
//        cout<<t2<<"n"<<endl;
//        for(int i=0;i<=t2;i++)cout<<Lcm[i]<<" ";
//        cout<<endl;
        if(t2==0)ansa=a,ansb=b;
        else dfs(0);
        if(ansa>ansb)swap(ansa,ansb);
        cout<<ansa<<" "<<ansb<<endl;
    }
    return 0;
}

• Re:想不到好的算法，只能枚举了 (4959) Manchery 2015-05-05 19:14:16