Project Euler 66

2016-08-28 09:09:11

ディオファントス方程式の一種、ペル方程式についての問題。

よくわからないので、yを1から試していく。

(defun diophantine-p (d)
  (labels ((f (y)
             (let ((x (+ 1 (* d y y))))
               (if (squarep x)
                   (floor (sqrt x))
                   (f (1+ y))))))
    (f 1)))

(defun euler-066-poor ()
  (apply #'max
         (mapcar #'diophantine-p
                 (loop for x from 2 upto 1000 unless (squarep x) collect x))))

これだと d=61 などで時間がかかる。

調べてみると、平方根の連分数展開とペル方程式の関係を基にした問題みたい。

方針としては、連分数展開をして、連分数の循環ごとに普通の分数化して、 ペル方程式を満たすかどうかをチェックすれば、解けそうなことがわかった。

連分数と Pell の方程式(PDF) 2.6.4 を参考に #64, #65 で使った関数を使う。

(defun init-contfrac (d)
  (reverse (cdr (continued-fraction d))))

(defun next-contfrac (cf)
  (append cf (cons (* 2 (car cf)) (cdr cf))))

(defun check-diophantine-equation (x y d)
  (= (* x x) (1+ (* d y y))))

(defun minimum-solution-of-diophantine-equation (d)
  (labels ((f (cf)
             (let* ((frac (nth-convergent cf))
                    (x (numerator frac))
                    (y (denominator frac)))
               (if (check-diophantine-equation x y d)
                   (cons d x)
                   (f (next-contfrac cf))))))
    (f (init-contfrac d))))

(defun euler-066 ()
  (caar (sort
         (mapcar #'minimum-solution-of-diophantine-equation
                 (loop for x from 2 upto 1000 unless (squarep x) collect x))
         #'> :key 'cdr)))

立ち読みした『はじめての数論』にも書いてあった。 「連分数のでんぐり返り世界」なんて楽しそうな章タイトル。