题目110：找出一种能够高效分析分析方程1/x + 1/y = 1/n的解的个数的算法

欧拉计划

Diophantine reciprocals II

In the following equation x, y, and n are positive integers.

It can be verified that when n = 1260 there are 113 distinct solutions and this is the least value of n for which the total number of distinct solutions exceeds one hundred.

What is the least value of n for which the number of distinct solutions exceeds four million?

NOTE: This problem is a much more difficult version of Problem 108 and as it is well beyond the limitations of a brute force approach it requires a clever implementation.

题目：

在如下方程中，x，y 和 n 都是正整数。

可以证明当 n = 1260 时，该方程一共有 113 个不同的解。而这也是使该方程不同解的个数超过 100 的最小的 n。

使得该方程不同解的个数超过四百万的最小的 n 是多少？

注意：该题目是题目108的一个更难的版本，并且已经超过了暴力手段的范围，所以你需要一个更聪明的实现方法。

代号-K

1. void getDiophantineReciprocals(int flag){
   
2.     int num = 1;
   
3.     int count;
   
4.     int i;
   
5.     ElementType numPow;
   
6.     while(true){
   
7.         count = 1;
   
8.         i = 2;
   
9.         num++;
   
10.         numPow = (ElementType) pow(num, 2);
    
11.         while(i <= num){
    
12.             if(numPow % i == 0){
    
13.                 count++;
    
14.             }
    
15.             i++;
    
16.         }
    
17.         if(count >= flag){
    
18.             printf("get it\n");
    
19.             printf("%d has %d\n", num,count);
    
20.             return;
    
21.         }
    
22.         printf("%d has %d\n",num,count);
    
23.     }
    
24. }
    
25. 
    
26. int main(int argc, char *argv[]){
    
27.     getDiophantineReciprocals(110);
    
28.     return 0;
    
29. }

复制代码

代号-K

下面的代码是利用质数进行求解的

1. void getDiophantinePrime(ElementType flag){
   
2.     ElementType prime[10000] = {0};
   
3.     QHash<ElementType, ElementType> ret;
   
4.     QHash<ElementType,ElementType>::iterator iter;
   
5.     prime[0] = 2;
   
6.     ElementType i = 1;
   
7.     ElementType j = 2;
   
8.     ElementType k = 0;
   
9.     bool primeFlag = true;
   
10.     ElementType value = 0;
    
11.     ElementType answer = 1;
    
12.     while(true){
    
13.         k = 0;
    
14.         j++;
    
15.         value = j;
    
16.         answer = 1;
    
17.         primeFlag = true;
    
18.         while(k < i){
    
19.             while(value % prime[k] == 0){
    
20.                 if(ret.find(prime[k]) != ret.end()){
    
21.                     ret.find(prime[k]).value() += 1;
    
22.                 }
    
23.                 else{
    
24.                     ret.insert(prime[k],1);
    
25.                 }
    
26.                 value /= prime[k];
    
27.                 primeFlag = false;
    
28.             }
    
29.             if(value > prime[k]) k++;
    
30.             else break;
    
31.         }
    
32.         if(primeFlag) prime[i++] = value;
    
33.         else{
    
34.             iter = ret.begin();
    
35.             while(iter != ret.end()){
    
36.                 answer *= (iter.value()*2+1);
    
37.                 iter++;
    
38.             }
    
39.             answer = (answer + 1)/2;
    
40.             printf("%lld has %lld\n",j, answer);
    
41.             if(answer >= flag){
    
42.                 printf("get it\n");
    
43.                 return;
    
44.             }
    
45.             ret.clear();
    
46.         }
    
47.     }
    
48. }
    
49. 
    
50. int main(int argc, char *argv[]){
    
51.     getDiophantinePrime(100);
    
52.     return 0;
    
53. }

复制代码

代号-K

还有一种简单的方法，可以直接进行计算出来
把一个大数分解成N个质数的乘积

