

 This page uses content from Wikipedia. The original article was at Tonelli-Shanks algorithm. The list of authors can be seen in the page history. As with Rosetta Code, the text of Wikipedia is available under the GNU FDL. (See links for details on variance)



In computational number theory, the Tonelli–Shanks algorithm is a technique for solving for x in a congruence of the form:


 x2 ≡ n (mod p)

where n is an integer which is a quadratic residue (mod p), p is an odd prime, and x,n ∈ Fp where Fp = {0, 1, ..., p - 1}.

It is used in cryptography techniques.



To apply the algorithm, we need the Legendre symbol:

The Legendre symbol (a | p) denotes the value of a(p-1)/2 (mod p).

 (a | p) ≡ 1    if a is a square (mod p)
 (a | p) ≡ -1    if a is not a square (mod p)
 (a | p) ≡ 0    if a ≡ 0 (mod p) 


Algorithm pseudo-code


All ≡ are taken to mean (mod p) unless stated otherwise.

 Input: p an odd prime, and an integer n .
 Step 0: Check that n is indeed a square: (n | p) must be ≡ 1 .
 Step 1: By factoring out powers of 2 from p - 1, find q and s such that p - 1 = q2s with q odd .
 If p ≡ 3 (mod 4) (i.e. s = 1), output the two solutions r ≡ ± n(p+1)/4 .
 Step 2: Select a non-square z such that (z | p) ≡ -1 and set c ≡ zq .
 Step 3: Set r ≡ n(q+1)/2, t ≡ nq, m = s .
 Step 4: Loop the following:
 If t ≡ 1, output r and p - r .
 Otherwise find, by repeated squaring, the lowest i, 0 < i < m , such that t2i ≡ 1 .
 Let b ≡ c2(m - i - 1), and set r ≡ rb, t ≡ tb2, c ≡ b2 and m = i .




Task

Implement the above algorithm.

Find solutions (if any) for

 n = 10 p = 13
 n = 56 p = 101
 n = 1030 p = 10009
 n = 1032 p = 10009
 n = 44402 p = 100049 
Extra credit

 n = 665820697 p = 1000000009 
 n = 881398088036 p = 1000000000039 
 n = 41660815127637347468140745042827704103445750172002 p = 10^50 + 577 


See also

 Modular exponentiation
 Cipolla's algorithm


' 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