POJ 2429 关于gcd和lcm的逆求

http://wenku.baidu.com/view/fbbed5a5f524ccbff12184af.html
讲两个算法的文档

题意：
给你两个数的gcd和lcm，求他们的值
如果有多组解，输出他们值和最小的那组解
题解：
对lcm/gcd这个数字质因数分解，为了让他满足gcd和lcm不变，每一种质因数只能存在于一个数字，所以先分解质因数，但对于较大的数是不是质数的判断要用miller_rabin，对一个较大数的质因数分解需要pollard_rho，返回一个大叔随机的质因数因子。
小的素数我们可以打表筛发处理，大叔检验它是不是素数只能使用miller_rabin

#include <iostream>
#include<stdio.h>
#include<string>
#include<map>
#include<vector>
#include<stdlib.h>
typedef long long ll;
using namespace std;

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

    }
    return res;
}

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;

}
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(m == 0){return n;}
    return get_gcd(m , n%m);
}

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

            }
        }

    }
}
bool isprime(ll n){
    if(n <= MAX_PRIME)return is_prime[n];
    else return miller_rabin(n,20);
}
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;
    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;
    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(){
    cin.tie(0);
    ios::sync_with_stdio(false);
    init_primes();
    ll a, b;
    while (cin >> a >> b)
    {
        pair<ll, ll> ans = solve(a, b);
        cout << ans.first << " " << ans.second << endl;
    }

}