题目135：求方程x^2 - y^2 - z^2=n的解的个数

欧拉计划

Same differences

Given the positive integers, x, y, and z, are consecutive terms of an arithmetic progression, the least value of the positive integer, n, for which the equation, x² − y² − z² = n, has exactly two solutions is n = 27:

34² − 27² − 20² = 12² − 9² − 6² = 27

It turns out that n = 1155 is the least value which has exactly ten solutions.

How many values of n less than one million have exactly ten distinct solutions?

题目：

给定等差数列中的三个数字 x，y 和 z，最小 x² − y² − z² = n 有两个解的最小的 n 是 27：

34² − 27² − 20² = 12² − 9² − 6² = 27

事实上 n = 1155 是使得上述方程具有 10 个解的最小值。

一百万以下有多少个数使得上述方程有 10 个不同的解？

jerryxjr1220

"""
(z+2d)^2-(z+d)^2-z^2=n
3d^2+2dz-z^2=n
(3d-z)*(d+z)=n
z<3d => d > z/3
"""
from numba import jit
import numpy as np 
@jit
def solve(target, limit):
        master = np.zeros(limit+1,dtype='int16')
        for z in range(1, limit):
                for d in range(z//3+1, limit):
                        if (3*d-z)*(d+z)>limit:
                                break
                        master[(3*d-z)*(d+z)] += 1
        count = 0
        for value in master:
                if value == target:
                        count += 1
        return count
print(solve(10, 1000000))

4989
[Finished in 5.2s]

guosl

/*
答案：4989
耗时：1.41674 (4核)
          6.37623 (单线程)
解题思路：令方程为：(x+d)^2 - x^2 - (x-d)^2 = n
=> (x+d)(3d-x)=n
所以对n进行因数分解为：n=m1 x m2
d=(m1+m2)/4,x=m1 - d
*/
#include <iostream>
#include <set>
#include <cmath>
#include <omp.h>
using namespace std;

int main(void)
{
  double t = omp_get_wtime();
  int nCount = 0;
#pragma omp parallel for reduction(+:nCount) schedule(dynamic,64)
  for (int n = 1155; n <= 1000000; ++n)
  {
    set<int> s;//记录所以整数解
    int d1 = (int)sqrt((double)n);
    for (int m1 = 1; m1 <= d1; ++m1) //枚举n的因数分解
    {
      if (n % m1 == 0)
      {
        int m2 = n / m1;
        if ((m1 + m2) % 4 != 0) //检查d是否是一个整数
          continue;
        int d = (m1 + m2) / 4;
        if (m1 - d > 0)
          s.insert(m1 - d);//得到一个整数解
        if (m2 - d > 0)
          s.insert(m2 - d);//得到另一个整数解
      }
    }
    if (s.size() == 10)
      ++nCount;
  }
  t = omp_get_wtime() - t;
  cout << nCount << endl << t << endl;
  return 0;
}