From: William Echard
Subject: Lisp code for calculating primes?
Date: 
Message-ID: <D2vt8C.98F@newshub.ccs.yorku.ca>
From: Micheal S. Hewett
Subject: Re: Lisp code for calculating primes?
Date: 
Message-ID: <3g1jq0$367@cascade.cs.utexas.edu>
Your message was blank, but I presume that you would like some
LISP code to find prime numbers.  Below is some code that
does just that.  It can find a 200-digit prime number in less
than 3 minutes on an RS/6000 in Allegro CL.  I wrote it while
playing around with the RSA encryption algorithm in class
last semester.  I'll place it in the public domain.

It can probably be optimized even further by adding more type declarations.
Note: finding a 200-digit prime number conses 95MB of garbage.

Call:

       (find-prime 10)

To find a 10-digit prime number.  You can easily modify the find-prime
function to test the primality of a given number.

Mike
(······@cs.utexas.edu)



;;; ********************************************************************
;;;
;;; prime.lisp   -  find prime numbers
;;;
;;; Micheal Hewett   11 October 1994
;;;                  University of Texas
;;;                  ······@cs
;;;
;;; Some info and algorithms are from "Elementary Number Theory and its
;;; Applications", Kenneth H. Rosen, Addison-Wesley, 1985.
;;;
;;; -------------------------------------------------
;;; Primality testing Functions
;;; -------------------------------------------------

