POJ2429_GCD & LCM Inverse【Miller Rabin素数测试】【Pollar Rho整数分解】

GCD & LCM Inverse

Time Limit: 2000MSMemory Limit: 65536K

Total Submissions: 9756Accepted: 1819

Description


Given two positive integers a and b, we can easily calculate the greatest common divisor (GCD) and the least common multiple (LCM) of a and b. But what about the inverse? That is: given GCD and LCM, finding a and b.
Input


The input contains multiple test cases, each of which contains two positive integers, the GCD and the LCM. You can assume that these two numbers are both less than 2^63.
Output


For each test case, output a and b in ascending order. If there are multiple solutions, output the pair with smallest a + b.
Sample Input


3 60
Sample Output


12 15
Source

POJ Achilles

题目大意：给你两个数a和b的最大公约数和最小公倍数，求a和b

（其中在满足条件的情况下，使a+b尽量小）

思路：最大公约数和最小公倍数的规模为2^63，暴力果断不行。

已知a*b = L(最小公倍数)*G(最大公约数)；

设p = L/a，q = L/b，s = L/G；

即p、q为a和b除去最大公约数的部分，且两者互质；

GCD（p，q） = 1，LCM（p，q） = p * q = L*L/(a*b) = L*L/(L*G) = L/G = s。

LCM（p，q） = s；

由上可得我们可由s求出a和b。此题就是让我们把s分解成两个互质数相乘的形式。

用Pollar Rho整数分解和Miller Rabin素数测试结合起来，将s的所有质因子分解出来。

因为GCD（p，q） = 1，所有相同的质数不能同时分到p和q中，应将相同的质数分开放。

这里我们把所有相同的质数当做一个整体。将这些数枚举相乘，找到最接近s的平方根且不

大于s的平方根的组合即为p，则q = s/p

最终a = L/p，b = L/q

例如 G L 为 2 120

s = 60 = 2 * 2 * 3 * 5 = 4 * 3 * 5。枚举进行组合，找到最接近根号60并不超过根号60的

值为5，即p = 5，则q = 60/5 = 12。最终a = 24，b = 10。

参考博文：http://blog.sina.com.cn/s/blog_69c3f0410100uac0.html

#include<stdio.h>#include<stdlib.h>#include<time.h>#include<math.h>#include<algorithm>using namespace std;#define MAX_VAL (pow(2.0,60))//miller_rabbin素性测试//__int64 mod_mul(__int64 x,__int64 y,__int64 mo)//{//    __int64 t;//    x %= mo;//    for(t = 0; y; x = (x<<1)%mo,y>>=1)//        if(y & 1)//            t = (t+x) %mo;////    return t;//}__int64 mod_mul(__int64 x,__int64 y,__int64 mo){    __int64 t,T,a,b,c,d,e,f,g,h,v,ans;    T = (__int64)(sqrt(double(mo)+0.5));    t = T*T - mo;    a = x / T;    b = x % T;    c = y / T;    d = y % T;    e = a*c / T;    f = a*c % T;    v = ((a*d+b*c)%mo + e*t) % mo;    g = v / T;    h = v % T;    ans = (((f+g)*t%mo + b*d)% mo + h*T)%mo;    while(ans < 0)        ans += mo;    return ans;}__int64 mod_exp(__int64 num,__int64 t,__int64 mo){    __int64 ret = 1, temp = num % mo;    for(; t; t >>=1,temp=mod_mul(temp,temp,mo))        if(t & 1)            ret = mod_mul(ret,temp,mo);    return ret;}bool miller_rabbin(__int64 n){    if(n == 2)        return true;    if(n < 2 || !(n&1))        return false;    int t = 0;    __int64 a,x,y,u = n-1;    while((u & 1) == 0)    {        t++;        u >>= 1;    }    for(int i = 0; i < 50; i++)    {        a = rand() % (n-1)+1;        x = mod_exp(a,u,n);        for(int j = 0; j < t; j++)        {            y = mod_mul(x,x,n);            if(y == 1 && x != 1 && x != n-1)                return false;            x = y;        }        if(x != 1)            return false;    }    return true;}//PollarRho大整数因子分解__int64 minFactor;__int64 gcd(__int64 a,__int64 b){    if(b == 0)        return a;    return gcd(b, a % b);}__int64 PollarRho(__int64 n, int c){    int i = 1;    srand(time(NULL));    __int64 x = rand() % n;    __int64 y = x;    int k = 2;    while(true)    {        i++;        x = (mod_exp(x,2,n) + c) % n;        __int64 d = gcd(y-x,n);        if(1 < d && d < n)            return d;        if(y == x)            return n;        if(i == k)        {            y = x;            k *= 2;        }    }}__int64 ans[1100],cnt;void getSmallest(__int64 n, int c){    if(n == 1)        return;    if(miller_rabbin(n))    {        ans[cnt++] = n;        return;    }    __int64 val = n;    while(val == n)        val = PollarRho(n,c--);    getSmallest(val,c);    getSmallest(n/val,c);}__int64 a,b,sq;void choose(__int64 s,__int64 val){    if(s >= cnt)    {        if(val > a && val <= sq)            a = val;            return;    }    choose(s+1,val);    choose(s+1,val*ans[s]);}int main(){    int T;    __int64 G,L;    while(~scanf("%I64d%I64d",&G,&L))    {        if(L == G)        {            printf("%I64d %I64d\n",G,L);            continue;        }        L /= G;        cnt = 0;        getSmallest(L,200);        sort(ans, ans+cnt);        int j = 0;        for(int 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 = (__int64)sqrt(L+0.0);        choose(0,1);        printf("%I64d %I64d\n",a*G,L/a*G);    }    return 0;}