题目链接:GCD & LCM Inverse
解题思路:这个给了a和b的GCD和LCM,让求出a和b,这里a和b有很多解,让输出a + b最小的一组。
gcd(a, b) = k
lcm(a, b) = m
a = pk, b = qk
a + b = (p + q)k
p * q * k = m
可以通过以上得出
p * q = m / k,gcd(p, q) = 1, 所以我们要求出q + p最小的。
我们可以知道当d是m / k的约数的时候并且gcd(d, m / k / d) = 1,那么我们这里可以看到,当p q的差距越来越小的时候p + q就越来越小,但是我们对于longlong没有sqrt,所以我们要另一种思路。
我们将m / k分解质因数,我们将相同的质因数相乘合成n个约数,这里为什么要合成,因为同一个约数既在p也在q的时候gcd(p, q) != 1,所以我们这样合成的都是两两互质的,所以搜索一遍就得出结果了。
这里涉及大整数分解,那就是Pollard rho算法,详细请见
#include<cstdio>
#include<cstring>
#include<iostream>
#include<cmath>
#include<ctime>
#include<algorithm>
#include<cstdlib>
#include<vector>
#define ll __int64
#define MAX 5500
const int S = 20;//素数测试误判率=1 / 2^S S越大误差率越小
using namespace std;
ll gcd(ll a, ll b){
return b == 0 ? a : gcd(b, a % b);
}
ll multi_mod(ll a, ll b, ll p){
ll ret = 0, q = a % p;
while(b){
if(b & 1){
ret = (ret + q) % p;
}
q = (q + q) % p;
b >>= 1;
}
return ret % p;
}
ll pow_mod(ll a, ll b, ll p){
ll ret = 1, q = a % p;
while(b){
if(b & 1){
ret = multi_mod(ret, q, p);
}
q = multi_mod(q, q, p);
b >>= 1;
}
return ret % p;
}
bool Miller_Rabin(ll p){
int i, j, k;
ll u, t, x, y;
if(p == 2) return true;
if(p % 2 == 0 || p == 1) return false;
u = p - 1, t = 0;
while(u % 2 == 0){
t++;
u >>= 1;
}
for(int i = 0; i < S; i++){
x = rand() % (p - 1) + 1;
x = pow_mod(x, u, p);
ll y = 0;
for(j = 0; j < t; j++){
y = multi_mod(x, x, p);
if(y == 1 && x != 1 && x != p - 1){
return false;
}
x = y;
}
if(y != 1) return false;
}
return true;
}
ll fact[MAX];
ll tot = 0;
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_mod(x, x, n) + c) % n;
ll d = gcd((y - x + n) % n, 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)){
fact[tot++] = n;
return;
}
ll d = n;
ll k = c;
while(d >= n){
d = pollard_rho(n, c--);
}
find(d, k);
find(n / d, k);
}
vector<ll> num;
int main(){
int i, j, k;
ll a, b, p;
ll n;
//srand(time(NULL));
while(cin >> a >> b){
ll gcd = a;
n = b / a;
tot = 0;
num.clear();
find(n, 120);
sort(fact, fact + tot);
for(i = 0; i < tot; i++){
if(i == 0 || fact[i] != fact[i - 1]){
num.push_back(fact[i]);
}
else{
num[num.size() - 1] *= fact[i];
}
}
p = n + 1;
a = 1, b = n;
ll tem, am;
for(i = 0; i < (1 << num.size()); i++){
int j = i;
tem = 1, am = 0;
while(j){
if(j & 1){
tem *= num[am];
}
j >>= 1;
am++;
}
ll tt = n / tem;
//cout << tt << " " << tem << endl;
if(tt + tem < p){
p = tt + tem;
a = min(tt, tem);
b = max(tt, tem);
}
}
cout << gcd * a << " " << gcd * b << endl;
}
return 0;
}