转载自: http://blog.sina.com.cn/s/blog_69c3f0410100uac0.html
http://www.hankcs.com/program/cpp/poj-2429-gcd-lcm-inverse.html
参考资料:http://www.cnblogs.com/kuangbin/archive/2012/08/19/2646396.html
题意
给出两个整数 m、n ( m =< n < 2^63),求出两个整数 a、b ,使他们的最大公约数为 m ,最小公倍数为 n ,如果有多组,取 a + b 值最小的那一组。
题解
由题意,则令 p = a / m , q = b / m , s = n / m,则有GCD(p , q) = 1 ,LCM(p , q) = p * q = s 。
由以上结论可知此题是让我们把 s 分解成两个互质数相乘的形式。我们可以把 pollard_rho 和 miller_rabin 算法相结合将 s 分解成素因子相乘的形式,由于GCD(p , q) = 1,所以多个相同的素数不能同时分到 p 和 q 里,因此,我们多个相同的素数当作一个整体处理。又由于 s 在 2^63 之内,不会超过14个不同的素数相乘,所以可以将这些数枚举相乘,找到最接近s的平方根且不大于 s 平方根的数即为 p ,q = s / p.
例如 m = 2, n = 120,则 s = 60 = 2 * 2 * 3 * 5 = 4 * 3 * 5 ,从 3 、4 、5中挑出任意个数相乘得到 3 , 4 , 5 , 12 , 15 , 20 , 60 ,其中最靠近根号s且不大于根号s的就是 5 ,即 p =5 , q = s / p = 12 ,所以最后结果为 a = p * m = 10 , b = q * m = 24 。
代码
版本1(c语言)
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
__int64 pr[100], f[100], t, ans;
__int64 mult(__int64 a, __int64 b, __int64 n) // a*b % n
{
__int64 s = 0;
while(b)
{
if(b & 1)
s = (s + a) % n;
a <<= 1;
b >>= 1;
if(a > n)
a = a % n;
}
return s;
}
__int64 miller_rabin(__int64 n, __int64 i) //素数测定,错误机率为1/4
{
__int64 j, k, b = 0, a = 1;
j = (n - 1) / (1 << i);
while(b % n == 0)
b = abs(rand()) % n;
while(j) //模指数运算
{
if(j & 1)
a = mult(a, b, n);
b = mult(b, b, n);
j >>= 1;
}
for(k = 1; k <= i; k++)
{
if(k < i - 1 && mult(a, a, n) == 1 && a != 1 && a != n - 1)
return 0;
a = mult(a, a, n);
}
if(a == 1)
return 1;
return 0;
}
__int64 prime(__int64 n)
{
//调用miller_rabin素数测试10次,正确概率为1 - (1/4)^10,几乎为1
__int64 m, k, i, s = 1;
if(n == 1)
return 0;
for(k = 0, m = n - 1; m % 2 == 0; m /= 2, k++)
;
for(i = 0, s = 1; i < 10; i++)
s *= miller_rabin(n, k);
return s;
}
__int64 gcd(__int64 a, __int64 b) //求最大公约数
{
return b ? gcd(b, a % b) : a;
}
__int64 pollard_rho(__int64 n) //pollard_rho随机算法求n的一个因子
{
__int64 i = 1, a, b, p, t = 2, c = 0;
a = b = abs(rand()) % (n - 1) + 1;
while(c == 0 || c == 2)
c = abs(rand()) % (n - 1) + 1;
do
{
i++;
p = gcd(n + b - a, n);
if(p > 1 && p < n)
return p;
if(i == t)
{
b = a;
t <<= 1;
}
a = (mult(a, a, n) + n - c) % n;
}while(a != b);
return n;
}
void add(__int64 n)//将素因子添加到素因子表中
{
__int64 i;
for(i = 0; i < t; i++)
if(pr[i] == n)
{
break;
}
if(i < t)
f[i]++; //第i个素因子个数加1
else //如果之前没有添加
{
pr[t] = n; //添加到表尾
f[t++] = 1; //设素因子个数初始化为1
}
}
void find(__int64 n)
{
__int64 p;
if(n < 2)
return;
if(prime(n))
add(n); //如果n为素数,添加到素因子表中
else //如果n不为素数,继续拆分n
{
p = pollard_rho(n);
find(p);//尝试拆分n的一个因子p
find(n / p); //尝试拆分n的一个因子n/p
}
}
void findx(__int64 i, __int64 x, __int64 q)
{
if(i == t)
return;
if(x > ans && x <= q)
ans = x; //如果因子x比ans更接近q,则更新ans
findx(i + 1, x, q);
x *= f[i];
if(x > ans && x <= q)
ans = x; //如果因子x比ans更接近q,则更新ans
findx(i + 1, x, q);
}
int main()
{
__int64 i, j, k, m, n, q;
srand((unsigned)time(NULL));
while(~scanf("%I64d%I64d", &m, &n))
{
t = 0;
n /= m; // n = s = n/m
find(n);//分解质因数
for(i = 0; i < t; i++) //将相同素因子作为一个整体处理
{
for(j = 1, k = pr[i]; j < f[i]; j++)
k *= pr[i];
f[i] = k;
}
q = (__int64)sqrt(n * 1.0);
ans = 1;
findx(0, 1, q); //枚举n的因子,找出小于q且最接近q的因子
printf("%I64d %I64d\n", ans * m, n / ans * m);
}
return 0;
}
版本2(c++):
#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[])
{
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;
}
return 0;
}