# HG changeset patch
# User Steve Losh <steve@stevelosh.com>
# Date 1460385647 0
# Node ID ef04c7b3d0b82611e81e8beea2d0f29482050606
# Parent  d4b651730553809ae36aa44c84fb8fba07228ea0
Problem 10, Miller-Rabin, etc

diff -r d4b651730553 -r ef04c7b3d0b8 .lispwords
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/.lispwords	Mon Apr 11 14:40:47 2016 +0000
@@ -0,0 +1,1 @@
+(1 repeat)
diff -r d4b651730553 -r ef04c7b3d0b8 Makefile
--- a/Makefile	Mon Apr 11 14:40:32 2016 +0000
+++ b/Makefile	Mon Apr 11 14:40:47 2016 +0000
@@ -1,7 +1,7 @@
 .PHONY: test
 
 quickutils.lisp: make-utilities.lisp
-	sbcl --noinform --load make-utilities.lisp  --eval '(quit)'
+	sbcl-rlwrap --noinform --load make-utilities.lisp  --eval '(quit)'
 
 test:
-	sbcl --noinform --load run-tests.lisp  --eval '(quit)'
+	sbcl-rlwrap --noinform --load run-tests.lisp  --eval '(quit)'
diff -r d4b651730553 -r ef04c7b3d0b8 euler.lisp
--- a/euler.lisp	Mon Apr 11 14:40:32 2016 +0000
+++ b/euler.lisp	Mon Apr 11 14:40:47 2016 +0000
@@ -145,6 +145,11 @@
                       :when (pythagorean-triplet-p a b c)
                       :do (return-from search (* a b c)))))))
 
+(defun problem-10 ()
+  (loop :for p :in (primes-below 2000000)
+        :sum p))
+
+
 ;;;; Tests
 (def-suite :euler)
 (in-suite :euler)
@@ -158,6 +163,7 @@
 (test p7 (is (= 104743 (problem-7))))
 (test p8 (is (= 23514624000 (problem-8))))
 (test p9 (is (= 31875000 (problem-9))))
+(test p10 (is (= 142913828922 (problem-10))))
 
 
 ; (run! :euler)
