题目157：7求解丢番图方程1/a+1/b= p/10^n

欧拉计划

Solving the diophantine equation ¹/_a+¹/_b= ^p/₁₀n

Consider the diophantine equation ¹/_a+¹/_b= ^p/₁₀n with a, b, p, n positive integers and a ≤ b.

For n=1 this equation has 20 solutions that are listed below:



How many solutions has this equation for 1 ≤ n ≤ 9?

题目：

考虑丢番图方程 ¹/_a+¹/_b= ^p/₁₀n ，a, b, p, n 都是正数，并且 a ≤ b

对于 n=1，则该方程有下面二十个解：

¹/₁+1/1=²⁰/₁₀            ¹/₁+¹/₂=¹⁵/₁₀           ¹/₁+¹/₅=¹²/₁₀            ¹/₁+¹/₁₀=¹¹/₁₀            ¹/₂+¹/₂=¹⁰/₁₀
¹/₂+¹/₅=⁷/₁₀              ¹/₂+¹/₁₀=⁶/₁₀            ¹/₃+¹/₆=⁵/₁₀            ¹/₃+¹/₁₅=⁴/₁₀              ¹/₄+¹/₄=⁵/₁₀
¹/₄+¹/₂₀=³/₁₀             ¹/₅+¹/₅=⁴/₁₀           ¹/₅+¹/₁₀=³/₁₀            ¹/₆+¹/₃₀=²/₁₀              ¹/₁₀+¹/₁₀=²/₁₀
¹/₁₁+¹/₁₁₀=¹/₁₀         ¹/₁₂+¹/₆₀=¹/₁₀          ¹/₁₄+¹/₃₅=¹/₁₀          ¹/₁₅+¹/₃₀=¹/₁₀          ¹/₂₀+¹/₂₀=¹/₁₀

请问，对于 1 ≤ n ≤ 9，该方程有多少个解？

guosl

先将p设为1，则原方程化为1/a + 1/b = 1/10^n。令m=10^n,a=m+x，b=m+y。则方程可以化为xy=m^2。对m^2进行因数分解可以求解，如果a,b有公因数p>1，则这个解可以化成1/a+1/b=p/10^n的解。

/*
答案：53490
*/
#include <iostream>
#include <cmath>
using namespace std;

const int N = 9;

unsigned long long mpow(int n, int m)
{
  unsigned long long p = 1;
  for (int i = 0; i < n; ++i)
    p *= 2;
  for (int i = 0; i < m; ++i)
    p *= 5;
  return p;
}

unsigned long long gcd(unsigned long long x, unsigned long long y)
{
  if (y == 0)
    return x;
  unsigned long long z = x % y;
  while (z != 0)
  {
    x = y;
    y = z;
    z = x % y;
  }
  return y;
}

unsigned long long getmp(unsigned long long k)
{
  if ((k % 2) == 0)
    return 2;
  unsigned long long u = (unsigned long long)sqrt((double)k);
  for (unsigned long long i = 3; i <= u; i += 2)
  {
    if (k % i == 0)
      return i;
  }
  return k;
}

int getdc(unsigned long long d)
{
  int nC = 1;
  while (d != 1)
  {
    unsigned long long p = getmp(d);
    int k = 0;
    while (d % p == 0)
    {
      ++k;
      d /= p;
    }
    nC *= k + 1;
  }
  return nC;
}

int main(void)
{
  int nCount = 0;
  unsigned long long M = 1;
  for (int n = 1; n <= N; ++n)
  {
    M *= 10;
    for (int i = 0; i <= 2 * n; ++i)
    {
      unsigned long long K1 = mpow(i, 0);
      if (K1 > M)
        break;
      for (int j = 0; j <= 2 * n; ++j)
      {
        unsigned long long K2 = K1 * mpow(0, j);
        if (K2 > M)
          break;
        unsigned long long a = M + K2, b = M + mpow(2 * n - i, 2 * n - j);
        unsigned long long d = gcd(a, b);
        nCount += getdc(d);
      }
    }
  }
  cout << nCount << endl;
  return 0;
}