Amicable chains

p95.lisp

@@ -0,0 +1,44 @@
+;;; Approach: Precompute all proper divisor sums and then check for
+;;; cycles.
+
+(defvar *divisor-sums*)
+
+(defun compute-divisor-sums (max)
+  "Precomputes a vector of all proper divisor sums up to MAX"
+  (ret (sums (make-array (1+ max) :initial-element 1))
+    (setf (aref sums 0) 0)
+    (setf (aref sums 1) 0)
+    (loop for n from 2 below max do
+      (loop for m upfrom (* 2 n) by n
+	    while (<= m max)
+	    do (incf (aref sums m) n)))))
+
+(defun sum-proper-divisors (n)
+  "Uses precomputed *DIVISOR-SUMS* to calculate the sum of proper divisors of N"
+  (aref *divisor-sums* n))
+
+(defun amicable-chain (n &optional (max 1000000))
+  "Returns the amicable chain starting with N (bailing out if any elements exceed MAX)"
+  (loop with first = n
+	with chain = (list first)
+	do (setf n (sum-proper-divisors n))
+	   (when (= n first)
+	     (return (nreverse chain)))
+	   (when (or (> n max)
+		     (<= n 1)
+		     (member n chain))
+	     (return nil))
+	   (push n chain)))
+
+(defun amicable-chains (&optional (max 1000000))
+  (let ((*divisor-sums* (compute-divisor-sums max)))
+    (loop with max-length = 0
+	  with best = nil
+	  for n from 1 upto max
+	  for chain = (amicable-chain n max)
+	  for length = (length chain)
+	  do (when (> length max-length)
+	       (setf max-length length)
+	       (setf best chain))
+	  finally (return (values (reduce #'min best)
+				  best)))))

shared.lisp

@@ -60,3 +60,7 @@
 (defun n-choose-r (n r)
   (/ (factorial n)
      (* (factorial r) (factorial (- n r)))))
+
+(defun divisiblep (n divisor)
+  "Returns non-NIL when N is divisible by DIVISOR"
+  (= 0 (mod n divisor)))