Tonelli-Shanks algorithm

Tonelli–Shanks algorithm

Tonelli-Shanks algorithm is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

In computational number theory, the Tonelli–Shanks algorithm is a technique for solving an equation of the form:

x² ≡ n (mod p)

─── where   n   is an integer which is a quadratic residue (mod p),   p   is an odd prime,   and

x,n   Є   Fp = {0, 1, ... p-1}

It is used in cryptography techniques.

To apply the algorithm, we need the Legendre symbol.

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

Algorithm pseudo-code

(copied from Wikipedia):

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. [Factors out powers of 2 from p-1] Define q -odd- and s such as p-1 = q * 2^s
  • if s = 1 , i.e p ≡ 3 (mod 4) , output the two solutions r ≡ +/- n^((p+1)/4) .
• Step 2. Select a non-square z such as (z | p) = -1 , and set c ≡ z^q .
• Step 3. Set r ≡ n ^((q+1)/2) , t ≡ n^q, m = s .
• Step 4. Loop.
  • if t ≡ 1 output r, p-r .
  • Otherwise find, by repeated squaring, the lowest i , 0 < i< m , such as t^(2^i) ≡ 1
  • Let b ≡ c^(2^(m-i-1)), and set r ≡ r*b, t ≡ t*b^2 , c ≡ b^2 and m = i.

Numerical Example

• n=10, p= 13.   See Wikipedia

Task

Implement the above.

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

EchoLisp

