GCD & LCM Inverse

题目：http://poj.org/problem?id=2429

题意：

给你两个数的gcd和lcm，[1, 2^63)。求a,b。使得a+b最小。

思路：

lcm = a * b / gcd 将lcm/gcd之后进行大数分解，形成a^x1 * b^x2 * c^x3…… 的形式，其中a,b,c为互不相同的质数。然后暴力枚举即可。

代码：

#include<map>
#include<set>
#include<queue>
#include<stack>
#include<cmath>
#include<cstdio>
#include<vector>
#include<string>
#include<fstream>
#include<cstring>
#include<ctype.h>
#include<iostream>
#include<algorithm>
#define INF (1<<30)
#define PI acos(-1.0)
#define mem(a, b) memset(a, b, sizeof(a))
#define rep(i, n) for (int i = 0; i < n; i++)
#define debug puts("===============")
typedef long long ll;
using namespace std;
const int S = 20; 

//计算 (a*b)%c.   a,b都是long long的数，直接相乘可能溢出的
//  a,b,c <2^63
ll mult_mod(ll a, ll b, ll c) {
    a %= c;
    b %= c;
    ll ret = 0;
    while(b) {
        if(b & 1) {
            ret += a;
            ret %= c;
        }
        a <<= 1;
        if(a >= c) a %= c;
        b >>= 1;
    }
    return ret;
}

//计算  x^n %c
ll pow_mod(ll x, ll n, ll mod) { //x^n%c
    if(n == 1)return x % mod;
    x %= mod;
    ll tmp = x;
    ll ret = 1;
    while(n) {
        if(n & 1) ret = mult_mod(ret, tmp, mod);
        tmp = mult_mod(tmp, tmp, mod);
        n >>= 1;
    }
    return ret;
}

//以a为基,n-1=x*2^t      a^(n-1)=1(mod n)  验证n是不是合数
//一定是合数返回true,不一定返回false
bool check(ll a, ll n, ll x, ll t) {
    ll ret = pow_mod(a, x, n);
    ll last = ret;
    for(int i = 1; i <= t; i++) {
        ret = mult_mod(ret, ret, n);
        if(ret == 1 && last != 1 && last != n - 1) return true; //合数
        last = ret;
    }
    if(ret != 1) return true;
    return false;
}

// Miller_Rabin()算法素数判定
//是素数返回true.(可能是伪素数，但概率极小)
//合数返回false;
bool Miller_Rabbin(ll n) {
    if(n < 2)return false;
    if(n == 2)return true;
    if((n & 1) == 0) return false;
    ll x = n - 1;
    ll t = 0;
    while((x & 1) == 0) {
        x >>= 1;
        t++;
    }
    for(int i = 0; i < S; i++) {
        ll a = rand() % (n - 1) + 1;
        if(check(a, n, x, t))
            return false;//合数
    }
    return true;
}

//************************************************
//pollard_rho 算法进行质因数分解
//************************************************
ll factor[100];//质因数分解结果（刚返回时是无序的）
int tot;//质因数的个数。数组小标从0开始

ll gcd(ll a, ll b) {
    if(a == 0) return 1;
    if(a < 0) return gcd(-a, b);
    while(b) {
        ll t = a % b;
        a = b;
        b = t;
    }
    return a;
}

ll Pollard_rho(ll x, ll c) {
    ll i = 1, k = 2;
    ll x0 = rand() % x;
    ll y = x0;
    while(1) {
        i++;
        x0 = (mult_mod(x0, x0, x) + c) % x;
        ll d = gcd(y - x0, x);
        if(d != 1 && d != x) return d;
        if(y == x0) return x;
        if(i == k) {
            y = x0;
            k += k;
        }
    }
}
//对n进行素因子分解
void findfac(ll n) {
    if(Miller_Rabbin(n)) { //素数
        factor[tot++] = n;
        return;
    }
    ll p = n;
    while(p >= n) p = Pollard_rho(p, rand() % (n - 1) + 1);
    findfac(p);
    findfac(n / p);
}

int vis[111] = {0};
ll g, lcm;
ll ans = (1LL << 61) , s;
int main () {
    while(~scanf("%lld%lld", &g, &lcm)) {
        lcm /= g;
        if (lcm == 1) {
            printf("%lld %lld\n", g, g);
            continue;
        }
        tot = 0;
        findfac(lcm);
        sort(factor, factor + tot);
        int cnt = 0;
        ll a[111];
        ll res = 1, now = factor[0];
        for (int i = 0; i < tot; i++) {
            if (factor[i] == now) {
                res = res * now;
            } else {
                a[cnt++] = res;
                now = factor[i];
                res = now;
            }
        }
        a[cnt++] = res;
        int t = (1 << cnt);
        ll A = 1, B = lcm;
        for (int i = 0; i < t; i++) {
            ll TA = 1, TB = 1;
            for (int j = 0; j < cnt; j++)
                if (i & (1 << j)) TA *= a[j];
                else TB *= a[j];
            if (TA + TB < A + B) {
                A = TA;
                B = TB;
            }
        }
        if (A > B) swap(A, B);
        printf("%lld %lld\n", A * g, B * g);
    }
    return 0;
}

POJ 2429 GCD & LCM Inverse （大数分解）,布布扣,bubuko.com