POJ 2429 -- GCD & LCM Inverse

http://poj.org/problem?id=2429

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

http://www.hankcs.com/program/cpp/poj-2429-gcd-lcm-inverse.html

分析：由于 G和L分别是a，b的最大公约数和最小公倍数，那么G一定整除L，否则就没有解了。

进一步可以知道L/G的素因子分解式中，同一项不会同时出现在a，b中，因为相同的在G中已经除掉了.

那么有：，所以我们现在只需要枚举L/G的约数，枚举的一个x在a中，那么另一个(L/G)/x就是在b中，但是我们知道

,,所以思路就是对L/G先分解，因为这个数很大，所以就用Pollard-rho分解法吧，分解之后，利用 dfs找因子，然后每找到一个就判断一次来记录a+b最小的，然后就分析完毕。

1. const int Times=10;  
2. const LL INF=(LL)1<<61;  
3. const int N=550;  
4.   
5. LL n,m,ct,cnt,mini,mina,minb,ans;  
6. LL fac[N],num[N];  
7.   
8. LL gcd(LL a,LL b)  
9. {  
10.     return b? gcd(b,a%b):a;  
11. }  
12.   
13. LL multi(LL a,LL b,LL m)  
14. {  
15.      LL ans=0;  
16.      while(b)  
17.      {  
18.          if(b&1)  
19.          {  
20.              ans=(ans+a)%m;  
21.              b--;  
22.          }  
23.          b>>=1;  
24.          a=(a+a)%m;  
25.      }  
26.      return ans;  
27. }  
28.   
29. LL quick_mod(LL a,LL b,LL m)  
30. {  
31.      LL ans=1;  
32.      a%=m;  
33.      while(b)  
34.      {  
35.          if(b&1)  
36.          {  
37.              ans=multi(ans,a,m);  
38.              b--;  
39.          }  
40.          b>>=1;  
41.          a=multi(a,a,m);  
42.      }  
43.      return ans;  
44. }  
45.   
46. bool Miller_Rabin(LL n)  
47. {  
48.     if(n==2) return true;  
49.     if(n<2||!(n&1)) return false;  
50.     LL a,m=n-1,x,y;  
51.     int k=0;  
52.     while((m&1)==0)  
53.     {  
54.         k++;  
55.         m>>=1;  
56.     }  
57.     for(int i=0;i
58.     {  
59.         a=rand()%(n-1)+1;  
60.         x=quick_mod(a,m,n);  
61.         for(int j=0;j
62.         {  
63.             y=multi(x,x,n);  
64.             if(y==1&&x!=1&&x!=n-1) return false;  
65.             x=y;  
66.         }  
67.         if(y!=1) return false;  
68.     }  
69.     return true;  
70. }  
71.   
72. LL Pollard_rho(LL n,LL c)  
73. {  
74.      LL x,y,d,i=1,k=2;  
75.      y=x=rand()%(n-1)+1;  
76.      while(true)  
77.      {  
78.          i++;  
79.          x=(multi(x,x,n)+c)%n;  
80.          d=gcd((y-x+n)%n,n);  
81.          if(1return d;  
82.          if(y==x) return n;  
83.          if(i==k)  
84.          {  
85.              y=x;  
86.              k<<=1;  
87.          }  
88.      }  
89. }  
90.   
91. void find(LL n,int c)  
92. {  
93.      if(n==1) return;  
94.      if(Miller_Rabin(n))  
95.      {  
96.          fac[ct++]=n;  
97.          return ;  
98.      }  
99.      LL p=n;  
100.      LL k=c;  
101.      while(p>=n) p=Pollard_rho(p,c--);  
102.      find(p,k);  
103.      find(n/p,k);  
104. }  
105.   
106. void dfs(LL c,LL value)  
107. {  
108.      LL s=1,a,b;  
109.      if(c==cnt)  
110.      {  
111.          a=value;  
112.          b=ans/a;  
113.          if(gcd(a,b)==1)  
114.          {  
115.              a*=n;  
116.              b*=n;  
117.              if(a+b
118.              {  
119.                  mini=a+b;  
120.                  mina=a;  
121.                  minb=b;  
122.              }  
123.          }  
124.          return ;  
125.      }  
126.      for(LL i=0;i<=num[c];i++)  
127.      {  
128.          if(s*value>mini) return;  
129.          dfs(c+1,s*value);  
130.          s*=fac[c];  
131.      }  
132. }  
133.   
134. int main()  
135. {  
136.     while(~scanf("%llu%llu",&n,&m))  
137.     {  
138.          if(n==m)  
139.          {  
140.              printf("%llu %llu\n",n,m);  
141.              continue;  
142.          }  
143.          mini=INF;  
144.          ct=cnt=0;  
145.          ans=m/n;  
146.          find(ans,120);  
147.          sort(fac,fac+ct);  
148.          num[0]=1;  
149.          int k=1;  
150.          for(LL i=1;i
151.          {  
152.              if(fac[i]==fac[i-1])  
153.                  ++num[k-1];  
154.              else  
155.              {  
156.        

http://blog.sina.com.cn/s/blog_c16aeb090101d2ab.html

逆最大公约数和最小公倍数：为了嘲讽 AOJ 0005 GCD and LCM，现在给出两个数的gcd和lcm，求原来的这两个数（限定两数之和最小）。

2.6 数学问题的解题窍门

作者显然高估了读者的数学修养，我对数论一点都不熟悉，这题光靠前面介绍的一点数论皮毛无从下手，还是看了人家的代码才知道要用Rabin-Miller强伪素数测试和Pollard r因数分解算法。这两个算法的详解放在最后面，不时回来复习一下。

基本思路是 

• (a / gcd) * (b / gcd) = lcm / gcd ，所以需要分解lcm / gcd 。将其分解为互质的两个数，如果这两个数之和最小，那么乘上gcd就是所求的答案。
• 但是题目数据太大，需要一个高效的素数检测算法，所以采用Rabin-Miller强伪素数测试
• 然后分解成质因子的n次方之积，从这些n次方中挑选一些作为x，剩下的作为y。枚举x和y的所有可能，找出最小值。

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;

}

// 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;

}

// Qualifier: Pollard 因数分解算法

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);

}

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;

}

关于两个算法的讲解，有一个超棒的文档。
Read full article from 2429 -- GCD & LCM Inverse