164 lines
3.7 KiB
Plaintext
164 lines
3.7 KiB
Plaintext
' version 11-04-2017
|
|
' compile with: fbc -s console
|
|
' maximum for p is 17 digits to be on the save side
|
|
|
|
' 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 Mod 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 Isprime(n As ULongInt, k As Long) As Long
|
|
' miller-rabin prime test
|
|
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
|
|
|
|
Function legendre_symbol (a As LongInt, p As LongInt) As LongInt
|
|
|
|
Dim As LongInt x = pow_mod(a, ((p -1) \ 2), p)
|
|
If p -1 = x Then
|
|
Return x - p
|
|
Else
|
|
Return x
|
|
End If
|
|
|
|
End Function
|
|
|
|
' ------=< MAIN >=------
|
|
|
|
Dim As LongInt b, c, i, k, m, n, p, q, r, s, t, z
|
|
|
|
For k = 1 To 7
|
|
Read n, p
|
|
Print "Find solution for n ="; n; " and p =";p
|
|
|
|
If legendre_symbol(n, p) <> 1 Then
|
|
Print n;" is not a quadratic residue"
|
|
Print
|
|
Continue For
|
|
End If
|
|
|
|
If p = 2 OrElse Isprime(p, 15) = FALSE Then
|
|
Print p;" is not a odd prime"
|
|
Print
|
|
Continue For
|
|
End If
|
|
|
|
s = 0 : q = p -1
|
|
Do
|
|
s += 1
|
|
q \= 2
|
|
Loop Until (q And 1) = 1
|
|
|
|
If s = 1 And (p Mod 4) = 3 Then
|
|
r = pow_mod(n, ((p +1) \ 4), p)
|
|
Print "Solution found:"; r; " and"; p - r
|
|
Print
|
|
Continue For
|
|
End If
|
|
|
|
z = 1
|
|
Do
|
|
z += 1
|
|
Loop Until legendre_symbol(z, p) = -1
|
|
c = pow_mod(z, q, p)
|
|
r = pow_mod(n, (q +1) \ 2, p)
|
|
t = pow_mod(n, q, p)
|
|
m = s
|
|
|
|
Do
|
|
i = 0
|
|
If (t Mod p) = 1 Then
|
|
Print "Solution found:"; r; " and"; p - r
|
|
Print
|
|
Continue For
|
|
End If
|
|
|
|
Do
|
|
i += 1
|
|
If i >= m Then Continue For
|
|
Loop Until pow_mod(t, 2 ^ i, p) = 1
|
|
b = pow_mod(c, (2 ^ (m - i -1)), p)
|
|
r = mul_mod(r, b, p)
|
|
c = mul_mod(b, b, p)
|
|
t = mul_mod(t, c, p)' t = t * b ^ 2
|
|
m = i
|
|
Loop
|
|
|
|
Next
|
|
|
|
Data 10, 13, 56, 101, 1030, 10009, 1032, 10009, 44402, 100049
|
|
Data 665820697, 1000000009, 881398088036, 1000000000039
|
|
|
|
' empty keyboard buffer
|
|
While InKey <> "" : Wend
|
|
Print : Print "hit any key to end program"
|
|
Sleep
|
|
End
|