ディオファントス方程式の一種、ペル方程式についての問題。
よくわからないので、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)))
立ち読みした『はじめての数論』にも書いてあった。 「連分数のでんぐり返り世界」なんて楽しそうな章タイトル。