POJ 2429 GCD & LCM Inverse 题解 《挑战程序设计竞赛》

POJ 2429 GCD & LCM Inverse

逆最大公约数和最小公倍数：为了嘲讽 AOJ 0005 GCD and LCM，现在给出两个数的gcd和lcm，求原来的这两个数（限定两数之和最小）。

2.6 数学问题的解题窍门

作者显然高估了读者的数学修养，我对数论一点都不熟悉，这题光靠前面介绍的一点数论皮毛无从下手，还是看了人家的代码才知道要用Rabin-Miller强伪素数测试和Pollard r因数分解算法。这两个算法的详解放在最后面，不时回来复习一下。

基本思路是 

• (a / gcd) * (b / gcd) = lcm / gcd ，所以需要分解lcm / gcd 。将其分解为互质的两个数，如果这两个数之和最小，那么乘上gcd就是所求的答案。

• 但是题目数据太大，需要一个高效的素数检测算法，所以采用Rabin-Miller强伪素数测试

• 然后分解成质因子的n次方之积，从这些n次方中挑选一些作为x，剩下的作为y。枚举x和y的所有可能，找出最小值。

#ifndef ONLINE_JUDGE
#pragma warning(disable : 4996)
#endif
#include <iostream>
#include <vector>
#include <map>
#include <cstdlib>
using namespace std;

typedef long long LL;

// return (a * b) % m
LL mod_mult(LL a, LL b, LL m) 
{
	LL res = 0;
	LL exp = a % m;
	while (b) 
	{
		if (b & 1) 
		{
			res += exp;
			if (res > m) res -= m;
		}
		exp <<= 1;
		if (exp > m) exp -= m;
		b >>= 1;
	}
	return res;
}

// return (a ^ b) % m
LL mod_exp(LL a, LL b, LL m) {
	LL res = 1;
	LL exp = a % m;
	while (b) 
	{
		if (b & 1) res = mod_mult(res, exp, m);
		exp = mod_mult(exp, exp, m);
		b >>= 1;
	}
	return res;
}

// 
//************************************
// Method:    miller_rabin
// FullName:  miller_rabin
// Access:    public 
// Returns:   bool 是否是素数
// Qualifier: Rabin-Miller强伪素数测试
// Parameter: LL n 待检测的数
// Parameter: LL times 
//************************************
bool miller_rabin(LL n, LL times) 
{
	if (n < 2) return false;
	if (n == 2) return true;
	if (!(n & 1)) return false;

	LL q = n - 1;
	int k = 0;
	while (q % 2 == 0) {
		k++;
		q >>= 1;
	}
	// n - 1 = 2^k * q (q是奇素数)
	// n是素数的话，一定满足下面条件
	// (i) a^q ≡ 1 (mod n)
	// (ii) a^q, a^2q,..., a^(k-1)q 中的某一个对n求模为-1
	//
	// 所以、当不满足(i)(ii)中的任何一个的时候，就有3/4的概率是合成数
	//
	for (int i = 0; i < times; ++i) 
	{
		LL a = rand() % (n - 1) + 1; // 从1,..,n-1随机挑一个数
		LL x = mod_exp(a, q, n);
		// 检查条件(i)
		if (x == 1) continue;
		// 检查条件(ii)
		bool found = false;
		for (int j = 0; j < k; j++) 
		{
			if (x == n - 1) 
			{
				found = true;
				break;
			}
			x = mod_mult(x, x, n);
		}
		if (found) continue;

		return false;
	}
	return true;
}

LL get_gcd(LL n, LL m) 
{
	if (n < m) swap(n, m);
	while (m) 
	{
		swap(n, m);
		m %= n;
	}
	return n;
}


