Does anyone know of a prime-number generator in CL that's available
under a liberal license? Failing that, GPL or an inexpensive
commercial license would be okay. Or do I need to go write my own?
--
/|_ .-----------------------.
,' .\ / | No to Imperialist war |
,--' _,' | Wage class war! |
/ / `-----------------------'
( -. |
| ) |
(`-. '--.)
`. )----'
Thomas F. Burdick wrote:
> Does anyone know of a prime-number generator in CL that's available
> under a liberal license? Failing that, GPL or an inexpensive
> commercial license would be okay. Or do I need to go write my own?
I s a Miller Rabin Pseudoprime generator enough?
ciao,
Jochen
--
http://www.dataheaven.de
Jochen Schmidt <···@dataheaven.de> writes:
> Thomas F. Burdick wrote:
>
> > Does anyone know of a prime-number generator in CL that's available
> > under a liberal license? Failing that, GPL or an inexpensive
> > commercial license would be okay. Or do I need to go write my own?
>
> I s a Miller Rabin Pseudoprime generator enough?
Oops, I should've mentioned that -- yes, this is for industrial use,
so I just need a controlled likelihood of the number being prime.
--
/|_ .-----------------------.
,' .\ / | No to Imperialist war |
,--' _,' | Wage class war! |
/ / `-----------------------'
( -. |
| ) |
(`-. '--.)
`. )----'
Search in sci.math.symbolic, I think, for
Copyright Gerald Roylance 1983, 1984, 1985, 1986
He posted code that I am running. It is fast, even for big numbers. His
post includes copyright permissions only for non-commercial use. You you
might be able to contact him for commercial permission.
The post is in a 46-article thread about a "Pseudoprimes too strong for
Maple". In fact, Maple's default number of iterations in the probability
test is (was?) 10. Upping the iterations fixed the false positives.
Roylance's test did not fail in any of the thread's cases for which
Maple failed. However, Roylance's iterations may have been higher. I use
30 and it's plenty fast enough for me.
Knuth Vol II 4.5.4 has a different algorithm but his comment on
probabilistic prime testing is interesting. His algorithm return a false
positive at most 1/4th of the time (and this is rare, but even in his
worst case...). Invoke it 25 times in a row and there is "less than a
chance in a quadrillion for one error in a billion primes."
Here is a sample of Roylance's code's speed. My function fastprimes
takes a starting number and an increment. Every number from start to (+
start inc) is tested. As I mentioned, I run each integer's test 30
times. I believe the starting number here is 123 digits.
CL-USER 314 > (time (fastprimes
109393923092309190109123902109129021901320981240981240912409124091241240909329032903290190120921901290129012901290210912091
1000))
Timing the evaluation of (FASTPRIMES
109393923092309190109123902109129021901320981240981240912409124091241240909329032903290190120921901290129012901290210912091
1000)
user time = 4.546
system time = 0.000
Elapsed time = 0:00:04
Allocation = 57075096 bytes standard / 14608 bytes fixlen
0 Page faults
(109393923092309190109123902109129021901320981240981240912409124091241240909329032903290190120921901290129012901290210912623
109393923092309190109123902109129021901320981240981240912409124091241240909329032903290190120921901290129012901290210912733)
Thomas F. Burdick wrote:
> Jochen Schmidt <···@dataheaven.de> writes:
>
>> Thomas F. Burdick wrote:
>>
>> > Does anyone know of a prime-number generator in CL that's available
>> > under a liberal license? Failing that, GPL or an inexpensive
>> > commercial license would be okay. Or do I need to go write my own?
>>
>> I s a Miller Rabin Pseudoprime generator enough?
>
> Oops, I should've mentioned that -- yes, this is for industrial use,
> so I just need a controlled likelihood of the number being prime.
You can try my implementation of the Miller Rabin algorithm
http://dataheaven.dnsalias.net/~neonsquare/crypto/math.tgz
ciao,
Jochen
--
http://www.dataheaven.de
Jochen Schmidt <···@dataheaven.de> writes:
> Thomas F. Burdick wrote:
>
> > Jochen Schmidt <···@dataheaven.de> writes:
> >
> >> Thomas F. Burdick wrote:
> >>
> >> > Does anyone know of a prime-number generator in CL that's available
> >> > under a liberal license? Failing that, GPL or an inexpensive
> >> > commercial license would be okay. Or do I need to go write my own?
> >>
> >> I s a Miller Rabin Pseudoprime generator enough?
> >
> > Oops, I should've mentioned that -- yes, this is for industrial use,
> > so I just need a controlled likelihood of the number being prime.
>
> You can try my implementation of the Miller Rabin algorithm
>
> http://dataheaven.dnsalias.net/~neonsquare/crypto/math.tgz
Wow, I didn't realize it would be so easy to implement. I'd read a
mathematical description of it years ago, and remember being confused
-- which was probably a result of the confusing math, not the
difficulty of implementation. But thank you, it's nice and speedy,
finding me all the primes I need and more in under 1 second on an old
SPARC :)
--
/|_ .-----------------------.
,' .\ / | No to Imperialist war |
,--' _,' | Wage class war! |
/ / `-----------------------'
( -. |
| ) |
(`-. '--.)
`. )----'
···@apocalypse.OCF.Berkeley.EDU (Thomas F. Burdick) writes:
> Does anyone know of a prime-number generator in CL that's available
> under a liberal license? Failing that, GPL or an inexpensive
> commercial license would be okay. Or do I need to go write my own?
The one in CMUCL is supposed to be good:
src/code/rand.lisp:
;;; **********************************************************************
;;;
;;; Functions for generating random numbers for CMU Common Lisp.
;;;
;;; Written by Raymond Toy. This is a replacement of the original
;;; rand.lisp in CMUCL. According to tests done by Peter VanEynde,
;;; the original generator fails some tests of randomness in
;;; Marsaglia's diehard test suite. This generator passes those
;;; tests.
;;;
;;; Also, the original generator lacks documentation. This remedies
;;; that, and includes references for the algorithms used. Hopefully
;;; better algorithms have been chosen that are at least as good and
;;; at least as fast as the original. In fact, the heart of the
;;; generator is geared towards generating floating point numbers
;;; which should be faster than the original. This should also give
;;; faster integer numbers since longer floating-point random numbers
;;; are generated. Tests have shown that this new generator is always
;;; a 10-20% faster and upto 3 times faster for double-floats.
or the much better:
src/code/rand-mt.lisp:
;;; **********************************************************************
;;;
;;; Support for the Mersenne Twister, MT19937, random number generator
;;; due to Matsumoto and Nishimura. This implementation has been
;;; placed in the public domain with permission from M. Matsumoto.
;;;
;;; Makoto Matsumoto and T. Nishimura, "Mersenne twister: A
;;; 623-dimensionally equidistributed uniform pseudorandom number
;;; generator.", ACM Transactions on Modeling and Computer Simulation,
;;; 1997, to appear.
Both are public domain.
Groetjes, Peter
--
It's logic Jim, but not as we know it. | ········@debian.org
"God, root, what is difference?" - Pitr| http://people.debian.org/~pvaneynd/
"God is more forgiving." - Dave Aronson| http://users.belgacom.net/pvaneynd/
Peter Van Eynde wrote:
> ···@apocalypse.OCF.Berkeley.EDU (Thomas F. Burdick) writes:
>
>> Does anyone know of a prime-number generator in CL that's available
>> under a liberal license? Failing that, GPL or an inexpensive
>> commercial license would be okay. Or do I need to go write my own?
>
> The one in CMUCL is supposed to be good:
[snipped references to CMUCLs rng]
There is a difference between random numbers and prime numbers ;-)
ciao,
Jochen
--
http://www.dataheaven.de
check clocc http://sourceforge.net/projects/clocc/
In the file .../clocc/src/cllib/math.lisp are a set of prime
number routines.
--
**************************************************
Dr Michael A. Koerber Micro$oft Free Zone
MIT/Lincoln Laboratory
Here's an implementation of the sieve that I did (it is distributed with the
Corman Lisp examples). I think it is pretty fast.
(defun sieve5 (n)
"Returns a list of all primes from 2 to n"
(declare (fixnum n) (optimize (speed 3) (safety 0)))
(let* ((a (make-array n :element-type 'bit :initial-element 0))
(result (list 2))
(root (isqrt n)))
(declare (fixnum root))
(do ((i 3 (the fixnum (+ i 2))))
((>= i n) (nreverse result))
(declare (fixnum i))
(progn (when (= (sbit a i) 0)
(push i result)
(if (< i root)
(do* ((inc (+ i i))
(j (* i i) (the fixnum (+ j inc))))
((>= j n))
(declare (fixnum j inc))
(setf (sbit a j) 1))))))))
Roger
----------------------------
On 27 Nov 2002 23:47:38 -0800, ···@apocalypse.OCF.Berkeley.EDU (Thomas F.
Burdick) wrote:
>Does anyone know of a prime-number generator in CL that's available
>under a liberal license? Failing that, GPL or an inexpensive
>commercial license would be okay. Or do I need to go write my own?
>
>--
> /|_ .-----------------------.
> ,' .\ / | No to Imperialist war |
> ,--' _,' | Wage class war! |
> / / `-----------------------'
> ( -. |
> | ) |
> (`-. '--.)
> `. )----'
·····@corman.net (Roger Corman) writes:
> Here's an implementation of the sieve that I did (it is distributed with the
> Corman Lisp examples). I think it is pretty fast.
>
> (defun sieve5 (n)
> "Returns a list of all primes from 2 to n"
> (declare (fixnum n) (optimize (speed 3) (safety 0)))
> (let* ((a (make-array n :element-type 'bit :initial-element 0))
> (result (list 2))
> (root (isqrt n)))
> (declare (fixnum root))
> (do ((i 3 (the fixnum (+ i 2))))
> ((>= i n) (nreverse result))
> (declare (fixnum i))
> (progn (when (= (sbit a i) 0)
> (push i result)
> (if (< i root)
This should be <=, not <, I think.
> (do* ((inc (+ i i))
> (j (* i i) (the fixnum (+ j inc))))
> ((>= j n))
> (declare (fixnum j inc))
> (setf (sbit a j) 1))))))))
Andras
>
> Roger
On 01 Dec 2002 03:44:47 +0100, ······@math.bme.hu (Simon Andr�s) wrote:
>·····@corman.net (Roger Corman) writes:
>> (defun sieve5 (n)
>> "Returns a list of all primes from 2 to n"
>> (declare (fixnum n) (optimize (speed 3) (safety 0)))
>> (let* ((a (make-array n :element-type 'bit :initial-element 0))
>> (result (list 2))
>> (root (isqrt n)))
>> (declare (fixnum root))
>> (do ((i 3 (the fixnum (+ i 2))))
>> ((>= i n) (nreverse result))
>> (declare (fixnum i))
>> (progn (when (= (sbit a i) 0)
>> (push i result)
>> (if (< i root)
>
>This should be <=, not <, I think.
>
>> (do* ((inc (+ i i))
>> (j (* i i) (the fixnum (+ j inc))))
>> ((>= j n))
>> (declare (fixnum j inc))
>> (setf (sbit a j) 1))))))))
No, marking multiples of the root is redundant. They will already have been
marked by previous loops of lower prime factors. All non-primes above the root
will include at least one prime factor less than the root.
For example:
(sieve1 25)
Example: You don't need to loop on 5, to mark off 10, 15 and 20, because they
are multiples of 2 or 3 respectively and have already been marked.
Note this algorithm only finds primes from 2 to n (exclusive) so if n is prime
it will not return it. I guess the doc string should be "Returns a list from
primes from 2 to (n-1)."
Roger
·····@corman.net (Roger Corman) writes:
> No, marking multiples of the root is redundant. They will already have been
> marked by previous loops of lower prime factors. All non-primes above the root
> will include at least one prime factor less than the root.
This would be true if root were (sqrt n), not (isqrt n). With isqrt,
however, if n is 26, then root is 5, and 25 is not marked.
CL-USER(19): (sieve5 26)
(2 3 5 7 11 13 17 19 23 25)
CL-USER(20):
More precisely:
Though the following is true
(forall q<n, q is odd and composite)(exists i < (sqrt n), i
prime)(exists k natural number) q = i^2 + 2ki
[take i to be the smallest prime factor of q], the same with isqrt in
place of sqrt is not. n=26, q=25 is a counterexample.
Andras
On 02 Dec 2002 13:24:43 +0100, ······@math.bme.hu (Simon Andr�s)
wrote:
>·····@corman.net (Roger Corman) writes:
>
>
>> No, marking multiples of the root is redundant. They will already have been
>> marked by previous loops of lower prime factors. All non-primes above the root
>> will include at least one prime factor less than the root.
>
>This would be true if root were (sqrt n), not (isqrt n). With isqrt,
>however, if n is 26, then root is 5, and 25 is not marked.
>
Yes, you are right of course. I thought of that last night after
posting it, realizing the bug. Your suggestion to simply change the
comparison to <= would fix this. It may be more efficient to just mark
(* root root) rather than do the whole loop.
Roger