0ac280dfa75f
primes -> euler
author | Steve Losh <steve@stevelosh.com> |
---|---|
date | Sun, 13 Aug 2017 00:06:58 -0400 |
parents | 60b451e2a6eb |
children | 6494ca97e78b |
branches/tags | (none) |
files | package.lisp sand.asd src/primes.lisp |
Changes
--- a/package.lisp Mon Aug 07 22:16:41 2017 -0400 +++ b/package.lisp Sun Aug 13 00:06:58 2017 -0400 @@ -7,16 +7,6 @@ (:export :average4)) -(defpackage :sand.primes - (:use - :cl - :losh - :iterate - :sand.quickutils - :sand.utils) - (:export - :primep)) - (defpackage :sand.random-numbers (:use :cl @@ -168,7 +158,6 @@ :cl :losh :iterate - :sand.primes :sand.quickutils :sand.utils) (:export
--- a/sand.asd Mon Aug 07 22:16:41 2017 -0400 +++ b/sand.asd Sun Aug 13 00:06:58 2017 -0400 @@ -48,7 +48,6 @@ (:module "src" :serial t :components ((:file "utils") - (:file "primes") (:file "graphs") (:file "graphviz") (:file "hanoi")
--- a/src/primes.lisp Mon Aug 07 22:16:41 2017 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,164 +0,0 @@ -(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))))))