【POJ2429】GCD & LCM Inverse-Pollard-rho分解+枚举

测试地址：GCD & LCM Inverse
题目大意：给定a和b的最大公因数（gcd）和最小公倍数（lcm），求a和b，其中a≤b，如果有多组解，输出使a+b最小的解。所有数≤263。
做法：这个题目应该使用Pollard-rho分解算法+枚举（或者有其他的方法？不管啦）。
我们来分析一下满足条件的a和b满足什么条件。显然a和b都含有gcd中含有的因子，又因为a×b=gcd×lcm，所以如果我们将a和b都除掉一个gcd，那么剩下来的因子就是lcmgcd的因子。所以我们可以对这个数进行质因数分解，然后将这些质因子分配到a和b上。
由于已经除掉了gcd，所以同一种质因子不能同时出现在两个数上，否则原来的gcd就不是gcd了，那么我们就可以把每一种质因子缩起来，然后枚举分配方案求出最优方案即可。可以证明，在题目所给的数据范围下，整数所含的质因子的种数不超过16（因为最小的16个质数乘起来已经超过了263），所以枚举方案的复杂度差不多是O(216)=O(65536)的。由于数据范围很大，质因数分解只能用Pollard-rho来做。需要注意的是，这一道题的数据范围已经超过了long long的存储范围，需要使用unsigned long long来存储才行。这就导致计算数中要是出现负数就会溢出，所以在Pollard-rho中的那一句gcd(y-x+n,n)就要写成gcd(n+y-x,n)（计算顺序问题）。这样我们就解决了这个问题。
以下是本人代码：

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <algorithm>
#define ll unsigned long long
using namespace std;
ll GCD,LCM,p[11]={2,3,5,7,11,13,17,19,23,29};
ll fac[110]={0},A,B,mins;
int tot;

ll mult(ll a,ll b,ll mod)
{
  ll s=a,sum=0;
  while(b)
  {
    if (b&1)
    {
      sum+=s;
      if (sum>=mod) sum-=mod;
    }
    b>>=1;
    s<<=1;
    if (s>=mod) s-=mod;
  }
  return sum;
}

ll power(ll a,ll b,ll mod)
{
  ll s=a,sum=1;
  while(b)
  {
    if (b&1) sum=mult(sum,s,mod);
    b>>=1;s=mult(s,s,mod);
  }
  return sum;
}

bool witness(ll n,ll a)
{
  ll p=power(a,n-1,n);
  if (p!=1) return 0;
  else
  {
    ll s=n-1;
    while(!(s%2)&&p==1)
    {
      s>>=1;
      p=power(a,s,n);
    }
    if (p==1||p==n-1) return 1;
    else return 0;
  }
}

bool miller_rabin(ll n)
{
  if (n<=29)
  {
    for(int i=0;i<=9;i++)
      if (p[i]==n) return 1;
    return 0;
  }
  for(int i=0;i<=9;i++)
    if (!witness(n,p[i])) return 0;
  return 1;
}

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

ll pollard_rho(ll n,ll c)
{
  ll x=rand()%n,y=x,d,i=1,k=2;
  while(1)
  {
    i++;
    x=(mult(x,x,n)+c)%n;
    d=gcd(n+y-x,n);
    if (d>1&&d<n) return d;
    if (y==x) return n;
    if (i==k) y=x,k<<=1;
  }
}

void find_factor(ll n)
{
  if (n==1) return;
  if (miller_rabin(n))
  {
    fac[++tot]=n;
    return;
  }
  ll p=n;
  while(p>=n) p=pollard_rho(p,rand()%(n-1)+1);
  find_factor(p);
  find_factor(n/p);
}

bool cmp(ll a,ll b) {return a<b;}

void work()
{
  int d=0;
  sort(fac+1,fac+tot+1,cmp);
  for(int i=1;i<=tot;i++)
  {
    if (fac[i]!=fac[i-1])
    {
      d++;
      if (i!=d) fac[d]=1;
    }
    if (i!=d) fac[d]*=fac[i];
  }
  tot=d;
  for(int i=0;i<(1<<tot);i++)
  {
    ll nowa=1,nowb=1,k=i;
    for(int j=1;j<=tot;j++)
    {
      if (k&1) nowb*=fac[j];
      else nowa*=fac[j];
      k>>=1;
    }
    if ((nowa+nowb)*GCD<mins) A=nowa*GCD,B=nowb*GCD,mins=(nowa+nowb)*GCD;
  }
}

int main()
{
  while(scanf("%lld%lld",&GCD,&LCM)!=EOF)
  {
    tot=0;A=B=1;mins=0;
    for(int i=0;i<=63;i++) mins+=1<<i;
    find_factor(LCM/GCD);
    work();
    if (A>B) {ll t=A;A=B;B=t;}
    printf("%lld %lld\n",A,B);
  }

  return 0;
}