POJ 2429 GCD & LCM Inverse Pollard Rho大数分解法

GCD & LCM Inverse

Time Limit: 2000MS		Memory Limit: 65536K
Total Submissions: 17281		Accepted: 3172

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

Source

POJ Achilles

题意很简单，就是告诉你两个数a, b的GCD和LCM，反过来让你求这两个数，因为可能有多组答案，所以你还要保证这两个数的和最小才行。
开始看题，这不就是个水题吗？一看数据范围，？？？
因为数据很大，而且涉及到取模，我们先把快速乘法取模和快速幂取模写好备用，都涉及到了二进制的思想，不清楚的可以看看网上大牛的博客。
因为我们要求a和b，而a*b / gcd = lcm, 两边同除gcd, 就得到了(a / gcd)*(b / gcd) = lcm / gcd，不难理解，a和b一定是互质的（虽然a和b不一定是质数），这个问题就转化成了求lcm / gcd的分解，并从它们里面找到a + b和最小并且互质的数。
思想讲完了，现在讲讲里面要用的算法。（搞懂这两个玄学算法我看了好久orz）


首先是Miller - Rabin素数测试算法，这里面涉及到这几个知识点：
费马小定理：假如p是质数，且(a, p) = 1，那么a ^ (p - 1)≡1(mod p)。即假如p是质数，且a, p互质，那么a的(p - 1)次方除以p的余数恒等于1。
Miller - Rabin测试：不断选取不超过n - 1的数，计算是否每次都有b ^ (n - 1) ≡ 1(mod n)，若每次都成立则n是素数，否则为合数。
二次探测定理：如果p是奇素数，对于0<x<p，若x ^ 2 mod p = 1 ，可推出 x = 1或p - 1。

假设我们要测试的数为n，在用这个方法的第一步，我们要特判一下这个数是不是1, 2和2的倍数，加快一些速度。然后保存一下n - 1（上面的知识点有的是在测试n - 1来判断n）。为了用二次探测定理，我们要先判断n - 1的二进制末尾有几个零（一个零就代表这个数可以除以一个2，放在幂上来说一个2就是一个平方，二次探测定理就是要用到x的平方的）让n - 1不断右移并保存一下0的个数，直到末尾为1为止。这里就相当于把n - 1转换了一下，变成了n - 1 = x*（2^t）, t代表零的个数，然后再用到Miller - Rabin测试，用随机数函数产生一个不超过n - 1的数a，先计算(a^x)mod n, 然后再用刚才数的零的个数把剩下的2补上，补2这个过程就相当于在联合二次探测定理判断。补全了后再判断最后算出来的结果。判断是不是素数具体在代码中实现，相信看完这个后不难理解。
还有就是Miller - Rabin测试是概率型的测试，并不是完全正确的，但一次成功的概率大概是四分之三，所以说在进行这个测试的时候我们要多产生几个随机数多测几次来提高准确率，这个题的话十多次就行了。
如果对这个方法还有所不懂，可以参见一下这位大佬的笔记https ://www.cnblogs.com/vongang/archive/2012/03/15/2398626.html

然后是Pollard Rho大数分解法。我感觉这个方法很玄学，涉及到概率论emmm，我可能讲不清为什么，如果你们想看看的话可以去这https ://wenku.baidu.com/view/3db5c7a6ad51f01dc381f156.html
首先我们先产生一个1到n的随机数x1，然后...再产生一个随机数x2，第二个随机数是用第一个算出来的伪随机数 :

x2 = （x1*x1 + c) % n。这个c随便选，主要是保证不同，还有这个x1换成2*x1, 3 * x1都行，主要是模拟一个随机过程。之后算x1 - x2的差值，如果差值小于0就加模数取模。用算出来这个不知道是多少的数（黑人问号）去与要拆分的那个数n求gcd，如果1<gcd&&gcd<n，那么这个数就被成功的随机拆成了两部分，gcd和n / gcd，然后用拆出来的这两个数继续拆；如果不成功就换个c值。这个方法在产生伪随机数多组过后可能形成一个环，又回到了原点，导致无限循环。这时我们要用一个方法，可以看做用这个数套用刚刚的伪随机数生产函数连续生产两次：相当于两个人在圆形操场一起跑步，一个人普通速度，一个人两倍速，当一定时间过后速度快的就会绕一圈追上那个速度慢的。这里也一样，如果会成环，那么在某个时刻产生的两个数会相等，这时直接退出，换个c的值继续拆分。

