(defun square-p (d)
(eq (isqrt d) (sqrt d)))
;;
;; ripped from http://www.cs.jhu.edu/~jason/software/fractions/fractions.lisp
;;
(defun true-sqrt-cf (d)
(assert (not (square-p d)))
(loop with (a b c) = (list 0 1 1)
with first-partial-quotient = (isqrt d)
with gcd
for partial-quotient = (truncate (+ a (isqrt (* b b d))) c)
collect partial-quotient
do (decf a (* c partial-quotient)) ;; subtract partial-quotient from x
(psetf a (* -1 a c) ;; take reciprocal via parallel assignment
b (* b c)
c (- (* b b d) (* a a)))
(setf gcd (gcd a b c) ;; put into lowest terms
a (/ a gcd)
b (/ b gcd)
c (/ c gcd))
until (= partial-quotient (* 2 first-partial-quotient)))) ;; about to repeat
;; answer problem 64
(defvar *count* 0)
(loop for i from 1 to 10000 do
(if (not (square-p i))
(if (evenp (length (true-sqrt-cf i)))
(incf *count*))))
(format t "~A~%" *count*)