POJ 2429 Miller-rabin素数判定 + pollard-rho质因子分解 + 埃氏筛法

题目不能说是很难，只是用到了许多数学上的知识（费马小定理，miller-radin，pollard-rho），还有一些算法上的知识DFS，辗转相除。

我也很菜，一个周末的时间都用在这个题目上了，但写了很多很多的注释，花费了大量的篇幅，浅谈了我对这些算法的拙见，希望能够帮助大家！

#include <iostream>
#include <vector>
#include <set>
#include <map>
#include <algorithm>
using namespace std;
// 无符号的64位整数，可以支持 18446744073709551615 (2^64 - 1 )
typedef unsigned long long ll;
// 我们把前30000的素数打个表
int primeLimit = 30000;
bool isPrime[30007];
// facPrime存放的是质数因子，例如对9000000000000000000来说，存放的是 2 2 2 2 ... 2 3 3 5 5 ...5 5 5 5 5，这些质因子相乘等于lcm/gcd
vector<ll> facPrime;
// facVector存放的是互质的因子，例如对9000000000000000000存放的是 2^18 3^9 5^18，即262144 9 3814697265625，这些因子相乘等于 lcm/gcd
vector<ll> facVector;
// 最终输出的结果
ll a, b;
// 米勒拉宾算法需要一个r和一个d，在下面会解释
ll r, d;
// 我们用埃氏筛法来找出 30000以内的素数
void sieve()
{
    for (int i = 1; i <= primeLimit; i++)
    {
        isPrime[i] = true;
    }
    // 0和1不是素数
    isPrime[0] = false;
    isPrime[1] = false;
    // 如果num不是素数，那么一定存在除了1之外，且小于等于根号num的数字i，使得 num是i的倍数，这就是埃氏筛法的原理
    for (int i = 2; i * i <= primeLimit; i++)
    {
        // 如果i不是素数，就跳过，这里跳过是为了加快效率，假如 i是18，那么18的倍数i一定是3的倍数，3的倍数在前面的循环已经设置了，那么就不用管了
        if (!isPrime[i])
        {
            continue;
        }
        // 如果i是素数，那么i的所有倍数一定都不是素数
        for (int j = i * 2; j <= primeLimit; j += i)
        {
            isPrime[j] = false;
        }
    }
}
// 生成一个小于p，且大于0的随机数
ll getLongRandom(ll p)
{
    int x = 0;
    // x需要大于0
    while (x == 0)
    {
        x = rand();
        // x需要是正数
        if (x < 0)
        {
            x = x * (-1);
        }
        // x需要小于p
        if (x >= p)
        {
            x = x % p;
        }
    }
    return (ll)x;
}
// 计算 (a * b) % mod，并且使得对于2^63以下的a,b，计算过程都不溢出，本质的原理就是乘法，只是我们自己算是10进制乘法，它是二进制乘法
ll mulMod(ll a, ll b, ll mod)
{
    ll res = 0;
    while (b)
    {
        // b&1就是b与1按位与的结果，当前仅当b的最低位是1时生效，如 十进制 11 & 1 = 二进制 1011 & 1 = 1
        if (b & 1)
        {
            res = (res + a) % mod;
        }
        // 左移就是*2，同时取余防止溢出，就相等于乘法中的进位
        a = (a << 1) % mod;
        // 右移就是/2，这样才能不断的计算最高位
        b = b >> 1;
    }
    return res;
}
// 用快速乘来优化的快速幂，快速幂可以去挑战程序设计里学一下，这里不再赘述
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;
}
// 根据大的奇数p，计算 (2^r)*d=p-1，其中d是一个奇数
void getRD(ll p)
{
    // 首先让 d = p-1，之后不断的除以2，直到d是一个奇数，同时每次除以2，r都+1，这样最终 (2^r)*d=p-1了
    d = p - 1;
    r = 0;
    // d转换为二进制最低位如果是1，d&1就是1，那么也就是一个奇数，反之，当d&1=0时，d也就是一个偶数
    // !(d&1) 就是 d&1 !=1 ，就是 d时一个偶数
    while (!(d & 1))
    {
        // d最左移1位，就是d/2
        d = d >> 1;
        r++;
    }
    // 例如 p=25时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
}
// 米勒拉宾算法判断素数
bool millerRabin(ll p)
{
    // 这里引入费马小定理， 当 p 时一个素数时，对任意的整数a，一定有 pow ( a , p ) % p == a % p
    // 当 1 <= a <= p - 1 我们对等式两边同时除以 a，得到 pow ( a , p - 1 ) % p == 1
    // 所以只要对于大整数 p 和随机数 a ，如果不满足这个规则，那么p一定不是素数
    ll a = getLongRandom(p);
    if (powMod(a, p - 1, p) != 1)
    {
        return false;
    }
    // 为了提高算法的准确性，我们把 pow ( a , p - 1 ) % p == 1 进行推导，这里用  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
    // 这里可以使用平方差公式进行推导
    // (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
    // 同时p是素数，如果这些项相乘 % p == 0，那么其中的某一项一定等于0，那么就推出了另外一个判断依据
    // 在  0 <= i < r 必定存在 a ^ ( (2 ^ i) * d) 等于 p 或等于 1，否则 p 一定不是素数
    // a ^ ( (2 ^ i) * d) = (a ^ d) ^ (2 ^ i)
    ll k = powMod(a, d, p);
    for (int i = 0; i < r; i++)
    {
        // 计算每一项
        ll parameter = powMod(k, powMod(2, i, p), p);
        if (parameter == 1 || parameter == p - 1)
        {
            return true;
        }
    }
    return false;
}
// 米勒拉宾具有概率性，我们连续判断20次，如果都满足，那么它一定是素数！
bool multiMillerRabin(ll p)
{
    // 对于固定的p，只需要计算1次，使得 (2 ^ ( r ) * d) == p - 1 的r和d
    getRD(p);
    for (int i = 0; i < 20; i++)
    {
        if (!millerRabin(p))
        {
            return false;
        }
    }
    return true;
}
// 质数判断方法，对于30000以内的小素数，直接用 “埃氏筛法” ，对于30000以上的，偶数直接可以确定不是质数，奇数的话，2 ^ 63 次幂以下，用米勒拉宾没问题
bool judgePrime(ll p)
{
    if (p <= primeLimit)
    {
        return isPrime[(int)p];
    }
    else if (!(p & 1))
    {
        // 举个例子，p = 30002 ，也就是二进制的 0111 0101 0011 0010 和 0000 0000 0000 0001 按位与一定是 0 ，!0就是1就是true
        return false;
    }
    else
    {
        // 直接用米勒拉宾判断20次
        return multiMillerRabin(p);
    }
}
// 辗转相除求出gcd
ll gcd(ll a, ll b)
{
    if (b == 0)
    {
        return a;
    }
    else
    {
        return gcd(b, a % b);
    }
}
// | a - b |  a,b相减的绝对值
ll absSub(ll a, ll b)
{
    if (a > b)
    {
        return a - b;
    }
    else
    {
        return b - a;
    }
}
// 使用ρ算法，来高效的分解质因子！
void pollardRho(ll p)
{
    // 如果p是素数，那就代表找到了一个质因子
    if (judgePrime(p) || p <= 1)
    {
        // facPrime存放的是质数因子，例如对9000000000000000000来说，存放的是 2 2 2 2 ... 2 3 3 5 5 ...5 5 5 5 5，这些质因子相乘等于lcm/gcd
        if (p >= 2)
        {
            facPrime.push_back(p);
            return;
        }
    }
    // 我还没有理解ρ算法特别深，但是可以浅谈下自己的拙见！如果错误的话，大家批评指正！
    // 它是产生随机数x，然后判断随机数x 和 p 是否存在 大于1的最大公约数 g
    // 存在的话就对 g 和 p / g 分别递归应用ρ算法，直到找到素数，是一个递归出口
    // 但是如果单纯的随机产生x，实在速度太慢
    // 这里引入一个生日悖论，不考虑闰年时，23个人中有两个人生日相同的概率超过 二分之一
    // 1 - ( 365.0 / 365.0 ) * ( 364.0 / 365.0 ) * ( 363.0 / 365.0 ) * ( 343.0 / 365.0 ) > 0.5
    // 所以我们也可以学一下这个生日悖论，产生一堆随机数 x1 x2 x3 x4 ... xn
    // 当n足够大时，那么其中的任意两个数 xi xj 的差与 p 存在大于1的最大公约数的概率会很大！
    // 同时如果循环的去算 xi - xj，也太慢了！所以让 xi 与 xi-1 存在一个关系，即 x2 = x1*x1 +c ,x3 = x2*x2+c,x4=x3*x3+c，(c和x1是小于p大于0的随机数）
    // 那么算出来的随机数就存在相互关系，产生的曲线是一个ρ型，当其中某一个 xj - xi 与 p 存在大于1 的最大公约数时
    // 那么对于 任意的 对于任意的 xk,k>j也一定有 xk - xi 与p存在大于1的最大公约数（这里是我的理解，我不会证明）
    // 同时曲线是一个 ρ 型，那么容易出现 循环，比如 x1 = 3 a = 2 p = 20
    // x1 = 3 ,x2 = (3*3+2)%20=11,x3 = (11*11 + 2) %20=3 ,x3又回到了x1，出现了si循环！
    // 所以我们让y以x两倍的速度去增长，即每次循环y算两次，x算一次，那么就类似于数学中的追击问题，定会有x和y相遇
    // 当x==y，我们重新生成随机的x1和c重新找因子！
    // 然后递归的时候，判断如果p是一个质数，就是递归出口！
    bool isFind = false;
    // 找不到就一直找，找到了直接出循环
    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)
            {
                // 找到了最大公约数，那么最大公约数就是因子，用 p 除以这个因子，可以得到另一个因子，再分别带进去循环，最终一定可以找到所有的质数因子
                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来判断因子间排列组合相乘的各种情况！找到最小的 a+b
void dfs(int i, ll current, ll p)
{
    if (i >= facVector.size())
    {
        return;
    }
    // 乘以这个因子的情况算一下
    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);
    // 不乘以这个因子的情况再算一下
    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))
    {
        // 不难看出，当 a b 的最大公约数是gcd，最小公倍数是lcm时
        // 一定有 ( lcm / gcd ) = ( a / gcd ) * ( b / gcd )
        // 同时也一定有( a / gcd ) * ( b / gcd )互质，如果它们不互质，那最大公约数不会是gcd，应该是它们的公约数 * gcd才对！
        // 因此我们要根据lcm和gcd去求解原数字，其实很简单
        // 定义 p = lcm/gcd ，把p变成质因子相乘的形式
        // 假设  p = 9000000000000000000 ，就把它变成 2 * 2 * 2 * 2 * 2 * 2 * 2 .. * 2 * 3 * 3 * 5 * 5 .... * 5
        // 设 (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)
        // 同时a1和b1还互质，不能存在公约数，那么a的质因子和b的质因子就不能有重复，有重复的话，gcd(a1,b1)无法等于1，a和b的gcd就与输入的gcd不符
        // 那么我们就推导一下思路：
        // 输入了 lcm 和 gcd，我们用 lcd / gcd 算出一个 p，然后把 p 推导成 p = a1 * b1 的形式，找出满足a1 和 b1 互质，且(a1+b1)最小的a1和b1
        // 然后输出 a1 * gcd 和 b1 * gcd就是结果
        // 同时当lcm与gcd相等时，不需要找，直接输出gcd就好，当且仅当a和b相等时，它们的gcd和lcm才相等，而且都等于a
        // 对于 9223372036854775807 规模的p， 把 p 推到成 p = a1 * b1 的过程并不容易，所以我们引入了 pollard-rho 算法，求出 p 的质因子
        // 把 p 写成一堆质数的次方，相乘的形式，然后 dfs，依次判断每一项乘上与不乘上的值，找出最小的 a1 + b1，最后输出答案
        // a1 * b1 = (2 ^ 18) * (3 ^ 2) * (5 ^ 18)，只要 a 和 b是右边三项中的任意1或多项 乘出来的，那么它们一定互质！
        ll p = lcm / gcd;
        if (p == 1)
        {
            // 当且仅当 a==b，的时候，lcm 与 gcd才相等，而且都等于a
            printf("%lld %lld\n", gcd, gcd);
        }
        else
        {
            // 质因子数组初始化
            if (facPrime.size() > 0)
            {
                facPrime.clear();
            }
            // 因子数组初始化
            if (facVector.size() > 0)
            {
                facVector.clear();
            }
            pollardRho(p);
            // 用来统计每个质因子出现数量的map
            map<ll, ll> countMap;
            // 用来对质因子去重的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;
}