75 lines
2.7 KiB
Common Lisp
75 lines
2.7 KiB
Common Lisp
;;; Approach: use continued fractions from Pell's equation, with the
|
|
;;; understanding that the error on the radicand is (EXPT (/ 1 Y) 2),
|
|
;;; so ensuring Y is greater than (EXPT 10 50) should be sufficient
|
|
|
|
;;; From problems 64 and 65
|
|
(defun compute-floor-of-square-root (n)
|
|
"Returns the floor of the square root of N"
|
|
(loop for last = nil then i
|
|
for i upfrom 1
|
|
when (> (* i i) n)
|
|
return last))
|
|
|
|
(defun find-square-root-period (n)
|
|
"Returns the period of the constants in the continued fraction which represents (SQRT N)"
|
|
(loop with states = (make-hash-table :test 'equal)
|
|
with a0 = (compute-floor-of-square-root n)
|
|
with a = a0
|
|
with m = 0
|
|
with d = 1
|
|
for constants = (list a0) then (cons a constants)
|
|
for index upfrom 1
|
|
for state = (list a m d)
|
|
for previous-index = (gethash state states)
|
|
do (when previous-index
|
|
(return (values (- index previous-index)
|
|
(reverse (rest constants)))))
|
|
(setf (gethash state states) index)
|
|
(setf m (- (* d a) m)) ; By subtracting A, inverting, and moving root to numerator
|
|
(setf d (/ (- n (* m m)) d)) ; Same process as previous
|
|
(setf a (floor (/ (+ a0 m) d))))) ; The constant is just the floor of the new fraction (... had to look this one up...)
|
|
|
|
(defun get-square-root-continued-fraction (n)
|
|
"Returns a circular list representing the continued fraction constants of (SQRT N)"
|
|
(multiple-value-bind (period constants) (find-square-root-period n)
|
|
(declare (ignore period))
|
|
(setf (rest (last constants)) (nthcdr 1 constants))
|
|
constants))
|
|
|
|
(defun compute-continued-fraction-until (constants iterations)
|
|
"Returns the continued fraction using CONSTANTS for ITERATIONS iterations"
|
|
(+ (first constants)
|
|
(if (> iterations 1)
|
|
(/ 1 (compute-continued-fraction-until (rest constants) (1- iterations)))
|
|
0)))
|
|
|
|
(defun compute-square-root (n digits)
|
|
"Computes the continued fraction until DIGITS fractional digits are exact"
|
|
(let ((continued-fraction (get-square-root-continued-fraction n))
|
|
(min (expt 10 (ceiling (/ (1- digits) 2)))))
|
|
(loop for i upfrom 200
|
|
for fraction = (compute-continued-fraction-until continued-fraction i)
|
|
until (> (denominator fraction) min)
|
|
finally (return fraction))))
|
|
|
|
(defun sum-digits (fraction n)
|
|
"Sums the first N digits of the decimal represenation of (positive) FRACTION"
|
|
(loop with sum = 0
|
|
for d = (truncate fraction)
|
|
repeat n
|
|
do (incf sum d)
|
|
(setf fraction (* (- fraction d) 10))
|
|
finally (return sum)))
|
|
|
|
(defun get-squares (max)
|
|
(loop for n upfrom 1
|
|
for square = (* n n)
|
|
until (> square max)
|
|
collect square))
|
|
|
|
(defun square-root-digital-expansion (&key (max 100) (digits 100))
|
|
(let ((squares (get-squares max)))
|
|
(loop for n from 2 upto max
|
|
unless (member n squares)
|
|
sum (sum-digits (compute-square-root n digits) digits))))
|