project euler 66

pe_66.clj

;; Problem 66  : 2011/11/29
;; "Elapsed time: 67.348225 msecs"

;; treat all num as
;; (sqrt X)+M
;; -----------
;;      N

(defn int-floor [x]
  (Math/round (Math/floor x)))

(defn inverse-rationalize [X M N]
  (let [new-N (/ (- X (* M M)) N)
	new-M (- M)]
    [new-M new-N]))
  
(defn make-proper [X M N]
  (let [nume (+ (Math/sqrt X) M)
	a (int-floor (/ nume N))
	new-M (- M (* a N))]
    [a new-M]))
	
(defn get-cf-list-loop
  "Get continued fraction parms of Sqrt X"
  ([X] (get-cf-list-loop X 0 1 nil))
  ([X M N top]
     (let [[a prop-M] (make-proper X M N)
	   [new-M new-N] (inverse-rationalize X prop-M N)]
       (lazy-seq
	(cons a
	      (if (= [new-M new-N] top)
		nil
		(get-cf-list-loop X
				  new-M
				  new-N
				  (if (empty? top)
				    [new-M new-N]
				    top))))))))

(defn calc-cf-partial-fraction [cf-list]
  (reduce #(+ (/ 1 %1) %2)
	  (drop 1 (reverse cf-list))))

(defn check-norm [X Y D]
  (= 1 (- (* X X) (* Y Y D))))

(defn my-numerator [n]
  (if (integer? n)
    n
    (numerator (rationalize n))))

(defn my-denominator [n]
  (if (integer? n)
    1
    (denominator (rationalize n))))

(defn solv-pell-eq [n]
  (let [fraction (calc-cf-partial-fraction (get-cf-list-loop n))
	num (my-numerator fraction)
	den (my-denominator fraction)]
    (if (check-norm num den n)
      [num den n]
      [(+ (* num num) (* den den  n)) (* 2 den num) n])))

(defn square? [n]
  (let [rt (Math/round (Math/sqrt n))]
    (= n (* rt rt))))

(defn pe-66 [limit]
  (reduce #(if (> (first %1) (first %2)) %1 %2)
	  (map solv-pell-eq
	       (filter (complement square?) (range 10 (inc limit))))))