題目不能說(shuō)是很難，只是用到了許多數(shù)學(xué)上的知識(shí)（費(fèi)馬小定理，miller-radin，pollard-rho），還有一些算法上的知識(shí)DFS，輾轉(zhuǎn)相除。

我也很菜，一個(gè)周末的時(shí)間都用在這個(gè)題目上了，但寫(xiě)了很多很多的注釋?zhuān)ㄙM(fèi)了大量的篇幅，淺談了我對(duì)這些算法的拙見(jiàn)，希望能夠幫助大家！

文章來(lái)源地址http://www.zghlxwxcb.cn/news/detail-648263.html

#include <iostream>
#include <vector>
#include <set>
#include <map>
#include <algorithm>
using namespace std;
// 無(wú)符號(hào)的64位整數(shù)，可以支持 18446744073709551615 (2^64 - 1 )
typedef unsigned long long ll;
// 我們把前30000的素?cái)?shù)打個(gè)表
int primeLimit = 30000;
bool isPrime[30007];
// facPrime存放的是質(zhì)數(shù)因子，例如對(duì)9000000000000000000來(lái)說(shuō)，存放的是 2 2 2 2 ... 2 3 3 5 5 ...5 5 5 5 5，這些質(zhì)因子相乘等于lcm/gcd
vector<ll> facPrime;
// facVector存放的是互質(zhì)的因子，例如對(duì)9000000000000000000存放的是 2^18 3^9 5^18，即262144 9 3814697265625，這些因子相乘等于 lcm/gcd
vector<ll> facVector;
// 最終輸出的結(jié)果
ll a, b;
// 米勒拉賓算法需要一個(gè)r和一個(gè)d，在下面會(huì)解釋
ll r, d;
// 我們用埃氏篩法來(lái)找出 30000以?xún)?nèi)的素?cái)?shù)
void sieve()
{
    for (int i = 1; i <= primeLimit; i++)
    {
        isPrime[i] = true;
    }
    // 0和1不是素?cái)?shù)
    isPrime[0] = false;
    isPrime[1] = false;
    // 如果num不是素?cái)?shù)，那么一定存在除了1之外，且小于等于根號(hào)num的數(shù)字i，使得 num是i的倍數(shù)，這就是埃氏篩法的原理
    for (int i = 2; i * i <= primeLimit; i++)
    {
        // 如果i不是素?cái)?shù)，就跳過(guò)，這里跳過(guò)是為了加快效率，假如 i是18，那么18的倍數(shù)i一定是3的倍數(shù)，3的倍數(shù)在前面的循環(huán)已經(jīng)設(shè)置了，那么就不用管了
        if (!isPrime[i])
        {
            continue;
        }
        // 如果i是素?cái)?shù)，那么i的所有倍數(shù)一定都不是素?cái)?shù)
        for (int j = i * 2; j <= primeLimit; j += i)
        {
            isPrime[j] = false;
        }
    }
}
// 生成一個(gè)小于p，且大于0的隨機(jī)數(shù)
ll getLongRandom(ll p)
{
    int x = 0;
    // x需要大于0
    while (x == 0)
    {
        x = rand();
        // x需要是正數(shù)
        if (x < 0)
        {
            x = x * (-1);
        }
        // x需要小于p
        if (x >= p)
        {
            x = x % p;
        }
    }
    return (ll)x;
}
// 計(jì)算 (a * b) % mod，并且使得對(duì)于2^63以下的a,b，計(jì)算過(guò)程都不溢出，本質(zhì)的原理就是乘法，只是我們自己算是10進(jìn)制乘法，它是二進(jìn)制乘法
ll mulMod(ll a, ll b, ll mod)
{
    ll res = 0;
    while (b)
    {
        // b&1就是b與1按位與的結(jié)果，當(dāng)前僅當(dāng)b的最低位是1時(shí)生效，如 十進(jìn)制 11 & 1 = 二進(jìn)制 1011 & 1 = 1
        if (b & 1)
        {
            res = (res + a) % mod;
        }
        // 左移就是*2，同時(shí)取余防止溢出，就相等于乘法中的進(jìn)位
        a = (a << 1) % mod;
        // 右移就是/2，這樣才能不斷的計(jì)算最高位
        b = b >> 1;
    }
    return res;
}
// 用快速乘來(lái)優(yōu)化的快速冪，快速冪可以去挑戰(zhàn)程序設(shè)計(jì)里學(xué)一下，這里不再贅述
ll powMod(ll a, ll b, ll mod)
{
    ll res = 1;
    while (b)
    {
        if (b & 1)
        {
            res = mulMod(res, a, mod);
        }
        a = mulMod(a, a, mod);
        b = b >> 1;
    }
    return res;
}
// 根據(jù)大的奇數(shù)p，計(jì)算 (2^r)*d=p-1，其中d是一個(gè)奇數(shù)
void getRD(ll p)
{
    // 首先讓 d = p-1，之后不斷的除以2，直到d是一個(gè)奇數(shù)，同時(shí)每次除以2，r都+1，這樣最終 (2^r)*d=p-1了
    d = p - 1;
    r = 0;
    // d轉(zhuǎn)換為二進(jìn)制最低位如果是1，d&1就是1，那么也就是一個(gè)奇數(shù)，反之，當(dāng)d&1=0時(shí)，d也就是一個(gè)偶數(shù)
    // !(d&1) 就是 d&1 !=1 ，就是 d時(shí)一個(gè)偶數(shù)
    while (!(d & 1))
    {
        // d最左移1位，就是d/2
        d = d >> 1;
        r++;
    }
    // 例如 p=25時(shí)d=24,r=0 -> d=12,r=1 -> d=6,r=2 -> d=3 r=3 最終求出 (2^3)*3 = 25 - 1，算出 d=3,r=3
}
// 米勒拉賓算法判斷素?cái)?shù)
bool millerRabin(ll p)
{
    // 這里引入費(fèi)馬小定理， 當(dāng) p 時(shí)一個(gè)素?cái)?shù)時(shí)，對(duì)任意的整數(shù)a，一定有 pow ( a , p ) % p == a % p
    // 當(dāng) 1 <= a <= p - 1 我們對(duì)等式兩邊同時(shí)除以 a，得到 pow ( a , p - 1 ) % p == 1
    // 所以只要對(duì)于大整數(shù) p 和隨機(jī)數(shù) a ，如果不滿(mǎn)足這個(gè)規(guī)則，那么p一定不是素?cái)?shù)
    ll a = getLongRandom(p);
    if (powMod(a, p - 1, p) != 1)
    {
        return false;
    }
    // 為了提高算法的準(zhǔn)確性，我們把 pow ( a , p - 1 ) % p == 1 進(jìn)行推導(dǎo)，這里用  a ^ (p-1) 代表a的 p - 1 次方
    //  a ^ ( p - 1 ) % p == 1
    // 定義 r d，使得 ( 2 ^ r ) * d = p - 1
    // a ^ ( ( 2 ^ r ) * d ) % p == 1
    // ( a ^ ( ( 2 ^ r ) * d ) - 1 ) % p == 0
    // 這里可以使用平方差公式進(jìn)行推導(dǎo)
    // (a ^ ( (2 ^ (r - 1) ) * d) - 1) * (a ^ ( (2 ^ (r - 1) ) * d) + 1 ) % p == 0
    // (a ^ ( (2 ^ (r - 2) ) * d) - 1) * (a ^ ( (2 ^ (r - 2) ) * d) + 1 ) * (a ^ ( (2 ^ (r - 1) ) * d) + 1) % p == 0
    // ....
    // (a ^ ( (2 ^ (r - r) ) * d) - 1) * (a ^ ( (2 ^ (r - r) ) * d) + 1) * ...  (a ^ ( (2 ^ (r - 1) ) * d) + 1) % p == 0
    // 同時(shí)p是素?cái)?shù)，如果這些項(xiàng)相乘 % p == 0，那么其中的某一項(xiàng)一定等于0，那么就推出了另外一個(gè)判斷依據(jù)
    // 在  0 <= i < r 必定存在 a ^ ( (2 ^ i) * d) 等于 p 或等于 1，否則 p 一定不是素?cái)?shù)
    // a ^ ( (2 ^ i) * d) = (a ^ d) ^ (2 ^ i)
    ll k = powMod(a, d, p);
    for (int i = 0; i < r; i++)
    {
        // 計(jì)算每一項(xiàng)
        ll parameter = powMod(k, powMod(2, i, p), p);
        if (parameter == 1 || parameter == p - 1)
        {
            return true;
        }
    }
    return false;
}
// 米勒拉賓具有概率性，我們連續(xù)判斷20次，如果都滿(mǎn)足，那么它一定是素?cái)?shù)！
bool multiMillerRabin(ll p)
{
    // 對(duì)于固定的p，只需要計(jì)算1次，使得 (2 ^ ( r ) * d) == p - 1 的r和d
    getRD(p);
    for (int i = 0; i < 20; i++)
    {
        if (!millerRabin(p))
        {
            return false;
        }
    }
    return true;
}
// 質(zhì)數(shù)判斷方法，對(duì)于30000以?xún)?nèi)的小素?cái)?shù)，直接用 “埃氏篩法” ，對(duì)于30000以上的，偶數(shù)直接可以確定不是質(zhì)數(shù)，奇數(shù)的話(huà)，2 ^ 63 次冪以下，用米勒拉賓沒(méi)問(wèn)題
bool judgePrime(ll p)
{
    if (p <= primeLimit)
    {
        return isPrime[(int)p];
    }
    else if (!(p & 1))
    {
        // 舉個(gè)例子，p = 30002 ，也就是二進(jìn)制的 0111 0101 0011 0010 和 0000 0000 0000 0001 按位與一定是 0 ，!0就是1就是true
        return false;
    }
    else
    {
        // 直接用米勒拉賓判斷20次
        return multiMillerRabin(p);
    }
}
// 輾轉(zhuǎn)相除求出gcd
ll gcd(ll a, ll b)
{
    if (b == 0)
    {
        return a;
    }
    else
    {
        return gcd(b, a % b);
    }
}
// | a - b |  a,b相減的絕對(duì)值
ll absSub(ll a, ll b)
{
    if (a > b)
    {
        return a - b;
    }
    else
    {
        return b - a;
    }
}
// 使用ρ算法，來(lái)高效的分解質(zhì)因子！
void pollardRho(ll p)
{
    // 如果p是素?cái)?shù)，那就代表找到了一個(gè)質(zhì)因子
    if (judgePrime(p) || p <= 1)
    {
        // facPrime存放的是質(zhì)數(shù)因子，例如對(duì)9000000000000000000來(lái)說(shuō)，存放的是 2 2 2 2 ... 2 3 3 5 5 ...5 5 5 5 5，這些質(zhì)因子相乘等于lcm/gcd
        if (p >= 2)
        {
            facPrime.push_back(p);
            return;
        }
    }
    // 我還沒(méi)有理解ρ算法特別深，但是可以淺談下自己的拙見(jiàn)！如果錯(cuò)誤的話(huà)，大家批評(píng)指正！
    // 它是產(chǎn)生隨機(jī)數(shù)x，然后判斷隨機(jī)數(shù)x 和 p 是否存在 大于1的最大公約數(shù) g
    // 存在的話(huà)就對(duì) g 和 p / g 分別遞歸應(yīng)用ρ算法，直到找到素?cái)?shù)，是一個(gè)遞歸出口
    // 但是如果單純的隨機(jī)產(chǎn)生x，實(shí)在速度太慢
    // 這里引入一個(gè)生日悖論，不考慮閏年時(shí)，23個(gè)人中有兩個(gè)人生日相同的概率超過(guò) 二分之一
    // 1 - ( 365.0 / 365.0 ) * ( 364.0 / 365.0 ) * ( 363.0 / 365.0 ) * ( 343.0 / 365.0 ) > 0.5
    // 所以我們也可以學(xué)一下這個(gè)生日悖論，產(chǎn)生一堆隨機(jī)數(shù) x1 x2 x3 x4 ... xn
    // 當(dāng)n足夠大時(shí)，那么其中的任意兩個(gè)數(shù) xi xj 的差與 p 存在大于1的最大公約數(shù)的概率會(huì)很大！
    // 同時(shí)如果循環(huán)的去算 xi - xj，也太慢了！所以讓 xi 與 xi-1 存在一個(gè)關(guān)系，即 x2 = x1*x1 +c ,x3 = x2*x2+c,x4=x3*x3+c，(c和x1是小于p大于0的隨機(jī)數(shù)）
    // 那么算出來(lái)的隨機(jī)數(shù)就存在相互關(guān)系，產(chǎn)生的曲線是一個(gè)ρ型，當(dāng)其中某一個(gè) xj - xi 與 p 存在大于1 的最大公約數(shù)時(shí)
    // 那么對(duì)于 任意的 對(duì)于任意的 xk,k>j也一定有 xk - xi 與p存在大于1的最大公約數(shù)（這里是我的理解，我不會(huì)證明）
    // 同時(shí)曲線是一個(gè) ρ 型，那么容易出現(xiàn) 循環(huán)，比如 x1 = 3 a = 2 p = 20
    // x1 = 3 ,x2 = (3*3+2)%20=11,x3 = (11*11 + 2) %20=3 ,x3又回到了x1，出現(xiàn)了si循環(huán)！
    // 所以我們讓y以x兩倍的速度去增長(zhǎng)，即每次循環(huán)y算兩次，x算一次，那么就類(lèi)似于數(shù)學(xué)中的追擊問(wèn)題，定會(huì)有x和y相遇
    // 當(dāng)x==y，我們重新生成隨機(jī)的x1和c重新找因子！
    // 然后遞歸的時(shí)候，判斷如果p是一個(gè)質(zhì)數(shù)，就是遞歸出口！
    bool isFind = false;
    // 找不到就一直找，找到了直接出循環(huán)
    while (!isFind)
    {
        ll x = getLongRandom(p);
        ll c = getLongRandom(p);
        ll y = x;
        y = (mulMod(y, y, p) + c) % p;
        while (x != y)
        {
            ll gcdNum = gcd(p, absSub(y, x));
            if (gcdNum > 1)
            {
                // 找到了最大公約數(shù)，那么最大公約數(shù)就是因子，用 p 除以這個(gè)因子，可以得到另一個(gè)因子，再分別帶進(jìn)去循環(huán)，最終一定可以找到所有的質(zhì)數(shù)因子
                pollardRho(gcdNum);
                pollardRho(p / gcdNum);
                isFind = true;
                break;
            }
            x = (mulMod(x, x, p) + c) % p;
            y = (mulMod(y, y, p) + c) % p;
            y = (mulMod(y, y, p) + c) % p;
        }
    }
}
// 用dfs來(lái)判斷因子間排列組合相乘的各種情況！找到最小的 a+b
void dfs(int i, ll current, ll p)
{
    if (i >= facVector.size())
    {
        return;
    }
    // 乘以這個(gè)因子的情況算一下
    ll currentA = current * facVector[i];
    ll currentB = p / currentA;
    if ((currentA + currentB) < (a + b))
    {
        a = currentA;
        b = currentB;
    }
    dfs(i + 1, current * facVector[i], p);
    // 不乘以這個(gè)因子的情況再算一下
    currentA = current;
    currentB = p / currentA;
    if ((currentA + currentB) < (a + b))
    {
        a = currentA;
        b = currentB;
    }
    dfs(i + 1, current, p);
}
int main()
{
    ll lcm, gcd;
    // 初始化埃氏篩法
    sieve();
    while (~scanf("%lld%lld", &gcd, &lcm))
    {
        // 不難看出，當(dāng) a b 的最大公約數(shù)是gcd，最小公倍數(shù)是lcm時(shí)
        // 一定有 ( lcm / gcd ) = ( a / gcd ) * ( b / gcd )
        // 同時(shí)也一定有( a / gcd ) * ( b / gcd )互質(zhì)，如果它們不互質(zhì)，那最大公約數(shù)不會(huì)是gcd，應(yīng)該是它們的公約數(shù) * gcd才對(duì)！
        // 因此我們要根據(jù)lcm和gcd去求解原數(shù)字，其實(shí)很簡(jiǎn)單
        // 定義 p = lcm/gcd ，把p變成質(zhì)因子相乘的形式
        // 假設(shè)  p = 9000000000000000000 ，就把它變成 2 * 2 * 2 * 2 * 2 * 2 * 2 .. * 2 * 3 * 3 * 5 * 5 .... * 5
        // 設(shè) (a/gcd)=a1,(b/gcd)=b1
        // 那么一定有 a1 * b1 = 2 * 2 * 2 * 2 * 2 * 2 * 2 .. * 2 * 3 * 3 * 5 * 5 .... * 5 = (2 ^ 18)*(3 ^ 2)*(5 ^ 18)
        // 同時(shí)a1和b1還互質(zhì)，不能存在公約數(shù)，那么a的質(zhì)因子和b的質(zhì)因子就不能有重復(fù)，有重復(fù)的話(huà)，gcd(a1,b1)無(wú)法等于1，a和b的gcd就與輸入的gcd不符
        // 那么我們就推導(dǎo)一下思路：
        // 輸入了 lcm 和 gcd，我們用 lcd / gcd 算出一個(gè) p，然后把 p 推導(dǎo)成 p = a1 * b1 的形式，找出滿(mǎn)足a1 和 b1 互質(zhì)，且(a1+b1)最小的a1和b1
        // 然后輸出 a1 * gcd 和 b1 * gcd就是結(jié)果
        // 同時(shí)當(dāng)lcm與gcd相等時(shí)，不需要找，直接輸出gcd就好，當(dāng)且僅當(dāng)a和b相等時(shí)，它們的gcd和lcm才相等，而且都等于a
        // 對(duì)于 9223372036854775807 規(guī)模的p， 把 p 推到成 p = a1 * b1 的過(guò)程并不容易，所以我們引入了 pollard-rho 算法，求出 p 的質(zhì)因子
        // 把 p 寫(xiě)成一堆質(zhì)數(shù)的次方，相乘的形式，然后 dfs，依次判斷每一項(xiàng)乘上與不乘上的值，找出最小的 a1 + b1，最后輸出答案
        // a1 * b1 = (2 ^ 18) * (3 ^ 2) * (5 ^ 18)，只要 a 和 b是右邊三項(xiàng)中的任意1或多項(xiàng) 乘出來(lái)的，那么它們一定互質(zhì)！
        ll p = lcm / gcd;
        if (p == 1)
        {
            // 當(dāng)且僅當(dāng) a==b，的時(shí)候，lcm 與 gcd才相等，而且都等于a
            printf("%lld %lld\n", gcd, gcd);
        }
        else
        {
            // 質(zhì)因子數(shù)組初始化
            if (facPrime.size() > 0)
            {
                facPrime.clear();
            }
            // 因子數(shù)組初始化
            if (facVector.size() > 0)
            {
                facVector.clear();
            }
            pollardRho(p);
            // 用來(lái)統(tǒng)計(jì)每個(gè)質(zhì)因子出現(xiàn)數(shù)量的map
            map<ll, ll> countMap;
            // 用來(lái)對(duì)質(zhì)因子去重的set
            set<ll> distinctSet;
            for (int i = 0; i < facPrime.size(); i++)
            {
                countMap[facPrime[i]] = 0;
                distinctSet.insert(facPrime[i]);
            }
            for (int i = 0; i < facPrime.size(); i++)
            {
                ll count = countMap[facPrime[i]];
                countMap[facPrime[i]] = count + 1;
            }
            for (set<ll>::iterator ite = distinctSet.begin(); ite != distinctSet.end(); ite++)
            {
                ll prime = *ite;
                prime = powMod(prime, countMap[prime], p);
                if (prime == 0)
                {
                    prime = p;
                }
                facVector.push_back(prime);
            }
            a = facVector[0];
            b = p / a;
            dfs(0, 1, p);
            if (a > b)
            {
                swap(a, b);
            }
            a = a * gcd;
            b = b * gcd;
            printf("%lld %lld\n", a, b);
        }
    }
    return 0;
}

到了這里，關(guān)于POJ 2429 Miller-rabin素?cái)?shù)判定 + pollard-rho質(zhì)因子分解 + 埃氏篩法的文章就介紹完了。如果您還想了解更多內(nèi)容，請(qǐng)?jiān)谟疑辖撬阉鱐OY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！