POJ 2429（GCD&LCM Inverse）

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

思路：

其实思路不复杂，题目也短小精炼，但是这个涉及到 Miller_Rabin 和 Pollard_rho 算法

 

直接上代码:

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<iostream>
#include<queue>
#include<stack>
#include<cmath>
#include<set>
#include<algorithm>
#include<vector>

using namespace std;
#define clc(a,b) memset(a,b,sizeof(a))
#define inf  (LL)1<<61
#define LL long long

const int Times = 10;
const int N = 5500;

LL n, m, ct, cnt;
LL minn, mina, minb, ans;
LL fac[N], num[N];

LL gcd(LL a, LL b)
{
    return b? gcd(b, a % b) : a;
}

LL multi(LL a, LL b, LL m)
{
    LL ans = 0;
    a %= m;
    while(b)
    {
        if(b & 1)
        {
            ans = (ans + a) % m;
            b--;
        }
        b >>= 1;
        a = (a + a) % m;
    }
    return ans;
}

LL quick_mod(LL a, LL b, LL m)
{
    LL ans = 1;
    a %= m;
    while(b)
    {
        if(b & 1)
        {
            ans = multi(ans, a, m);
            b--;
        }
        b >>= 1;
        a = multi(a, a, m);
    }
    return ans;
}

bool Miller_Rabin(LL n)
{
    if(n == 2) return true;
    if(n < 2 || !(n & 1)) return false;
    LL m = n - 1;
    int k = 0;
    while((m & 1) == 0)
    {
        k++;
        m >>= 1;
    }
    for(int i=0; i<Times; i++)
    {
        LL a = rand() % (n - 1) + 1;
        LL x = quick_mod(a, m, n);
        LL y = 0;
        for(int j=0; j<k; j++)
        {
            y = multi(x, x, n);
            if(y == 1 && x != 1 && x != n - 1) return false;
            x = y;
        }
        if(y != 1) return false;
    }
    return true;
}

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

void find(LL n, int c)
{
    if(n == 1) return;
    if(Miller_Rabin(n))
    {
        fac[ct++] = n;
        return ;
    }
    LL p = n;
    LL k = c;
    while(p >= n) p = pollard_rho(p, c--);
    find(p, k);
    find(n / p, k);
}

void dfs(LL dept, LL tem=1)
{
    if(dept == cnt)
    {
        LL a = tem;
        LL b = ans / a;
        if(gcd(a, b) == 1)
        {
            a *= n;
            b *= n;
            if(a + b < minn)
            {
                minn = a + b;
                mina = a;
                minb = b;
            }
        }
        return ;
    }
    for(int i=0; i<=num[dept]; i++)
    {
        if(tem > minn) return;
        dfs(dept + 1, tem);
        tem *= fac[dept];
    }
}

int main()
{
    while(~scanf("%llu %llu", &n, &m))
    {
        if(n == m)
        {
            printf("%llu %llu
",n,m);
            continue;
        }
        minn = inf;
        ct = cnt = 0;
        ans = m / n;
        find(ans, 120);
        sort(fac, fac + ct);
        num[0] = 1;
        int k = 1;
        for(int i=1; i<ct; i++)
        {
            if(fac[i] == fac[i-1])
                ++num[k-1];
            else
            {
                num[k] = 1;
                fac[k++] = fac[i];
            }
        }
        cnt = k;
        dfs(0, 1);
        if(mina > minb) swap(mina, minb);
        printf("%llu %llu
",mina, minb);
    }
    return 0;
}

　　

补一个 扩展欧几里得算法

 分析：

ax1+by1=gcd(a,b)

        =gcd(b,a%b)=bx2+a%by2

        =bx2+(a-(a/b)*b)y2

        =bx2+ay2-(a/b)*by2

        =b(x2-(a/b)*y2)+ay2

　　

待定系数:

x1=y2

y1=x2-(a/b)*y2

// 说明：a/b是取整除法

代码：

int ext_gcd(int a,int b,int& x,int& y){

　　int t,ret;

　　if (!b){

　　　　x=1,y=0;

　　　　return a;

　　}

　　ret=ext_gcd(b,a%b,x,y);

　　t=x,x=y,y=t-a/b*y;

　　return ret;

}