Author: Tord Romstad
Date: 03:28:46 05/28/05
Go up one level in this thread
On May 28, 2005 at 05:24:51, Joseph Tadeusz wrote:
>On May 27, 2005 at 11:12:52, Vincent Diepeveen wrote:
>
>>A retry :
>>
>>// Diepeveen test to find mersenne primes
>>// test whether p is prime p = 2^n - 1
>>// when given n is prime
>>
>>t = 3;
>>for( j = 1 ; j < n ; j++ )
>> t = t*t; // mod p
>> if( t == p-3 ) printf("prime\n"); else printf("composite\n");
>>
>>3
>>9
>>19
>>20
>>28 <==> p-3 = 31 - 3 = 28
>>
>>Vincent
>
>Does this also work for primes larger than 31?
Yes, it works until the numbers get so big that you get an
integer overflow (which happens very quickly if you use
C with 32-bit or 64-bit ints). The algorithm in itself is
correct.
If you want to test with bigger numbers, try to find some
bignum artithmetics library for C (the GMP library is supposed
to be good, I think) or use a programming language with
built-in support for arbitrary-precision integer arithmetic,
like Common Lisp. Below you will find Common Lisp code
for the classical Lucas-Lehmer test as well as Vincent's variant:
(defun prime-p (n)
"Tests whether N is a prime number."
(when (>= n 2)
(loop for i from 2 to (floor (sqrt n))
when (= 0 (mod n i)) return nil
finally (return t))))
(defun mersenne-prime-p (p)
"Tests whether 2^P - 1 is a prime, assuming that P is prime."
(loop with n = (- (expt 2 p) 1)
for i from 2 to p
for s = 4 then (mod (- (* s s) 2) n)
finally (return (= 0 s))))
(defun vincent-mersenne-prime-p (p)
"Tests whether 2^P - 1 is a prime, assuming that P is prime."
(loop with n = (- (expt 2 p) 1)
for j from 1 to p
for s = 3 then (mod (* s s) n)
finally (return (= (- n 3) s))))
Finding all Mersenne primes below 2^1000 - 1 can now be done like this:
CL-USER> (loop for x from 3 to 1000
when (and (prime-p x) (mersenne-prime-p x))
collect x)
(3 5 7 13 17 19 31 61 89 107 127 521 607)
Let's have a look at the biggest number, 2^607 - 1:
CL-USER> (- (expt 2 607) 1)
531137992816767098689588206552468627329593117727031923199444138200403559860852242739162502265229285668889329486246501015346579337652707239409519978766587351943831270835393219031728127
Tord
This page took 0 seconds to execute
Last modified: Thu, 15 Apr 21 08:11:13 -0700
Current Computer Chess Club Forums at Talkchess. This site by Sean Mintz.