POJ 2429 GCD & LCM Inverse （挑战程序设计竞赛，Miller_Rabin算法， pollard_rho 算法）

挑战程序设计竞赛， 135页，素数测试， 大数因数分解
题目意思：
给出两个数 ， gcd 和 lcm 分别表示 a和b的最大公约数和最小公倍数。
求a和b. 如果有多组，输出 a + b 最小的一组。

本题要点：
1、Miller_Rabin算法：
https://www.cnblogs.com/fzl194/p/9046117.html
2、pollard_rho 算法：
https://www.cnblogs.com/fzl194/p/9047710.html
3、 对 b / a 做因式分解后，在有 dfs 搜出若干个因子的总乘积， 找出答案。
4、 a == b时候，做特判， 否则会 RE

#include <cstdio>
#include <cstring>
#include <ctime>
#include <iostream>
#include <cstdlib>
#include <algorithm>
using namespace std;
typedef long long ll;
const int repeat = 20, MaxN = 5500;
ll fac[MaxN], num[MaxN], fac_all[MaxN];
ll ct, cnt, n;

ll gcd(ll a, ll b)
{
	return b? gcd(b, a % b) : a;
}

ll multi(ll a, ll b, ll m)	//求 a*b % m
{
	ll ans = 0;
	a %= m;
	while(b)
	{
		if(b & 1)
			ans = (ans + a) % m;
		b /= 2;
		a = (a + a) % m;
	}
	return ans;
}

ll pow(ll a, ll b, ll m)	//求 a^b % m
{
	ll ans = 1;
	a %= m;
	while(b)
	{
		if(b & 1)
			ans = multi(a, ans, m);
		b /= 2;
		a = multi(a, a, m);
	}
	ans %= m;
	return ans;
}

bool Miller_Rabin(ll n)
{
	if(2 == n || 3 == n)
		return true;
	if(n % 2 == 0 || 1 == n)	//偶数和1
		return false;
	ll d = n - 1;	// n - 1 分解为 2^s * d
	int s = 0;
	while(!(d & 1))
	{
		++s, d >>= 1;
	}
	//srand((unsigned)time(NULL));
	for(int i = 0; i < repeat; ++i)	//重复 repeat 次
	{
		ll a = rand() % (n - 3) + 2;	//选取一个随机数 [2, n - 1)
		ll x = pow(a, d, n);
		ll y = 0;
		for(int j = 0; j < s; ++j)
		{
			y = multi(x, x, n);
			if(1 == y && x != 1 && x != n - 1)
				return false;
			x = y;
		}
		if(y != 1)	//费马定理， 不满足费马定理的，就不是素数
			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 = gcd((y - x + n) % n, n);
		if(1 < d && d < n)
			return d;	// d 是n的因子
		if(x == y)
			return n;	//找到循环， 返回 n，重新来
		if(i == k)
		{
			y = x;
			k <<= 1;
		}
	}
}

void find(ll n, int c)
{	
	if(1 == n)
		return;
	if(Miller_Rabin(n))
	{
		fac[ct++] = n;
		return;
	}
	ll p = n, k = c;
	while(p >= n)
		p = pollard_rho(p, c--);
	find(p, k);
	find(n / p, k);
}

ll ans_a, ans_b, min_sum;

void dfs(int index, long long now)	//now 表示当前的额乘积
{
	if(index == cnt)
	{
		if(min_sum > now + n / now)
		{
			ans_a = now, ans_b = n / now;	
			min_sum = now + n / now;
		}
		return;
	}
	dfs(index + 1, now);
	dfs(index + 1, now * fac_all[index]);
}

int main()
{
	ll a, b;
	while(scanf("%lld%lld", &a, &b) != EOF)
	{
		if(a == b)
		{	
			printf("%lld %lld\n", a, b);
			continue;
		}
		n = b / a;
		ct = 0;
		find(n, 20);
		sort(fac, fac + ct);
		num[0] = 1;
		int k = 1; 
		for(int i = 1; i < ct; ++i)
		{
			if(fac[i] == fac[i - 1])
			{
				++num[k - 1];
			}else{
				num[k] = 1;
				fac[k++] = fac[i];
			}
		}
		cnt = k;
		for(int i = 0; i < cnt; ++i)
		{
			fac_all[i] = 1;
			for(int j = 0; j < num[i]; ++j)
			{
				fac_all[i] *= fac[i];
			}
		}
		min_sum = b / a * 2;
		dfs(0, 1);
		if(ans_a > ans_b)
		{
			swap(ans_a, ans_b);
		}
		printf("%lld %lld\n", ans_a * a, ans_b * a);
	}
	return 0;
}

/*
3 60
*/

/*
12 15
*/