方法讲完了，现在进入正题。我们的目标是要拆lcm / gcd，在这个拆的过程中判断一下拆出来的数是不是素数，是素数的话就保存下来，不是就继续拆。最后对拆出来的因子进行dfs找出a + b最小的组合。注意当输入的gcd和lcm相等的时候直接输出。

萌新第一次写，有不懂或者有错误请提出！

代码：

#include<iostream>
#include<cstring>
#include<algorithm>
#include<math.h>
#include<cstdio>
#include<vector>
#include<queue>
#include<set>
#include<map>
#include<string>
#include<stack>
using namespace std;
typedef long long int LL;
const LL inf = (LL)1 << 61;
const LL maxn = 600;
LL fac[maxn], rfac[maxn];//fac保存因子，rfac保存相同的因子有多少
LL ansa, ansb, gcd, lcm, cnt, rcnt, minn, ab;
LL multi(LL a, LL b, LL mod)
{
	LL ans = 0;
	a %= mod;
	while (b)
	{
		if (b & 1)
			ans = (ans + a) % mod;
		a = (a + a) % mod;
		b >>= 1;
	}
	return ans;
}
LL qmod(LL a, LL b, LL mod)
{
	LL ans = 1;
	a %= mod;
	while (b)
	{
		if (b & 1)
			ans = multi(ans,a,mod);
		a = multi(a, a, mod);
		b >>= 1;
	}
	return ans;
}
LL fgcd(LL a, LL b)
{
	if (b == 0)
		return a;
	else
		return fgcd(b, a%b);
}
bool miller(LL n)
{
	if (n == 2)
		return true;
	if (n < 2 || !(n & 1))
		return false;
	LL m = n - 1;
	LL c2 = 0;
	while (!(m & 1))
	{
		c2++;
		m >>= 1;
	}
	for (int i = 0; i < 14; i++)//测试次数，十多次就行了，时间短的话可以少选几次
	{
		LL a = rand() % (n - 1) + 1;
		if (a == n)
			continue;
		LL x = qmod(a, m, n);
		LL y = 0;
		for (int j = 0; j < c2; j++)
		{
			y = multi(x, x, n);
			if (y == 1 && x != 1 && x != n - 1)  //如果p是奇素数，对于0<x<p，若x^2 mod p =1 ，可推出 x=1或p-1。
				return false;
			x = y;
		}	
		if (y != 1)//y为模出来的最后结果
			return false;
	}
	return true;
}
LL pollard_rho(LL n, LL c)
{
	LL i = 1, k = 2;
	LL x = rand() % (n - 1) + 1;
	LL y=x;
	while (true)
	{
		i++;
		x = (multi(x, x, n) + c) % n;
		LL d = fgcd((y - x + n) % n, n);
		if (1 < d&&d < n)
			return d;
		if (x == y)
			return n;
		if (i == k)//两倍速操作
		{
			y = x;
			k <<= 1;
		}
	}
}

void find(LL n, LL c)
{
	if (n == 1)
		return;
	if (miller(n))
	{
		fac[cnt++] = n;
		return;
	}
	LL p = n;
	LL k = c;
	while (p >= n)
		p = pollard_rho(p, c--);
	find(p, k);
	find(n / p, k);
}
void dfs(LL deeep, LL num)
{
	if (rcnt == deeep)
	{
		LL a = num;
		LL b = ab / a;
		if (fgcd(a, b) == 1)
		{
			a *= gcd;
			b *= gcd;
			if (a + b < minn)
			{
				minn = a + b;
				ansa = a;
				ansb = b;
			}
		}
		return;
	}
	for (int i = 0; i <= rfac[deeep]; i++)
	{
		if (num > minn)
			return;
		dfs(deeep + 1, num);
		num *= fac[deeep];
	}
}
int main()
{
	while (~scanf("%lld %lld", &gcd, &lcm))
	{
		if (gcd == lcm)
		{
			printf("%lld %lld\n", gcd, lcm);
			continue;
		}
		minn = inf;
		ab = lcm / gcd;
		cnt = 0;
		find(ab, 120);
		sort(fac, fac + cnt);
		rfac[0] = 1;
		int k = 1;
		for (int i = 1; i<cnt; i++)
		{
			if (fac[i] == fac[i - 1])
				rfac[k - 1]++;
			else
			{
				rfac[k] = 1;
				fac[k++] = fac[i];
			}
		}
		rcnt = k;
		dfs(0, 1);
		if (ansa>ansb)
			swap(ansa, ansb);
		printf("%lld %lld\n", ansa, ansb);
	}
	return 0;
}