Problem 66

Project Euler

Problem 66
各 D <= 1000 (D は平方数でない整数) に対して、

x^2 - Dy^2 = 1 (x, y は整数)

を満たす最小の x を求め、それが最大となる D を求めよ、という問題。
最初は brute force で解こうとしたけど、D = 61 で詰まった。全然整数解が見つからない。
そこで、別の方法を探した。以下、その解法。


どうやら、この手の不定方程式 (ディオファントス方程式) はペル方程式と呼ばれているらしい。
ペル方程式 - Wikipedia
解き方は
http://www004.upp.so-net.ne.jp/s_honma/pell/pell.htm
を参考にした。
リンク先のをまとめて、最終的な結論を書くと、

という漸化式があって、

g_i = g_{i+1} または h_i = h_{i+1}

となる i を見つければ、

h_i = h_{i+1} のときは x = x_i x_{i+1} + D y_i y_{i+1}     y = x_i y_{i+1} + x_{i+1} y_i
g_i = g_{i+1} のときは x = x_i^2 + D y_i^2     y = 2 x_i y_i

が、ペル方程式 x^2 - Dy^2 = \pm 1 の解。
しかし問題は +1 の場合のみなので、 x^2 < Dy^2 の場合はその次の解である

x' = x^2 + Dy^2
y' = 2xy

が求める解となる。

#include <math.h>
#include <gmp.h>

#define M 1000

/* Solve Pell Equation */

int main(void)
{
  mpz_t g, gp, h, hp, k, kp, k0, y, yp, ypp, x, xp, ansx, ansy, maxx;
  mpz_init(g);
  mpz_init(gp);
  mpz_init(h);
  mpz_init(hp);
  mpz_init(k);
  mpz_init(kp);
  mpz_init(k0);
  mpz_init(y);
  mpz_init(yp);
  mpz_init(ypp);
  mpz_init(x);
  mpz_init(xp);
  mpz_init(ansx);
  mpz_init(ansy);
  mpz_init_set_ui(maxx, 0);
  unsigned long D, maxd = 0;
  for (D = 2; D <= M; D++) {
    unsigned long d = lround(floor(sqrt(D)));
    if (d*d == D) {
      /* square number */
      continue;
    }
    mpz_set_ui(gp, 0);
    mpz_set_ui(hp, 1);
    mpz_set_ui(k0, d);
    mpz_set(kp, k0);
    mpz_set_ui(yp, 0);
    mpz_set_ui(ypp, 1);
    mpz_set_ui(xp, 1);

    for(;;) {
      //g = -gp + kp*hp;
      mpz_mul(g, kp, hp); mpz_sub(g, g, gp);
      //h = (D - g*g)/hp;
      mpz_set_ui(h, D); mpz_submul(h, g, g);  mpz_divexact(h, h, hp);
      //k = (k0 + g)/h;
      mpz_add(k, k0, g);  mpz_fdiv_q(k, k, h);
      //y = ypp + kp*yp;
      mpz_mul(y, kp, yp); mpz_add(y, y, ypp);
      //x = g*y + h*yp;
      mpz_mul(x, g, y);  mpz_addmul(x, h, yp);
      if (mpz_cmp(h, hp) == 0 || mpz_cmp(gp, g) == 0) {
        break;
      }
      mpz_set(gp, g);
      mpz_set(hp, h);
      mpz_set(kp, k);
      mpz_set(ypp, yp);
      mpz_set(yp, y);
      mpz_set(xp, x);
    }

    if (mpz_cmp(h, hp) == 0) {
      //ansx = (xp*x + D*yp*y)/h;
      mpz_mul_ui(ansx, yp, D);  mpz_mul(ansx, ansx, y);  mpz_addmul(ansx, xp, x);  mpz_divexact(ansx, ansx, h);
      //ansy = (xp*y + x*yp)/h;
      mpz_mul(ansy, xp, y); mpz_addmul(ansy, x, yp);  mpz_divexact(ansy, ansy, h);
    }
    else {
      //ansx = (xp*xp + D*yp*yp)/hp;
      mpz_mul_ui(ansx, yp, D);  mpz_mul(ansx, ansx, yp);  mpz_addmul(ansx, xp, xp); mpz_divexact(ansx, ansx, hp);
      //ansy = 2*xp*yp/hp;
      mpz_mul_ui(ansy, xp, 2);  mpz_mul(ansy, ansy, yp);  mpz_divexact(ansy, ansy, hp);
    }

    /* use x and y for temporary */
    //x = ansx*ansx
    mpz_mul(x, ansx, ansx);
    //y = D*ansy*ansy
    mpz_mul_ui(y, ansy, D); mpz_mul(y, y, ansy);
    if (mpz_cmp(x, y) < 0) {
      /* use g as temporary */
      mpz_set(g, ansx);
      //ansx = ansx*ansx + D*ansy*ansy;
      mpz_mul_ui(ansx, ansy, D);  mpz_mul(ansx, ansx, ansy);  mpz_addmul(ansx, g, g);
      //ansy = 2*ansx*ansy;
      mpz_mul(ansy, ansy, g); mpz_mul_ui(ansy, ansy, 2);
    }
    
    if (mpz_cmp(ansx, maxx) > 0) {
      mpz_set(maxx, ansx);
      maxd = D;
    }
  }
  gmp_printf("when D=%lu, minimal x=%Zu\n", maxd, maxx);
  
  mpz_clear(g);
  mpz_clear(gp);
  mpz_clear(h);
  mpz_clear(hp);
  mpz_clear(k);
  mpz_clear(kp);
  mpz_clear(k0);
  mpz_clear(y);
  mpz_clear(yp);
  mpz_clear(ypp);
  mpz_clear(x);
  mpz_clear(xp);
  mpz_clear(ansx);
  mpz_clear(ansy);
  mpz_clear(maxx);
  return 0;
}