<lang scheme> (require 'bigint)

test equality mod p

(define-syntax-rule (mod= a b p) (zero? (% (- a b) p)))

assign mod p

(define-syntax-rule (mod:≡ s v p) (set! s (% v p)))

(define (Legendre a p) (powmod a (/ (1- p) 2) p))

(define (Tonelli n p)

   (unless (= 1 (Legendre n p)) (error "not a square (mod p)" (list n p)))
   (define q (1- p))
   (define s 0)

(while (even? q) (/= q 2) (++ s)) (if (= s 1) (powmod n (/ (1+ p) 4) p) (begin (define z (for ((z (in-range 2 p))) #:break (= (1- p) (Legendre z p)) => z ))

(define c (powmod z q p)) (define r (powmod n (/ (1+ q) 2) p)) (define t (powmod n q p)) (define m s) (define t2 0) (while #t #:break (mod= 1 t p) => r (mod:≡ t2 (* t t) p) (define i (for ((i (in-range 1 m))) #:break (mod= t2 1 p) => i (mod:≡ t2 (* t2 t2) p))) (define b (powmod c (expt 2 (- m i 1)) p)) (mod:≡ r (* r b) p) (mod:≡ c (* b b) p) (mod:≡ t (* t c) p) (set! m i))))) </lang>

Output:

(define ttest 
	`((10 13) (56 101) (1030 10009) (44402 100049)  
	(665820697 1000000009) 
	(881398088036  1000000000039)
	(41660815127637347468140745042827704103445750172002  ,(+ 1e50 577))))  
	     	
(define (task ttest)
	(for ((test ttest))
		(define n (first test))
		(define p (second test))
		(define r (Tonelli n p))
		(assert (mod= (* r r) n p))
		(printf "n = %d p = %d" n p)
		(printf "\t  roots : %d %d"  r (- p r))))

(task ttest)
n = 10 p = 13
  roots : 7 6
n = 56 p = 101
  roots : 37 64
n = 1030 p = 10009
  roots : 1632 8377
n = 44402 p = 100049
  roots : 30468 69581
n = 665820697 p = 1000000009
  roots : 378633312 621366697
n = 881398088036 p = 1000000000039
  roots : 791399408049 208600591990
n = 41660815127637347468140745042827704103445750172002 
p = 100000000000000000000000000000000000000000000000577
  roots : 32102985369940620849741983987300038903725266634508    
  67897014630059379150258016012699961096274733366069
(Tonelli 1032 10009)
❌ error: not a square (mod p) (1032 10009)

FreeBASIC

LongInt version

<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.

1. Ifndef TRUE

   #Define FALSE 0
   #Define TRUE Not FALSE

1. 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>

Output:

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

GMP version

Library: GMP

<lang freebasic>' version 12-04-2017 ' compile with: fbc -s console

1. Include Once "gmp.bi"

Data "10", "13", "56", "101", "1030", "10009", "1032", "10009" Data "44402", "100049", "665820697", "1000000009" Data "881398088036", "1000000000039" Data "41660815127637347468140745042827704103445750172002" ' p = 10^50 + 577

' ------=< MAIN >=------

Dim As uLong k Dim As ZString Ptr zstr Dim As String n_str, p_str

Dim As Mpz_ptr b, c, i, m, n, p, q, r, s, t, z, tmp b = Allocate(Len(__Mpz_struct)) : Mpz_init(b) c = Allocate(Len(__Mpz_struct)) : Mpz_init(c) i = Allocate(Len(__Mpz_struct)) : Mpz_init(i) m = Allocate(Len(__Mpz_struct)) : Mpz_init(m) n = Allocate(Len(__Mpz_struct)) : Mpz_init(n) p = Allocate(Len(__Mpz_struct)) : Mpz_init(p) q = Allocate(Len(__Mpz_struct)) : Mpz_init(q) r = Allocate(Len(__Mpz_struct)) : Mpz_init(r) s = Allocate(Len(__Mpz_struct)) : Mpz_init(s) t = Allocate(Len(__Mpz_struct)) : Mpz_init(t) z = Allocate(Len(__Mpz_struct)) : Mpz_init(z) tmp = Allocate(Len(__Mpz_struct)) : Mpz_init(tmp)

For k = 1 To 8

   Read n_str
   Mpz_set_str(n, n_str, 10)
   If k < 8 Then
       Read p_str
       Mpz_set_str(p, p_str, 10)
   Else
       p_str = "10^50 + 577"
       Mpz_set_str(p, "1" + String(50, "0"), 10)
       Mpz_add_ui(p, p, 577)
   End If

   Print "Find solution for n = "; n_str; " and p = "; p_str

   If Mpz_legendre(n, p) <> 1 Then
       Print n_str; " is not a quadratic residue"
       Print
       Continue For
   End If

   If Mpz_tstbit(p, 0) = 0 OrElse Mpz_probab_prime_p(p, 20) = 0 Then
       Print p_str; "is not a odd prime"
       Print
       Continue For
   End If

   Mpz_set_ui(s, 0) : Mpz_set(q, p) : Mpz_sub_ui(q, q, 1) ' q = p -1
   Do
       Mpz_add_ui(s, s, 1)
       Mpz_fdiv_q_2exp(q, q, 1)
   Loop Until Mpz_tstbit(q, 0) = 1

   If Mpz_cmp_ui(s, 1) = 0 Then
       If Mpz_tstbit(p, 1) = 1 Then
           Mpz_add_ui(tmp, p, 1)
           Mpz_fdiv_q_2exp(tmp, tmp, 2)         ' tmp = p +1 \ 4
           Mpz_powm(r, n, tmp, p)
           zstr = Mpz_get_str(0, 10, r)
           Print "Solution found: "; *zstr;
           Mpz_sub(r, p, r)
           zstr = Mpz_get_str(0, 10, r)
           Print " and "; *zstr
           Print
           Continue For
       End If
   End If

   Mpz_set_ui(z, 1)
   Do
       Mpz_add_ui(z, z, 1)
   Loop Until Mpz_legendre(z, p) = -1
   Mpz_powm(c, z, q, p)
   Mpz_add_ui(tmp, q, 1)
   Mpz_fdiv_q_2exp(tmp, tmp, 1)
   Mpz_powm(r, n, tmp, p)
   Mpz_powm(t, n, q, p)
   Mpz_set(m, s)

   Do
       Mpz_set_ui(i, 0)
       Mpz_mod(tmp, t, p)
       If Mpz_cmp_ui(tmp, 1) = 0 Then
           zstr = Mpz_get_str(0, 10, r)
           Print "Solution found: "; *zstr;
           Mpz_sub(r, p, r)
           zstr = Mpz_get_str(0, 10, r)
           Print " and "; *zstr
           Print
           Continue For
       End If

       Mpz_set_ui(q, 1)
       Do
           Mpz_add_ui(i, i, 1)
           If Mpz_cmp(i, m) >= 0 Then
               Continue For
           end if
           Mpz_mul_ui(q, q, 2)                  ' q = 2^i
           Mpz_powm(tmp, t, q, p)
       Loop Until Mpz_cmp_ui(tmp, 1) = 0

       Mpz_set_ui(q, 2)
       Mpz_sub(tmp, m, i) : Mpz_sub_ui(tmp, tmp, 1) : Mpz_powm(tmp, q, tmp, p)
       Mpz_powm(b, c, tmp, p)
       Mpz_mul(r, r, b) : Mpz_mod(r, r, p)
       Mpz_mul(tmp, b, b) : Mpz_mod(c, tmp, p)
       Mpz_mul(tmp, t, c) : Mpz_mod(t, tmp, p)
       Mpz_set(m, i)
   Loop

Next

Mpz_clear(b) : Mpz_clear(c) : Mpz_clear(i) : Mpz_clear(m) Mpz_clear(n) : Mpz_clear(p) : Mpz_clear(q) : Mpz_clear(r) Mpz_clear(s) : Mpz_clear(t) : Mpz_clear(z) : Mpz_clear(tmp)

' empty keyboard buffer While InKey <> "" : Wend Print : Print "hit any key to end program" Sleep End</lang>

Output:

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

Find solution for n = 41660815127637347468140745042827704103445750172002 and p = 10^50 + 577
Solution found: 32102985369940620849741983987300038903725266634508 and 67897014630059379150258016012699961096274733366069

Go

int

Implementation following Wikipedia, using similar variable names, and using the int type for simplicity. <lang go>package main

import "fmt"

// Arguments n, p as described in WP // If Legendre symbol != 1, ok return is false. Otherwise ok return is true, // R1 is WP return value R and for convenience R2 is p-R1. func ts(n, p int) (R1, R2 int, ok bool) {

   // a^e mod p
   powModP := func(a, e int) int {
       s := 1
       for ; e > 0; e-- {
           s = s * a % p
       }
       return s
   }
   // Legendre symbol, returns 1, 0, or -1 mod p -- that's 1, 0, or p-1.
   ls := func(a int) int {
       return powModP(a, (p-1)/2)
   }
   // argument validation
   if ls(n) != 1 {
       return 0, 0, false
   }
   // WP step 1, factor out powers two.
   // variables Q, S named as at WP.
   Q := p - 1
   S := 0
   for Q&1 == 0 {
       S++
       Q >>= 1
   }
   // WP step 1, direct solution
   if S == 1 {
       R1 = powModP(n, (p+1)/4)
       return R1, p - R1, true
   }
   // WP step 2, select z, assign c
   z := 2
   for ; ls(z) != p-1; z++ {
   }
   c := powModP(z, Q)
   // WP step 3, assign R, t, M
   R := powModP(n, (Q+1)/2)
   t := powModP(n, Q)
   M := S
   // WP step 4, loop
   for {
       // WP step 4.1, termination condition
       if t == 1 {
           return R, p - R, true
       }
       // WP step 4.2, find lowest i...
       i := 0
       for z := t; z != 1 && i < M-1; {
           z = z * z % p
           i++
       }
       // WP step 4.3, using a variable b, assign new values of R, t, c, M
       b := c
       for e := M - i - 1; e > 0; e-- {
           b = b * b % p
       }
       R = R * b % p
       c = b * b % p // more convenient to compute c before t
       t = t * c % p
       M = i
   }

}

func main() {

   fmt.Println(ts(10, 13))
   fmt.Println(ts(56, 101))
   fmt.Println(ts(1030, 10009))
   fmt.Println(ts(1032, 10009))
   fmt.Println(ts(44402, 100049))

}</lang>

Output:

7 6 true
37 64 true
1632 8377 true
0 0 false
30468 69581 true

big.Int

For the extra credit, we use big.Int from the math/big package of the Go standard library. While the method call syntax is not as easy on the eyes as operator syntax, the package provides modular exponentiation and even the Legendre symbol as the Jacobi function. <lang go>package main

import (

   "fmt"
   "math/big"

)

func ts(n, p big.Int) (R1, R2 big.Int, ok bool) {

   if big.Jacobi(&n, &p) != 1 {
       return
   }
   var one, Q big.Int
   one.SetInt64(1)
   Q.Sub(&p, &one)
   S := 0
   for Q.Bit(0) == 0 {
       S++
       Q.Rsh(&Q, 1)
   }
   if S == 1 {
       R1.Exp(&n, R1.Rsh(R1.Add(&p, &one), 2), &p)
       R2.Sub(&p, &R1)
       return R1, R2, true
   }
   var z, c big.Int
   for z.SetInt64(2); big.Jacobi(&z, &p) != -1; z.Add(&z, &one) {
   }
   c.Exp(&z, &Q, &p)
   var R, t big.Int
   R.Exp(&n, R.Rsh(R.Add(&Q, &one), 1), &p)
   t.Exp(&n, &Q, &p)
   M := S
   for {
       if t.Cmp(&one) == 0 {
           R2.Sub(&p, &R)
           return R, R2, true
       }
       i := 0
       // reuse z as a scratch variable
       for z.Set(&t); z.Cmp(&one) != 0 && i < M-1; {
           z.Mod(z.Mul(&z, &z), &p)
           i++
       }
       // and instead of a new scratch variable b, continue using z
       z.Set(&c)
       for e := M - i - 1; e > 0; e-- {
           z.Mod(z.Mul(&z, &z), &p)
       }
       R.Mod(R.Mul(&R, &z), &p)
       c.Mod(c.Mul(&z, &z), &p)
       t.Mod(t.Mul(&t, &c), &p)
       M = i
   }

}

func main() {

   var n, p big.Int
   n.SetInt64(665820697)
   p.SetInt64(1000000009)
   R1, R2, ok := ts(n, p)
   fmt.Println(&R1, &R2, ok)

   n.SetInt64(881398088036)
   p.SetInt64(1000000000039)
   R1, R2, ok = ts(n, p)
   fmt.Println(&R1, &R2, ok)
   n.SetString("41660815127637347468140745042827704103445750172002", 10)
   p.SetString("100000000000000000000000000000000000000000000000577", 10)
   R1, R2, ok = ts(n, p)
   fmt.Println(&R1)
   fmt.Println(&R2)

}</lang>

Output:

378633312 621366697 true
791399408049 208600591990 true
32102985369940620849741983987300038903725266634508
67897014630059379150258016012699961096274733366069

Library

It gets better; the library has a ModSqrt function that uses Tonelli-Shanks internally. Output is same as above. <lang go>package main

import (

   "fmt"
   "math/big"

)

func main() {

   var n, p, R1, R2 big.Int
   n.SetInt64(665820697)
   p.SetInt64(1000000009)
   R1.ModSqrt(&n, &p)
   R2.Sub(&p, &R1)
   fmt.Println(&R1, &R2)

   n.SetInt64(881398088036)
   p.SetInt64(1000000000039)
   R1.ModSqrt(&n, &p)
   R2.Sub(&p, &R1)
   fmt.Println(&R1, &R2)

   n.SetString("41660815127637347468140745042827704103445750172002", 10)
   p.SetString("100000000000000000000000000000000000000000000000577", 10)
   R1.ModSqrt(&n, &p)
   R2.Sub(&p, &R1)
   fmt.Println(&R1)
   fmt.Println(&R2)

}</lang>

J

Implementation:

<lang J>leg=: dyad define

 x (y&|)@^ (y-1)%2

)

tosh=:dyad define

 assert. 1=1 p: y [ 'y must be prime'
 assert. 1=x leg y [ 'x must be square mod y'
 pow=. y&|@^
 if. 1=m=. {.1 q: y-1 do.
   r=. x pow (y+1)%4 
 else.
   z=. 1x while. 1>: z leg y do. z=.z+1 end.
   c=. z pow q=. (y-1)%2^m
   r=. x pow (q+1)%2
   t=. x pow q
   while. t~:1 do.
     n=. t
     i=. 0
     whilst. 1~:n do.
       n=. n pow 2
       i=. i+1
     end.
     r=. y|r*b=. c pow 2^m-i+1
     m=. i
     t=. y|t*c=. b pow 2
   end.
 end.
 y|(,-)r

)</lang>

Task examples:

<lang J> 10 tosh 13 7 6

  56 tosh 101

37 64

  1030 tosh 10009

1632 8377

  1032 tosh 10009

|assertion failure: tosh | 1=x leg y['x must be square mod y'

  44402 tosh 100049

30468 69581

  665820697x tosh 1000000009x

378633312 621366697

  881398088036 tosh 1000000000039x

791399408049 208600591990

  41660815127637347468140745042827704103445750172002x tosh (10^50x)+577

32102985369940620849741983987300038903725266634508 67897014630059379150258016012699961096274733366069</lang>

Kotlin

Translation of: Go

<lang scala>// version 1.1.3

import java.math.BigInteger

data class Solution(val root1: BigInteger, val root2: BigInteger, val exists: Boolean)

val bigZero = BigInteger.ZERO val bigOne = BigInteger.ONE val bigTwo = BigInteger.valueOf(2L) val bigFour = BigInteger.valueOf(4L) val bigTen = BigInteger.TEN

fun ts(n: Long, p: Long) = ts(BigInteger.valueOf(n), BigInteger.valueOf(p))

fun ts(n: BigInteger, p: BigInteger): Solution {

   fun powModP(a: BigInteger, e: BigInteger) = a.modPow(e, p)

   fun ls(a: BigInteger) = powModP(a, (p - bigOne) / bigTwo)

   if (ls(n) != bigOne) return Solution(bigZero, bigZero, false)
   var q = p - bigOne
   var ss = bigZero
   while (q.and(bigOne) == bigZero) {
       ss = ss + bigOne
       q = q.shiftRight(1)
   }

   if (ss == bigOne) {
       val r1 = powModP(n, (p + bigOne) / bigFour)
       return Solution(r1, p - r1, true)
   }

   var z = bigTwo
   while (ls(z) != p - bigOne) z = z + bigOne
   var c = powModP(z, q)
   var r = powModP(n, (q + bigOne) / bigTwo)
   var t = powModP(n, q)
   var m = ss

   while (true) {
       if (t == bigOne) return Solution(r, p - r, true)
       var i = bigZero
       var zz = t
       while (zz != bigOne && i < m - bigOne) {
           zz  = zz * zz % p
           i = i + bigOne
       }
       var b = c
       var e = m - i - bigOne
       while (e > bigZero) {
           b = b * b % p
           e = e - bigOne
       }
       r = r * b % p
       c = b * b % p
       t = t * c % p
       m = i
   }

}

fun main(args: Array<String>) {

   val pairs = listOf<Pair<Long, Long>>(
       10L to 13L, 
       56L to 101L, 
       1030L to 10009L,
       1032L to 10009L,
       44402L to 100049L,
       665820697L to 1000000009L,
       881398088036L to 1000000000039L
   )

   for (pair in pairs) {
       val (n, p) = pair
       val (root1, root2, exists) = ts(n, p)
       println("n = $n")
       println("p = $p")
       if (exists) {
           println("root1 = $root1")
           println("root2 = $root2")
       }
       else println("No solution exists")
       println()
   }

   val bn = BigInteger("41660815127637347468140745042827704103445750172002")
   val bp = bigTen.pow(50) + BigInteger.valueOf(577L)
   val (broot1, broot2, bexists) = ts(bn, bp)
   println("n = $bn")
   println("p = $bp")
   if (bexists) {
       println("root1 = $broot1")
       println("root2 = $broot2")
   }
   else println("No solution exists")    

}</lang>

Output:

n = 10
p = 13
root1 = 7
root2 = 6

n = 56
p = 101
root1 = 37
root2 = 64

n = 1030
p = 10009
root1 = 1632
root2 = 8377

n = 1032
p = 10009
No solution exists

n = 44402
p = 100049
root1 = 30468
root2 = 69581

n = 665820697
p = 1000000009
root1 = 378633312
root2 = 621366697

n = 881398088036
p = 1000000000039
root1 = 791399408049
root2 = 208600591990

n = 41660815127637347468140745042827704103445750172002
p = 100000000000000000000000000000000000000000000000577
root1 = 32102985369940620849741983987300038903725266634508
root2 = 67897014630059379150258016012699961096274733366069

Perl 6

Works with: Rakudo version 2016.10

Translation of the Wikipedia pseudocode, heavily influenced by Sidef and Python.

<lang perl6># Legendre operator (𝑛│𝑝) sub infix:<│> (Int \𝑛, Int \𝑝 where 𝑝.is-prime && (𝑝 != 2)) {

   given 𝑛.expmod( (𝑝-1) div 2, 𝑝 ) {
       when 0  {  0 }
       when 1  {  1 }
       default { -1 }
   }

}

sub tonelli-shanks ( \𝑛, \𝑝 where (𝑛│𝑝) > 0 ) {

   my $𝑄 = 𝑝 - 1;
   my $𝑆 = 0;
   $𝑄 +>= 1 and $𝑆++ while $𝑄 %% 2;
   return 𝑛.expmod((𝑝+1) div 4, 𝑝) if $𝑆 == 1;
   my $𝑐 = ((2..𝑝).first: (*│𝑝) < 0).expmod($𝑄, 𝑝);
   my $𝑅 = 𝑛.expmod( ($𝑄+1) +> 1, 𝑝 );
   my $𝑡 = 𝑛.expmod( $𝑄, 𝑝 );
   while ($𝑡-1) % 𝑝 {
       my $b;
       my $𝑡2 = $𝑡² % 𝑝;
       for 1 .. $𝑆 {
           if ($𝑡2-1) %% 𝑝 {
               $b = $𝑐.expmod(1 +< ($𝑆-1-$_), 𝑝);
               $𝑆 = $_;
               last;
           }
           $𝑡2 = $𝑡2² % 𝑝;
       }
       $𝑅 = ($𝑅 * $b) % 𝑝;
       $𝑐 = $b² % 𝑝;
       $𝑡 = ($𝑡 * $𝑐) % 𝑝;
   }
   $𝑅;

}

my @tests = (

   (10, 13),
   (56, 101),
   (1030, 10009),
   (1032, 10009),
   (44402, 100049),
   (665820697, 1000000009),
   (881398088036, 1000000000039),
   (41660815127637347468140745042827704103445750172002,
     100000000000000000000000000000000000000000000000577)

);

for @tests -> ($n, $p) {
   try my $t = tonelli-shanks($n, $p);
   say "No solution for ({$n}, {$p})." and next if !$t or ($t² - $n) % $p;
   say "Roots of $n are ($t, {$p-$t}) mod $p";

}</lang>

Output:

Roots of 10 are (7, 6) mod 13
Roots of 56 are (37, 64) mod 101
Roots of 1030 are (1632, 8377) mod 10009
No solution for (1032, 10009).
Roots of 44402 are (30468, 69581) mod 100049
Roots of 665820697 are (378633312, 621366697) mod 1000000009
Roots of 881398088036 are (791399408049, 208600591990) mod 1000000000039
Roots of 41660815127637347468140745042827704103445750172002 are (32102985369940620849741983987300038903725266634508, 67897014630059379150258016012699961096274733366069) mod 100000000000000000000000000000000000000000000000577

PicoLisp

Translation of: Go

<lang PicoLisp># from @lib/rsa.l (de **Mod (X Y N)

  (let M 1
     (loop
        (when (bit? 1 Y)
           (setq M (% (* M X) N)) )
        (T (=0 (setq Y (>> 1 Y)))
           M )
        (setq X (% (* X X) N)) ) ) )

(de legendre (N P)

  (**Mod N (/ (dec P) 2) P) )

(de ts (N P)

  (and
     (=1 (legendre N P))
     (let
        (Q (dec P)
           S 0
           Z 0
           C 0
           R 0
           D 0
           M 0
           B 0
           I 0 )
        (until (bit? 1 Q)
           (setq Q (>> 1 Q))
           (inc 'S) )
        (if (=1 S)
           (list
              (setq @@ (**Mod N (/ (inc P) 4) P))
              (- P @@) )
           (setq Z 2)
           (until (= (legendre Z P) (dec P))
              (inc 'Z) )
           (setq
              C (**Mod Z Q P)
              R (**Mod N (/ (inc Q) 2) P)
              D (**Mod N Q P)
              M S )
           (until (=1 D)
              (zero I)
              (for
                 (Z
                    D
                    (and (<> Z 1) (< I (dec M)))
                    (setq Z (% (* Z Z) P)) )
                 (inc 'I) )
              (setq B C)
              (for
                 (Z
                    (- M I 1)
                    (> Z 0) (dec Z) )
                 (setq B (% (* B B) P)) )
              (setq
                 R (% (* R B) P)
                 C (% (* B B) P)
                 D (% (* D C) P)
                 M I ) )
           (list R (- P R)) ) ) ) )

(println (ts 10 13)) (println (ts 56 101)) (println (ts 1030 10009)) (println (ts 1032 10009)) (println (ts 44402 100049)) (println (ts 665820697 1000000009)) (println (ts 881398088036 1000000000039)) (println (ts 41660815127637347468140745042827704103445750172002 (+ (** 10 50) 577)))</lang>

Output:

(7 6)
(37 64)
(1632 8377)
NIL
(30468 69581)
(378633312 621366697)
(791399408049 208600591990)
(32102985369940620849741983987300038903725266634508 67897014630059379150258016012699961096274733366069)

Python

Translation of: EchoLisp

Works with: Python version 3

<lang python>def legendre(a, p):

   return pow(a, (p - 1) // 2, p)

def tonelli(n, p):

   assert legendre(n, p) == 1, "not a square (mod p)"
   q = p - 1
   s = 0
   while q % 2 == 0:
       q //= 2
       s += 1
   if s == 1:
       return pow(n, (p + 1) // 4, p)
   for z in range(2, p):
       if p - 1 == legendre(z, p):
           break
   c = pow(z, q, p)
   r = pow(n, (q + 1) // 2, p)
   t = pow(n, q, p)
   m = s
   t2 = 0
   while (t - 1) % p != 0:
       t2 = (t * t) % p
       for i in range(1, m):
           if (t2 - 1) % p == 0:
               break
           t2 = (t2 * t2) % p
       b = pow(c, 1 << (m - i - 1), p)
       r = (r * b) % p
       c = (b * b) % p
       t = (t * c) % p
       m = i
   return r

if __name__ == '__main__':

   ttest = [(10, 13), (56, 101), (1030, 10009), (44402, 100049),

(665820697, 1000000009), (881398088036, 1000000000039),

            (41660815127637347468140745042827704103445750172002, 10**50 + 577)]
   for n, p in ttest:
       r = tonelli(n, p)
       assert (r * r - n) % p == 0
       print("n = %d p = %d" % (n, p))
       print("\t  roots : %d %d" % (r, p - r))</lang>

Output:

n = 10 p = 13
	  roots : 7 6
n = 56 p = 101
	  roots : 37 64
n = 1030 p = 10009
	  roots : 1632 8377
n = 44402 p = 100049
	  roots : 30468 69581
n = 665820697 p = 1000000009
	  roots : 378633312 621366697
n = 881398088036 p = 1000000000039
	  roots : 791399408049 208600591990
n = 41660815127637347468140745042827704103445750172002 p = 100000000000000000000000000000000000000000000000577
	  roots : 32102985369940620849741983987300038903725266634508 67897014630059379150258016012699961096274733366069

Racket

Translation of: EchoLisp

<lang racket>#lang racket

(require math/number-theory)

(define (Legendre a p)

 (modexpt a (quotient (sub1 p) 2)))

(define (Tonelli n p (err (λ (n p) (error "not a square (mod p)" (list n p)))))

 (with-modulus p
   (unless (= 1 (Legendre n p)) (err n p))

   (define-values (q s)
     (let even?-q-loop ((q (sub1 p)) (s 0))
       (if (even? q)
           (even?-q-loop (quotient q 2) (add1 s))
           (values q s))))
   
   (cond
     [(= s 1)
      (modexpt n (/ (add1 p) 4))]
     [else
      (define z (for/first ((z (in-range 2 p)) #:when (= (sub1 p) (Legendre z p))) z)) 
      (let loop ((c (modexpt z q))
                 (r (modexpt n (quotient (add1 q) 2)))
                 (t (modexpt n q))
                 (m s))
        (cond
          [(mod= 1 t)
           r]
          [else
           (define-values (t2 m′) (for/fold ((t2 (modsqr t)) (i 1))
                                            ((j (in-range 1 m)) #:final (mod= t2 1))
                                    (values (modsqr t2) j)))
           (define b (modexpt c (expt 2 (- m m′ 1))))
           (define c′ (modsqr b))
           (loop c′ (mod* r b) (mod* t c′) m′)]))])))

(module+ test

 (require rackunit)

 (define ttest 
   `((10 13)
     (56 101)
     (1030 10009)
     (44402 100049)  
     (665820697 1000000009) 
     (881398088036  1000000000039)
     (41660815127637347468140745042827704103445750172002
      ,(+ #e1e50 577))))  

 (define (task ttest)
   (for ((test ttest))
     (define n (first test))
     (define p (second test))
     (define r (Tonelli n p))
     (printf "n = ~a p = ~a~%  roots : ~a ~a~%" n p r (- p r))))

 (task ttest)

 (check-exn exn:fail? (λ () (Tonelli 1032 1009))))</lang>

Output:

n = 10 p = 13
  roots : 7 6
n = 56 p = 101
  roots : 37 64
n = 1030 p = 10009
  roots : 1632 8377
n = 44402 p = 100049
  roots : 30468 69581
n = 665820697 p = 1000000009
  roots : 378633312 621366697
n = 881398088036 p = 1000000000039
  roots : 791399408049 208600591990
n = 41660815127637347468140745042827704103445750172002 p = 100000000000000000000000000000000000000000000000577
  roots : 32102985369940620849741983987300038903725266634508 67897014630059379150258016012699961096274733366069

Sidef

Translation of: Python

<lang ruby>func tonelli(n, p) {

   legendre(n, p) == 1 || die "not a square (mod p)"
   var q = p-1
   var s = valuation(q, 2)
   s == 1 ? return(powmod(n, (p + 1) >> 2, p)) : (q >>= s)
   var c = powmod(2 ..^ p -> first {|z| legendre(z, p) == -1}, q, p)
   var r = powmod(n, (q + 1) >> 1, p)
   var t = powmod(n, q, p)
   var m = s
   var t2 = 0
   while (!p.divides(t - 1)) {
       t2 = ((t * t) % p)
       var b
       for i in (1 ..^ m) {
           if (p.divides(t2 - 1)) {
               b = powmod(c, 1 << (m - i - 1), p)
               m = i
               break
           }
           t2 = ((t2 * t2) % p)
       }

       r = ((r * b) % p)
       c = ((b * b) % p)
       t = ((t * c) % p)
   }
   return r

}

var tests = [

   [10, 13], [56, 101], [1030, 10009], [44402, 100049],
   [665820697, 1000000009], [881398088036, 1000000000039],
   [41660815127637347468140745042827704103445750172002, 10**50 + 577],

]

for n,p in tests {

   var r = tonelli(n, p)
   assert((r*r - n) % p == 0)
   say "Roots of #{n} are (#{r}, #{p-r}) mod #{p}"

}</lang>

Output:

Roots of 10 are (7, 6) mod 13
Roots of 56 are (37, 64) mod 101
Roots of 1030 are (1632, 8377) mod 10009
Roots of 44402 are (30468, 69581) mod 100049
Roots of 665820697 are (378633312, 621366697) mod 1000000009
Roots of 881398088036 are (791399408049, 208600591990) mod 1000000000039
Roots of 41660815127637347468140745042827704103445750172002 are (32102985369940620849741983987300038903725266634508, 67897014630059379150258016012699961096274733366069) mod 100000000000000000000000000000000000000000000000577

zkl

Translation of: EchoLisp

<lang zkl>var BN=Import("zklBigNum"); fcn modEq(a,b,p) { (a-b)%p==0 } fcn legendre(a,p){ a.powm((p - 1)/2,p) }

fcn tonelli(n,p){ //(BigInt,Int|BigInt)

  _assert_(legendre(n,p)==1, "not a square (mod p)"+vm.arglist);
  q,s:=p-1,0;
  while(q.isEven){ q/=2; s+=1; }
  if(s==1) return(n.powm((p+1)/4,p));
  z:=[BN(2)..p].filter1('wrap(z){ legendre(z,p)==(p-1) });
  c,r,t,m,t2:=z.powm(q,p), n.powm((q+1)/2,p), n.powm(q,p), s, 0;
  while(not modEq(t,1,p)){
     t2=(t*t)%p;
     i:=1; while(not modEq(t2,1,p)){ i+=1; t2=(t2*t2)%p; } // assert(i<m)
     b:=c.powm(BN(1).shiftLeft(m-i-1), p);
     r,c,t,m = (r*b)%p, (b*b)%p, (t*c)%p, i;
  }
  r

}</lang> <lang zkl>ttest:=T(T(10,13), T(56,101), T(1030,10009), T(44402,100049),

  T(665820697,1000000009), T(881398088036,1000000000039),
  T("41660815127637347468140745042827704103445750172002", BN(10).pow(50) + 577),
  T(1032,10009) );

foreach n,p in (ttest){ n=BN(n);

  r:=tonelli(n,p);
  assert((r*r-n)%p == 0,"(r*r-n)%p == 0 : %s,%s,%s-->%s".fmt(r,n,p,(r*r-n)%p));
  println("n=%d p=%d".fmt(n,p));
  println("   roots: %d %d".fmt(r, p-r));

}</lang>

Output:

n=10 p=13
   roots: 7 6
n=56 p=101
   roots: 37 64
n=1030 p=10009
   roots: 1632 8377
n=44402 p=100049
   roots: 30468 69581
n=665820697 p=1000000009
   roots: 378633312 621366697
n=881398088036 p=1000000000039
   roots: 791399408049 208600591990
n=41660815127637347468140745042827704103445750172002 p=100000000000000000000000000000000000000000000000577
   roots: 32102985369940620849741983987300038903725266634508 67897014630059379150258016012699961096274733366069
VM#1 caught this unhandled exception:
   AssertionError : not a square (mod p)L(1032,10009)
Stack trace for VM#1 ():
   bbb.assert addr:13  args(2) reg(0) 
   bbb.tonelli addr:29  args(2) reg(10) R
...