101 lines
2.4 KiB
Text
101 lines
2.4 KiB
Text
' version 29-11-2016
|
|
' compile with: fbc -s console
|
|
|
|
' TRUE/FALSE are built-in constants since FreeBASIC 1.04
|
|
' But we have to define them for older versions.
|
|
#Ifndef TRUE
|
|
#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 pow_mod(b As ULongInt, power As ULongInt, modulus As ULongInt) As ULongInt
|
|
' returns b ^ power mod modulus
|
|
Dim As ULongInt x = 1
|
|
|
|
While power > 0
|
|
If (power And 1) = 1 Then
|
|
' x = (x * b) Mod modulus
|
|
x = mul_mod(x, b, modulus)
|
|
End If
|
|
' b = (b * b) Mod modulus
|
|
b = mul_mod(b, b, modulus)
|
|
power = power Shr 1
|
|
Wend
|
|
|
|
Return x
|
|
|
|
End Function
|
|
|
|
Function miller_rabin_test(n As ULongInt, k As Integer) As Byte
|
|
|
|
If n > 9223372036854775808ull Then ' limit 2^63, pow_mod/mul_mod can't handle bigger numbers
|
|
Print "number is to big, program will end"
|
|
Sleep
|
|
End
|
|
End If
|
|
|
|
' 2 is a prime, if n is smaller then 2 or n is even then n = composite
|
|
If n = 2 Then Return TRUE
|
|
If (n < 2) OrElse ((n And 1) = 0) Then Return FALSE
|
|
|
|
Dim As ULongInt a, x, n_one = n - 1, d = n_one
|
|
Dim As UInteger s
|
|
|
|
While (d And 1) = 0
|
|
d = d Shr 1
|
|
s = s + 1
|
|
Wend
|
|
|
|
While k > 0
|
|
k = k - 1
|
|
a = Int(Rnd * (n -2)) +2 ' 2 <= a < n
|
|
x = pow_mod(a, d, n)
|
|
If (x = 1) Or (x = n_one) Then Continue While
|
|
For r As Integer = 1 To s -1
|
|
x = pow_mod(x, 2, n)
|
|
If x = 1 Then Return FALSE
|
|
If x = n_one Then Continue While
|
|
Next
|
|
If x <> n_one Then Return FALSE
|
|
Wend
|
|
Return TRUE
|
|
|
|
End Function
|
|
' ------=< MAIN >=------
|
|
|
|
Randomize Timer
|
|
|
|
Dim As Integer total
|
|
Dim As ULongInt y, limit = 2^63-1
|
|
|
|
For y = limit - 1000 To limit
|
|
If miller_rabin_test(y, 5) = TRUE Then
|
|
total = total + 1
|
|
Print y,
|
|
End If
|
|
Next
|
|
|
|
Print : Print
|
|
Print total; " primes between "; limit - 1000; " and "; y -1
|
|
|
|
' empty keyboard buffer
|
|
While Inkey <> "" : Wend
|
|
Print : Print "hit any key to end program"
|
|
Sleep
|
|
End
|