POJ：2429-GCD & LCM Inverse（素数判断神题）（Millar-Rabin素性判断和Pollard-rho因子分解）

原题链接：http://poj.org/problem?id=2429

GCD & LCM Inverse

Time Limit: 2000MS Memory Limit: 65536K
Total Submissions: 17639 Accepted: 3237

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

---

解题心得：

1. 题意就是给你LCM和GCD，要求你输出两个数a和b，a和b要满足题目所给的LCM和GCD同时要使a+b最小。
2. 之前做了给你a和b叫你输出LCM和GCD，所以在看到这个题的时候以为难度不大，但是由于数据量真的太大，所以去网上看了一下大神们的代码，结果就花费了两天的时间去学习前备知识Millar-Rabin强伪素数判断和Pollard-rho素数分解，然后再慢慢码代码。
3. 关于Millar-Rabin是个什么，可以去看一下这篇文章《素数和素性判断》，上面说得十分清楚，然后看完这个前备知识之后剩下的就很好解决了，还有就是感觉有一点点玄学的Pollard-rho素数分解。
4. 然后就可以看这个题的思路了，给你的是LCM和GCD，可以设得到的答案为a和b，可以肯定的是（a/gcd）*（b/gcd）= lcm/gcd。那么（lcm/gcd）就可以看作两个互素数的乘积（a/gcd和b/gcd一定是两个互素的数），这个时候lcm和gcd都是知道的，就可以用lcm/gcd，假设c=lcm/gcd，那么就可以将c分解成多个素数的乘积，既然素数的乘积，那么在分解的时候必须要判断当前数是不是素数，使用素筛范围不可能太大，就要用到millar-rabin判断素数，如果是素数就停止分解，如果不是，首先可以筛选的素数范围内得到就直接得到一些素数因子，否则使用Pollard-rho算法分解出因子，然后递归，最后找到所有的素数因子之后再选择，这就只能枚举，可以使用状压。
5. 大概思路就是这样的，实在不是很明白的可以先看看代码，了解整体结构。但是还有需要注意的地方就是数据范围太大，不可以两个数直接相乘，需要用快数乘法，化乘法为加法，原理和快速幂一样，这样一边加一边mod才不会超过long long范围。

---

#include <algorithm>
#include <stdio.h>
#include <vector>
#include <cstring>
#include <map>
#include <iostream>
using namespace std;
typedef long long ll;
const int maxn = 2e5;

bool prim[maxn];
vector <ll> prim_num;

ll __gcd(ll n,ll m) {
    if(m == 0)
        return n;
    return __gcd(m,n%m);
}

ll mod_mult(ll x,ll q,ll mod) {//快速乘法
    ll res = 0;
    while(q) {
        if(q&1) {
            res = res + x;
            if(res > mod) res -= mod;
        }
        x <<= 1;
        if(x > mod) x -= mod;
        q >>= 1;
    }
    return res;
}

ll mod_exp(ll x,ll q,ll mod) {//快速幂
    ll res = 1;
    while(q) {
        if(q&1)
            res = mod_mult(res, x, mod);//注意不能直接相乘
        x = mod_mult(x,x,mod);
        q >>= 1;
    }
    return res;
}

bool millar_rabin(ll n) {
    ll q = n-1;
    if(n < 2 || !(n&1))
        return false;
    if(n == 2)
        return true;
    ll k = 0;
    while(q%2 == 0) {
        k++;
        q >>= 1;
    }

    for(int i=0;i<20;i++) {//这里重复判断20次，将错误的机率减小
        ll temp = rand()%(n-1) + 1;
        ll x = mod_exp(temp,q,n);
        if(x == 1)
            continue;
        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;
}

void get_prim() {//素数筛选
    prim[1] = prim[0] = true;
    for(int i=2;i<maxn;i++) {
        if(prim[i])
            continue;
        prim_num.push_back(i);
        for(int j=i*2;j<maxn;j+=i)
            prim[j] = true;
    }
}

bool is_prim(ll x) {//判断是不是素数
    if(x < maxn)
        return !prim[x];
    return millar_rabin(x);
}

ll pollard_rho(ll n,ll c) {
    ll x = 2;
    ll y = 2;
    ll d = 1;

    while(d == 1 || d == n) {//分解到为1和n的因子则继续分解
        x = mod_mult(x,x,n) + c;
        y = mod_mult(y,y,n) + c;
        y = mod_mult(y,y,n) + c;
        d = __gcd(n,x>y?x-y:y-x);
    }
    return d;
}

void factorize(map<ll,ll>& maps,ll n) {
    if(is_prim(n))//如果n已经是素数了，那么就不用分解了
        maps[n]++;
    else {
        for(int i=0;i<prim_num.size();i++) {//可以在筛选出来的素数中找到的直接找
            while(n%prim_num[i] == 0) {
                maps[prim_num[i]]++;
                n /= prim_num[i];
            }
        }
        if(n != 1) {
            if(is_prim(n))//如果在前面已经分解完了那就结束，不然继续分解超过筛选出素数范围
                maps[n]++;
            else {
                ll d = pollard_rho(n,1);//如果分解出一个因子是d，那么另一个因子一定是n/d
                factorize(maps, d);//递归下去
                factorize(maps,n/d);
            }
        }
    }
}

pair <ll,ll> Solve(ll n,ll m) {
    ll c = m/n;
    map <ll,ll> maps;//用于储存分解出来的素数因子
    factorize(maps,c);//分解c

    vector <ll> mult_factors;
    map <ll,ll> :: iterator iter1;
    for(iter1=maps.begin();iter1!=maps.end();iter1++) {
        ll mul = 1;
        while(iter1->second){
            mul *= iter1->first;
            iter1->second--;
        }
        mult_factors.push_back(mul);
    }

    ll Max = 1e18,ans_x = 1,ans_y = c;
    for(ll i=0;i<(1 << mult_factors.size());i++) {//枚举所有可能的因子，找到最小的那一个，总共有2的n次方种可能性，n是分解出来的因子数，然后用状压，枚举
        ll x = 1;
        for(ll j=0;j<mult_factors.size();j++) {
            if(i & (1 << j))
                x *= mult_factors[j];
        }
        ll y = c/x;//如果一个因子为x，则另一个必定为c/x
        if(x < y && x + y < Max) {
            Max = x + y;
            ans_x = x;
            ans_y = y;
        }
    }
    pair<ll,ll> p;
    p.first = ans_x*n;//记得要乘上之前除去的GCD
    p.second = ans_y*n;
    return p;
}

int main() {
    ll gcd,lcm;
    get_prim();
    while(scanf("%lld%lld",&gcd,&lcm) != EOF) {
        pair<ll,ll> ans = Solve(gcd,lcm);
        printf("%lld %lld
",ans.first,ans.second);
    }
    return 0;
}