题目110：找出一种能够高效分析分析方程1/x + 1/y = 1/n的解的个数的算法
Diophantine reciprocals IIIn 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的一个更难的版本，并且已经超过了暴力手段的范围，所以你需要一个更聪明的实现方法。
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;
} 下面的代码是利用质数进行求解的
void getDiophantinePrime(ElementType flag){
ElementType prime = {0};
QHash<ElementType, ElementType> ret;
QHash<ElementType,ElementType>::iterator iter;
prime = 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 == 0){
if(ret.find(prime) != ret.end()){
ret.find(prime).value() += 1;
}
else{
ret.insert(prime,1);
}
value /= prime;
primeFlag = false;
}
if(value > prime) k++;
else break;
}
if(primeFlag) prime = 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;
} 还有一种简单的方法，可以直接进行计算出来
把一个大数分解成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种可能
聪明的人应该知道其中的规律了吧{:5_90:}
那么超过四佰万的数呢，应该怎么计算呢
首先假设(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 于 2022-2-7 11:29 编辑
本题的数学理论的结论请参考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 = { 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;//当前的计算值
int x = 1, fs1 = fs*(2 * x + 1);//当前z1^2的因子个数
while (z1 < M)//进行搜索
{
dfs(k + 1, fs1, z1);
++x;
z1 *= p;
fs1 = fs*(2 * x + 1);
}
}
int main(void)
{
double t = omp_get_wtime();
for (int i = 0; i < 15; ++i) //先求出一个初步的估值
M *= p;
dfs(0, 1, 1);//深搜结果
t = omp_get_wtime() - t;
cout << M << endl << t << endl;//输出结果
return 0;
}
页:
[1]