# HG changeset patch # User Steve Losh # Date 1473027637 0 # Node ID dc2c9931b634c9c5480e1d8f9c7ff5430da931dd # Parent 778814a3ff72633c3366c5b1ebdba9effa3fbd2f SICP streams diff -r 778814a3ff72 -r dc2c9931b634 package.lisp --- 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 diff -r 778814a3ff72 -r dc2c9931b634 sand.asd --- 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 diff -r 778814a3ff72 -r dc2c9931b634 src/primes.lisp --- /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)))))) diff -r 778814a3ff72 -r dc2c9931b634 src/streams.lisp --- /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) + + +