题目不能说是很难，只是用到了许多数学上的知识（费马小定理，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;
}
![[HDLBits] Exams/m2014 q4d](https://img-blog.csdnimg.cn/img_convert/b8d5dcb031b10e9c3ac7f81d36e10fa6.png)