1. 
   
2. 例如：
   
3. 510510 = 2*3*5*7*11*13*17
   
4. 那么共有
   
5. ((3^7)+1)/2 = 1094种可能
   
6. 
   
7. 再例如：
   
8. 60 = 2^2*3*5
   
9. 那么共有
   
10. ((5*3*3)+1)/2 = 23种可能
    
11. 120 = 2^3*3*5
    
12. 那么共有
    
13. ((7*3*3)+1)/2 = 32种可能
    

复制代码

聪明的人应该知道其中的规律了吧

那么超过四佰万的数呢，应该怎么计算呢

1. 首先假设(x+1)/2 = n(n是大于四百万的数)
   
2. x = 2*n - 1;
   
3. 令 m = log3(x)
   
4. 则 这个数 X = 2*3*5*7*11*.....（大约m个质数）
   
5. 这个X 就差不多是答案了，可以再进行优化
   
6. X = 2^2*3*5*7.....(大约m-1个质数)
   
7. 不断优化，可以求出最小的值

复制代码

guosl

本题的数学理论的结论请参考108题中我的论述，这里就不再啰嗦了。
做本题还需要知道的一件事是：一个数x的素因数分解为p1^k1 ... pn^kn，则这个数的因子个数为：(k1+1)*...*(kn+1)。
由于本题中涉及的数值太大，我们使用了mpir大数值库（具体参考：http://www.mpir.org/）。
对于值M=2*3*...*47（即前15个素数的乘积）就是一个满足条件的数。但可能不是最小的值，但是根据这个值我们可以得出以下的结论（数学的证明略）：
1. 这个最小值小于等于M。
2. 如果令N为欲求的最小值那么在这个N的素因数分解中一定不会出现值比47还大的素数。
3. N的素因数分解中一定是从2开始的连续若干个素数。
为此，我们只需用深度搜索即可。

1. /*
   
2. 答案：9350130049860600
   
3. 耗时：0.108634秒
   
4. */
   
5. #include <iostream>
   
6. #include <algorithm>
   
7. #include <mpirxx.h>
   
8. #include <omp.h>
   
9. using namespace std;
   
10. 
    
11. const int N = 8000000;//为n^2的因子数，同时也是n^2的分解数的2倍，为因子数搜索达到要求的值
    
12. const int p[15] = { 2,3,5,7,11,13,17,19,23,29,31,37,41,43,47 };
    
13. mpz_class M = 1; //搜索的上界也是最后的结果
    
14. 
    
15. void dfs(int k, //搜索的深度
    
16.   int fs,       //当前值z^2的因子数
    
17.   mpz_class z)  //当前的值
    
18. {
    
19.   if (k >= 15)//迭代深度到底了
    
20.     return;
    
21.   if (z > M)//迭代的计算值超过了当前的值，没有再做下去的意义了
    
22.     return;
    
23.   if (fs > N)//因子数达到了要求
    
24.   {
    
25.     M = min(M, z);//保留最小值
    
26.     return;
    
27.   }
    
28.   mpz_class z1 = z*p[k];//当前的计算值
    
29.   int x = 1, fs1 = fs*(2 * x + 1);//当前z1^2的因子个数
    
30.   while (z1 < M)//进行搜索
    
31.   {
    
32.     dfs(k + 1, fs1, z1);
    
33.     ++x;
    
34.     z1 *= p[k];
    
35.     fs1 = fs*(2 * x + 1);
    
36.   }
    
37. }
    
38. 
    
39. int main(void)
    
40. {
    
41.   double t = omp_get_wtime();
    
42.   for (int i = 0; i < 15; ++i) //先求出一个初步的估值
    
43.     M *= p[i];
    
44.   dfs(0, 1, 1);//深搜结果
    
45.   t = omp_get_wtime() - t;
    
46.   cout << M << endl << t << endl;//输出结果
    
47.   return 0;
    
48. }

复制代码