Miller–Rabin primality test: Difference between revisions
added FreeBASIC
(→{{header|Fortran}}: Split.) |
(added FreeBASIC) |
||
Line 1,526:
Thus, there is no avoiding a special MODEXP function, even for small test numbers.
=={{header|FreeBASIC}}==
Using the task pseudo code
===Up to 2^63-1===
<lang freebasic>' 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</lang>
{{out}}
<pre>9223372036854774893 9223372036854774917 9223372036854774937
9223372036854774959 9223372036854775057 9223372036854775073
9223372036854775097 9223372036854775139 9223372036854775159
9223372036854775181 9223372036854775259 9223372036854775279
9223372036854775291 9223372036854775337 9223372036854775351
9223372036854775399 9223372036854775417 9223372036854775421
9223372036854775433 9223372036854775507 9223372036854775549
9223372036854775643 9223372036854775783
23 primes between 9223372036854774808 and 9223372036854775808</pre>
===Using Big Integer library===
{{libheader|GMP}}
<lang freebasic>' 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
#Include Once "gmp.bi"
#Macro big_int(a)
Dim As Mpz_ptr a = Allocate( Len( __mpz_struct))
Mpz_init(a)
#EndMacro
Dim Shared As __gmp_randstate_struct rnd_seed
Function miller_rabin(big_n As Mpz_ptr, num_of_tests As ULong) As Byte
#Macro cleanup
mpz_clear(n_1) : mpz_clear(a) : mpz_clear(d)
mpz_clear(n_2) : mpz_clear(x)
#EndMacro
If mpz_cmp_ui(big_n, 1) <= 0 Then
Print "Number not allowed"
Sleep 5000
End If
If mpz_cmp_ui(big_n, 2) = 0 OrElse mpz_cmp_ui(big_n, 3) = 0 Then
Return TRUE ' 2 = prime , 3 = prime
End If
If mpz_tstbit(big_n, 0) = 0 Then Return FALSE ' even number, no prime
Dim As ULong r, s
big_int(n_1) : big_int(n_2) : big_int(a) : big_int(d) : big_int(x)
mpz_sub_ui(n_1, big_n, 1) : mpz_sub_ui(n_2, big_n, 2) : mpz_set(d, n_1)
While mpz_tstbit(d, 0) = 0
mpz_fdiv_q_2exp(d, d, 1)
s = s +1
Wend
While num_of_tests > 0
num_of_tests -= 1
mpz_urandomm(a, @rnd_seed, n_2)
mpz_add_ui(a, a, 2)
mpz_powm(x, a, d, big_n)
If mpz_cmp_ui(x, 1) = 0 Or mpz_cmp(x, n_1) = 0 Then Continue While
For r = 1 To s -1
mpz_powm_ui(x, x, 2, big_n)
If mpz_cmp_ui(x, 1) = 0 Then
cleanup
Return FALSE
End If
If mpz_cmp(x, n_1) = 0 Then Continue While
Next
If mpz_cmp(x, n_1) <> 0 Then
cleanup
Return FALSE
End If
Wend
cleanup
Return TRUE
End Function
' ------=< MAIN >=------
Dim As ZString Ptr gmp_str : gmp_str = Allocate(10000000)
big_int(big_n)
Randomize Timer
gmp_randinit_mt(@rnd_seed)
gmp_randseed_ui(@rnd_seed, Rnd * &HFFFFFFFF) ' seed the random number generator
Dim As Long x
For x = 2 To 100
mpz_set_ui(big_n, x)
If miller_rabin(big_n, 5) = TRUE Then
Print Using "####"; x;
End If
Next
Print : Print
For x = 2 To 3300
mpz_set_ui(big_n, 1)
mpz_mul_2exp(big_n, big_n, x)
mpz_sub_ui(big_n,big_n,1)
If miller_rabin(big_n, 5) = TRUE Then
gmp_str = Mpz_get_str(0, 10, big_n)
Print "2^";Str(x);"-1 = prime"
End If
Next
Print
gmp_randclear(@rnd_seed)
mpz_clear(big_n)
DeAllocate(gmp_str)
' empty keyboard buffer
While InKey <> "" : Wend
Print : Print "hit any key to end program"
Sleep
End</lang>
{{out}}
<pre> 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
2^2-1 = prime
2^3-1 = prime
2^5-1 = prime
2^7-1 = prime
2^13-1 = prime
2^17-1 = prime
2^19-1 = prime
2^31-1 = prime
2^61-1 = prime
2^89-1 = prime
2^107-1 = prime
2^127-1 = prime
2^521-1 = prime
2^607-1 = prime
2^1279-1 = prime
2^2203-1 = prime
2^2281-1 = prime
2^3217-1 = prime</pre>
=={{header|FunL}}==
| |||