# Project Euler Problem 66 Problem 66 asks us to find integer solutions of a quadratic Diophantine equation.

My first, exceptionally naive approach was to write:

(use '[clojure.contrib.math :only (sqrt)])

(defn square? [n] (= n (* (int (sqrt n)) (int (sqrt n)))))

(defn minimal-soln [D]
(first
(for [x (iterate inc 1)
y (range 1 x)
:when (= 1 (- (* x x) (* D y y)))]
x)))

(defn euler-66 []
(reduce max (map minimal-soln (remove square? (range 1 1001)))))

(time (euler-66))  ;; Kabooom!


It may be entertaining for you to work out why the universe will probably end before the above computation will…

A more efficient method for finding the solution of Pell’s Equation is to use Chakravala’s Method and iterate until an integer solution is found.

(use '[clojure.contrib.math :only (sqrt abs)])

(defn square? [n] (= n (* (int (sqrt n)) (int (sqrt n)))))

(defn argmin [f coll] (reduce #(if (< (f %1) (f %2)) %1 %2) coll))

(defn argmax [f coll] (reduce #(if (> (f %1) (f %2)) %1 %2) coll))

(defn minimal-soln
"Returns the minimal solution to Pell's Equation using Chakravala's method."
[D]
(loop [a (bigint (sqrt D))
b 1
k (- (* a a) D)]
(if (= k 1)
a
(let [m (argmin #(abs (- (* % %) D))
(filter #(integer? (/ (+ a (* % b)) k))
(range (* 2 (int (sqrt D))))))]
(recur (bigint (/ (+ (* a m) (* D b)) (abs k)))
(/ (+ a (* b m)) (abs k))
(/ (- (* m m ) D) k))))))

(defn euler-66 []
(argmax minimal-soln (remove square? (range 1 1001))))

(time (euler-66)) ;; "Elapsed time: 1818.137004 msecs"


The bounds on this problem probably can probably be tightened and the computation of m streamlined further, but I’m going to leave that for more interested readers.