//************************************
// Method:    pollard_rho
// FullName:  pollard_rho
// Access:    public 
// Returns:   LL
// Qualifier: Pollard 因数分解算法
// Parameter: LL n
// Parameter: int c
//************************************
LL pollard_rho(LL n, int c) 
{
	LL x = 2;
	LL y = 2;
	LL d = 1;
	while (d == 1) 
	{
		x = mod_mult(x, x, n) + c;
		y = mod_mult(y, y, n) + c;
		y = mod_mult(y, y, n) + c;
		d = get_gcd((x - y >= 0 ? x - y : y - x), n);
	}
	if (d == n) return pollard_rho(n, c + 1);
	return d;
}

#define MAX_PRIME 200000
vector<int> primes;
vector<bool> is_prime;

// 先生成MAX_PRIME内的素数表
void init_primes() 
{
	is_prime = vector<bool>(MAX_PRIME + 1, true);
	is_prime[0] = is_prime[1] = false;
	for (int i = 2; i <= MAX_PRIME; ++i) 
	{
		if (is_prime[i]) 
		{
			primes.push_back(i);
			for (int j = i * 2; j <= MAX_PRIME; j += i) 
			{
				is_prime[j] = false;
			}
		}
	}
}

// 判断是否是素数，优先查表，如果n很大用Rabin-Miller强伪素数测试
bool isPrime(LL n) 
{
	if (n <= MAX_PRIME) return is_prime[n];
	else return miller_rabin(n, 20);
}

// 分解成素因子，小数用素数表，大数用Pollard 因数分解算法
void factorize(LL n, map<LL, int>& factors) 
{
	if (isPrime(n)) 
	{
		factors[n]++;
	}
	else 
	{
		for (int i = 0; i < primes.size(); ++i) 
		{
			int p = primes[i];
			while (n % p == 0) 
			{
				factors[p]++;
				n /= p;
			}
		}
		if (n != 1) 
		{
			if (isPrime(n)) 
			{
				factors[n]++;
			}
			else 
			{
				LL d = pollard_rho(n, 1);
				factorize(d, factors);
				factorize(n / d, factors);
			}
		}
	}
}

pair<LL, LL> solve(LL a, LL b) 
{
	LL c = b / a;
	map<LL, int> factors;
	factorize(c, factors);

	vector<LL> mult_factors;	// 每个质因子的n次方，比如2^2和5^1
	for (map<LL, int>::iterator it = factors.begin(); it != factors.end(); it++) 
	{
		LL mul = 1;
		while (it->second) 
		{
			mul *= it->first;
			it->second--;
		}
		mult_factors.push_back(mul);
	}

	LL best_sum = 1e18, best_x = 1, best_y = c;
	// 这是一个取数的过程，一共 2^size 种情况
	for (int mask = 0; mask < (1 << mult_factors.size()); ++mask)
	{
		LL x = 1;
		for (int i = 0; i < mult_factors.size(); ++i)
		{
			if (mask & (1 << i)) x *= mult_factors[i];
		}
		LL y = c / x;
		if (x < y && x + y < best_sum) 
		{
			best_sum = x + y;
			best_x = x;
			best_y = y;
		}
	}
	return make_pair(best_x * a, best_y * a);
}

///////////////////////////SubMain//////////////////////////////////
int main(int argc, char *argv[])
{
#ifndef ONLINE_JUDGE
	freopen("in.txt", "r", stdin);
	freopen("out.txt", "w", stdout);
#endif
	cin.tie(0);
	ios::sync_with_stdio(false);

	init_primes();
	LL a, b;
	while (cin >> a >> b) 
	{
		pair<LL, LL> ans = solve(a, b);
		cout << ans.first << " " << ans.second << endl;
	}
#ifndef ONLINE_JUDGE
	fclose(stdin);
	fclose(stdout);
	system("out.txt");
#endif
	return 0;
}
///////////////////////////End Sub//////////////////////////////////

关于两个算法的讲解，有一个超棒的文档。

 知识共享署名-非商业性使用-相同方式共享：码农场 » POJ 2429 GCD & LCM Inverse 题解 《挑战程序设计竞赛》