dc2c9931b634

SICP streams
[view raw] [browse files]
author Steve Losh <steve@stevelosh.com>
date Sun, 04 Sep 2016 22:20:37 +0000
parents 778814a3ff72
children ba6e58d41e02
branches/tags (none)
files package.lisp sand.asd src/primes.lisp src/streams.lisp

Changes

--- a/package.lisp	Tue Aug 30 23:54:54 2016 +0000
+++ b/package.lisp	Sun Sep 04 22:20:37 2016 +0000
@@ -13,6 +13,18 @@
   (:shadowing-import-from #:cl-arrows
     #:->))
 
+(defpackage #:sand.primes
+  (:use
+    #:cl
+    #:cl-arrows
+    #:losh
+    #:iterate
+    #:sand.graphviz
+    #:sand.quickutils
+    #:sand.utils)
+  (:export
+    #:primep))
+
 (defpackage #:sand.random-numbers
   (:use
     #:cl
@@ -126,6 +138,18 @@
   (:export
     ))
 
+(defpackage #:sand.streams
+  (:use
+    #:cl
+    #:cl-arrows
+    #:losh
+    #:iterate
+    #:sand.primes
+    #:sand.quickutils
+    #:sand.utils)
+  (:export
+    ))
+
 #+sbcl
 (defpackage #:sand.ffi
   (:use
--- a/sand.asd	Tue Aug 30 23:54:54 2016 +0000
+++ b/sand.asd	Sun Sep 04 22:20:37 2016 +0000
@@ -36,6 +36,7 @@
    (:module "src"
     :serial t
     :components ((:file "utils")
+                 (:file "primes")
                  (:file "graphviz")
                  (:file "random-numbers")
                  (:file "ascii")
@@ -44,6 +45,7 @@
                  #+sbcl (:file "ffi")
                  (:file "binary-decision-diagrams")
                  (:file "huffman-trees")
+                 (:file "streams")
                  (:file "color-difference")
                  (:module "terrain"
                   :serial t
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/primes.lisp	Sun Sep 04 22:20:37 2016 +0000
@@ -0,0 +1,164 @@
+(in-package #:sand.primes)
+
+(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/
+  (let ((result (list)))
+    (while (evenp n) ; handle 2, the only even prime factor
+      (push 2 result)
+      (setf n (/ n 2)))
+    (loop :for i :from 3 :to (sqrt n) :by 2 ; handle odd (prime) divisors
+          :do (while (dividesp n i)
+                (push i result)
+                (setf n (/ n i))))
+    (when (> n 2) ; final check in case we ended up with a prime
+      (push n result))
+    (nreverse result)))
+
+
+(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, `n` might be prime.  If `nil` is returned it is definitely
+  composite.
+
+  "
+  (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.
+
+  "
+  ;; https://rosettacode.org/wiki/Miller%E2%80%93Rabin_primality_test#Common_Lisp
+  (cond
+    ((< 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 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 (primep i)
+                     :collect i)))))
+
+(defun primes-upto (n)
+  "Return the prime numbers less than or equal to `n`."
+  (primes-below (1+ n)))
+
+
+(defun nth-prime (n)
+  "Return the `n`th prime number."
+  (if (= n 1)
+    2
+    (loop :with seen = 1
+          :for i :from 3 :by 2
+          :when (primep i)
+          :do (progn
+                (incf seen)
+                (when (= seen n)
+                  (return i))))))
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/streams.lisp	Sun Sep 04 22:20:37 2016 +0000
@@ -0,0 +1,98 @@
+(in-package #:sand.streams)
+
+
+;;;; Streams from SICP
+
+;;;; Delay/Force --------------------------------------------------------------
+(defun memoize (thunk)
+  (let ((done nil) (result nil))
+    (lambda ()
+      (if done
+        result
+        (setf done t
+              result (funcall thunk))))))
+
+(defmacro delay (&body body)
+  `(memoize (lambda () ,@body)))
+
+(defun force (delay)
+  (funcall delay))
+
+
+;;;; Basic Streams ------------------------------------------------------------
+(defun scar (stream)
+  (car stream))
+
+(defun scdr (stream)
+  (if (cdr stream)
+    (force (cdr stream))
+    nil))
+
+(defmacro scons (car cdr)
+  `(cons ,car (delay ,cdr)))
+
+(defun snullp (stream)
+  (null stream))
+
+(defun empty-stream ()
+  nil)
+
+
+;;;; Stream Operations --------------------------------------------------------
+(defun stream-nth (n stream)
+  (if (zerop n)
+    (scar stream)
+    (stream-nth (1- n) (scdr stream))))
+
+(defun stream-map (function stream)
+  (if (snullp stream)
+    (empty-stream)
+    (scons (funcall function (scar stream))
+           (stream-map function (scdr stream)))))
+
+(defun stream-do (function stream)
+  (if (snullp stream)
+    (values)
+    (progn (funcall function (scar stream))
+           (stream-do function (scdr stream)))))
+
+(defun print-stream (stream)
+  (stream-do #'pr stream))
+
+(defun stream-range (low high)
+  (if (> low high)
+    (empty-stream)
+    (scons low (stream-range (1+ low) high))))
+
+(defun stream-filter (function stream)
+  (if (snullp stream)
+    (empty-stream)
+    (if (funcall function (scar stream))
+      (scons (scar stream) (stream-filter function (scdr stream)))
+      (stream-filter function (scdr stream)))))
+
+(defun stream-take (n stream)
+  (if (or (snullp stream) (zerop n))
+    (empty-stream)
+    (scons (scar stream) (stream-take (1- n) (scdr stream)))))
+
+(defun stream-drop (n stream)
+  (if (or (snullp stream) (zerop n))
+    (scdr stream)
+    (stream-drop (1- n) (scdr stream))))
+
+
+;;;; Scratch ------------------------------------------------------------------
+(trace primep)
+
+(->> (stream-range 10000 1000000)
+  (stream-filter #'primep)
+  (stream-take 2)
+  print-stream)
+
+(stream-nth 2 (stream-filter #'primep (stream-range 1 10)))
+
+(untrace primep)
+
+
+