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

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;
}