Mathematics:GCD & LCM Inverse(POJ 2429)

　　题目大意，不翻译了，就是上面链接的意思。

　　具体思路就是要根据数论来，设a和b的GCD（最大公约数）和LCM（最小公倍数），则a/GCD*b/GCD=LCM/GCD，我们只用枚举LCM/GCD的所有质因数就可以了，然后把相应的质因数乘以GCD即可得出答案。

　　找素数很简单，用Miller_Rabin求素数的方法，可以多求几次提高正确率，原理就是用的费马定理：如果P是素数，则A^(p-1)mod P恒等于1，为了绕过Carmichael数，采用费马小定理：如果n是素数，则存在x(0<x<n)，(x*x)mod n 要么是1要么是n-1，否则，x就是合数。

　　另外就是要把因数分解了，这里暴力解法貌似可以，但是那种方法太笨了（也就是枚举），我们采用一个O(N^1/4)的算法Pollard_Rho快速因数分解，具体证明点我，我们把结论用上就可以了，具体将在代码中呈现。

　　因数分解以后，剩下就只用DFS就可以了，注意要把所有重复的质因数先成起来，那样我们找解的时候就可以保证我们找到的解都是互质的

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

　　　　http://www.cnblogs.com/chenxiwenruo/p/3567526.html

  1 #include <iostream>
  2 #include <functional>
  3 #include <algorithm>
  4 #define MAX_N 1000
  5 
  6 using namespace std;
  7 
  8 long long gcd(long long, long long);
  9 bool Miller_Rabin(const long long);
 10 long long witness(long long, long long, long long);
 11 long long Pollard_Rho(long long,long long);
 12 long long Multi_Mod(long long, long long);
 13 void Find_Factors(long long, int *const, long long);
 14 void DFS(long long, long long, int,const int);
 15 static long long factors[MAX_N],factors_sum[MAX_N],a, b, min_sum;
 16 
 17 int main(void)
 18 {
 19     long long n, GCD, LCM;
 20 
 21     while (~scanf("%lld %lld", &GCD, &LCM))
 22     {
 23         int len = 0, len_factors_sum = 0;
 24         n = LCM / GCD; 
 25         if (Miller_Rabin(n))
 26             printf("%lld %lld
", GCD, n*GCD);
 27         else if (LCM == GCD)
 28             printf("%lld %lld
", GCD, GCD);
 29         else
 30         {
 31             Find_Factors(n, &len, 120);//120是经验值
 32             sort(factors, factors + len);
 33             factors_sum[0] = factors[0];
 34             for (int i = 1; i < len; i++)//把相同的质数全部成起来，那么当DFS的时候我们只用找这些乘积就可以了（保证a,b一定互素）
 35             {
 36                 if (factors[i] == factors[i - 1])
 37                     factors_sum[len_factors_sum] *= factors[i];
 38                 else
 39                     factors_sum[++len_factors_sum] = factors[i];
 40             }
 41             a = factors[0]; b = n / a;
 42             min_sum = b + a;
 43             DFS(1, 1, 0, len_factors_sum + 1);//找解，DFS枚举就可以了
 44             if (a < b)
 45                 printf("%lld %lld
", a*GCD, b*GCD);
 46             else
 47                 printf("%lld %lld
", b*GCD, a*GCD);
 48         }
 49     }
 50     return 0;
 51 }
 52 
 53 void Find_Factors(long long n,int *const len,long long times)
 54 {
 55     if (n == 1)
 56         return;
 57     else if (Miller_Rabin(n))
 58     {
 59         factors[(*len)++] = n;//分解质因数
 60         return;
 61     }
 62     else
 63     {
 64         long long p = n;
 65         long long k = times;
 66         while (p >= n)
 67             p = Pollard_Rho(n, k--);
 68         Find_Factors(p, len, times);
 69         Find_Factors(n / p, len, times);
 70     }
 71 }
 72 
 73 bool Miller_Rabin(const long long n)
 74 {
 75     if (n == 2)
 76         return true;
 77     if (n == 1 || n < 0)
 78         return false;
 79     else
 80     {
 81         for (int i = 0; i < 5; i++)//Miller_Rabin测试方法+费马定理，叠5次减少出错几率
 82             if (!(witness((long long)rand() % (n - 1) + 1, n - 1, n) == 1))
 83                 return false;
 84         return true;
 85     }
 86 }
 87 
 88 long long witness(long long coe, long long level, long long n)
 89 {
 90     long long x, y;
 91     if (level == 0)
 92         return 1;//到达最后一层，开始后序遍历
 93 
 94     x = witness(coe, level >> 1, n);//level以幂次递减
 95     if (x == 0)
 96         return 0;//如果x出的结果是0，那么n一定是一个合数
 97     
 98     y = (x*x) % n;
 99     if (y == 1 && x != 1 && x != n - 1)
100         return 0;//费马小定理，如果一个数是素数，则x*x对n的模一定是1或者是n-1,如果不是，则是合数
101     if (level % 2  == 1)
102         y = (coe*y) % n;//和幂运算的道理是一样的
103 
104     return y;
105 }
106 
107 long long Pollard_Rho(long long n,long long c)
108 {
109     long long x, y, k = 2, d;
110     y = x = rand() % (n - 1) + 1;//y和x的初始值都是定任意一个常数，然后直到找到非平常因子为止
111 
112     for (int i = 1;; i++)
113     {
114         x = (Multi_Mod(x, n) + c) % n;//算f(x),f(x)的定义见多项式乘法f(x)=x^2+c
115         d = gcd((y - x + n) % n, n);//计算|y-x|与n的最大公因数，当y==x时，返回n，说明在这个c下无法产生非平常因子
116         if (1 < d && d < n)
117             return d;//如果得出d，那么d就是因数之一（不一定是质数，要继续判断）
118         else if (y == x)
119             return n;
120         else if (i == k)//brent判据，目的就是找到在偶数周期内找到gcd(x(k)-x(i/2))
121         {
122             y = x;
123             k <<= 1;
124         }
125     }
126     return n;
127 }
128 
129 long long gcd(long long a, long long b)
130 {
131     if (b == 0) return a;
132     return gcd(b, a%b);
133 }
134 
135 long long Multi_Mod(long long x, long long mod)//Pollard_Rho用到的多项式算法y=(x^2 + c)mod n
136 {
137     long long ans = 0, p = x;
138 
139     while (p)
140     {
141         if (p & 1)
142             ans = (ans + x) % mod;
143         x = (x << 1) % mod;//记得取模
144         p >>= 1;
145     }
146     return ans;
147 }
148 
149 void DFS(long long tmpa, long long tmpb, int level, const int len_factors_sum)
150 {
151     if (level == len_factors_sum)
152     {
153         if (tmpa + tmpb < min_sum)
154         {
155             a = tmpa; b = tmpb;
156             min_sum = tmpa + tmpb;
157         }
158     }
159     else
160     {
161         DFS(tmpa*factors_sum[level], tmpb, level + 1, len_factors_sum);
162         DFS(tmpa, tmpb*factors_sum[level], level + 1, len_factors_sum);
163     }
164 }