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.