Lucas-Lehmer Test for Mersenne Primes

[Zurtex]

I'm trying to come up with a program that works out Mersenne prime numbers, I've successfully managed to get one to work that brute force attacks the problem and gets me up to 2^2281 - 1 as about the highest prime I can get with the computer still running fine.

I've been reading up and it seems one of the best ways is to use the Lucas-Lehmer Test (http://mathworld.wolfram.com/Lucas-LehmerTest.html) however I am fairly new to MATLAB and don't quite seem to get it to work. Could anyone give me some points where I am going wrong please.

(Below is code I have used to test if 2^p - 1 is prime or not)
(z records the biggest value of p such that 2^p - 1 is prime)

p=sym(2)

for p=2:100
    if issymprime(p)==sym('true')
        s = 4
        for s=4:(p-2)
            s = mod(s^2 - 2,2^p -1)
        end
        if s==mod(0,2^p-1)
            z=p
        end
    end
end

 

[Zurtex]

O.K just noticed a huge mistake in my programming, but can I leave the thread here for a sec in case it is still horribly wrong?

 

[Zurtex]

Really sorry for the triple post but I've got it to work with a slight problem:

p=sym(2)

for p=2:40
    if issymprime(p)==sym('true')
        s = 1
        r = 4
        for s=1:(p-2)
            r = mod(r^2 - 2,2^p -1)
        end
        if r==mod(0,2^p-1)
            z=p
        end
    end
end

That works, up to about 2^19 - 1, however after that it doesn't seem to pick anymore Mersenne prime numbers up. Could someone tell me what to do please.

 

[shmoe]

Hi, I think your problem is with how MATLAB stores data, it uses double precision variables, which are 64-bits and accurate to 15 decimal places or so. When computing your r's for p=31, at one stage you calculate 1416317953^2-2, which has 18 digits, before reducing mod 2^31-1. This is where your program bites it.

You can probably have it compute 1416317953^2-2 in pieces, reducing mod 2^31-1 each time. There is also an integer type, but I believe it only goes to 32 bits, not enough for you.

I don't have much experience with large integer operations in matlab, so I don't know if there's a clever way around this. On the other hand, Maple is quite capable of handling things like this.

 

[Zurtex]

Thanks, I thought that was my problem and managed to come up with a maple solution. My lecturer has already made some functions for us that can handle symbolic input including a round() function which seems to work great and does so by using a couple of maple functions

Anyway it isn't quite as fast as I hoped but I can probably change that over time when I get to know MATLAB a bit better.

 

[T.Rex]

Look at GIMPS

Hi,
You should have a look at the GIMPS project: http://www.mersenne.org/ .
You'll find useful information and a very fast program on ia32/Windows machines.
You could also use the GLucas program for Unix machines. It is multi-threaded.
About the theory, look at Ribenboim's book: "The little book of bigger primes" or at Williams' book: "Edouard Lucas and Primality Testing".
Does Mathlab use FFT for multiplying big numbers ?
Regards,
Tony

 

[Zurtex]

Thanks, with time I could work a much faster algorithm but I had to submit what I had today.

Unfortunatly I have little understanding of how MatLab works at the moment, but I will look into it.