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).




 
GET "libhdr"
 
LET M(n) = (1 << n) - 1
 
LET isMersennePrime(p) =
    p < 3 -> p = 2,
    VALOF {
        LET n = M(p)
        LET s = 4
        FOR i = 1 TO p-2 DO {
            muldiv(s, s, n) // ignore quotient; remainder is in result2
            s := result2 - 2
            s := s + (n & s < 0)
        }
        RESULTIS s = 0
    }
 
LET start() = VALOF {
    LET primes = #x28208A20A08A28AC // bitmask of primes upto 63
 
    writes("These Mersenne numbers are prime: ")
    FOR k = 0 TO 63 DO
        IF (primes & 1 << k) ~= 0 & isMersennePrime(k) THEN
            writef("M%d  ", k)
 
    wrch('*n')
    RESULTIS 0
}
 