diff -r d4b651730553 -r ef04c7b3d0b8 package.lisp
--- a/package.lisp	Mon Apr 11 14:40:32 2016 +0000
+++ b/package.lisp	Mon Apr 11 14:40:47 2016 +0000
@@ -2,6 +2,7 @@
   (:use #:cl #:euler.quickutils)
   (:export
     #:random-exclusive
+    #:repeat
     #:dividesp))
 
 (defpackage #:euler
diff -r d4b651730553 -r ef04c7b3d0b8 primes.lisp
--- a/primes.lisp	Mon Apr 11 14:40:32 2016 +0000
+++ b/primes.lisp	Mon Apr 11 14:40:47 2016 +0000
@@ -1,5 +1,11 @@
 (in-package #:euler)
 
+(define-constant carmichael-numbers ; from oeis
+  '(561 1105 1729 2465 2821 6601 8911 10585 15841 29341 41041 46657 52633
+    62745 63973 75361 101101 115921 126217 162401 172081 188461 252601 278545
+    294409 314821 334153 340561 399001 410041 449065 488881 512461)
+  :test #'equal)
+
 (defun prime-factorization (n)
   "Return the prime factors of `n`."
   ;; from http://www.geeksforgeeks.org/print-all-prime-factors-of-a-given-number/
@@ -16,57 +22,142 @@
     (nreverse result)))
 
 
-(defun fermat-prime-p (p &optional (tests 10))
-  "Return whether `p` might be prime.
+(defun expmod (base exp m)
+  "Return base^exp % m quickly."
+  ;; From SICP and
+  ;; https://rosettacode.org/wiki/Miller%E2%80%93Rabin_primality_test#Common_Lisp
+  ;;
+  ;; We want to avoid bignums as much as possible.  This computes (base^exp % m)
+  ;; without having to deal with huge numbers by taking advantage of the fact
+  ;; that:
+  ;;
+  ;;     (x * y) % m
+  ;;
+  ;; is equivalent to:
+  ;;
+  ;;     ((x % m) * (y % m)) % m
+  ;;
+  ;; So for the cases where `exp` is even, we can split base^exp into an x and
+  ;; y both equal to base^(exp/2) and use the above trick to handle them
+  ;; separately.  Even better, we can just compute it once and square it.
+  ;;
+  ;; We also make it tail recursive by keeping a running accumulator:
+  ;;
+  ;;    base^exp * acc
+  (labels
+      ((recur (base exp acc)
+         (cond
+           ((zerop exp) acc)
+           ((evenp exp)
+            (recur (rem (square base) m)
+                   (/ exp 2)
+                   acc))
+           (t
+            (recur base
+                   (1- exp)
+                   (rem (* base acc) m))))))
+    (recur base exp 1)))
+
+
+(defun fermat-prime-p (n &optional (tests 10))
+  "Return whether `n` might be prime.
 
   Checks `tests` times.
 
-  If `t` is returned, `p` might be prime.  If `nil` is returned it is definitely
+  If `t` is returned, `n` might be prime.  If `nil` is returned it is definitely
   composite.
 
   "
-  (if (or (zerop p) (= 1 p))
-    nil
-    (flet ((fermat-check (a)
-             (= 1 (mod (expt a (1- p)) p))))
-      (loop :repeat tests
-            :when (not (fermat-check (random-exclusive 0 p)))
-            :do (return nil)
-            :finally (return t)))))
+  (flet ((fermat-check (a)
+           (= (expmod a n n) a)))
+    (loop :repeat tests
+          :when (not (fermat-check (random-exclusive 0 n)))
+          :do (return nil)
+          :finally (return t))))
+
+(defun factor-out (n factor)
+  "Factor the all the `factor`s out of `n`.
+
+  Turns `n` into:
+
+    factor^e * d
+
+  where `d` is no longer divisible by `n`, and returns `e` and `d`.
+
+  "
+  (loop :for d = n :then (/ d factor)
+        :for e = 0 :then (1+ e)
+        :while (dividesp d factor)
+        :finally (return (values e d))))
+
+(defun miller-rabin-prime-p (n &optional (k 10))
+  "Return whether `n` might be prime.
+
+  If `t` is returned, `n` is probably prime.
+  If `nil` is returned, `n` is definitely composite.
 
-(defun brute-force-prime-p (p)
-  "Return (slowly) whether `p` is prime."
+  "
+  ;; https://rosettacode.org/wiki/Miller%E2%80%93Rabin_primality_test#Common_Lisp
   (cond
-    ((or (= p 0) (= p 1)) nil)
-    ((= p 2) t)
-    ((evenp p) nil)
-    (t (loop :for divisor :from 3 :to (sqrt p)
-          :when (dividesp p divisor)
+    ((< n 2) nil)
+    ((< n 4) t)
+    ((evenp n) nil)
+    (t (multiple-value-bind (r d)
+           (factor-out (1- n) 2)
+         (flet ((strong-liar-p (a)
+                  (let ((x (expmod a d n)))
+                    (or (= x 1)
+                        (loop :repeat r
+                              :for y = x :then (expmod y 2 n)
+                              :when (= y (1- n))
+                              :do (return t))))))
+           (loop :repeat k
+                 :for a = (random-exclusive 1 (1- n))
+                 :always (strong-liar-p a)))))))
+
+(defun brute-force-prime-p (n)
+  "Return (slowly) whether `n` is prime."
+  (cond
+    ((or (= n 0) (= n 1)) nil)
+    ((= n 2) t)
+    ((evenp n) nil)
+    (t (loop :for divisor :from 3 :to (sqrt n)
+          :when (dividesp n divisor)
           :do (return nil)
           :finally (return t)))))
 
+(defun primep (n)
+  "Return (less slowly) whether `n` is prime."
+  (cond
+    ;; short-circuit a few edge/common cases
+    ((< n 2) nil)
+    ((= n 2) t)
+    ((evenp n) nil)
+    ((< n 100000) (brute-force-prime-p n))
+    (t (miller-rabin-prime-p n))))
+
 
 (defun primes-below (n)
-  "Return (slowly) the prime numbers less than `n`."
+  "Return the prime numbers less than `n`."
   (cond
     ((<= n 2) (list))
     ((= n 3) (list 2))
     (t (cons 2 (loop :for i :from 3 :by 2 :below n
-                     :when (brute-force-prime-p i)
+                     :when (primep i)
                      :collect i)))))
 
 (defun primes-upto (n)
-  "Return (slowly) the prime numbers less than or equal to `n`."
+  "Return the prime numbers less than or equal to `n`."
   (primes-below (1+ n)))
 
 
 (defun nth-prime (n)
-  "Return (slowly) the `n`th prime number."
+  "Return the `n`th prime number."
   (if (= n 1)
     2
     (loop :with seen = 1
           :for i :from 3 :by 2
-          :when (brute-force-prime-p i)
+          :when (primep i)
           :do (progn
                 (incf seen)
                 (when (= seen n)
diff -r d4b651730553 -r ef04c7b3d0b8 utils.lisp
--- a/utils.lisp	Mon Apr 11 14:40:32 2016 +0000
+++ b/utils.lisp	Mon Apr 11 14:40:47 2016 +0000
@@ -8,3 +8,8 @@
   "Return whether `n` is evenly divisible by `divisor`."
   (zerop (mod n divisor)))
 
+
+(defmacro repeat (n &body body)
+  "Repeat `body` `n` times."
+  `(dotimes (,(gensym) ,n)
+     ,@body))