POJ2429_GCD&LCM Inverse【Miller Rabin素数テスト】【Pollar Rho整数分解】

GCD & LCM Inverse
Time Limit: 2000MS Memory Limit: 65536K
Total Submissions: 9756Accepted: 1819
Description Given two positive integers a and b, we can easily calculate the greatest common divisor (GCD) and the least common multiple (LCM) of a and b. But what about the inverse? That is: given GCD and LCM, finding a and b. Input The input contains multiple test cases, each of which contains two positive integers, the GCD and the LCM. You can assume that these two numbers are both less than 2^63. Output For each test case, output a and b in ascending order. If there are multiple solutions, output the pair with smallest a + b. Sample Input 3 60 Sample Output 12 15 Source
POJ Achilles
あなたに2つの数aとbの最大公約数と最小公倍数をあげて、aとbを求めます
(ただし条件を満たす場合はa+bをできるだけ小さくする)
構想:最大公約数と最小公倍数の規模は2^63で、暴力は断固としてだめだ.
a*b=L(最小公倍数)*G(最大公約数);
p=L/a、q=L/b、s=L/Gとする.
すなわち、p、qは、aおよびbが最大公約数を除去する部分であり、両者が互いに質的である.
GCD(p,q) = 1,LCM(p,q) = p * q = L*L/(a*b) = L*L/(L*G) = L/G = s.
LCM(p,q) = s;
上から得ることができて、私たちはsからaとbを求めることができます.この問題はsを2つの互質数に乗算する形式に分解することである.
Pollar Rho整数分解とMiller Rabin素数試験を組み合わせてsの全質量因子を分解した.
GCD(p,q)=1のため,すべての同じ素数をpとqに同時に分けることができず,同じ素数分を開放すべきである.
ここではすべての同じ質量数を全体としています.これらの数を列挙して乗算し、sに最も近い平方根を見つけ、
sより大きい平方根の組み合わせがpであればq=s/p
最終a=L/p,b=L/q
例えばG Lは2 120
s = 60 = 2 * 2 * 3 * 5 = 4 * 3 * 5.列挙して組み合わせて、最も近いルート番号60がルート番号60を超えないものを見つけます.
値は5、すなわちp=5であり、q=60/5=12である.最終的にはa＝24,b＝10であった.
参考文献:http://blog.sina.com.cn/s/blog_69c3f0410100uac0.html

#include<stdio.h>
#include<stdlib.h>
#include<time.h>
#include<math.h>
#include<algorithm>
using namespace std;
#define MAX_VAL (pow(2.0,60))
//miller_rabbin 
//__int64 mod_mul(__int64 x,__int64 y,__int64 mo)
//{
//    __int64 t;
//    x %= mo;
//    for(t = 0; y; x = (x<<1)%mo,y>>=1)
//        if(y & 1)
//            t = (t+x) %mo;
//
//    return t;
//}

__int64 mod_mul(__int64 x,__int64 y,__int64 mo)
{
    __int64 t,T,a,b,c,d,e,f,g,h,v,ans;
    T = (__int64)(sqrt(double(mo)+0.5));

    t = T*T - mo;
    a = x / T;
    b = x % T;
    c = y / T;
    d = y % T;
    e = a*c / T;
    f = a*c % T;
    v = ((a*d+b*c)%mo + e*t) % mo;
    g = v / T;
    h = v % T;
    ans = (((f+g)*t%mo + b*d)% mo + h*T)%mo;
    while(ans < 0)
        ans += mo;
    return ans;
}

__int64 mod_exp(__int64 num,__int64 t,__int64 mo)
{
    __int64 ret = 1, temp = num % mo;
    for(; t; t >>=1,temp=mod_mul(temp,temp,mo))
        if(t & 1)
            ret = mod_mul(ret,temp,mo);

    return ret;
}

bool miller_rabbin(__int64 n)
{
    if(n == 2)
        return true;
    if(n < 2 || !(n&1))
        return false;
    int t = 0;
    __int64 a,x,y,u = n-1;
    while((u & 1) == 0)
    {
        t++;
        u >>= 1;
    }
    for(int i = 0; i < 50; i++)
    {
        a = rand() % (n-1)+1;
        x = mod_exp(a,u,n);
        for(int j = 0; j < t; j++)
        {
            y = mod_mul(x,x,n);
            if(y == 1 && x != 1 && x != n-1)
                return false;
            x = y;
        }
        if(x != 1)
            return false;
    }
    return true;
}
//PollarRho 
__int64 minFactor;
__int64 gcd(__int64 a,__int64 b)
{
    if(b == 0)
        return a;
    return gcd(b, a % b);
}

__int64 PollarRho(__int64 n, int c)
{
    int i = 1;
    srand(time(NULL));
    __int64 x = rand() % n;
    __int64 y = x;
    int k = 2;
    while(true)
    {
        i++;
        x = (mod_exp(x,2,n) + c) % n;
        __int64 d = gcd(y-x,n);
        if(1 < d && d < n)
            return d;
        if(y == x)
            return n;
        if(i == k)
        {
            y = x;
            k *= 2;
        }
    }
}
__int64 ans[1100],cnt;
void getSmallest(__int64 n, int c)
{
    if(n == 1)
        return;
    if(miller_rabbin(n))
    {
        ans[cnt++] = n;
        return;
    }
    __int64 val = n;
    while(val == n)
        val = PollarRho(n,c--);
    getSmallest(val,c);
    getSmallest(n/val,c);
}
__int64 a,b,sq;
void choose(__int64 s,__int64 val)
{
    if(s >= cnt)
    {
        if(val > a && val <= sq)
            a = val;
            return;
    }
    choose(s+1,val);
    choose(s+1,val*ans[s]);
}

int main()
{
    int T;
    __int64 G,L;
    while(~scanf("%I64d%I64d",&G,&L))
    {
        if(L == G)
        {
            printf("%I64d %I64d
",G,L);
            continue;
        }
        L /= G;
        cnt = 0;
        getSmallest(L,200);
        sort(ans, ans+cnt);
        int j = 0;
        for(int i = 1; i < cnt; i++)
        {
            while(ans[i-1] == ans[i] && i < cnt)
                ans[j] *= ans[i++];
            if ( i < cnt )
                ans[++j] = ans[i];
        }

        cnt = j+1;
        a = 1;
        sq = (__int64)sqrt(L+0.0);
        choose(0,1);
        printf("%I64d %I64d
",a*G,L/a*G);
    }
    return 0;
}