POJ 2429 - GCD & LCM Inverse(数论)

题目:

http://poj.org/problem?id=2429

题意:

已知两个数的最大公约数和最小公倍数,求出两个数(和最小).

思路:

因为题目的数据特别大,所以要用Rabin-Miller强伪素数测试和Pollard因数分解算法.

基本思路是:

(a / gcd) * (b / gcd) = lcm / gcd ，所以需要分解lcm / gcd 。将其分解为互质的两个数，如果这两个数之和最小，那么乘上gcd就是所求的答案。

论文链接http://wenku.baidu.com/view/fbbed5a5f524ccbff12184af.html 抄了代码=.=

CODE:

#include <iostream>
#include <vector>
#include <map>
#include <cstdlib>
using namespace std;

typedef long long ll;

// return (a * b) % m
ll mod_mult(ll a, ll b, ll m)
{
    ll res = 0;
    ll exp = a % m;
    while (b)
    {
        if (b & 1)
        {
            res += exp;
            if (res > m) res -= m;
        }
        exp <<= 1;
        if (exp > m) exp -= m;
        b >>= 1;
    }
    return res;
}

// return (a ^ b) % m
ll mod_exp(ll a, ll b, ll m) {
    ll res = 1;
    ll exp = a % m;
    while (b)
    {
        if (b & 1) res = mod_mult(res, exp, m);
        exp = mod_mult(exp, exp, m);
        b >>= 1;
    }
    return res;
}
// Qualifier: Rabin-Miller强伪素数测试
bool miller_rabin(ll n, ll times)
{
    if (n < 2) return false;
    if (n == 2) return true;
    if (!(n & 1)) return false;

    ll q = n - 1;
    int k = 0;
    while (q % 2 == 0) {
        k++;
        q >>= 1;
    }
    // n - 1 = 2^k * q (q是奇素数)
    // n是素数的话，一定满足下面条件
    // (i) a^q ≡ 1 (mod n)
    // (ii) a^q, a^2q,..., a^(k-1)q 中的任何一个对n求模都为-1
    //
    // 所以、当不满足(i)(ii)中的任何一个的时候，就有3/4的概率是合成数
    //
    for (int i = 0; i < times; ++i)
    {
        ll a = rand() % (n - 1) + 1; // 从1,..,n-1随机挑一个数
        ll x = mod_exp(a, q, n);
        // 检查条件(i)
        if (x == 1) continue;
        // 检查条件(ii)
        bool found = false;
        for (int j = 0; j < k; j++)
        {
            if (x == n - 1)
            {
                found = true;
                break;
            }
            x = mod_mult(x, x, n);
        }
        if (found) continue;

        return false;
    }
    return true;
}

ll get_gcd(ll n, ll m)
{
    if (n < m) swap(n, m);
    while (m)
    {
        swap(n, m);
        m %= n;
    }
    return n;
}
// Qualifier: Pollard 因数分解算法
ll pollard_rho(ll n, int c)
{
    ll x = 2;
    ll y = 2;
    ll d = 1;
    while (d == 1)
    {
        x = mod_mult(x, x, n) + c;
        y = mod_mult(y, y, n) + c;
        y = mod_mult(y, y, n) + c;
        d = get_gcd((x - y >= 0 ? x - y : y - x), n);
    }
    if (d == n) return pollard_rho(n, c + 1);
    return d;
}

#define MAX_PRIME 200000
vector<int> primes;
vector<bool> is_prime;

// 先生成MAX_PRIME内的素数表
void init_primes()
{
    is_prime = vector<bool>(MAX_PRIME + 1, true);
    is_prime[0] = is_prime[1] = false;
    for (int i = 2; i <= MAX_PRIME; ++i)
    {
        if (is_prime[i])
        {
            primes.push_back(i);
            for (int j = i * 2; j <= MAX_PRIME; j += i)
            {
                is_prime[j] = false;
            }
        }
    }
}

// 判断是否是素数，优先查表，如果n很大用Rabin-Miller强伪素数测试
bool isPrime(ll n)
{
    if (n <= MAX_PRIME) return is_prime[n];
    else return miller_rabin(n, 20);
}

// 分解成素因子，小数用素数表，大数用Pollard 因数分解算法
void factorize(ll n, map<ll, int>& factors)
{
    if (isPrime(n))
    {
        factors[n]++;
    }
    else
    {
        for (int i = 0; i < primes.size(); ++i)
        {
            int p = primes[i];
            while (n % p == 0)
            {
                factors[p]++;
                n /= p;
            }
        }
        if (n != 1)
        {
            if (isPrime(n))
            {
                factors[n]++;
            }
            else
            {
                ll d = pollard_rho(n, 1);
                factorize(d, factors);
                factorize(n / d, factors);
            }
        }
    }
}

pair<ll, ll> solve(ll a, ll b)
{
    ll c = b / a;
    map<ll, int> factors;
    factorize(c, factors);

    vector<ll> mult_factors; // 每个质因子的n次方，比如2^2和5^1
    for (map<ll, int>::iterator it = factors.begin(); it != factors.end(); it++)
    {
        ll mul = 1;
        while (it->second)
        {
            mul *= it->first;
            it->second--;
        }
        mult_factors.push_back(mul);
    }

    ll best_sum = 1e18, best_x = 1, best_y = c;
    // 这是一个取数的过程，一共 2^size 种情况
    for (int mask = 0; mask < (1 << mult_factors.size()); ++mask)
    {
        ll x = 1;
        for (int i = 0; i < mult_factors.size(); ++i)
        {
            if (mask & (1 << i)) x *= mult_factors[i];
        }
        ll y = c / x;
        if (x < y && x + y < best_sum)
        {
            best_sum = x + y;
            best_x = x;
            best_y = y;
        }
    }
    return make_pair(best_x * a, best_y * a);
}

int main()
{
    init_primes();
    ll a, b;
    while (cin >> a >> b)
    {
        pair<ll, ll> ans = solve(a, b);
        cout << ans.first << " " << ans.second << endl;
    }
    return 0;
}