0ac280dfa75f

primes -> euler
[view raw] [browse files]
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))))))