POJ 2429 -- miller-rabin素数测试，Pollard_rho素因子分解

转载自： 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;
}