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




' version 18-09-2015
' compile with: fbc -s console
 
#Ifndef TRUE        ' define true and false for older freebasic versions
    #Define FALSE 0
    #Define TRUE Not FALSE
#EndIf
 
Function mul_mod(a As ULongInt, b As ULongInt, modulus As ULongInt) As ULongInt
    ' returns a * b mod modulus
 
    Dim As ULongInt x , y = a ' a mod modulus, but a is already smaller then modulus
 
    While b > 0
        If (b And 1) = 1 Then
            x = (x + y) Mod modulus
        End If
        y = (y Shl 1) Mod modulus
        b = b Shr 1
    Wend
    Return x
 
End Function
 
Function LLT(p As UInteger) As Integer
 
    Dim As ULongInt s = 4, m = 1
    m = m Shl p : m = m - 1       ' m = 2 ^ p - 1
 
    For i As Integer = 2 To p - 1
        s = mul_mod(s, s, m) - 2
    Next
 
    If s = 0 Then Return TRUE Else Return FALSE
 
End Function
 
' ------=< MAIN >=------
 
Dim As UInteger p
 
Print
' M2 can not be tested, we start with 3
for p = 3 To 63
    If LLT(p) = TRUE Then Print " M";Str(p);
Next
 
Print
' empty keyboard buffer
While Inkey <> "" : Wend
Print : Print "hit any key to end program"
Sleep
End