;$Id: PRIMES.LSP 0.5006 2000/12/31 10:59:00 rurban Rel $ -*-AutoLISP-*-
;;; Time-stamp: <2000-12-31 11:31:42 rurban>
;;; Copyright (c) 1998 by Reini Urban
;;; Available for free at http://xarch.tu-graz.ac.at/autocad/stdlib/
;;;
;;; Permission to use, copy, modify and distribute this software and its
;;; documentation for any purpose is hereby granted without fee, provided
;;; 1) that the above copyright notice appear in all copies,
;;; 2) that the copyright notice, this permission notice and the pointer
;;; where to download the source code for free appear in the
;;; supporting documentation of source code distributions,
;;; 3) that the name of Reini Urban not be used in advertising or
;;; publicity pertaining to distribution of the software and
;;; 4) that modifications without changing the defined function and
;;; symbol names may not be published, distributed nor copied
;;; without specific, written prior permission.
;;;
;;; No Warranty
;;; Reini Urban makes no representations about the suitability of this
;;; software for any purpose, without even the implied warranty of
;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. It is provided
;;; "as is" without express or implied warranty.
;;; See the full Disclaimer for all detailed warranty exclusions.
;;; --------------------------------------------------------------------
;;; Primes
;;;
;;; Note:
;;; This was seperated from STDMATH because it is only rarely used.
;;; It is not loaded with STDLIB and should be loaded with
;;; (std-require 'PRIMES)
;;;
;;; STATUS:
;;; No overflow checks are yet implemented for STD-PRIMES.
;;; No probability testing on larger numbers. Only 100% methods.
;;; STD-PRIMEP and STD-PRIME-FACTORIZATION are fast enough.
;;; STD-PRIMES not at all.
;;; Todo:
;;; (STD-PROBABLE-PRIME-P n) with small factor and fermat test
;;; as in rsaref/source/prime.c
;;;
;;; Algorithms:
;;; Eratosthenes's sieve with two global tables of primes numbers.
;;; std-%prime-sieve does no other consing but adding found
;;; primes to the global list.
;;; The initial list could be larger, so far it has only all
;;; primes up to 100. The 168 numbers up to 997 can be considered.
;;; Factorization by division
;;; std-%primep stores the last found divisor if it fails which
;;; is used by std-prime-factorization.
;;; Complexity:
;;; Prime testing is fastest, getting all numbers is slow but required
;;; up the root of the tested number for factorization and testing.
;;; Factorization is slow for large numbers.
;;; see test/PRIMETST.LSP for benchmarks and sample code.
;;; Numbers larger than *max-longint* cannot be calculated precisely
;;; without bignums.
;;; The reasonable ranges on todays hardware are
;;; <40000 for STD-PRIMES (all primes) ~10min
;;; practical performance limit
;;; <10^15 for STD-PRIMEP (test): ~20ms
;;; upper numerical bound
;;; <3.2^10 for STD-PRIME-FACTORIZATION ~140ms
;;; upper bound if primes<=40000
;;; Probabilistic methods for large numbers are not implemented,
;;; (STD-PROBABLE-PRIME-P n) and are not reasonable because of the lack
;;; of bignums and the needed accuracy. e.g. (std-primep 1e32) => T
;;; $Log: PRIMES.LSP $
;;; Revision 0.5006 2000/12/31 10:59:00 rurban
;;; Release, see Changes
;;;
;;; Revision 0.5004 2000/09/20 12:50:35 rurban
;;; 0.5004 release, see Changes
;;;
;;; Revision 0.5003 2000/07/11 15:59:52 rurban
;;; 0.5003 release, see Changes
;;;
;;; 2000-07-07 12:02:18 rurban
;;; added std-%(un)protect-assign wrapper
;;; Revision 0.5002 2000/06/19 08:31:18 rurban
;;; 0.5002 release, see Changes
;;;
;;; Revision 0.5000 2000/05/14 17:41:36 rurban
;;; 0.5 release
;;;
;;; 2000-04-26 18:19:49 rurban
;;; renumbered MSG ID's
;|
(STD-PRIMEP n) ; is it a prime number?
(STD-PRIMES n) ; list of prime numbers up to n
(STD-PRIME-FACTORIZATION n) ; unique list of multiplicants for n
Used globals:
*LASTDIV* the last found divisor in primality testing
*PRIMES* list of yet calculated primes (grows)
*SIEVE* part of *primes* up to the square-root of n
|;
;;; internal functions start with std-%
;;; ===================================================================73
;|#+ VL|;
(if std-%unprotect-assign
(std-%unprotect-assign
'(STD-PRIME-FACTORIZATION
STD-PRIMES
STD-PRIMEP
STD-PROBABLE-PRIME-P
)))
;|END #+ VL|;
;;; Globals:
(if (not *PRIMES*) (setq *PRIMES*
'(97 89 83 79 73 71 67 61 59 53 47 43 41 37 31 29 23 19 17 13 11 7
5 3 2)))
(if (not *SIEVE*) (setq *SIEVE* *PRIMES*))
;;; -------------------------------------------------------------------73
;;; Primes
;;; Factorization into prime divisors of descending order
;;; (STD-PRIME-FACTORIZATION 7) => (7)
;;; (STD-PRIME-FACTORIZATION 8) => (2 2 2)
;;; (STD-PRIME-FACTORIZATION 1118) => (43 13 2)
;;; (STD-PRIME-FACTORIZATION 1000000) => (5 5 5 5 5 5 2 2 2 2 2 2)
(defun STD-PRIME-FACTORIZATION (n / x)
(while (not (std-%prime-sieve n))
(setq x (cons *lastdiv* x)
n (/ n *lastdiv*)))
(if (= n 1) x (cons n x))
)
;;; Return descending list of prime numbers up to n without 1
;;; 6 => (5 3 2)
;;; no overflow checks needed.
;;; Using Eratosthenes's sieve
;;; (STD-PRIMES 2500): 2403 ms,(STD-PRIMES 10000): 19088 ms,
;;; (STD-PRIMES 1000000000007) : stopped after 23hrs.
(defun STD-PRIMES (n / i)
(if (zerop (rem n 2)) (setq n (1- n)))
(if (<= n (car *primes*))
(std-%primes-upto n)
(progn
(setq i (car *primes*))
(while (<= (setq i (+ i 2)) n)
(if (std-primep i)
(setq *primes* (cons i *primes*))))
*primes*)))
;;; Is the given number a prime number?
;;; Caches previous results to speed up consecutive calls on larger
;;; numbers. We don't apply probability tests for larger numbers yet.
;;; (STD-PRIMEP 10000005): 0 ms, (STD-PRIMEP 10000007): 361 ms
(defun STD-PRIMEP (n)
(if (<= n (car *primes*))
(not (null (member n *primes*)))
(std-%prime-sieve n)
)
)
;;; Test the sieve up to root
;;; no overflow checks needed. stores the failing divisor in *lastdiv*
;;; uses two global lists* *primes* and *sieves*, both contain decsending
;;; primes, but *sieves* is truncated up to the root of n to speed it up.
(defun std-%prime-sieve (n / i p root)
(setq root (fix (sqrt n)))
;; keep a list of divisors to be tested. this is the tail of
;; *primes* from root. keep it global for std-primes
(setq i (car *sieve*))
(setq root (max 2 (if (zerop (rem root 2)) (1- root) root)))
(if (or (not i)
(< i root))
;; *sieve* too small, start with all we get from *primes*
(setq *sieve* (std-%primes-upto root)))
(if (> i root) ; nil is like 0 with <
;; we already have enough divisors in our sieve, just test it.
;; *sieve* too large, use only the relevant tail
;; otherwise we would have to destroy our costly created *sieve*
(std-%primep n (std-%primes-upto root))
(progn
;; test the current sieve
(setq p (std-%primep n *sieve*))
(setq i (car *sieve*))
;;fill the sieve until a divisor is found or the seed exceeds root
(while (and p (<= (setq i (+ i 2)) root))
(if (std-%primep i *sieve*)
(progn
(setq *sieve* (cons i *sieve*))
;; not really needed, only to speed up the next calls:
(if (> i (car *primes*))
(setq *primes* (cons i *primes*)))
(if (zerop (rem n i))
(setq p nil)
)
)
)
)
p
)
)
)
;;; Is i not dividable by one of nums? (relative prime?)
;;; Try all divisors but break when one returns zero
;;; Sideeffect: store *lastdiv* for the factorization
;;; Note: of course we could also use (= 1 (std-gcd (cons i div)))
;;; but this is more efficient.
(defun std-%primep (i div)
;; could maybe improved to omit the div test by adding
;; a unique TRUE element to the end and test for it at the end.
(setq div (reverse div))
(while (and div (not (zerop (rem i (car div)))))
(setq div (cdr div)))
(setq *lastdiv* (car div))
(null div))
;;; return descending primes upto max. n from the cached list
;;; same as (std-member-if '(lambda (x) (<= x n)) *primes*))
;;; (17 13 11 7 5 3 2), 5 => (5 3 2)
(defun std-%primes-upto (n / lst)
(setq lst *primes*) ;return only the relevant tail of primes
(while (< n (car lst)) (setq lst (cdr lst)))
lst
)
(defun STD-PROBABLE-PRIME-P (n)
(std-error (list "STD-PROBABLE-PRIME-P"
(std-msg ;|MSG18|;" not yet implemented"))))
;;; -------------------------------------------------------------------73
;|#+ VL|;
(if std-%protect-assign
(std-%protect-assign
'(STD-PRIME-FACTORIZATION
STD-PRIMES
STD-PRIMEP
STD-PROBABLE-PRIME-P
)))
;|END #+ VL|;
;;; Module dependencies - None
(std-provide "PRIMES")