Square root digits
This commit is contained in:
parent
9f2fc1d75e
commit
bf29830bdd
|
@ -0,0 +1,74 @@
|
|||
;;; 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))))
|
Loading…
Reference in New Issue