Lucas-Lehmer Test:

for 



p


{\displaystyle p}

 an odd prime, the Mersenne number 




2

p


−
1


{\displaystyle 2^{p}-1}

 is prime if and only if 




2

p


−
1


{\displaystyle 2^{p}-1}

 divides 



S
(
p
−
1
)


{\displaystyle S(p-1)}

 where 



S
(
n
+
1
)
=
(
S
(
n
)

)

2


−
2


{\displaystyle S(n+1)=(S(n))^{2}-2}

, and 



S
(
1
)
=
4


{\displaystyle S(1)=4}

.



Task

Calculate all Mersenne primes up to the implementation's
maximum precision, or the 47th Mersenne prime   (whichever comes first).




function [mNumber,mersennesPrime] = mersennePrimes()
 
    function isPrime = lucasLehmerTest(thePrime)
 
        llResidue = 4;
        mersennesPrime = (2^thePrime)-1;
 
        for i = ( 1:thePrime-2 )
            llResidue = mod( ((llResidue^2) - 2),mersennesPrime );        
        end         
 
        isPrime = (llResidue == 0);
 
    end
 
    %Because IEEE764 Double is the highest precision number we can
    %represent in MATLAB, the highest Mersenne Number we can test is 2^52.
    %In addition, because we have this cap, we can only test up to the
    %number 30 for Mersenne Primeness. When we input 31 into the
    %Lucas-Lehmer test, during the computation of the residue, the
    %algorithm multiplies two numbers together the result of which is
    %greater than 2^53. Because we require every digit to be significant,
    %this leads to an error. The Lucas-Lehmer test should say that M31 is a
    %Mersenne Prime, but because of the rounding error in calculating the
    %residues caused by floating-point arithmetic, it does not. So M30 is
    %the largest number we test.
 
    mNumber = (3:30);
 
    [isPrime] = arrayfun(@lucasLehmerTest,mNumber);
 
    mNumber = [2 mNumber(isPrime)];
    mersennesPrime = (2.^mNumber) - 1;
 
end