POJ 2429 - GCD & LCM Inverse(数论)

题目:

http://poj.org/problem?id=2429

题意:

已知两个数的最大公约数和最小公倍数,求出两个数(和最小).

思路:

因为题目的数据特别大,所以要用Rabin-Miller强伪素数测试和Pollard因数分解算法.

基本思路是:

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

论文链接http://wenku.baidu.com/view/fbbed5a5f524ccbff12184af.html 抄了代码=.=

CODE:

#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;
}
// Qualifier: Rabin-Miller强伪素数测试
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;
}
// Qualifier: Pollard 因数分解算法
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
;
else return miller_rabin(n, 20);
}

// 分解成素因子，小数用素数表，大数用Pollard 因数分解算法
void factorize(ll n, map<ll, int>& factors)
{
if (isPrime(n))
{
factors
++;
}
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
++;
}
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);
}

int main()
{
init_primes();
ll a, b;
while (cin >> a >> b)
{
pair<ll, ll> ans = solve(a, b);
cout << ans.first << " " << ans.second << endl;
}
return 0;
}