Hi,
I'm still writing on a version of primep that will prove primality for
all inputs of reasonable size (whatever that means).
Because of events not really under my control, I am progressing a lot
slower than I hoped, but I have a fairly functional probable prime
test and a rigorous test for numbers smaller than 341550071728321.
I wanted to improve on it before I publish it, but since people are
talking about primep a lot recently maybe it's better to release the
code as is.
Here is it:
;;; compute x^n mod m. No error checking on the arguments!
(defun %expt-mod (x n m)
(if (zerop n)
1
(do ((f (- (integer-length n) 2) (- f 1))
(y x))
((< f 0) y)
(setq y (mod (* y y) m))
(when (logbitp f n)
(setq y (mod (* y x) m))))))
;;; if p is provably prime return t else nil.
;;; for the sufficiency of miller rabin tests look at literature cited in
;;; The Prime Pages http://www.utm.edu/research/primes/
(defun primep (n)
"If n is prime return true, else nil. Proves primality! n must be a positive integer."
(when (<= n 1) (return-from primep nil)) ; 0 is not a prime and 1 is a unity,
; also not a prime
(dolist (p '(2 3 5 7 11)) ; weed out a few small primes
(when (zerop (mod n p)) ; actually almost 80% of all numbers
(return-from primep (eql n p))))
;; now let's do a few miller-rabin tests
;; these are sufficient to prove primality up to 341550071728321
;; and will weed out many other composites
(let* ((n-1 (1- n)) ; split up n-1 = 2^s * r with r odd
(s (loop for s upfrom 1 ;n-1 is always even here since n is odd
until (logbitp s n-1)
finally (return s)))
(r (ash n-1 (- s))))
(unless (and (%miller-rabin-test n 2 r s) (%miller-rabin-test n 3 r s))
(return-from primep nil))
(when (< n 1373653) ; miller-rabin tests to base 2 and 3 are sufficient
(return-from primep t)) ; to prove primality below 1373653
(unless (%miller-rabin-test n 5 r s)
(return-from primep nil))
(when (< n 25326001) ; miller-rabin tests to base 2, 3 and 5 are sufficient
(return-from primep t)) ; to prove primality below 25326001
(unless (and (%miller-rabin-test n 7 r s) (%miller-rabin-test n 11 r s))
(return-from primep nil))
(when (< n 2152302898747) ; miller-rabin tests to base 2, 3, 5, 7 an 11
(return-from primep t)) ; are sufficient for primality below 2152302898747
(unless (%miller-rabin-test n 13 r s)
(return-from primep nil))
(when (< n 3474739660383) ; miller-rabin tests to base 2,3,5,7,11 and 13 suffice
(return-from primep t)) ; to prove primality below 3474739660383
(unless (%miller-rabin-test n 17 r s)
(return-from primep nil))
(when (< n 341550071728321) ; miller-rabin tests to base 2,3,5,7,11,13 and 17
(return-from primep t)) ; suffice to prove primality below 341550071728321
;;do a few more miller-rabin tests with random bases
(loop with limit = (min most-positive-fixnum (- n 21))
repeat 10
for a fixnum = (+ (random limit) 20)
unless (%miller-rabin-test n a r s) do (return-from primep nil))
)
;; now we should turn to lucas pseudoprimality tests and then to the big guns
(%prove-prime n)
)
(defun %prove-prime (n)
"prove or disprove that n is prime"
'dont-know)
;;; Return nil if n is a composite number and t if it is a pseudo-prime to base a
(defun %miller-rabin-test (n a r s)
"do a miller-rabin test on n to base a. n-1 = 2^s * r with r odd."
(let ((y (%expt-mod a r n))
(n-1 (1- n)))
(if (and (/= y 1) (/= y n-1))
(loop for j from 1
while (and (<= j (1- s)) (/= y n-1))
do (setf y (mod (* y y) n))
never (= y 1)
finally (return (= y n-1)))
t)))
(defun gcd-extended (a b)
"given integers a and b compute d, u, v such that au+bv = d = (a,b)"
(when (not (and (integerp a) (integerp b)))
(error "GCD-EXTENDED: both arguments must be integers"))
(when (zerop b) (values a 1 0))
(do ((u 1) (d a) (v1 0) (v3 b) (t1 0))
((zerop v3) (if (plusp d)
(values d u (/ (- d (* a u)) b))
(values (- d) (- u) (/ (- (* a u) d) b))))
(declare (integer u d v1 v3 t1))
(multiple-value-bind (q t3) (floor d v3)
(setq t1 (- u (* q v1))
u v1
d v3
v1 t1
v3 t3))))
(defun mod-inv (x m)
"given integers x and m, compute the inverse of x mod m"
(unless (and (integerp x) (integerp m))
(error "MOD-INV: both arguments must be integers"))
(multiple-value-bind (d u)
(gcd-extended x m)
(if (eql d 1)
(mod u m)
(error "MOD-INV: arguments must be relatively prime"))))
(defun expt-mod (x n m)
"given integers x n m, compute x^n mod m"
(unless (and (integerp x) (integerp n) (integerp m))
(error "EXPTMOD: all arguments must be integers"))
(cond ((zerop n) 1)
((minusp n) (%expt-mod (mod-inv x m) (- n) m))
(t (%expt-mod (mod x m) n m))))
As it currently stands, the code returns 'dont-know if primality or
compositeness can't be proven with miller rabin.
I hope this helps, and if I ever get more time I promise to first
implement the lucas psudo-primality test and then an ecm test (or an
equivalent).
'Andreas
--
Wherever I lay my .emacs, there's my $HOME.