解题报告 之 SOJ2106 GCD & LCM Inverse

解题报告 之 SOJ2106 GCD & LCM Inverse

Description

We are all familiar with the Euclid's algorithm. Given two positive integers a, b. we can easily find the greatest common divisor (GCD) and least common multiple (LCM) of a, b through this method. But what about the inverse? That is: given GCD and LCM of a and b, can you find a, b ?

Input

The input will contain multiple test cases, each of which contains two positive integers: GCD and LCM. You can safely assume these integers will be less than 2^63-1.

Output

For each test case, output a, 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

LangLang @ Achilles Team

题目大意：给出你两个数的gcd和lcm，求这两个数a，b（a<=b），如果有多组解则输出a+b最小的一组。

分析：乍看这道题似乎很easy，其实是大坑。首先明确a*b==gcd*lcm。然后因为有：

a=p1^a1*p2^a2*p3^a3……

b=p1^b1*p2^b2*p3^b3……

那么可见gcd是a和b都要有的。那么lcm中剩下的每个质因子pk要么全部分给a，要么全部分给b，因为如果将质因子pk为给a一部分，分给b一部分，那么gcd势必可以继续扩大。我们求出所有分配方法更新最小的a+b即可。

所以我们大致步骤先计算ex=lcm/gcd，表示剩下可以分配的质因子的积；再将ex质因子分解；再搜索每种分配方案，更新a、b、minab即可得到答案。其中质因子分解必须采用pollard_rho算法，因为很大。

（如果还不会p_r算法的请参考另一篇博文 http://blog.youkuaiyun.com/maxichu/article/details/45459533）

因为我的分解大数时偷懒用了map，导致dfs时很不方便，大家最好注意。然后dfs注意剪枝，剪枝条件是目前的分配下都有一个因子超过当前的minab。

上代码：

#include<iostream>
#include<algorithm>
#include<cstring>
#include<vector>
#include<map>
using namespace std;
typedef long long ll;

const int times = 20;
map<ll, int> m;
int cnt = 0;
ll mina, minb, minab;
ll extra;

ll Random( ll n )
{
	return ((double)rand() / RAND_MAX*n + 0.5);
}

long long q_mul( ll a, ll b, ll c )
{
	int res = 0;
	while(b)
	{
		if(b & 1)
			res =(res+ a)%c;
		b /= 2;
		a=(a+a)%c;
	}
	return res;
}

long long q_pow( ll a, ll b, ll c )
{
	int res = 1;
	while(b)
	{
		if(b & 1)
			res = (res*a)%c;
		b /= 2;
		a=(a*a)%c;
	}
	return res;
}

bool witness( ll a, ll n )
{
	ll tem = n - 1;
	int j = 0;
	while(tem % 2 == 0)
	{
		tem /= 2;	
		j++;
	}

	ll x = q_pow( a, tem, n );
	if(x == 1 || x == n - 1)	return true;
	while(j--)
	{
		x = q_mul( x, x, n );
		if(x == n - 1) return true;
	}
	return false;
}

bool miller_rabin( ll n )
{
	if(n == 2)
		return true;
	if(n < 2 || n % 2 == 0)
		return false;

	for(int i = 1; i <= times; i++)
	{
		ll a = Random( n-2 )+1;
		if(!witness( a, n ))	return false;
	}
	return true;
}

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

ll pollard_rho( ll n, ll c )
{
	int x, y, d, i = 1, k = 2;
	x = Random( n - 1 ) + 1;
	y = x;
	while(1)
	{
		i++;
		x = (q_mul( x, x, n ) + c) % n;
		d = gcd( y - x, n );
		if(1 < d&&d < n)
		{
			return d;
		}
		if(y == x)	return n;

		if(i == k)
		{
			y = x;
			k <<= 1;
		}
	}
}

void find( ll n, ll c )
{
	if(n == 1) return;
	if(miller_rabin( n ))
	{
		m[n]++;
		cnt++;
		return;
	}
	ll p = n;
	while(p >= n)
	{
		p = pollard_rho( p, c-- );
	}
	find( p, c );
	find( n / p, c );
}

void dfs( ll dep, ll a, map<long long, int>::iterator c )
{

	if(a > minab) return; //剪枝
	if(dep == m.size())//已经分配完
	{
		ll b = extra / a;	
		if(a + b < minab)
		{
			mina = a;
			minb = b;
			minab = a + b;
		}
		return;
	}
	ll na = a*q_pow( c->first, c->second, extra+1 );
	map<long long, int>::iterator tc = c;
	tc++;
	dfs( dep + 1, na, tc );	//当前因子分配给a
	dfs( dep + 1, a, tc );		//当前因子不分配给a
	
}

int main()
{
	long long g, l;
	while(cin >> g>>l)
	{
		m.clear();
		cnt=0;
		minab = 0x7fffffffffffffff ;  // 初始化 ll 最大值
		extra = l / g;
		find( extra, 1245 );			// 质因子分解，1245是随意输的数
		dfs( 0, 1, m.begin() );		//dfs(0,1,m.begin()) 表示当前分配了0个数，其中a已经分配成为了1，现在要分配的数m.begin()
		if(mina > minb) swap( mina, minb ); //要求a<=b
		cout << (mina*g)<<" "<<(minb*g)<< endl; //别忘了乘上gcd
	}
	return 0;
}

数论加油中，，表示poyal是我的弱项啊。。。肿么办