형태의 디오판토스 방정식 풀기
의 번째 근사분수(convergents)를 라 하면, 펠 방정식의 해 중 를 최소화하는 해 은 를 만족시키는 가 된다. 는 다음과 같이 구할 수 있다.
첫 번째 근사분수는 이므로 는 각각 다음과 같이 초기화할 수 있다.
두 번째 근사분수는 이므로, 를 구하는 식에 을 대입해 계산하면 의 값은 다음과 같아야 한다.
이를 바탕으로, 를 계산해 이 되는 경우를 구하면 가 펠 방정식의 해가 된다.
문제 64에서 의 연분수꼴을 구하는 함수를 이미 구현했으므로 그대로 가져다 쓸 수 있다. 펠 방정식의 기본 해를 구하는 함수는 다음과 같이 작성할 수 있다. 루프를 돌면서 처음 찾은 해가 를 최소화하는 해가 된다.
(ns p066
(:require [util :refer [perfect-square?]]
[p064 :refer [expand-continued-fraction]]))
(defn find-fundamental-solution
"find the fundamental solution of pell's equation, x^2 - dy^2 = 1"
[d]
(let [cf (expand-continued-fraction d)
as (lazy-cat cf (cycle (rest cf)))]
(loop [h2 0, h1 1,
k2 1, k1 0,
as as, n 0]
(if (and (>= n 1) (= 1 (-' (*' h1 h1) (*' d k1 k1))))
[h1 k1]
(recur h1 (+' (*' (first as) h1) h2)
k1 (+' (*' (first as) k1) k2)
(rest as) (inc n))))))
이제 를 2부터 1000까지 변경하면서 펠 방정식에서 를 최소화하는 기본해를 구한 다음, 이 중 값이 가장 큰 경우의 를 구하면 된다. 완전제곱수인 경우에는 을 제외한 다른 정수해가 존재하지 않으므로 이를 제외해야 한다.
(defn solve []
(->> (range 2 (inc 1000))
(filter (complement perfect-square?))
(map (fn [d] [d (first (find-fundamental-solution d))]))
(apply (partial max-key second))))
실행 결과는 다음과 같다.
p066=> (time (solve)) "Elapsed time: 21.273114 msecs" 6?1