poj2429(miller_robin算法和pollard分解质因数)

来源:互联网 发布:淘宝店铺数据查询 编辑:程序博客网 时间:2024/05/22 12:27
/*translation:给出两个数的gcd,lcm的值,求两个数分别是多少?有多解时输出两数和最小的一组。solution:miller-robin和pollard素因子分解算法这道题可以参考poj1811(模板),根据规律a*b == gcd(a,b)*lcm(a,b)很容易想出分解质因数然后凑出a,b来暴力比较求出答案。但是由于数据范围太大。所以常规的素数分解算法不行。套上pollard算法的模板即可。然后dfs求出最靠近sqrt(y)但是不超过y的数字,这个数字就是所求的第一个数。注意一开始x,y相乘可能溢出,根据x,y的特性,可以用分解y/x的质因数来代替。因为y*x和y/x只有次数不同。note:# 注意原来pollard模板不能对1分解质因数,会RE。在原来的基础上在一开始加上了对n==1的情况的处理即可date:2016.10.26*/#include <iostream>#include <cstring>#include <cstdio>#include <cmath>#include <cstdlib>#include <vector>#include <ctime>#include <utility>using namespace std;const int max_prime = 1e6;const int s = 20;typedef long long ll;vector<ll> factor;vector<ll> prime_factor;ll x, y;//x:gcd, y:lcmll multi_mod(ll a, ll b, ll n){    a %= n;    b %= n;    ll res = 0;    while(b)    {        if(b & 1)        {            res += a;            if(res >= n) res -= n;        }        a <<= 1;        if(a >= n) a -= n;        b >>= 1;    }    return res;}ll pow_mod(ll x, ll n, ll mod){    if(n == 1)  return x % mod;    x %= mod;    ll tmp = x;    ll res = 1;    while(n > 0)    {        if(n & 1)   res = multi_mod(res, tmp, mod);        tmp = multi_mod(tmp, tmp, mod);        n >>= 1;    }    return res;}bool witness(ll a, ll n, ll u, ll t){    ll res = pow_mod(a, u, n);    ll last = res;    for(int i = 0; i < t; i++)    {        res = multi_mod(res, res, n);        if(res == 1 && last != 1 && last != n-1)    return true;        last = res;    }    if(res != 1)    return true;    else        return false;}bool robin_miller(ll n) //判断是否素数{    if(n < 2)    return false;    if(n == 2)  return true;    if(!(n & 1))    return false;    ll u = n - 1, t = 0;    while(!(u & 1)) u >>= 1, t++;    if(t >= 1 && (u & 1) == 1)    {        for(int i = 0; i < s; i++)        {            ll a = rand() % (n-1) + 1;            if(witness(a, n, u, t)) return false;   //不是素数        }    }    return true;    //是素数}ll gcd(ll a, ll b){    if(a == 0)  return 1;    if(a < 0)    return gcd(-a, b);    while(b)    {        ll t = a % b;        a = b;        b = t;    }    return a;}ll pollard_rho(ll x, ll c){    ll i = 1, x0 = rand() % x;    ll y = x0;    ll k = 2;    for(;;)    {        i++;        x0 = (multi_mod(x0, x0, x) + c) % x;        ll d = gcd(y - x0, x);  //这里传给gcd的参数可能负数,注意分别讨论        if(d != 1 && d != x) return d;  //这里表明找到了一个因子        if(y == x0) return x;        if(i == k)        {            y = x0;            k += k;        }    }}void find_fac(ll n){    if(n == 1)    return;    if(robin_miller(n))    {        prime_factor.push_back(n);        return;    }    ll p = n;    while(p >= n) p = pollard_rho(p, rand() % (n-1) + 1);    find_fac(p);    find_fac(n / p);}ll dfs(int p, int len, ll res, ll m){if(res > m)return 0;if(p == len || res == m)return res;ll tmp1 = dfs(p+1, len, res, m), tmp2;if(res * factor[p] + 0.0 > m)return tmp1;else{tmp2 = dfs(p+1, len, res*factor[p], m);return tmp1 > tmp2 ? tmp1 : tmp2;}}int main(){//freopen("in.txt", "r", stdin);while(cin >> x >> y){srand(time(NULL));prime_factor.clear();factor.clear();y /= x;find_fac(y);ll tmp = y;for(int i = 0; i < prime_factor.size(); i++){ll res = 1;while(tmp % prime_factor[i] == 0)res *= prime_factor[i], tmp /= prime_factor[i];if(res != 1)factor.push_back(res);}int len = factor.size();ll m = (ll)sqrt(y + 0.0);ll a = dfs(0, len, 1, m);ll b = y / a;printf("%lld %lld\n", a * x, b * x);}return 0;}

0 0
原创粉丝点击