Mathematics:GCD & LCM Inverse(POJ 2429)

　　　　　　　　　　　　　　　　


　　　　　　　　　　根据最大公约数和最小公倍数求原来的两个数

　　题目大意，不翻译了，就是上面链接的意思。

　　具体思路就是要根据数论来，设a和b的GCD（最大公约数）和LCM（最小公倍数），则a/GCD*b/GCD=LCM/GCD，我们只用枚举LCM/GCD的所有质因数就可以了，然后把相应的质因数乘以GCD即可得出答案。

　　找素数很简单，用Miller_Rabin求素数的方法，可以多求几次提高正确率，原理就是用的费马定理：如果P是素数，则A^(p-1)mod P恒等于1，为了绕过Carmichael数，采用费马小定理：如果n是素数，则存在x(0<x<n)，(x*x)mod n 要么是1要么是n-1，否则，x就是合数。

　　另外就是要把因数分解了，这里暴力解法貌似可以，但是那种方法太笨了（也就是枚举），我们采用一个O(N^1/4)的算法Pollard_Rho快速因数分解，具体证明点我，我们把结论用上就可以了，具体将在代码中呈现。

　　因数分解以后，剩下就只用DFS就可以了，注意要把所有重复的质因数先成起来，那样我们找解的时候就可以保证我们找到的解都是互质的

　　参考http://www.hankcs.com/program/cpp/poj-2429-gcd-lcm-inverse.html

　　　　/article/6437000.html

#include <iostream>
#include <functional>
#include <algorithm>
#define MAX_N 1000

using namespace std;

long long gcd(long long, long long);
bool Miller_Rabin(const long long);
long long witness(long long, long long, long long);
long long Pollard_Rho(long long,long long);
long long Multi_Mod(long long, long long);
void Find_Factors(long long, int *const, long long);
void DFS(long long, long long, int,const int);
static long long factors[MAX_N],factors_sum[MAX_N],a, b, min_sum;

int main(void)
{
long long n, GCD, LCM;

while (~scanf("%lld %lld", &GCD, &LCM))
{
int len = 0, len_factors_sum = 0;
n = LCM / GCD;
if (Miller_Rabin(n))
printf("%lld %lld\n", GCD, n*GCD);
else if (LCM == GCD)
printf("%lld %lld\n", GCD, GCD);
else
{
Find_Factors(n, &len, 120);//120是经验值
sort(factors, factors + len);
factors_sum[0] = factors[0];
for (int i = 1; i < len; i++)//把相同的质数全部成起来，那么当DFS的时候我们只用找这些乘积就可以了（保证a,b一定互素）
{
if (factors[i] == factors[i - 1])
factors_sum[len_factors_sum] *= factors[i];
else
factors_sum[++len_factors_sum] = factors[i];
}
a = factors[0]; b = n / a;
min_sum = b + a;
DFS(1, 1, 0, len_factors_sum + 1);//找解，DFS枚举就可以了
if (a < b)
printf("%lld %lld\n", a*GCD, b*GCD);
else
printf("%lld %lld\n", b*GCD, a*GCD);
}
}
return 0;
}

void Find_Factors(long long n,int *const len,long long times)
{
if (n == 1)
return;
else if (Miller_Rabin(n))
{
factors[(*len)++] = n;//分解质因数
return;
}
else
{
long long p = n;
long long k = times;
while (p >= n)
p = Pollard_Rho(n, k--);
Find_Factors(p, len, times);
Find_Factors(n / p, len, times);
}
}

bool Miller_Rabin(const long long n)
{
if (n == 2)
return true;
if (n == 1 || n < 0)
return false;
else
{
for (int i = 0; i < 5; i++)//Miller_Rabin测试方法+费马定理，叠5次减少出错几率
if (!(witness((long long)rand() % (n - 1) + 1, n - 1, n) == 1))
return false;
return true;
}
}

long long witness(long long coe, long long level, long long n)
{
long long x, y;
if (level == 0)
return 1;//到达最后一层，开始后序遍历

x = witness(coe, level >> 1, n);//level以幂次递减
if (x == 0)
return 0;//如果x出的结果是0，那么n一定是一个合数

y = (x*x) % n;
if (y == 1 && x != 1 && x != n - 1)
return 0;//费马小定理，如果一个数是素数，则x*x对n的模一定是1或者是n-1,如果不是，则是合数
if (level % 2  == 1)
y = (coe*y) % n;//和幂运算的道理是一样的

return y;
}

long long Pollard_Rho(long long n,long long c)
{
long long x, y, k = 2, d;
y = x = rand() % (n - 1) + 1;//y和x的初始值都是定任意一个常数，然后直到找到非平常因子为止

for (int i = 1;; i++)
{
x = (Multi_Mod(x, n) + c) % n;//算f(x),f(x)的定义见多项式乘法f(x)=x^2+c
d = gcd((y - x + n) % n, n);//计算|y-x|与n的最大公因数，当y==x时，返回n，说明在这个c下无法产生非平常因子
if (1 < d && d < n)
return d;//如果得出d，那么d就是因数之一（不一定是质数，要继续判断）
else if (y == x)
return n;
else if (i == k)//brent判据，目的就是找到在偶数周期内找到gcd(x(k)-x(i/2))
{
y = x;
k <<= 1;
}
}
return n;
}

long long gcd(long long a, long long b)
{
if (b == 0) return a;
return gcd(b, a%b);
}

long long Multi_Mod(long long x, long long mod)//Pollard_Rho用到的多项式算法y=(x^2 + c)mod n
{
long long ans = 0, p = x;

while (p)
{
if (p & 1)
ans = (ans + x) % mod;
x = (x << 1) % mod;//记得取模
p >>= 1;
}
return ans;
}

void DFS(long long tmpa, long long tmpb, int level, const int len_factors_sum)
{
if (level == len_factors_sum)
{
if (tmpa + tmpb < min_sum)
{
a = tmpa; b = tmpb;
min_sum = tmpa + tmpb;
}
}
else
{
DFS(tmpa*factors_sum[level], tmpb, level + 1, len_factors_sum);
DFS(tmpa, tmpb*factors_sum[level], level + 1, len_factors_sum);
}
}