Project Euler

[ホーム](https://redmine.noppi.jp) - [[Wiki|Project Euler]]

# [[Problem 66]]

## Diophantine Equation

Consider quadratic Diophantine equations of the form:

$$x^2 - Dy^2 = 1$$

For example, when $D=13$, the minimal solution in $x$ is $649^2 - 13 \times 180^2 = 1$.

It can be assumed that there are no solutions in positive integers when $D$ is square.

<p>By finding minimal solutions in $x$ for $D = \{2, 3, 5, 6, 7\}$, we obtain the following:</p>

\begin{align}

3^2 - 2 \times 2^2 &= 1\\

2^2 - 3 \times 1^2 &= 1\\

{\color{red}{\mathbf 9}}^2 - 5 \times 4^2 &= 1\\

5^2 - 6 \times 2^2 &= 1\\

8^2 - 7 \times 3^2 &= 1

\end{align}

Hence, by considering minimal solutions in $x$ for $D \le 7$, the largest $x$ is obtained when $D=5$.

Find the value of $D \le 1000$ in minimal solutions of $x$ for which the largest value of $x$ is obtained.

## ディオファントス方程式

次の形式の, 2次のディオファントス方程式を考えよう:

$$x^2 - Dy^2 = 1$$

たとえば D=13 のとき, x を最小にする解は 6492 - 131802 = 1 である.

D が平方数(square)のとき, 正整数のなかに解は存在しないと考えられる.

<p>D = {2, 3, 5, 6, 7} に対して x を最小にする解は次のようになる:

\begin{align}

3^2 - 2 \times 2^2 &= 1\\

2^2 - 3 \times 1^2 &= 1\\

{\color{red}{\mathbf 9}}^2 - 5 \times 4^2 &= 1\\

5^2 - 6 \times 2^2 &= 1\\

8^2 - 7 \times 3^2 &= 1

\end{align}

</p>

したがって, D ≤ 7 に対して x を最小にする解を考えると, D=5 のとき x は最大である.

D ≤ 1000 に対する x を最小にする解で, x が最大になるような D の値を見つけよ.

```scheme

;;; こっちが正解

(import (scheme base)

        (gauche base)

        (util match)

        (scheme inexact))

;;; a + (√b - c) / d → (a b c d)

;;;

;;; (√b - c) / d を有理化する →

;;; (d * (√b + c)) / (b - c^2)

;;;

;;; d / (b - c^2) を既約分数にして 1/q と置く →

;;; (√b + c) / q

(define (next-fraction lis)

  (assume (= (length lis) 4))

  (match-let1 (a b c d)

              lis

              (let-values ([(isqrt-b _) (exact-integer-sqrt b)])

                (if (<= d (- (sqrt b) c))

                  `(,(+ a (* d isqrt-b)) ,b ,(+ c isqrt-b) ,d)

                  (let* ([q (/ (- b (* c c))

                               d)]

                         [next-a (div (+ isqrt-b c)

                                      q)]

                         [next-b b]

                         [next-c (- (* q next-a) c)]

                         [next-d q])

                    `(,next-a ,next-b ,next-c ,next-d))))))

(define (continued-fraction-list num)

  (assume (exact-integer? num))

  (assume (<= 2 num))

  (let* ([first-fraction (next-fraction `(0 ,num 0 1))]

         [first-int (car first-fraction)]

         [end-fraction (next-fraction first-fraction)]

         [result `(,(car end-fraction) ,first-int)])

    (let loop ([current-fraction (next-fraction end-fraction)]

               [result result])

      (if (equal? current-fraction end-fraction)

        (reverse result)

        (loop (next-fraction current-fraction)

              (cons (car current-fraction) result))))))

(define (cf-calc lis)

  (case (length lis)

    [(0) (errorf "cf-calcの引数が異常です : ~s~%" lis)]

    [(1) `(,(car lis), 1)]

    [else

      (let loop ([pn-1 (car lis)]

                 [pn (+ (* (car lis)

                           (cadr lis))

                        1)]

                 [qn-1 1]

                 [qn (cadr lis)]

                 [rest (cddr lis)])

        (if (null? rest)

          `(,pn ,qn)

          (loop pn

                (+ (* (car rest) pn)

                   pn-1)

                qn

                (+ (* (car rest) qn)

                   qn-1)

                (cdr rest))))]))

(define (diophantine d)

  (let ([cf-list (continued-fraction-list d)])

    (if (odd? (length cf-list))

      (append (cf-calc (drop-right cf-list 1))

              d)

      (let* ([temp-list-1 (cdr cf-list)]

             [temp-list-2 (append cf-list temp-list-1)]

             [temp-list (drop-right temp-list-2 1)])

        (append (cf-calc temp-list) d)))))

(define (non-square-num-list num)

  (filter (^n

            (let-values ([(_ b) (exact-integer-sqrt n)])

              (not (zero? b))))

          (iota (- num 1) 2)))

(define answer-66

  (cddr

    (fold (^[lis result]

            (if (< (car result) (car lis))

              lis

              result))

          `(0 0 . 0) 

          (map (cut diophantine <>)

               (non-square-num-list 1000)))))

(format #t "66: ~d~%" answer-66)

```

```scheme

;;; 素直なアルゴリズムで書いたら止まらなくなったので

;;; このやり方は使えない

(import (scheme base)

        (gauche base))

(define (d-list)

  (filter (^n

            (let-values ([(_ b) (exact-integer-sqrt n)])

              (not (zero? b))))

          (iota 999 2)))

(define (diophantine d)

  (let loop ([y 1])

    (let ([dy2 (* d (square y))])

      (let-values ([(x-1 dif) (exact-integer-sqrt dy2)])

        (if (zero? dif)

          (loop (+ y 1))

          (let ([x (+ x-1 1)])

            (if (= (- (square x)

                      dy2)

                   1)

              `(,x ,y . ,d)

              (loop (+ y 1)))))))))

;;; D = 109 の時に死んだので連分数展開を使用した方法に切り替える必要がある

;;; ちなみに、D = 109 の時の x と y の値は次の通りらしい…

;;; x = 158070671986249, y = 15140424455100

;;; (see https://ja.wikipedia.org/wiki/%E3%83%9A%E3%83%AB%E6%96%B9%E7%A8%8B%E5%BC%8F)

(define (diophantine-1000)

  (map (cut diophantine <>)

       (d-list)))

(define answer-66)

(format #t "66: ~d~%" answer-66)

```