POJ 2429 GCD & LCM Inverse （整数分解，由gcd+lcm求a,b)

题意：给你两个数a,b的最大公约数和最小公倍数，求a,b。（有多组a,b的情况下取a+b最小的）

题解：令 c = a * b / gcd(a,b)，对 c 因式分解

假如c = p1^k1 * p2^k2 * p3^k3

令 d1 = p1^k1, d2 = p2^k2, d3 = p3^k3

然后从d1,d2,d3中选取某几项，使得它们的积 s 最接近sqrt(c)，且<=sqrt(c)

那么 a = s * gcd(a,b)， b = c / s * gcd(a,b)

#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<ctime>
#include<algorithm>
using namespace std;

#define lint __int64
lint ans[1000], cnt;
lint G, L, a, b, sq;

lint gcd(lint a,lint b)
{
if ( b == 0 )
return a;
return gcd ( b, a % b );
}

lint mod_mult ( lint a, lint b, lint n )
{
lint ret = 0;
a = a % n;
while ( b >= 1 )
{
if ( b & 1 )
{
ret += a;
if ( ret >= n ) ret -= n;
}
a = a << 1;
if ( a >= n ) a -= n;
b = b >> 1;
}
return ret;
}

lint mod_exp ( lint a, lint b, lint n )
{
lint ret = 1;
a = a % n;
while ( b >= 1 )
{
if ( b & 1 )
ret =  mod_mult(ret,a,n);
a = mod_mult(a,a,n);
b = b >> 1;
}
return ret;
}

bool witness ( lint a, lint n )
{
int i, t = 0;
lint m = n - 1, x, y;
while ( m % 2 == 0 ) { m >>= 1; t++; }
x = mod_exp (a, m, n);

for ( i = 1; i <= t; i++ )
{
y = mod_exp ( x, 2, n );
if( y==1 && x!=1 && x!=n-1 )
return true;
x = y;
}
if ( y != 1 ) return true;
return false;
}

bool miller_rabin ( lint n, int times = 10 )
{
if ( n == 2 ) return true;
if ( n == 1 || n % 2 == 0 ) return false;

srand ( time(NULL) );
for ( int i = 1; i <= times; i++ )
{
lint a = rand() % (n-1) + 1;
if ( witness(a,n) ) return false;
}
return true;
}

lint rho ( lint n, int c )
{
lint i, k, x, y, d;
srand ( time(NULL) );
i = 1;  k = 2;
y = x = rand() % n;
while ( true )
{
i++;
x = (mod_mult(x,x,n)+c) % n;
d = gcd ( y - x, n );
if ( d > 1 && d < n ) return d;
if ( y == x ) break;
if ( i == k ) { y = x; k *= 2; }
}
return n;
}

void pollard ( lint n, int c )
{
if ( n == 1 )  return;
if ( miller_rabin(n) ) { ans[cnt++] = n; return; }
lint m = n;
while ( m >= n )
m = rho ( m, c-- );
pollard ( m, c );
pollard ( n / m, c );
}

void choose ( lint s, lint val )
{
if ( s >= cnt )
{
if ( val > a && val <= sq )
a = val;
return;
}
choose ( s + 1, val );
choose ( s + 1, val * ans[s] );
}

int main()
{
while ( scanf("%I64d%I64d",&G,&L) != EOF )
{
if ( L == G )
{
printf("%I64d %I64d\n",L,G);
continue;
}
L /= G;
cnt = 0;
pollard ( L, 107 );
sort ( ans, ans + cnt );
int i, j = 0;
for ( i = 1; i < cnt; i++ )
{
while ( ans[i-1] == ans[i] && i < cnt )
ans[j] *= ans[i++];
if ( i < cnt ) ans[++j] = ans[i];
}

cnt = j + 1; a = 1;
sq = (lint)sqrt(L+0.0);
choose ( 0, 1 );
printf("%I64d %I64d\n", a*G, L/a*G);
}
return 0;
}