PS:图片由作者用wps制作,使用请注明链接,O(∩_∩)O谢谢!
如对Miller-Rabin有疑问,请参考作者其他博客:Miller-Rabin素性测试算法详解
代码如下:
#include <algorithm>
#include <iostream>
#include <cstring>
#include <cstdlib>
#include <cstdio>
#include <cmath>
using namespace std;
typedef long long int ll;
const int MAX_SIZE = 1000;
const ll INF = 0x3f3f3f3f3f;
ll fac[MAX_SIZE], ct, cnt;
ll m, n; // GCD : m, LCM : n
ll GCD(ll x, ll y)
{
return !y ? x : GCD(y, x % y);
}
typedef long long int ll;
ll mod_mul(ll a, ll b, ll mod)
{
ll res = 0;
while (b)
{
if (b & 1)
res = (res + a) % mod;
a = (a + a) % mod;
b >>= 1;
}
return res;
}
ll mod_pow(ll a, ll n, ll mod)
{
ll res = 1;
while (n)
{
if (n & 1)
res = mod_mul(res, a, mod);
a = mod_mul(a, a, mod);
n >>= 1;
}
return res;
}
// Miller-Rabin随机算法检测n是否为素数
bool Miller_Rabin(ll n)
{
if (n == 2)
return true;
if (n < 2 || !(n & 1))
return false;
ll m = n - 1, k = 0;
while (!(m & 1))
{
k++;
m >>= 1;
}
for (int i = 1; i <= 20; i++) // 20为Miller-Rabin测试的迭代次数
{
ll a = rand() % (n - 1) + 1;
ll x = mod_pow(a, m, n);
ll y;
for (int j = 1; j <= k; j++)
{
y = mod_mul(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 i = 1, k = 2;
ll x = rand() % (n - 1) + 1, y = x;
while (1)
{
i++;
x = (mod_mul(x, x, n) + c) % n;
ll d = GCD((y - x + n) % n, n);
if (d != 1 && d != n) // if (d > 1 && 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) // 1不处理,因为1的质因数分解方式有无穷种
return;
if (Miller_Rabin(n))
{
fac[ct++] = n;
return;
}
ll p = n;
ll k = c;
while (p >= n)
p = Pollard_Rho(p, c--);
find(p, k);
find(n / p, k);
}
ll ans, s;
void dfs(ll cur, ll val)
{
if (val > ans)
ans = val;
for (ll i = cur; i < cnt; i++)
{
ll t = fac[i] * val;
if (t <= s)
dfs(i + 1, t);
}
}
/*
void dfs(ll cur, ll val)
{
if (val > ans)
ans = val;
if (cur == cnt)
return;
if (val * fac[cur] <= s)
dfs(cur + 1, val * fac[cur]);
dfs(cur + 1, val);
}
*/
int main()
{
//freopen("test.txt", "r", stdin);
while (~scanf("%I64d%I64d", &m, &n))
{
ct = 0;
find(n / m, 11); // 因式分解后不一定是递增的,有可能出现13,3,3的情况
sort(fac, fac + ct);
cnt = 0;
int t = fac[0];
for (int i = 1; i < ct; i++)
{
if (fac[i] == fac[i - 1])
t *= fac[i];
else
{
fac[cnt++] = t;
t = fac[i];
}
}
fac[cnt++] = t;
ans = -INF;
s = (ll)sqrt(n * 1.0 / m);
dfs(0, 1); // 代码中的两种dfs均是可行的
printf("%I64d %I64d\n", ans * m, n / ans);
}
return 0;
}