832 lines
28 KiB
Common Lisp
832 lines
28 KiB
Common Lisp
;;; -*- Package: USER; Mode: LISP; Syntax: Common-lisp -*-
|
|
;;(in-package "USER")
|
|
|
|
;; Author: Jason Eisner <jason@cs.jhu.edu>, 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* "~/.scratch/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))
|
|
|