当前位置: 代码迷 >> 综合 >> POJ 2429 GCD LCM Inverse Pollard_rho大数因子分解
  详细解决方案

POJ 2429 GCD LCM Inverse Pollard_rho大数因子分解

热度:83   发布时间:2024-01-13 17:27:02.0

这题的话。假设原来的数是a和b,

给出GCD和LCM以后

我们只需要求一下

LCM/GCD 的因子  就可以求出 a/gcd 和 b/gcd 

然后求出的这两个东西之间不能有公约数

所以一个数取因子的时候要把某个因子全拿走

用的是Pollard_rho进行分解这种比较大的数。

复杂度的话,假设分解出一个因子是p,用的时间是O(sqrt(p))

总体来说比你直接sqrt(n)分解快多了


#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <algorithm>
#include <ctime>
#include <map>
#define MAXN 10
#define INF 1LL<<61
#define C 16381
using namespace std;
typedef long long LL;
LL gcd(LL a, LL b)
{return b ? gcd(b, a % b) : a;
}LL multi(LL a, LL b, LL m)
{LL ans = 0;while(b){if(b & 1){ans = (ans + a) % m;b--;}b >>= 1;a = (a + a) % m;}return ans;
}LL quick_mod(LL a, LL b, LL m)
{LL ans = 1;a %= m;while(b){if(b & 1){ans = multi(ans, a, m);b--;}b >>= 1;a = multi(a, a, m);}return ans;
}bool Miller_Rabin(LL n)
{if(n == 2) return true;if(n < 2 || !(n & 1)) return false;LL a, m = n - 1, x, y;int k = 0;while((m & 1) == 0){k++;m >>= 1;}for(int i = 0; i < MAXN; i++){a = rand() % (n - 1) + 1;x = quick_mod(a, m, n);for(int j = 0; j < k; j++){y = multi(x, x, n);if(y == 1 && x != 1 && x != n - 1) return false;x = y;}if(y != 1) return false;}return true;
}LL Pollard_rho(LL n, LL c)
{LL x, y, d, i = 1, k = 2;y = x = rand() % (n - 1) + 1;while(true){i++;x = (multi(x, x, n) + c) % n;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;}}
}
LL fac[555];
int cnt;
map<LL, int>mp;
void find(LL n, int c)
{if(n == 1) return;if(Miller_Rabin(n)){mp[n]++;return ;}LL p = n;LL k = c;while(p >= n) p = Pollard_rho(p, c--);find(p, k);find(n / p,k);
}int main()
{LL a, b;while(scanf("%I64d%I64d", &a, &b) != EOF){cnt = 0;LL tmp = b / a;if(a == b){printf("%I64d %I64d\n", a, b);continue;}mp.clear();find(tmp, C);cnt = 0;for(map<LL, int>:: iterator it = mp.begin(); it != mp.end(); it++){LL t = 1;for(int j = 0; j < it -> second; j++)t = t * (it -> first);fac[cnt++] = t;}//for(int i = 0; i < cnt; i++)// printf("%I64d\n", fac[i]);LL ansa = INF, ansb = INF;for(int i = 0; i < (1 << cnt); i++){LL ta = 1;for(int j = 0; j < cnt; j++)if(i & (1 << j))ta *= fac[j];if(ansa + ansb > ta + b / a / ta)ansa = ta, ansb = b / a / ta;}if(ansa > ansb) swap(ansa, ansb);printf("%I64d %I64d\n", ansa * a, ansb * a);}return 0;
}