挑战程序设计解题报告 2.6.1辗转相除法

挑战程序设计解题报告 2.6.1辗转相除法

1.POJ 2429

        这道理的意思就是告诉你GCD和LCM让你求原数，但是有个问题，很明显有很多情况，所以题目要求输出满足a+b最小的情况，怎样才能最小呢，已知GCD(a, b) == n, lcm(a, b) == m,那么a * b == n / m（由GCD和LCM的定义可知），问题变成了找到两个数a,b使得a * b == n / m、且a+b最小，即f(x) = a * b = a  +  (n / m) / a，这个函数是对勾函数，最小值是a等于分子的平方根，如图：

 

        但是盲目去找是不行的，所以我们需要将n/m分解成素因子，由于同一个素因子可能有很多，所以需要先合并在一起，然后从这些数字及其乘积的组合中去寻找最接近上图最低值的点，首先是大整数的分解需要用到Pollard_rho，然后是大素数判断Miller_Rabin，这些都是模板，假设我们现在得到数字以后，只需要使用挑战程序设计这本书的31页的穷遏搜索即可，也就是从一堆数字中选几个数字乘在一起使得它的值接近给定值。

        例如2 120这组数据，分解形成 2 2 3 5，合并形成3 4 5，搜索的时候检查3 4 5 12 15 20 60最终发现5更接近sqrt（120 / 2），所以结果就是 5 和 12 即为 10 24

代码如下：

const int S = 10;    //米勒罗宾素数测试

map<LL, LL>v;    //记录并合并多余的素因子
vector<LL>ans;    //记录最终需要进行搜索的数字
LL mult_mod(LL a, LL b, LL c)    //计算a * b % c
{
    a %= c;
    b %= c;
    LL ret = 0;
    while(b)
    {
        if(b & 1)
        {
            ret += a;
            ret %= c;
        }
        a <<= 1;
        if(a >= c)a %= c;
        b >>= 1;
    }
    return ret;
}

LL pow_mod(LL x, LL n, LL mod)  //计算a^b%c
{
    if(n == 1)
        return x % mod;
    x %= mod;
    LL tmp = x;
    LL ret = 1;
    while(n)
    {
        if(n & 1) ret = mult_mod(ret, tmp, mod);
        tmp = mult_mod(tmp, tmp, mod);
        n >>= 1;
    }
    return ret;
}

bool Miller_Rabin(LL n)     //模板 米勒罗宾大素数测试
{
    if(n == 2)
        return true;
    if(n < 2 || !(n & 1))
        return false;
    int t = 0;
    const int S = 10;

    LL a, x, y, u = n - 1;
    while((u & 1) == 0)
    {
        t++;
        u >>= 1;
    }
    for(int i = 0; i < S; i++)
    {
        a = rand() % (n - 1) + 1;
        x = pow_mod(a, u, n);
        for(int j = 0; j < t; j++)
        {
            y = mult_mod(x, x, n);
            if(y == 1 && x != 1 && x != n - 1)
                return false;
            x = y;
        }
        if(x != 1)
            return false;
    }
    return true;
}

LL gcd(LL a, LL b)    //gcd
{
    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)    //得到其中一个因子，配合findfac一起找到所有的因子
{
    LL i = 1, k = 2;
    LL x0 = rand() % x;
    LL y = x0;
    while(1)
    {
        i++;
        x0 = (mult_mod(x0, x0, x) + c) % x;
        LL d = gcd(y - x0, x);
        if(d != 1 && d != x) return d;
        if(y == x0) return x;
        if(i == k)
        {
            y = x0;
            k += k;
        }
    }
}

void findfac(LL n)   
{
    if(n == 1)
        return;
    if(Miller_Rabin(n))
    {
        if(v[n] == 0)    //记录并合并素因子
            v[n] = n;
        else
            v[n] *= n;
        return;
    }
    LL p = n;
    while(p >= n)
        p = Pollard_rho(p, rand() % (n - 1) + 1);
    findfac(p);    //利用二分的思想找其余的因子
    findfac(n / p);    //
}

void bfs(LL i, LL sum, LL& key, LL p)
{
    if(abs(key - p) > abs(sum - p))    //用key保存最接近p的值
        key = sum;
    if(i == ans.size())
        return;     //递归出口
    bfs(i + 1, sum, key, p);    //不选i+1这个数字 
    bfs(i + 1, sum * ans[i], key, p);    //选上i+1这个数字
}

int main()
{
    LL N, M;
    while(cin >> N >> M)
    {
        v.clear();
        ans.clear();
        LL m = M / N;
        findfac(m);   //找到因子
        for(map<LL, LL>::iterator it = v.begin(); it != v.end(); ++it)
        {
            ans.push_back(it->second);    //提取到vector中
        }
        LL a = 1;
        bfs(0, 1, a, sqrt(m));    //找到最接近M / N的值
        LL b = m / a;
        if(a < b)
            cout << a * N << " " << b * N << endl;
        else
            cout << b * N << " " << a * N << endl;
    }
    return 0;
}

        这道题为了解决素数检测，因式分解。花了很多时间，虽然是模板但是这两个算法的说明，我将会在另一篇文章中指出。