(in-package 'user)

(defparameter *max-tests* 50
  "How many times do you want to call Miller's test for primality?
A value of 100 gives virtual certainty.  p is (1/4)^n")


(defun factor-out-twos (possible-prime)
  "Returns two values: s and o, where ((o * 2^s) + 1)= possible-prime."

  ;;; Find odd factor 'o', such that (2^s)*o = possible-prime.
  ;;; This is 't' in Rosen's book (p. 156).

  (let* ((o          (1- possible-prime))
	 (s          0)
	 (next-o     0)
	 (remainder  0)
	)
    (loop
     (multiple-value-setq
      (next-o remainder)
      (floor o 2))

     (when (not (zerop remainder))
       (return))

     (incf s)
     (setq o next-o)
     )

    (values s o)     ;; Return values
    )
  )  

(defun square (x) (* x x))


(defun expmod (base exponent modulus)
  "Returns (b^e) mod m"

  (declare (type integer base exponent modulus))

  ;;This is defined in Abelson & Sussman (1985).

  (cond ((zerop exponent) 1)
     
	((evenp exponent)
	 (rem (square (expmod base (/ exponent 2) modulus)) modulus))

	(T  (rem (* base (expmod base (- exponent 1) modulus)) modulus))
	)
  )


(defun millers-test (possible-prime base s o)
  "Returns T if the prime passes Miller's test for the base.
S is the largest number such that ((2^s)*o + 1) = possible-prime."

  (declare (type integer possible-prime base s o))

  (let* ((pp-1       (1- possible-prime))
	 (mod-value  (expmod base o possible-prime))
	 )
    (declare (type integer pp-1 mod-value))

    (or (= 1 mod-value)  ;; Test: b^o = 1 mod p

	(dotimes (j s NIL)                    ;;b^(o*2^j) = -1 mod p
	  (when (= pp-1 mod-value)
	    (return T))
	  (setq mod-value (mod (* mod-value mod-value)
			       possible-prime)) ; Test: b^(no) => b^(2no)
	  )
	)
    )
  )
	

(defun find-prime (number-of-digits)
  "Returns a prime number greater than 10^(n-1).
This functions uses Rabin's probabilistic primality test."

  (do* ((p        (1+ (expt 10 (1- number-of-digits)))   (+ p 2))
	(base     0)
	(s        0)
	(o        0)
	)
    (nil nil)

    (multiple-value-setq
     (s o)
     (factor-out-twos p)
     )

    ;; Check *max-tests* random bases smaller than the prime.
    ;; If it passes Miller's test for 100 of them, it is almost
    ;; certainly a prime (p < 10^-60).
    ;; 
    (when
      (dotimes (i *max-tests* T)
	(setq base (random p))
	(unless (millers-test p base s o)
	  (return NIL))
	)
      (return p)
      )
    )
  )
From: Scott D. Anderson
Subject: Re: Lisp code for calculating primes?
Date: 
Message-ID: <ANDERSON.95Jan25173600@earhart.cs.umass.edu>
In article <··········@newshub.ccs.yorku.ca> ·······@trocaz.yorku.ca
(William Echard) seems to want Lisp code to calculate primes.

The following code may do what you want. 

;;;; -*- Mode:Common-Lisp; Package:user; Base:10 -*-

;;;; **************************************************************************
;;;; **************************************************************************
;;;; *                                                                        *
;;;; *                     Simple Prime Number Functions                      *
;;;; *                                                                        *
;;;; **************************************************************************
;;;; **************************************************************************
;;;
;;; Written by: Scott D. Anderson
;;;             Experimental Knowledge Systems Laboratory
;;;             Paul R. Cohen, Principal Investigator
;;;             David L. Westbrook, Systems Manager
;;;             David M. Hart, Laboratory Manager
;;;             Department of Computer Science
;;;             University of Massachusetts
;;;             Amherst, Massachusetts 01003.
;;;
;;; * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
;;;
;;;  03-09-94 File Created.  (anderson)
;;;
;;; * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
;;; --*--

(in-package #+CLTL2 USER #-CLTL2 'USER)

;;; --*--
;;; ***************************************************************************

(defvar *sieve-of-eratosthenes* nil
  "A bitvector for representing whether a number is prime.  1 means composite, 0 means prime.")

(defun sieve-of-eratosthenes (max)
  "Fills *sieve-of-eratosthenes* up to `max.' Then a simple aref tells whether a number is prime."
  (declare (optimize speed (safety 0)))
  (check-type max fixnum)
  (cond ((null *sieve-of-eratosthenes*)
	 (setf *sieve-of-eratosthenes*
	       (make-array (1+ max) :element-type 'bit :initial-element 0)))
	((<= (length *sieve-of-eratosthenes*) max)
	 (let ((old *sieve-of-eratosthenes*))
	   (setf *sieve-of-eratosthenes*
		 (make-array (1+ max) :element-type 'bit :initial-element 0))
	   (replace *sieve-of-eratosthenes* old)))
	(t
	 (return-from sieve-of-eratosthenes)))
  (let ((sieve *sieve-of-eratosthenes*))
    (let ((prime 2)
	  (max-p (isqrt max))
	  (ops   0))
      (declare (type fixnum prime max-p ops))
      (do () (nil)
	(do* ((p (lsh prime 1) (+ prime p)))
	     ((< max p))
	  (declare (type fixnum p))
	  (incf ops)
	  (setf (aref sieve p) 1))
	(do () (nil)
	  (incf prime)
	  (when (< max-p prime)
	    (return-from sieve-of-eratosthenes ops))
	  (when (zerop (aref sieve prime))
	    (return)))))))

(defun primep (x)
  (check-type x fixnum)
  (sieve-of-eratosthenes x)
  (zerop (aref *sieve-of-eratosthenes* x)))

(defun print-primes (high)
  (check-type high fixnum)
  (sieve-of-eratosthenes high)
  (loop for i from 2 to high
	when (zerop (aref *sieve-of-eratosthenes* i)) do
	(format t "~d~%" i))
  (values))

(defun test-prime-functions ()
  (setf *sieve-of-eratosthenes* nil)
  (sieve-of-eratosthenes 50)
  (setf *sieve-of-eratosthenes* nil)
  (print-primes 50)
  (format t "~%")
  (print-primes 100)
  (format t "~2%primep(~d) = ~s~%" 101 (primep 101)))

;;; ***************************************************************************
;;; EOF
comp.lang.lisp / 1995 / jan