--- 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)
+
+
+