;;; -*- Package: USER; Mode: LISP; Syntax: Common-lisp -*- (in-package "USER") ;; Author: Jason Eisner , December 1993 ;; ;; Use this code at your own risk; please give appropriate credit. ;; See primes.pdf for explanation and discussion. ;; ;; Because Common Lisp has bignum support built in, it is a handy ;; language for experimenting with very large prime numbers. ;;; ***************** ;;; User macros to reload this file. ;;; ***************** (defparameter *workfile* "~/lisp/primes") (defmacro L (&optional (filename *workfile*)) `(load ,filename)) (defmacro CL (&optional (filename *workfile*)) `(progn (compile-file ,filename) (load ,filename))) ;;; *************************** ;;; Useful macros & operations. ;;; *************************** ;; Logical implication. (defmacro implies (antecedent consequent) `(or (not ,antecedent) ,consequent)) ;; Accesses the nth value (counting from 0) ;; of a form that returns multiple values. ;; ;; This macro is supposed to be built into Common ;; Lisp, but it was left out of this implementation. (defmacro nth-value (n form) `(nth ,n (multiple-value-list ,form))) ;;; ***************************** ;;; Basic mathematical operations ;;; ***************************** (defmacro divides (a b) ; returns T if a | b, NIL otherwise `(zerop (mod ,b ,a))) (defun square (x) (* x x)) ;; returns a list of bits constituting the binary expansion of n, ;; from most to least significant. T indicates a 1 bit, NIL a 0 ;; bit. Example: (binary-expansion 11) ==> (T NIL T T). (defun binary-expansion (n) (loop with answer = nil if (evenp n) do (push nil answer) else do (push t answer) and do (decf n) do (setf n (/ n 2)) until (zerop n) finally (return answer))) ;; Quickly finds b^e (mod m). (defun exptmod (b e m) (if (zerop e) 1 (mod (if (evenp e) (square (exptmod b (/ e 2) m)) (* b (exptmod b (1- e) m))) m))) ;; A bit faster than exptmod. It's recursive rather than iterative, and ;; the caller must provide a precomputed binary expansion of the exponent. ;; ;; (This is advantageous because we may have to reuse the same exponent ;; many times.) (defun exptmod-fast (b binary-e m) (loop with prod = 1 for bit in binary-e do (setf prod (square prod)) if bit do (setf prod (* b prod)) do (setf prod (mod prod m)) finally (return prod))) ;; Integer part of the cube root---avoids floating-point goofups. ;; (See remarks at pi-Meissel.) (defun icuberoot (n) (floor (expt (+ 0.5 n) 1/3))) ;; Integration by the trapezoidal rule. (defun integrate (fn lower upper dx) ;; correct dx so that we have a whole # of intervals (let* ((n (round (/ (abs (- upper lower)) dx))) (dx (/ (- upper lower) n))) (* dx (+ (loop for i from 1 to (- n 1) for x from (+ lower dx) by dx sum (funcall fn x)) (* 1/2 (funcall fn lower)) (* 1/2 (funcall fn upper)))))) ;;; ***************** ;;; Primality Testing ;;; ***************** ;; Almost (but not quite) the stupidest primality test conceivable. (defun prime?-slow (n) (loop for i from 2 to (isqrt n) never (divides i n))) ;; Returns a list of all primes in [2,n], in increasing order. ;; Uses the stupid method prime?-slow to test numbers individually. ;; Doesn't amortize the work at all. (defun primes-tested-slow (n) (loop for i from 2 to n when (prime?-slow i) collect i)) (defun pi-tested-slow (n) (length (primes-tested-slow n))) ;; Returns an array holding the primes in [2,n], in increasing order. ;; ;; The array has a fill pointer which registers its "active ;; length." The active elements of the array are taken to be 0 ;; up to but not including the fill pointer. (The built-in length ;; function respects fill pointers.) ;; ;; An unproved theorem from Part II Number Theory says that ;; pi(n) is about n/log(n). We initially allow for 1.2 times this ;; many elements, but we use an adjustable array that will automatically ;; get bigger if necessary. ;; ;; Notes: ;; 1. Takes no square roots. (a significant linear speedup) ;; 2. Stops as soon as it knows a candidate is composite ;; (divides by primes in the order 2, 3, 5, ... sqrt(candidate)). (defun primes-tested (n) (let ((primes (make-array (floor (* 1.2 (/ n (log n)))) :fill-pointer 0 ; no primes in array yet :adjustable t)) ; can get more space if necessary (n-divisors 0) ; number of primes being used as divisors (i-biggest-testable 0)) ; largest integer testable using only ; the first n-divisors primes as divisors (loop for i from 2 to n ;; First be very careful about the case where we're not sure ;; we have enough divisors to test i for primality. ;; ;; We can certainly test integers up to p^2 by using only ;; primes up to p. ;; ;; Indeed, we can do slightly better most of the time: ;; we can test integers up to q^2-1, where q is the least ;; prime exceeding p. (if one has been found!) ;; ;; The array holds all primes in [2,i-1], which is ;; always enough to test i ... so if n-divisors = (length primes), ;; we can proceed with a clear conscience, even if ;; i > i-biggest-testable. do (if (and (> i i-biggest-testable) (< n-divisors (length primes))) (progn (incf n-divisors) (setf i-biggest-testable (if (< n-divisors (length primes)) (1- (square (svref primes n-divisors))) (square (svref primes (1- n-divisors))))))) ;; Now see if i is a prime, and if so, add it to the array. (if (loop for index from 0 to (1- n-divisors) never (divides (svref primes index) i)) ;; stops ASAP (vector-push-extend i primes))) primes)) ;; Finds pi(n) by using primes-tested. ;; ;; We only use primes-tested to generate the primes up to sqrt(n). ;; (We must test larger numbers for primality, of course, but we don't need ;; to store them!) (defun pi-tested (n) (let* ((sqrt[n] (max 2 (isqrt n))) (divisors (primes-tested sqrt[n]))) (+ (length divisors) (loop for i from (1+ sqrt[n]) to n count (loop for index from 0 to (1- (length divisors)) never (divides (svref divisors index) i)))))) ;;; ******************************* ;;; Probabilistic primality testing ;;; ******************************* (defvar max-error-prob) (defvar spot-checks) (setf max-error-prob 1E-10) (setf spot-checks (ceiling (- (/ (log max-error-prob) (log 4))))) (defun Fermat-pseudoprime? (n base) (= 1 (exptmod base (1- n) n))) ;; A fast version: the caller must supply the binary-expansion of n-1. (defmacro Fermat-pseudoprime?-fast (n base binary-e) `(= 1 (exptmod-fast ,base ,binary-e ,n))) ;; Test whether n is a strong pseudoprime to the given base. ;; n-1 = t*(2^s), for t odd. ;; binary-t is the binary-expansion of t. (defun strong-pseudoprime? (n base binary-t s) ;; where n-1 = t*(2^s) (let ((base^t (exptmod-fast base binary-t n))) (or (= 1 base^t) (loop with prod = base^t for r from 0 to (1- s) thereis (= (1- n) prod) do (setf prod (mod (square prod) n)))))) ;; Pick an integer from [1, n) that is coprime to n. (defun random-base (n) (loop with base do (setf base (1+ (random (1- n)))) until (= 1 (gcd base n)) finally (return base))) ;; The Miller-Rabin test. Relies on the fact that if n is prime, ;; it is an strong pseudoprime to every base in [1, n), whereas if ;; it's composite, it is a strong pseudoprime to at most 1/4 of these ;; bases. ;; ;; We do enough random spot checks to bring the probability of an ;; error below max-error-prob. (defun prime?-probably (n) (loop with tee = (1- n) with s = 0 with binary-t initially (loop while (evenp tee) do (setf tee (/ tee 2)) do (incf s)) ;; now n-1 = 2^s * t with t odd initially (setf binary-t (binary-expansion tee)) for i from 1 to spot-checks always (strong-pseudoprime? n (random-base n) binary-t s))) ;; Perfect testing of large primes, if the Generalized Riemann Hypothesis is true. ;; We only have to test against bases up to log n. (defun prime?-GRH (n) (loop with tee = (1- n) with s = 0 with binary-t initially (loop while (evenp tee) do (setf tee (/ tee 2)) do (incf s)) ;; now n-1 = 2^s * t with t odd initially (setf binary-t (binary-expansion tee)) for base from 1 to (floor (log n)) when (= 1 (gcd base n)) always (strong-pseudoprime? n base binary-t s))) (defun pi-tested-probably (n) (loop for i from 2 to n count (prime?-probably i))) (defun pi-tested-GRH (n) (loop for i from 2 to n count (prime?-probably i))) ;; First prime after n. To speed this up, we can look ;; in an arithmetic progression starting at n -- for example, ;; at the odd numbers after n. (defun next-prime (n &optional (step 1)) (assert (= (gcd n step) 1)) (loop for i from n by step when (prime?-probably i) do (return i))) ;;; **************** ;;; Sieving Methods. ;;; **************** ;; Sieve the numbers from 1 to n, and return the resulting primes in a list. ;; ;; As a second value, we return an array with elements from 1 to n; the ;; primes are the elements with non-nil entries. (defun primes-sieved (n) (loop with sieve = (make-array (1+ n) :initial-element T) ; all elements initially flagged prime with divisor = 2 while (<= (square divisor) n) ; for all prime divisors < sqrt(n) ;; Strike out multiples of the divisor. do (loop for i from (+ divisor divisor) to n by divisor do (setf (svref sieve i) nil)) ;; Find the next non-prime divisor, if any. do (loop do (incf divisor) until (or (svref sieve divisor) (> divisor n))) ;; At the end, collect up the remaining primes and return them. ;; Return the sieve itself as a second value. finally (setf (svref sieve 0) nil ; these aren't primes either (svref sieve 1) nil) (return (values (loop for i from 2 to n when (svref sieve i) collect i) sieve)))) (defun pi-sieved (n) (length (primes-sieved n))) (defmacro link (i j) `(setf (svref next ,i) ,j (svref prev ,j) ,i)) ;; A LINEAR algorithm. ;; ;; We keep a doubly linked list of the integers that might be ;; prime. As multiples are struck out, they're removed from ;; this list. ;; ;; The list is implemented as a pair of arrays, prev and next. ;; If i is an integer on the list, and 1 <= k <= sqrt(n), then ;; ;; Prev[i] = the next smallest list element (or 1 if none such) ;; Next[i] = the next largest list element (or n+1 if none such) ;; S[k] = the largest list element <= n/k ;; Sinv[i] = {k: S(k) = i} (stored as a list) ;; ;; Note: 2 is always the smallest list element; S[1] is always the largest ;; list element. Either fact would enable us to enumerate the primes in ;; the sieve (we don't currently). (defun pi-sieved-linear (n) (loop with list-size = (1- n) ;; all nos. in [2,n] are initially thought prime with sqrt-n = (isqrt n) with prev = (make-array (+ n 2)) with next = (make-array (+ n 2)) with Sinv = (make-array (+ n 1)) with S = (make-array (+ sqrt-n 1)) initially (loop for i from 1 to n ; set up linked list do (link i (1+ i))) (loop for k from 1 to sqrt-n ; set up S and Sinv for S[k] = (floor (/ n k)) do (setf (svref S k) S[k]) (push k (svref Sinv S[k]))) with divisor = 2 while (<= divisor sqrt-n) ; for all prime divisors <= sqrt(n) ;; Multiplier starts at S[divisor] -- the biggest potential prime that, ;; when multiplied by divisor, still yields a number <= n. do (loop with multiplier = (svref S divisor) while (>= multiplier divisor) ; needn't go smaller for composite = (* multiplier divisor) for predecessor = (svref prev composite) for successor = (svref next composite) ;; now strike out the composite number. do (link predecessor successor) (decf list-size) (loop for k in (svref Sinv composite) do (setf (svref S k) predecessor) (push k (svref Sinv predecessor))) ;; move on to the next multiplier. do (setf multiplier (svref prev multiplier))) do (setf divisor (svref next divisor)) finally (return list-size))) ;; A variant of pi-sieved-linear, where we start with ;; all multiples of 2 and 3 already struck out. ;; ;; This optimization doesn't change the complexity of ;; the algorithm -- it just saves some cycles, by ;; eliminating two passes and by shortening the ;; initialization of the arrays. (defun pi-sieved-linear-fast (n) (loop with list-size = (+ 1 (- n (floor (/ n 2)) (floor (/ n 3))) (floor (/ n 6))) with sqrt-n = (isqrt n) with prev = (make-array (+ n 7)) with next = (make-array (+ n 7)) with Sinv = (make-array (+ n 1)) with S = (make-array (+ sqrt-n 1)) initially (link 2 3) (link 3 5) (link 5 7) (loop with i = 7 while (<= i n) for j = (+ i 6) do (link i (+ i 4)) (link (+ i 4) j) do (setf i j)) (loop for k from 1 to sqrt-n ; set up S and Sinv for S[k] = (floor (/ n k)) ; might not be in the list do (loop until (svref prev S[k]) do (decf S[k])) (setf (svref S k) S[k]) (push k (svref Sinv S[k]))) with divisor = 5 ;; ** From here on, the code is exactly the same as in pi-sieved-linear. ** while (<= divisor sqrt-n) ; for all prime divisors <= sqrt(n) ;; Multiplier starts at S[divisor] -- the biggest potential prime that, ;; when multiplied by divisor, still yields a number <= n. do (loop with multiplier = (svref S divisor) while (>= multiplier divisor) ; needn't go smaller for composite = (* multiplier divisor) for predecessor = (svref prev composite) for successor = (svref next composite) ;; now strike out the composite number. do (link predecessor successor) (decf list-size) (loop for k in (svref Sinv composite) do (setf (svref S k) predecessor) (push k (svref Sinv predecessor))) ;; move on to the next multiplier. do (setf multiplier (svref prev multiplier))) do (setf divisor (svref next divisor)) finally (return list-size))) ;;; ****************** ;;; Legendre's formula ;;; ****************** ;; Computes phi(n, a), i.e., the number of integers in [1,n] ;; not divisible by any of the first a primes. ;; Primelist is a list of the first a primes, in REVERSE order. (defun phi (n a primelist) (cond ((zerop a) (floor n)) ((< n 1) 0) (t (- (phi n (1- a) (rest primelist)) (phi (/ n (first primelist)) (1- a) (rest primelist)))))) ;; Legendre's formula. (defun pi-Legendre (n) (let* ((sqrt[n] (isqrt n)) (pr (primes-sieved sqrt[n])) (a (length pr))) (+ a -1 (phi n a (reverse pr))))) ;; Returns an array holding the values of phi(t, k) for ALL 0 <= t <= m, ;; where m is the product of the first k primes. ;; ;; The first k primes must be provided. We use a sieving method to ;; determine the values efficiently, more or less as suggested in Riesel. (defun phi-array (m primes) (let ((sieve (make-array (1+ m) :initial-element t))) ;; Cross out the integers that are divisible by at least one ;; of the given primes. (loop for p in primes do (loop for i from 0 to m by p do (setf (svref sieve i) nil))) ;; Now change the entries of the sieve from logical flags to ;; numbers: at each entry, we count the number of previous entries ;; that are indivisible by all the primes, i.e., still have t entries. (loop for i from 0 to m count (svref sieve i) into phi do (setf (svref sieve i) phi)) ;; Return the result. sieve)) (defvar *k* nil) ;; Value of k from our last call to pi-Legendre-fast (defvar *phi-array* nil) ;; The phi-array we built then -- we may be able to ;; reuse it! (defvar *m* nil) ;; The value of m we used then ;; Set the above global variables to be appropriate to a new ;; value of k. If the value of k hasn't changed, leave the ;; old values alone (they're still appropriate.) (defun initialize-phi-array (k) (unless (eq k *k*) (let ((primes (loop ;; the first k primes (k is very small) for i from 2 when (prime?-slow i) collect i and count T into count until (= count k)))) (setf *k* k *m* (apply #'* primes) *phi-array* (phi-array *m* primes)))) (values)) ;; Given a phi-array of the correct form, computes phi quickly. ;; (See documentation at phi-array.) ;; ;; The optimization as described in the project handout uses ;; the Euler phi function of m. Since we know m's factors, this ;; would be easy to find. However, since it's exactly phi(m,k), ;; it was just as convenient to compute it as part of the phi-array. (defun phi-fast (n a primelist &optional (k *k*) (m *m*) (phi-array *phi-array*)) (cond ((= a k) (multiple-value-bind (s tee) (floor n m) ;; so n = s*m + t, 0 <= t < m (+ (* s (svref phi-array m)) (svref phi-array tee)))) ((zerop a) ;; another termination condition, in case we're (floor n)) ;; initially called with (a < k) ((< n 1) 0) (t (- (phi-fast n (1- a) (rest primelist) k m phi-array) (phi-fast (floor n (first primelist)) (1- a) (rest primelist)))))) (defun pi-Legendre-fast (n &optional (k 6)) ;; k=0 gives vanilla pi-Legendre (let* ((sqrt[n] (isqrt n)) (primes (primes-sieved sqrt[n])) (a (length primes))) ;; Recompute the global variables unless they're appropriate from the last call. (initialize-phi-array k) (+ a -1 (phi-fast n a (reverse primes))))) ;;; ***************** ;;; Meissel's formula ;;; ***************** ;; By scanning through a sieve, in linear time, creates an array ;; of length n with a very special form: ;; If n is prime, the nth element is pi(n). ;; If n is not prime, the nth element is a negative value giving ;; the previous prime -- or 0 if none. (defun pi-array (n) (loop with array = (nth-value 1 (primes-sieved n)) with current-pi = 0 with last-prime = 0 for i from 0 to n if (svref array i) do (incf current-pi) (setf (svref array i) current-pi last-prime i ) else do (setf (svref array i) (- last-prime)) finally (return array))) ;; Finds pi(n) by looking in a pi-array. (defun lookup-pi (n pi-array) (let ((entry (svref pi-array n))) (if (> entry 0) entry (svref pi-array (- entry))))) ;; Finds the greatest prime <= n, by looking ;; in a pi-array. Returns 0 if none such. (defun lookup-prime-floor (n pi-array) (let ((entry (svref pi-array n))) (if (> entry 0) n (- entry)))) ;; Computes pi from Meissel's formula. ;; ;; Here pi-array has the form produced by the pi-array routine above. ;; It must have length AT LEAST sqrt(x). ;; ;; Important: We use an icuberoot function I've defined ;; to get the integer part of the cube root. (This is analagous ;; to isqrt.) If we don't do this, our answer is occasionally ;; off by one. For example, pi(343) would be miscomputed, because ;; this computer thinks the cube root of 343 is 6.999999999999999. ;; So would any pi(n) that is computed in terms of pi(343). (defun pi-Meissel (x &optional (pi-array (pi-array (isqrt x)))) (if (<= x 1) 0 (let* ((cuberoot (icuberoot x)) (sqrt (isqrt x)) (c (lookup-pi cuberoot pi-array)) (b (lookup-pi sqrt pi-array))) (initialize-phi-array 6) ;; this call usually does nothing (+ (phi-fast x c (loop ; find the first c primes, in reverse order with p = (lookup-prime-floor cuberoot pi-array) until (zerop p) collect p do (setf p (lookup-prime-floor (1- p) pi-array)))) (* 1/2 (+ b c -2) (- b c -1)) (- (loop with p = (lookup-prime-floor sqrt pi-array) while (> p cuberoot) sum (pi-Meissel (floor x p) pi-array) do (setf p (lookup-prime-floor (1- p) pi-array)))))))) ;;; ******************* ;;; Approximating pi(x) ;;; ******************* ;; Computes Li(x) from its formal definition. (defun Li-integrand-slow (s) (/ (log s))) (defun Li-slow (x &optional (ds 1.0)) (+ 1.045 (integrate #'Li-integrand-slow 2 x ds))) ;; Faster computation using a change of variable -- in the integration ;; above, we're using narrow intervals on an integrand that gets flatter ;; and flatter. (defun Li-integrand (u) (/ (exp u) u)) (defun Li (x &optional (du 0.001)) (+ 1.045 (integrate #'Li-integrand (log 2) (log x) du))) ;; The approximation to pi(x) given by the prime number theorem. (defun log-approx (n) (/ n (log n))) ;;; *********************** ;;; Tabulating our results. ;;; *********************** ;; Calls the given function on n = 10, 20, 30, ... 90, 100, 200, ... ;; up to upperbound (if supplied), throwing away any results. ;; ;; We use this to print a table. (defun drive-tabulation (function &optional upperbound (factor 10)) (loop with interval = factor with arg = interval while (implies upperbound (<= arg upperbound)) for total-calls from 0 do (funcall function arg) (if (and (divides (1- factor) total-calls) (> total-calls 0)) (setf interval (* interval factor))) (incf arg interval)) (values)) ;; Calls all the functions on the argument, and collect the results ;; into a list. (defun eval-fns (fns arg) (loop for fn in fns collect (funcall fn arg))) ;; Prints the print names of the list elements as column headers, splitting ;; them as readably as possible across several lines if necessary. (defun print-headers (headers) (loop with done? = nil with names = (loop for object in headers collect (string-capitalize (format nil "~A" object))) for print = (loop initially (setf done? t) for i from 0 to (1- (length names)) for name = (elt names i) for length = (length name) for cut-point = (if (<= length 10) length (1+ (or (position #\- name :end 10 :from-end t) 9))) for print-part = (subseq name 0 cut-point) for remaining = (subseq name cut-point) ;; to end collect print-part do (setf (elt names i) remaining) (when (string/= remaining "") (setf done? nil))) do (print-table-line " " print) until done?) (values)) ;; Prints the list as a table line. (defun print-table-line (row-label list) (format t "~&~A~{~,12T~A~}" row-label list)) ;; Tabulates the values of a function for 10, 20, 30, ... 90, 100, 200, ... ;; Here fns may be either a single function or a list of functions. ;; Their print names will be used as the headers, unless a list of other ;; headers (strings) is provided as a keyword argument. ;; ;; There are some other optional keyword arguments. The multiplicative factor ;; may be specified, and so may an upper bound that says how far to tabulate. ;; Finally, if the list upperbounds has an nth element, it serves as a ;; further upperbound on the nth function. (defun tabulate (fns &key upperbounds (upperbound (when upperbounds (apply #'max upperbounds))) (factor 10) headers) ;; If there's only a single function, make it a list of one function. ;; (Careful, because some function specifiers look like lists.) (unless (and (listp fns) (not (eq (car fns) 'function)) (not (eq (car fns) 'lambda))) (setf fns (list fns))) (print-headers (or headers fns)) (drive-tabulation #'(lambda (arg) (print-table-line arg (loop for fn in fns for i from 0 for upperbound = (nth i upperbounds) when (implies upperbound (<= arg upperbound)) collect (funcall fn arg) else collect " "))) upperbound factor)) ;; Tabulates all our methods as far as appropriate. (defun tabulate-pi-methods () (tabulate '(pi-tested-slow pi-tested pi-tested-probably pi-tested-GRH pi-sieved pi-sieved-linear pi-sieved-linear-fast pi-Legendre pi-Legendre-fast pi-Meissel) :upperbounds '( 100000 100000 10000 10000 500000 100000 100000 1000000 10000000 10000000))) ;; Tabulates approximations to pi. (defun tabulate-approximations () (print-headers '(" pi" " x/log x" " li" " log-ratio" " li-ratio" " log-diff" " li-diff")) (drive-tabulation #'(lambda (n) (format t "~&~A~{~,12T~11,3F~}" n (let ((pi-true (pi-Meissel n)) (pi-log (log-approx n)) (pi-li (Li n))) (list pi-true pi-log pi-li (/ pi-true pi-log) (/ pi-true pi-li) (- pi-true pi-log) (- pi-true pi-li))))) 10000000))