题目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

void getDiophantineReciprocals(int flag){
    int num = 1;
    int count;
    int i;
    ElementType numPow;
    while(true){
        count = 1;
        i = 2;
        num++;
        numPow = (ElementType) pow(num, 2);
        while(i <= num){
            if(numPow % i == 0){
                count++;
            }
            i++;
        }
        if(count >= flag){
            printf("get it\n");
            printf("%d has %d\n", num,count);
            return;
        }
        printf("%d has %d\n",num,count);
    }
}

int main(int argc, char *argv[]){
    getDiophantineReciprocals(110);
    return 0;
}

代号-K

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

void getDiophantinePrime(ElementType flag){
    ElementType prime[10000] = {0};
    QHash<ElementType, ElementType> ret;
    QHash<ElementType,ElementType>::iterator iter;
    prime[0] = 2;
    ElementType i = 1;
    ElementType j = 2;
    ElementType k = 0;
    bool primeFlag = true;
    ElementType value = 0;
    ElementType answer = 1;
    while(true){
        k = 0;
        j++;
        value = j;
        answer = 1;
        primeFlag = true;
        while(k < i){
            while(value % prime[k] == 0){
                if(ret.find(prime[k]) != ret.end()){
                    ret.find(prime[k]).value() += 1;
                }
                else{
                    ret.insert(prime[k],1);
                }
                value /= prime[k];
                primeFlag = false;
            }
            if(value > prime[k]) k++;
            else break;
        }
        if(primeFlag) prime[i++] = value;
        else{
            iter = ret.begin();
            while(iter != ret.end()){
                answer *= (iter.value()*2+1);
                iter++;
            }
            answer = (answer + 1)/2;
            printf("%lld has %lld\n",j, answer);
            if(answer >= flag){
                printf("get it\n");
                return;
            }
            ret.clear();
        }
    }
}

int main(int argc, char *argv[]){
    getDiophantinePrime(100);
    return 0;
}

代号-K

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

例如：
510510 = 2*3*5*7*11*13*17
那么共有
((3^7)+1)/2 = 1094种可能

再例如：
60 = 2^2*3*5
那么共有
((5*3*3)+1)/2 = 23种可能
120 = 2^3*3*5
那么共有
((7*3*3)+1)/2 = 32种可能

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

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

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

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开始的连续若干个素数。
为此，我们只需用深度搜索即可。

/*
答案：9350130049860600
耗时：0.108634秒
*/
#include <iostream>
#include <algorithm>
#include <mpirxx.h>
#include <omp.h>
using namespace std;

const int N = 8000000;//为n^2的因子数，同时也是n^2的分解数的2倍，为因子数搜索达到要求的值
const int p[15] = { 2,3,5,7,11,13,17,19,23,29,31,37,41,43,47 };
mpz_class M = 1; //搜索的上界也是最后的结果

void dfs(int k, //搜索的深度
  int fs,       //当前值z^2的因子数
  mpz_class z)  //当前的值
{
  if (k >= 15)//迭代深度到底了
    return;
  if (z > M)//迭代的计算值超过了当前的值，没有再做下去的意义了
    return;
  if (fs > N)//因子数达到了要求
  {
    M = min(M, z);//保留最小值
    return;
  }
  mpz_class z1 = z*p[k];//当前的计算值
  int x = 1, fs1 = fs*(2 * x + 1);//当前z1^2的因子个数
  while (z1 < M)//进行搜索
  {
    dfs(k + 1, fs1, z1);
    ++x;
    z1 *= p[k];
    fs1 = fs*(2 * x + 1);
  }
}

int main(void)
{
  double t = omp_get_wtime();
  for (int i = 0; i < 15; ++i) //先求出一个初步的估值
    M *= p[i];
  dfs(0, 1, 1);//深搜结果
  t = omp_get_wtime() - t;
  cout << M << endl << t << endl;//输出结果
  return 0;
}