Tonelli-Shanks algorithm: Difference between revisions
added FreeBASIC
Thundergnat (talk | contribs) m (→{{header|Perl 6}}: Minor code simplifications, add missing failing test) |
(added FreeBASIC) |
||
Line 151:
❌ error: not a square (mod p) (1032 10009)
</pre>
=={{header|FreeBASIC}}==
<lang FreeBASIC>' 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</lang>
{{out}}
<pre>Find solution for n = 10 and p = 13
Solution found: 7 and 6
Find solution for n = 56 and p = 101
Solution found: 37 and 64
Find solution for n = 1030 and p = 10009
Solution found: 1632 and 8377
Find solution for n = 1032 and p = 10009
1032 is not a quadratic residue
Find solution for n = 44402 and p = 100049
Solution found: 30468 and 69581
Find solution for n = 665820697 and p = 1000000009
Solution found: 378633312 and 621366697
Find solution for n = 881398088036 and p = 1000000000039
Solution found: 791399408049 and 208600591990</pre>
=={{header|Go}}==
| |||