From bf29830bdd188be4c39cc16d2a455ae19ddff81f Mon Sep 17 00:00:00 2001
From: scms <scms@tilde.team>
Date: Fri, 1 Mar 2024 15:14:39 -0800
Subject: [PATCH] Square root digits

---
 p80.lisp | 74 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 74 insertions(+)
 create mode 100644 p80.lisp

diff --git a/p80.lisp b/p80.lisp
new file mode 100644
index 0000000..c84f1dc
--- /dev/null
+++ b/p80.lisp
@@ -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))))