Tonelli-Shanks algorithm

Tonelli-Shanks algorithm
You are encouraged to solve this task according to the task description, using any language you may know.

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:

x² ≡ n (mod p)

where n is an integer which is a quadratic residue (mod p), p is an odd prime, and x,n ∈ F_p where F_p = {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 = q2^s 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 ≡ z^q .
• Step 3: Set r ≡ n^(q+1)/2, t ≡ n^q, 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 t²ⁱ ≡ 1 .
  • Let b ≡ c^{2(m - i - 1)}, and set r ≡ rb, t ≡ tb², c ≡ b² 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

11l

Translation of: Python

F legendre(a, p)
   R pow(a, (p - 1) I/ 2, p)

F tonelli(n, p)
   assert(legendre(n, p) == 1, ‘not a square (mod p)’)
   V q = p - 1
   V s = 0
   L q % 2 == 0
      q I/= 2
      s++
   I s == 1
      R pow(n, (p + 1) I/ 4, p)
   V z = 2
   L
      I p - 1 == legendre(z, p)
         L.break
      z++
   V c = pow(z, q, p)
   V r = pow(n, (q + 1) I/ 2, p)
   V t = pow(n, q, p)
   V m = s
   V t2 = BigInt(0)
   L (t - 1) % p != 0
      t2 = (t * t) % p
      V i = 1
      L(ii) 1 .< m
         I (t2 - 1) % p == 0
            i = ii
            L.break
         t2 = (t2 * t2) % p
      V b = pow(c, Int64(1 << (m - i - 1)), p)
      r = (r * b) % p
      c = (b * b) % p
      t = (t * c) % p
      m = i
   R r

V ttest = [(BigInt(10), BigInt(13)), (BigInt(56), BigInt(101)), (BigInt(1030), BigInt(10009)), (BigInt(44402), BigInt(100049)),
           (BigInt(665820697), BigInt(1000000009)), (BigInt(881398088036), BigInt(1000000000039)),
           (BigInt(‘41660815127637347468140745042827704103445750172002’), BigInt(10) ^ 50 + 577)]
L(n, p) ttest
   V r = tonelli(n, p)
   assert((r * r - n) % p == 0)
   print(‘n = #. p = #.’.format(n, p))
   print("\t  roots : #. #.".format(r, p - r))

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

AArch64 Assembly

Works with: as version Raspberry Pi 3B version Buster 64 bits
or android 64 bits with application Termux

/* ARM assembly AARCH64 Raspberry PI 3B or android 64 bits */
/* program tonshan64.s   */

/*******************************************/
/* Constantes file                         */
/*******************************************/
/* for this file see task include a file in language AArch64 assembly*/
.include "../includeConstantesARM64.inc"

/*******************************************/
/* Initialized data                        */
/*******************************************/
.data
szMessStartPgm:          .asciz "Program 64 bits start \n"
szMessEndPgm:            .asciz "Program normal end.\n"
szMessError:             .asciz "\033[31mError  !!!\n"
szMessErrGen:            .asciz "Error end program.\n"
szMessOverflow:          .asciz "Overflow function modulo.\n"
szMessNoSolution:        .asciz "No solution.\n"
szCarriageReturn:        .asciz "\n"

/* datas message display */
szMessEntry:             .asciz "Number : @ modulo : @ ==> "
szMessResult:            .asciz "Racine 1 : @ Racine 2 : @  \n"


qNumberN:            .quad 44402
qNumberP:            .quad 100049
/*******************************************/
/* UnInitialized data                      */
/*******************************************/
.bss 
.align 4
sZoneConv:               .skip 24
/*******************************************/
/*  code section                           */
/*******************************************/
.text
.global main 
main:                               // program start
    ldr x0,qAdrszMessStartPgm       // display start message
    bl affichageMess
    
    mov x0,10
    mov x1,13 
    bl displayEntry
    bl computeTonSha
    bl displayResult
    
    mov x0,56
    mov x1,101
    bl displayEntry
    bl computeTonSha
    bl displayResult
    
    mov x0,1030
    mov x1,10009
    bl displayEntry
    bl computeTonSha
    bl displayResult
    
    mov x0,1032
    mov x1,10009
    bl displayEntry
    bl computeTonSha
    bcs 1f
    bl displayResult
1:
    ldr x4,qAdrqNumberN
    ldr x0,[x4]
    ldr x4,qAdrqNumberP
    ldr x1,[x4]
    bl displayEntry
    bl computeTonSha
    bl displayResult

    ldr x0,qAdrszMessEndPgm         // display end message
    bl affichageMess
    b 100f
99:                                 // display error message 
    ldr x0,qAdrszMessError
    bl affichageMess
100:                                // standard end of the program
    mov x0, #0                      // return code
    mov x8, #EXIT                   // request to exit program
    svc 0                           // perform system call
qAdrszMessStartPgm:        .quad szMessStartPgm
qAdrszMessEndPgm:          .quad szMessEndPgm
qAdrszMessError:           .quad szMessError
qAdrszMessNoSolution:      .quad szMessNoSolution
qAdrszCarriageReturn:      .quad szCarriageReturn
qAdrqNumberN:              .quad qNumberN
qAdrqNumberP:              .quad qNumberP

qAdrszMessResult:          .quad szMessResult
qAdrsZoneConv:             .quad sZoneConv

/******************************************************************/
/*     algorithm Tonelli–Shanks                                   */ 
/******************************************************************/
/* x0 contains number */
/* x1 contains modulus */
/* x0 return root 1 */
/* x1 return root 2 */
computeTonSha:
    stp x10,lr,[sp,-16]!    // save  registres
    stp x2,x3,[sp,-16]!     // save  registres
    stp x4,x5,[sp,-16]!     // save  registres
    stp x6,x7,[sp,-16]!     // save  registres
    stp x8,x9,[sp,-16]!     // save  registres
    stp x11,x12,[sp,-16]!   // save  registres
    mov x9,x0               // save number
    mov x10,x1              // save modulo p
    mov x2,x10
    sub x1,x2,1
    lsr x1,x1,1
    bl moduloPuR64
    bcs 100f                // error ?
    cmp x0,#1
    bne 20f
    sub x5,x10,1
    mov x6,#1               // s
1:
    lsr x5,x5,#1            // div by 2
    tst x5,1                // even ?
    cinc x6,x6,eq           // yes count
    beq 1b                  // and loop
                            // x5 = q
    cmp x6,#1               // s = 1 ?
    bne 3f
    add x1,x10,1            // compute root 1
    lsr x1,x1,#2            // p + 1 / 4
    mov x0,x9               // n
    mov x2,x10              // p
    bl moduloPuR64
    bcs 100f                // error ?
    neg x1,x0               // compute root 2 = - root 1
    b 100f                  // and end
3:
    mov x7,#3               // z  
4:
    mov x0,x7
    mov x2,x10              //  p
    sub x1,x2,1 
    lsr x1,x1,1             // power = p - 1 / 2
    bl moduloPuR64
    bcs 100f                // error ?
    cmp x0,#1
    cinc x7,x7,eq           // si égal à 1
    cinc x7,x7,eq
    beq 4b
    cmp x0,0          
    cinc x7,x7,eq           // si egal à 0
    cinc x7,x7,eq
    beq 4b
    mov x0,x7               // z
    mov x1,x5               // q
    mov x2,x10              // p
    bl moduloPuR64
    bcs 100f                // error ?
    mov x12,x0              // c = z pow q mod p

    add x1,x5,1             // = q +1
    lsr x1,x1,1             // div 2
    mov x0,x9               // n
    mov x2,x10              // p
    bl moduloPuR64
    mov x4,x0               // r =  n puis (q+1)/2 mod p
    
    mov x0,x9               // n
    mov x1,x5               // = q
    mov x2,x10              // p
    bl moduloPuR64
    bcs 100f                // error ?
    mov x5,x0               // reuse r5 = t = n pow q mod p

8:                          // begin loop
    cmp x5,1
    beq 10f
    mov x0,x5               // t
    mov x1,x6               // m
    mov x2,x10              // p
    bl searchI              // search i for t puis 2 puis i = 1 mod p
    cmp x0,-1               // not find -> no solution
    beq 20f
    mov x9,x0               // i
    sub x8,x6,x0            // compute b
    sub x8,x8,1             // m - i - 1
    mov x1,1
    lsl x1,x1,x8
    mov x0,x12
    mov x2,x10              // p
    bl moduloPuR64
    bcs 100f                // error ?
    mov x7,x0               // b = c puis 2 puis 2 puis m-i-1  à verifier
    
    mul x0,x7,x4            // r = r * b mod p
    umulh x1,x7,x4
    mov x2,x10
    bl divisionReg128U
    mov x4,x3               // r mod p  
    mul x0,x7,x7
    umulh x1,x7,x7
    mov x2,x10
    bl divisionReg128U
    mov x12,x3              // c mod p  

    mul x0,x5,x12
    umulh x1,x5,x12
    mov x2,x10
    bl divisionReg128U
    mov x5,x3               // t mod p  
     
    mov x6,x9               // m = i
    b 8b
9:

10:
    mov x0,x4               // x0 return root 1
    sub x1,x10,x0           //  x1 return root 2
    cmn x0,0                // carry à zero roots ok
   b 100f
20:
    ldr x0,qAdrszMessNoSolution
    bl affichageMess
    
    mov x0,0
    mov x1,0
    cmp x0,0               // carry to 1 No solution
100:
    ldp x11,x12,[sp],16 
    ldp x8,x9,[sp],16 
    ldp x6,x7,[sp],16 
    ldp x4,x5,[sp],16 
    ldp x2,x3,[sp],16 
    ldp x10,lr,[sp],16      // restaur des  2 registres
    ret                     // retour adresse lr x30

/******************************************************************/
/*     search i                                               */ 
/******************************************************************/
// x0 contains t
// x1 contains maxi
// x2 contains modulo
searchI:
    stp x1,lr,[sp,-16]!    // save  registres
    stp x2,x3,[sp,-16]!    // save  registres
    stp x4,x5,[sp,-16]!    // save  registres
    stp x6,x7,[sp,-16]!    // save  registres
    mov x4,x0              // t
    mov x6,x1              // m
    mov x3,1               // i
1:
    mov x5,1
    lsl x5,x5,x3           // compute 2 power i

    mov x0,x4
    mov x1,x5
    bl moduloPuR64         // compute t pow 2 pow i mod p
    cmp x0,1               // = 1 ?
    beq 3f                 // yes it is ok
    add x3,x3,1            // next i
    cmp x3,x6
    blt 1b                 // loop 
    mov x0,-1              // not find
    b 100f
3:
    mov x0,x3              // return i 
100:
    ldp x6,x7,[sp],16 
    ldp x4,x5,[sp],16 
    ldp x2,x3,[sp],16 
    ldp x1,lr,[sp],16         // restaur des  2 registres
    ret                       // retour adresse lr x30
/******************************************************************/
/*     display numbers                                             */ 
/******************************************************************/
/* x0 contains number */
/* x1 contains modulo */
displayEntry:
    stp x0,lr,[sp,-16]!        // save  registres
    stp x1,x2,[sp,-16]!        // save  registres
    mov x2,x1                  // root 2
    ldr x1,qAdrsZoneConv       // convert root 1 in r0
    bl conversion10S           // convert ascii string
    ldr x0,qAdrszMessEntry
    ldr x1,qAdrsZoneConv
    bl strInsertAtCharInc      // and put in message
    mov x3,x0
    mov x0,x2                  // racine 2
    ldr x1,qAdrsZoneConv
    bl conversion10S           // convert ascii string
    mov x0,x3
    ldr x1,qAdrsZoneConv
    bl strInsertAtCharInc      // and put in message
    bl affichageMess
100:
    ldp x1,x2,[sp],16 
    ldp x0,lr,[sp],16          // restaur des  2 registres
    ret                        // retour adresse lr x30
qAdrszMessEntry:   .quad szMessEntry
/******************************************************************/
/*     display roots                                               */ 
/******************************************************************/
/* x0 contains root 1 */
/* x1 contains root 2 */
displayResult:
    stp x1,lr,[sp,-16]!        // save  registres
    stp x2,x3,[sp,-16]!        // save  registres
    mov x2,x1                  // root 2
    ldr x1,qAdrsZoneConv       // convert root 1 in r0
    bl conversion10S           // convert ascii string
    ldr x0,qAdrszMessResult
    ldr x1,qAdrsZoneConv
    bl strInsertAtCharInc      // and put in message
    mov x3,x0
    mov x0,x2                  // racine 2
    ldr x1,qAdrsZoneConv
    bl conversion10S           // convert ascii string
    mov x0,x3
    ldr x1,qAdrsZoneConv
    bl strInsertAtCharInc      // and put in message
    bl affichageMess
100:
    ldp x2,x3,[sp],16 
    ldp x1,lr,[sp],16          // restaur des  2 registres
    ret                        // retour adresse lr x30
/**************************************************************/
/********************************************************/
/*   Calcul modulo de b puissance e modulo m  */
/*    Exemple 4 puissance 13 modulo 497 = 445         */
/********************************************************/
/* x0  nombre  */
/* x1 exposant */
/* x2 modulo   */
moduloPuR64:
    stp x1,lr,[sp,-16]!        // save  registres
    stp x3,x4,[sp,-16]!        // save  registres
    stp x5,x6,[sp,-16]!        // save  registres
    stp x7,x8,[sp,-16]!        // save  registres
    stp x9,x10,[sp,-16]!        // save  registres
    cbz x0,100f
    cbz x1,100f
    mov x8,x0
    mov x7,x1
    mov x6,1                   // resultat
    udiv x4,x8,x2
    msub x9,x4,x2,x8           // contient le reste
1:
    tst x7,1
    beq 2f
    mul x4,x9,x6
    umulh x5,x9,x6
    mov x6,x4
    mov x0,x6
    mov x1,x5
    bl divisionReg128U
    cbnz x1,99f                // overflow
    mov x6,x3
2:
    mul x8,x9,x9
    umulh x5,x9,x9
    mov x0,x8
    mov x1,x5
    bl divisionReg128U
    cbnz x1,99f                // overflow
    mov x9,x3
    lsr x7,x7,1
    cbnz x7,1b
    cmn x0,0                   // carry à zero pas d'erreur  
    mov x0,x6                  // result
    b 100f
99:
    ldr x0,qAdrszMessOverflow
    bl  affichageMess
    cmp x0,0                   // carry à un car erreur
    mov x0,-1                  // code erreur

100:
    ldp x9,x10,[sp],16          // restaur des  2 registres
    ldp x7,x8,[sp],16          // restaur des  2 registres
    ldp x5,x6,[sp],16          // restaur des  2 registres
    ldp x3,x4,[sp],16          // restaur des  2 registres
    ldp x1,lr,[sp],16          // restaur des  2 registres
    ret                        // retour adresse lr x30
qAdrszMessOverflow:         .quad  szMessOverflow
/***************************************************/
/*   division d un nombre de 128 bits par un nombre de 64 bits */
/***************************************************/
/* x0 contient partie basse dividende */
/* x1 contient partie haute dividente */
/* x2 contient le diviseur */
/* x0 retourne partie basse quotient */
/* x1 retourne partie haute quotient */
/* x3 retourne le reste */
divisionReg128U:
    stp x6,lr,[sp,-16]!        // save  registres
    stp x4,x5,[sp,-16]!        // save  registres
    mov x5,#0                  // raz du reste R
    mov x3,#128                // compteur de boucle
    mov x4,#0                  // dernier bit
1:    
    lsl x5,x5,#1               // on decale le reste de 1
    tst x1,1<<63               // test du bit le plus à gauche
    lsl x1,x1,#1               // on decale la partie haute du quotient de 1
    beq 2f
    orr  x5,x5,#1              // et on le pousse dans le reste R
2:
    tst x0,1<<63
    lsl x0,x0,#1               // puis on decale la partie basse 
    beq 3f
    orr x1,x1,#1               // et on pousse le bit de gauche dans la partie haute
3:
    orr x0,x0,x4               // position du dernier bit du quotient
    mov x4,#0                  // raz du bit
    cmp x5,x2
    blt 4f
    sub x5,x5,x2                // on enleve le diviseur du reste
    mov x4,#1                   // dernier bit à 1
4:
                               // et boucle
    subs x3,x3,#1
    bgt 1b    
    lsl x1,x1,#1               // on decale le quotient de 1
    tst x0,1<<63
    lsl x0,x0,#1              // puis on decale la partie basse 
    beq 5f
    orr x1,x1,#1
5:
    orr x0,x0,x4                  // position du dernier bit du quotient
    mov x3,x5
100:
    ldp x4,x5,[sp],16          // restaur des  2 registres
    ldp x6,lr,[sp],16          // restaur des  2 registres
    ret                        // retour adresse lr x30

/********************************************************/
/*        File Include fonctions                        */
/********************************************************/
/* for this file see task include a file in language AArch64 assembly */
.include "../includeARM64.inc"

Output:

Program 64 bits start
Number : +10 modulo : +13 ==> Racine 1 : +7 Racine 2 : +6
Number : +56 modulo : +101 ==> Racine 1 : +37 Racine 2 : +64
Number : +1030 modulo : +10009 ==> Racine 1 : +1632 Racine 2 : +8377
Number : +1032 modulo : +10009 ==> No solution.
Number : +44402 modulo : +100049 ==> Racine 1 : +30468 Racine 2 : +69581
Program normal end.

ARM Assembly

Works with: as version Raspberry Pi
or android 32 bits with application Termux

/* ARM assembly Raspberry PI  or android 32 bits */
/* program tonshan.s   */

/* REMARK 1 : this program use routines in a include file 
   see task Include a file language arm assembly 
   for the routine affichageMess conversion10 
   see at end of this program the instruction include */
/* for constantes see task include a file in arm assembly */
/************************************/
/* Constantes                       */
/************************************/
.include "../constantes.inc"

/*******************************************/
/* Initialized data                        */
/*******************************************/
.data
szMessStartPgm:          .asciz "Program 32 bits start \n"
szMessEndPgm:            .asciz "Program normal end.\n"
szMessError:             .asciz "\033[31mError  !!!\n"
szMessErrGen:            .asciz "Error end program.\n"
szMessOverflow:          .asciz "Overflow function modulo.\n"
szMessNoSolution:        .asciz "No solution.\n"
szCarriageReturn:        .asciz "\n"

/* datas message display */
szMessEntry:             .asciz "Number : @ modulo : @ ==> "
szMessResult:            .asciz "Racine 1 : @ Racine 2 : @  \n"

iNumberN:                .int 1030
iNumberP:                .int 10009

iNumberN1:               .int 1032
iNumberP1:               .int 10009

iNumberN2:               .int 44402
iNumberP2:               .int 100049

/*******************************************/
/* UnInitialized data                      */
/*******************************************/
.bss 
.align 4
sZoneConv:               .skip 24
/*******************************************/
/*  code section                           */
/*******************************************/
.text
.global main 
main:                               // program start
    ldr r0,iAdrszMessStartPgm       // display start message
    bl affichageMess
    
    mov r0,#10
    mov r1,#13 
    bl displayEntry                 // display entry number
    bl computeTonSha                // compute roots
    bl displayResult                // display roots
    
    mov r0,#56
    mov r1,#101
    bl displayEntry
    bl computeTonSha
    bl displayResult
    

    ldr r4,iAdriNumberN
    ldr r0,[r4]
    ldr r4,iAdriNumberP
    ldr r1,[r4]
    bl displayEntry
    bl computeTonSha
    bl displayResult
    
    ldr r4,iAdriNumberN1
    ldr r0,[r4]
    ldr r4,iAdriNumberP1
    ldr r1,[r4]
    bl displayEntry
    bl computeTonSha
    bcs 1f
    bl displayResult
1:
    ldr r4,iAdriNumberN2
    ldr r0,[r4]
    ldr r4,iAdriNumberP2
    ldr r1,[r4]
    bl displayEntry
    bl computeTonSha
    bl displayResult

    ldr r0,iAdrszMessEndPgm         // display end message
    bl affichageMess
    b 100f
99:                                 // display error message 
    ldr r0,iAdrszMessError
    bl affichageMess
100:                                // standard end of the program
    mov r0, #0                      // return code
    mov r7, #EXIT                   // request to exit program
    svc 0                           // perform system call
iAdrszMessStartPgm:        .int szMessStartPgm
iAdrszMessEndPgm:          .int szMessEndPgm
iAdrszMessError:           .int szMessError
iAdrszMessNoSolution:      .int szMessNoSolution
iAdrszCarriageReturn:      .int szCarriageReturn
iAdriNumberN:              .int iNumberN
iAdriNumberP:              .int iNumberP
iAdriNumberN1:             .int iNumberN1
iAdriNumberP1:             .int iNumberP1
iAdriNumberN2:             .int iNumberN2
iAdriNumberP2:             .int iNumberP2

iAdrszMessResult:          .int szMessResult
iAdrsZoneConv:             .int sZoneConv

/******************************************************************/
/*     algorithm Tonelli–Shanks                                   */ 
/******************************************************************/
/* r0 contains number */
/* r1 contains modulus */
/* r0 return root 1 */
/* r1 return root 2 */
computeTonSha:
    push {r2-r12,lr}

    mov r9,r0               // save number
    mov r10,r1              // save modulo p
    mov r2,r10
    sub r1,r2,#1
    lsr r1,r1,#1
    bl moduloPuR32
    cmp r0,#1
    bne 20f
    sub r5,r10,#1
    mov r6,#1               // s
1:
    lsr r5,r5,#1            // div by 2
    tst r5,#1                // even ?
    addeq r6,#1
    beq 1b                  // and loop
                            // r5 = q
    cmp r6,#1               // s = 1 ?
    bne 3f
    add r1,r10,#1            // compute root 1
    lsr r1,r1,#2            // p + 1 / 4
    mov r0,r9               // n
    mov r2,r10              // p
    bl moduloPuR32
    neg r1,r0               // compute root 2 = - root 1
    b 100f                  // and end
3:
    mov r7,#3               // z  
4:
    mov r0,r7
    mov r2,r10              //  p
    sub r1,r2,#1 
    lsr r1,r1,#1             // power = p - 1 / 2
    bl moduloPuR32
    cmp r0,#1
    addeq r7,#2
    beq 4b
    cmp r0,#0
    addeq r7,#2
    beq 4b
    mov r0,r7               // z
    mov r1,r5               // q
    mov r2,r10              // p
    bl moduloPuR32
    mov r12,r0              // c = z pow q mod p

    add r1,r5,#1             // = q +1
    lsr r1,r1,#1             // div 2
    mov r0,r9               // n
    mov r2,r10              // p
    bl moduloPuR32
    mov r4,r0               // r =  n puis (q+1)/2 mod p
    
    mov r0,r9               // n
    mov r1,r5               // = q
    mov r2,r10              // p
    bl moduloPuR32
    mov r5,r0               // reuse r5 = t = n pow q mod p

8:                          // begin loop
    cmp r5,#1
    beq 10f
    mov r0,r5               // t
    mov r1,r6               // m
    mov r2,r10              // p
    bl searchI              // search i for t puis 2 puis i = 1 mod p
    cmp r0,#-1               // not find -> no solution
    beq 20f
    mov r9,r0               // i
    sub r8,r6,r0            // compute b
    sub r8,r8,#1             // m - i - 1
    mov r1,#1
    lsl r1,r1,r8
    mov r0,r12
    mov r2,r10              // p
    bl moduloPuR32
    mov r7,r0               // b = c puis 2 puis 2 puis m-i-1  à verifier
    
    umull r0,r1,r7,r4            // r = r * b mod p
    mov r2,r10
    bl division32R
    mov r4,r2               // r mod p  
    umull r0,r1,r7,r7
    mov r2,r10
    bl division32R
    mov r12,r2              // c mod p  

    umull r0,r1,r5,r12
    mov r2,r10
    bl division32R
    mov r5,r2               // t mod p  
     
    mov r6,r9               // m = i
    b 8b
9:

10:
    mov r0,r4               // r0 return root 1
    sub r1,r10,r0           //  r1 return root 2
    cmn r0,#0               // carry à zero roots ok
    b 100f
20:
    ldr r0,iAdrszMessNoSolution
    bl affichageMess
    
    mov r0,#0
    mov r1,#0
    cmp r0,#0               // carry to 1 No solution
100:
    pop {r2-r12,lr}         // restaur registers
    bx lr                   // return
/******************************************************************/
/*     search i                                               */ 
/******************************************************************/
// r0 contains t
// r1 contains maxi
// r2 contains modulo
// r0 return i
searchI:
    push {r1-r6,lr}

    mov r4,r0               // t
    mov r6,r1               // m
    mov r3,#1               // i
1:
    mov r5,#1
    lsl r5,r5,r3            // compute 2 power i

    mov r0,r4
    mov r1,r5
    bl moduloPuR32          // compute t pow 2 pow i mod p
    cmp r0,#1               // = 1 ?
    beq 3f                  // yes it is ok
    add r3,r3,#1            // next i
    cmp r3,r6
    blt 1b                  // loop 
    mov r0,#-1              // not find
    b 100f
3:
    mov r0,r3              // return i 
100:
    pop {r1-r6,lr}         // restaur registers
    bx lr                  // return
/******************************************************************/
/*     display numbers                                             */ 
/******************************************************************/
/* r0 contains number */
/* r1 contains modulo */
displayEntry:
    push {r0-r3,lr}
    mov r2,r1                  // root 2
    ldr r1,iAdrsZoneConv       // convert root 1 in r0
    bl conversion10S           // convert ascii string
    ldr r0,iAdrszMessEntry
    ldr r1,iAdrsZoneConv
    bl strInsertAtCharInc      // and put in message
    mov r3,r0
    mov r0,r2                  // racine 2
    ldr r1,iAdrsZoneConv
    bl conversion10S           // convert ascii string
    mov r0,r3
    ldr r1,iAdrsZoneConv
    bl strInsertAtCharInc      // and put in message
    bl affichageMess
100:
    pop {r0-r3,lr}             // restaur registers
    bx lr                      // return
iAdrszMessEntry:   .int szMessEntry
/******************************************************************/
/*     display roots                                               */ 
/******************************************************************/
/* r0 contains root 1 */
/* r1 contains root 2 */
displayResult:
    push {r1-r3,lr}
    mov r2,r1                  // root 2
    ldr r1,iAdrsZoneConv       // convert root 1 in r0
    bl conversion10S           // convert ascii string
    ldr r0,iAdrszMessResult
    ldr r1,iAdrsZoneConv
    bl strInsertAtCharInc      // and put in message
    mov r3,r0
    mov r0,r2                  // racine 2
    ldr r1,iAdrsZoneConv
    bl conversion10S           // convert ascii string
    mov r0,r3
    ldr r1,iAdrsZoneConv
    bl strInsertAtCharInc      // and put in message
    bl affichageMess
100:

    pop {r1-r3,lr}             // restaur registers
    bx lr                      // return
/********************************************************/
/*   Calcul modulo de b puissance e modulo m  */
/*    Exemple 4 puissance 13 modulo 497 = 445         */
/*                                             */
/********************************************************/
/* r0  nombre  */
/* r1 exposant */
/* r2 modulo   */
/* r0 return result  */
moduloPuR32:
    push {r1-r7,lr}    @ save registers  
    cmp r0,#0          @ verif <> zero 
    beq 90f
    cmp r1,#0          @ verif <> zero 
    moveq r0,#0
    beq 90f
    cmp r2,#0          @ verif <> zero 
    moveq r0,#0
    beq 90f            @ 
1:
    mov r4,r2          @ save modulo
    mov r5,r1          @ save exposant 
    mov r6,r0          @ save base
    mov r3,#1          @ start result

    mov r1,#0          @ division de r0,r1 par r2
    bl division32R
    mov r6,r2          @ base <- remainder
2:
    tst r5,#1          @  exposant even or odd
    beq 3f
    umull r0,r1,r6,r3
    mov r2,r4
    bl division32R
    mov r3,r2          @ result <- remainder
3:
    umull r0,r1,r6,r6
    mov r2,r4
    bl division32R
    mov r6,r2          @ base <- remainder

    lsr r5,#1          @ left shift 1 bit
    cmp r5,#0          @ end ?
    bne 2b
    mov r0,r3
90:
    cmn r0,#0          @ no error
100:                   @ fin standard de la fonction
    pop {r1-r7,lr}     @ restaur des registres
    bx lr              @ retour de la fonction en utilisant lr    

/***************************************************/
/*   division number 64 bits in 2 registers by number 32 bits */
/***************************************************/
/* r0 contains lower part dividende   */
/* r1 contains upper part dividende   */
/* r2 contains divisor   */
/* r0 return lower part quotient    */
/* r1 return upper part quotient    */
/* r2 return remainder               */
division32R:
    push {r3-r9,lr}    @ save registers
    mov r6,#0          @ init upper upper part remainder  !!
    mov r7,r1          @ init upper part remainder with upper part dividende
    mov r8,r0          @ init lower part remainder with lower part dividende
    mov r9,#0          @ upper part quotient 
    mov r4,#0          @ lower part quotient
    mov r5,#32         @ bits number
1:                     @ begin loop
    lsl r6,#1          @ shift upper upper part remainder
    lsls r7,#1         @ shift upper  part remainder
    orrcs r6,#1        
    lsls r8,#1         @ shift lower  part remainder
    orrcs r7,#1
    lsls r4,#1         @ shift lower part quotient
    lsl r9,#1          @ shift upper part quotient
    orrcs r9,#1
                       @ divisor sustract  upper  part remainder
    subs r7,r2
    sbcs  r6,#0        @ and substract carry
    bmi 2f             @ négative ?
    
                       @ positive or equal
    orr r4,#1          @ 1 -> right bit quotient
    b 3f
2:                     @ negative 
    orr r4,#0          @ 0 -> right bit quotient
    adds r7,r2         @ and restaur remainder
    adc  r6,#0 
3:
    subs r5,#1         @ decrement bit size 
    bgt 1b             @ end ?
    mov r0,r4          @ lower part quotient
    mov r1,r9          @ upper part quotient
    mov r2,r7          @ remainder
100:                   @ function end
    pop {r3-r9,lr}     @ restaur registers
    bx lr  
/***************************************************/
/*      ROUTINES INCLUDE                           */
/***************************************************/
.include "../affichage.inc"

Output:

Program 32 bits start
Number :         +10 modulo :         +13 ==> Racine 1 :          +7 Racine 2 :          +6
Number :         +56 modulo :        +101 ==> Racine 1 :         +37 Racine 2 :         +64
Number :       +1030 modulo :      +10009 ==> Racine 1 :       +1632 Racine 2 :       +8377
Number :       +1032 modulo :      +10009 ==> No solution.
Number :      +44402 modulo :     +100049 ==> Racine 1 :      +30468 Racine 2 :      +69581
Program normal end.

C

Version 1

Translation of: C#

#include <stdbool.h>
#include <stdint.h>
#include <stdio.h>

uint64_t modpow(uint64_t a, uint64_t b, uint64_t n) {
    uint64_t x = 1, y = a;
    while (b > 0) {
        if (b % 2 == 1) {
            x = (x * y) % n; // multiplying with base
        }
        y = (y * y) % n; // squaring the base
        b /= 2;
    }
    return x % n;
}

struct Solution {
    uint64_t root1, root2;
    bool exists;
};

struct Solution makeSolution(uint64_t root1, uint64_t root2, bool exists) {
    struct Solution sol;
    sol.root1 = root1;
    sol.root2 = root2;
    sol.exists = exists;
    return sol;
}

struct Solution ts(uint64_t n, uint64_t p) {
    uint64_t q = p - 1;
    uint64_t ss = 0;
    uint64_t z = 2;
    uint64_t c, r, t, m;

    if (modpow(n, (p - 1) / 2, p) != 1) {
        return makeSolution(0, 0, false);
    }

    while ((q & 1) == 0) {
        ss += 1;
        q >>= 1;
    }

    if (ss == 1) {
        uint64_t r1 = modpow(n, (p + 1) / 4, p);
        return makeSolution(r1, p - r1, true);
    }

    while (modpow(z, (p - 1) / 2, p) != p - 1) {
        z++;
    }

    c = modpow(z, q, p);
    r = modpow(n, (q + 1) / 2, p);
    t = modpow(n, q, p);
    m = ss;

    while (true) {
        uint64_t i = 0, zz = t;
        uint64_t b = c, e;
        if (t == 1) {
            return makeSolution(r, p - r, true);
        }
        while (zz != 1 && i < (m - 1)) {
            zz = zz * zz % p;
            i++;
        }
        e = m - i - 1;
        while (e > 0) {
            b = b * b % p;
            e--;
        }
        r = r * b % p;
        c = b * b % p;
        t = t * c % p;
        m = i;
    }
}

void test(uint64_t n, uint64_t p) {
    struct Solution sol = ts(n, p);
    printf("n = %llu\n", n);
    printf("p = %llu\n", p);
    if (sol.exists) {
        printf("root1 = %llu\n", sol.root1);
        printf("root2 = %llu\n", sol.root2);
    } else {
        printf("No solution exists\n");
    }
    printf("\n");
}

int main() {
    test(10, 13);
    test(56, 101);
    test(1030, 10009);
    test(1032, 10009);
    test(44402, 100049);

    return 0;
}

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

Version 2

// return (a * b) % mod, avoiding overflow errors while doing modular multiplication.
static unsigned multiplication_modulo(unsigned a, unsigned b, const unsigned mod) {
	unsigned res = 0, tmp;
	for (b %= mod; a; a & 1 ? b >= mod - res ? res -= mod : 0, res += b : 0, a >>= 1, (tmp = b) >= mod - b ? tmp -= mod : 0, b += tmp);
	return res % mod;
}

// return (n ^ exp) % mod
static unsigned mod_pow(unsigned n, unsigned exp, const unsigned mod) {
	unsigned res = 1;
	for (n %= mod; exp; exp & 1 ? res = multiplication_modulo(res, n, mod) : 0, n = multiplication_modulo(n, n, mod), exp >>= 1);
	return res;
}

static unsigned tonelli_shanks_1(const unsigned n, const unsigned mod) {
	// return root such that (root * root) % mod congruent to n % mod.
	// return 0 if no solution to the congruence exists.
	// mod is assumed odd prime.
	const unsigned a = n % mod;
	unsigned res, b, c, d, e, f, g, h;
	if (mod_pow(a, (mod - 1) >> 1, mod) != 1)
		res = 0;
	else
		switch (mod & 7) {
			case 3 : case 7 :
				res = mod_pow(a, (mod + 1) >> 2, mod);
				break;
			case 5 :
				res = mod_pow(a, (mod + 3) >> 3, mod);
				if (multiplication_modulo(res, res, mod) != a){
					b = mod_pow(2, (mod - 1) >> 2, mod);
					res = multiplication_modulo(res, b, mod);
				}
				break;
			default :
				if (a == 1)
					res = 1;
				else {
					for (c = mod - 1, d = 2; d < mod && mod_pow(d, c >> 1, mod) != c; ++d);
					for (e = 0; !(c & 1); ++e, c >>= 1);
					f = mod_pow(a, c, mod);
					b = mod_pow(d, c, mod);
					for (h = 0, g = 0; h < e; h++) {
						d = mod_pow(b, g, mod);
						d = multiplication_modulo(d, f, mod);
						d = mod_pow(d, 1 << (e - 1 - h), mod);
						if (d == mod - 1)
							g += 1 << h;
					}
					f = mod_pow(a, (c + 1) >> 1, mod);
					b = mod_pow(b, g >> 1, mod);
					res = multiplication_modulo(f, b, mod);
				}
		}
	return res;
}

// return root such that (root * root) % mod congruent to n % mod.
// return 0 (the default value of a) if no solution to the congruence exists.
static unsigned tonelli_shanks_2(unsigned n, const unsigned mod) {
	unsigned a = 0, b = mod - 1, c, d = b, e = 0, f = 2, g;
	if (mod_pow(n, b >> 1, mod) == 1) {
		for (; !(d & 1); ++e, d >>= 1);
		if (e == 1)
			a = mod_pow(n, (mod + 1) >> 2, mod);
		else {
			for (; b != mod_pow(f, b >> 1, mod); ++f);
			for (b = mod_pow(f, d, mod), a = mod_pow(n, (d + 1) >> 1, mod), c = mod_pow(n, d, mod), g = e; c != 1; g = d) {
				for (d = 0, e = c, --g; e != 1 && d < g; ++d)
					e = multiplication_modulo(e, e, mod);
				for (f = b, n = g - d; n--;)
					f = multiplication_modulo(f, f, mod);
				a = multiplication_modulo(a, f, mod);
				b = multiplication_modulo(f, f, mod);
				c = multiplication_modulo(c, b, mod);
			}
		}
	}
	return a;
}

#include <assert.h>
int main() {
	unsigned n, mod, root ; /* root_2 = mod - root */

	n = 27875, mod = 26371, root = tonelli_shanks_1(n, mod);
	assert(root == 14320); // 14320 * 14320  mod  26371 = 1504     and   1504 =    27875 mod 26371

	n = 1111111111, mod = 1111111121, root = tonelli_shanks_1(n, mod);
	assert(root == 88664850);

	n = 5258, mod = 3851, root = tonelli_shanks_1(n, mod);
	assert(root == 0); // no solution to the congruence exists.
}

A is assumed odd prime, the algorithm requires O(log A + r * r) multiplications modulo A, where r is the power of 2 dividing A − 1.

C#

Translation of: Java

using System;
using System.Collections.Generic;
using System.Numerics;

namespace TonelliShanks {
    class Solution {
        private readonly BigInteger root1, root2;
        private readonly bool exists;

        public Solution(BigInteger root1, BigInteger root2, bool exists) {
            this.root1 = root1;
            this.root2 = root2;
            this.exists = exists;
        }

        public BigInteger Root1() {
            return root1;
        }

        public BigInteger Root2() {
            return root2;
        }

        public bool Exists() {
            return exists;
        }
    }

    class Program {
        static Solution Ts(BigInteger n, BigInteger p) {
            if (BigInteger.ModPow(n, (p - 1) / 2, p) != 1) {
                return new Solution(0, 0, false);
            }

            BigInteger q = p - 1;
            BigInteger ss = 0;
            while ((q & 1) == 0) {
                ss = ss + 1;
                q = q >> 1;
            }

            if (ss == 1) {
                BigInteger r1 = BigInteger.ModPow(n, (p + 1) / 4, p);
                return new Solution(r1, p - r1, true);
            }

            BigInteger z = 2;
            while (BigInteger.ModPow(z, (p - 1) / 2, p) != p - 1) {
                z = z + 1;
            }
            BigInteger c = BigInteger.ModPow(z, q, p);
            BigInteger r = BigInteger.ModPow(n, (q + 1) / 2, p);
            BigInteger t = BigInteger.ModPow(n, q, p);
            BigInteger m = ss;

            while (true) {
                if (t == 1) {
                    return new Solution(r, p - r, true);
                }
                BigInteger i = 0;
                BigInteger zz = t;
                while (zz != 1 && i < (m - 1)) {
                    zz = zz * zz % p;
                    i = i + 1;
                }
                BigInteger b = c;
                BigInteger e = m - i - 1;
                while (e > 0) {
                    b = b * b % p;
                    e = e - 1;
                }
                r = r * b % p;
                c = b * b % p;
                t = t * c % p;
                m = i;
            }
        }

        static void Main(string[] args) {
            List<Tuple<long, long>> pairs = new List<Tuple<long, long>>() {
                new Tuple<long, long>(10, 13),
                new Tuple<long, long>(56, 101),
                new Tuple<long, long>(1030, 10009),
                new Tuple<long, long>(1032, 10009),
                new Tuple<long, long>(44402, 100049),
                new Tuple<long, long>(665820697, 1000000009),
                new Tuple<long, long>(881398088036, 1000000000039),
            };

            foreach (var pair in pairs) {
                Solution sol = Ts(pair.Item1, pair.Item2);
                Console.WriteLine("n = {0}", pair.Item1);
                Console.WriteLine("p = {0}", pair.Item2);
                if (sol.Exists()) {
                    Console.WriteLine("root1 = {0}", sol.Root1());
                    Console.WriteLine("root2 = {0}", sol.Root2());
                } else {
                    Console.WriteLine("No solution exists");
                }
                Console.WriteLine();
            }

            BigInteger bn = BigInteger.Parse("41660815127637347468140745042827704103445750172002");
            BigInteger bp = BigInteger.Pow(10, 50) + 577;
            Solution bsol = Ts(bn, bp);
            Console.WriteLine("n = {0}", bn);
            Console.WriteLine("p = {0}", bp);
            if (bsol.Exists()) {
                Console.WriteLine("root1 = {0}", bsol.Root1());
                Console.WriteLine("root2 = {0}", bsol.Root2());
            } else {
                Console.WriteLine("No solution exists");
            }
        }
    }
}

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

C++

#include <cstdint>
#include <iostream>
#include <vector>

struct Pair {
	uint64_t n;
	uint64_t p;
};

struct Solution {
	uint64_t root1;
	uint64_t root2;
	bool is_square;
};

uint64_t multiply_modulus(uint64_t a, uint64_t b, const uint64_t& modulus) {
    a %= modulus; b %= modulus;
    if ( b < a ) {
    	uint64_t temp = a; a = b; b = temp;
    }

    uint64_t result = 0;
    while ( a > 0 ) {
        if ( a % 2 == 1 ) {
        	result = ( result + b ) % modulus;
        };
        b = ( b << 1 ) % modulus;
        a >>= 1;
    }
    return result;
}

uint64_t power_modulus(uint64_t base, uint64_t exponent, const uint64_t& modulus) {
	if ( modulus == 1 ) {
		return 0;
	}

	base %= modulus;
	uint64_t result = 1;
	while ( exponent > 0 ) {
		if ( ( exponent & 1 ) == 1 ) {
			result = multiply_modulus(result, base, modulus);
		}
		base = multiply_modulus(base, base, modulus);
		exponent >>= 1;
	}
	return result;
}

uint64_t legendre(const uint64_t& a, const uint64_t& p) {
    return power_modulus(a, ( p - 1 ) / 2, p);
}

Solution tonelli_shanks(const uint64_t& n, const uint64_t& p) {
	if ( legendre(n, p) != 1 ) {
		return Solution(0, 0, false);
	}

	// Factor out powers of 2 from p - 1
    uint64_t q = p - 1;
    uint64_t s = 0;
    while ( q % 2 == 0 ) {
        q /= 2;
        s += 1;
    }

    if ( s == 1 ) {
	    uint64_t result = power_modulus(n, ( p + 1 ) / 4, p);
	    return Solution(result, p - result, true);
    }

    // Find a non-square z such as ( z | p ) = -1
    uint64_t z = 2;
	while ( legendre(z, p) != p - 1 ) {
		z += 1;
    }

    uint64_t c = power_modulus(z, q, p);
    uint64_t t = power_modulus(n, q, p);
    uint64_t m = s;
    uint64_t result = power_modulus(n, ( q + 1 ) >> 1, p);

    while ( t != 1 ) {
        uint64_t i = 1;
        z = multiply_modulus(t, t, p);
        while ( z != 1 && i < m - 1 ) {
            i += 1;
            z = multiply_modulus(z, z, p);
        }
        uint64_t b = power_modulus(c, 1 << ( m - i - 1 ), p);
        c = multiply_modulus(b, b, p);
        t = multiply_modulus(t, c, p);
        m = i;
        result = multiply_modulus(result, b, p);
    }
    return Solution(result, p - result, true);
}

int main() {
	const std::vector<Pair> tests = { Pair(10, 13), Pair(56, 101), Pair(1030, 1009), Pair(1032, 1009),
		Pair(44402, 100049), Pair(665820697, 1000000009), Pair(881398088036, 1000000000039) };

	for ( const Pair& test : tests ) {
		Solution solution = tonelli_shanks(test.n, test.p);
		std::cout << "n = " << test.n << ", p = " << test.p;
		if ( solution.is_square == true ) {
			std::cout << " has solutions: " << solution.root1 << " and " << solution.root2 << std::endl << std::endl;
		} else {
			std::cout << " has no solutions because n is not a square modulo p" << std::endl << std::endl;
		}
	}
}

Output:

n = 10, p = 13 has solutions: 7 and 6

n = 56, p = 101 has solutions: 37 and 64

n = 1030, p = 1009 has solutions: 651 and 358

n = 1032, p = 1009 has no solutions because n is not a square modulo p

n = 44402, p = 100049 has solutions: 30468 and 69581

n = 665820697, p = 1000000009 has solutions: 378633312 and 621366697

n = 881398088036, p = 1000000000039 has solutions: 791399408049 and 208600591990

Clojure

(defn find-first
 " Finds first element of collection that satisifies predicate function pred "
  [pred coll]
  (first (filter pred coll)))

(defn modpow
  " b^e mod m (using Java which solves some cases the pure clojure method has to be modified to tackle--i.e. with large b & e and
    calculation simplications when gcd(b, m) == 1 and gcd(e, m) == 1) "
  [b e m]
  (.modPow (biginteger b) (biginteger e) (biginteger m)))

(defn legendre [a p]
  (modpow a (quot (dec p) 2) p)
)

(defn tonelli [n p]
  " Following Wikipedia https://en.wikipedia.org/wiki/Tonelli%E2%80%93Shanks_algorithm "
  (assert (= (legendre n p) 1) "not a square (mod p)")
  (loop [q (dec p)                                                  ; Step 1 in Wikipedia
         s 0]
    (if (zero? (rem q 2))
      (recur (quot q 2) (inc s))
      (if (= s 1)
        (modpow n (quot (inc p) 4) p)
        (let [z (find-first #(= (dec p) (legendre % p)) (range 2 p))] ; Step 2 in Wikipedia
          (loop [
                 M s
                 c (modpow z q p)
                 t (modpow n q p)
                 R (modpow n (quot (inc q) 2) p)]
            (if (= t 1)
              R
              (let [i (long (find-first #(= 1 (modpow t (bit-shift-left 1 %) p)) (range 1 M))) ; Step 3
                    b (modpow c (bit-shift-left 1 (- M i 1)) p)
                    M i
                    c (modpow b 2 p)
                    t (rem (* t c) p)
                    R (rem (* R b) p)]
                (recur M c t R)
                )
              )
            )
          )
        )
      )
    )
  )


; Testing--using Python examples
(doseq [[n p]  [[10, 13], [56, 101], [1030, 10009], [44402, 100049],
                [665820697, 1000000009], [881398088036, 1000000000039],
                [41660815127637347468140745042827704103445750172002, 100000000000000000000000000000000000000000000000577]]
        :let [r (tonelli n p)]]
  (println (format "n: %5d p: %d \n\troots: %5d %5d" (biginteger n) (biginteger p) (biginteger r) (biginteger (- p r)))))

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

D

Translation of: Kotlin

import std.bigint;
import std.stdio;
import std.typecons;

alias Pair = Tuple!(long, "n", long, "p");

enum BIGZERO = BigInt("0");
enum BIGONE = BigInt("1");
enum BIGTWO = BigInt("2");
enum BIGTEN = BigInt("10");

struct Solution {
    BigInt root1, root2;
    bool exists;
}

/// https://en.wikipedia.org/wiki/Modular_exponentiation#Right-to-left_binary_method
BigInt modPow(BigInt b, BigInt e, BigInt n) {
    if (n == 1) return BIGZERO;
    BigInt result = 1;
    b = b % n;
    while (e > 0) {
        if (e % 2 == 1) {
            result = (result * b) % n;
        }
        e >>= 1;
        b = (b*b) % n;
    }
    return result;
}

Solution ts(long n, long p) {
    return ts(BigInt(n), BigInt(p));
}

Solution ts(BigInt n, BigInt p) {
    auto powMod(BigInt a, BigInt e) {
        return a.modPow(e, p);
    }

    auto ls(BigInt a) {
        return powMod(a, (p-1)/2);
    }

    if (ls(n) != 1) return Solution(BIGZERO, BIGZERO, false);
    auto q = p - 1;
    auto ss = BIGZERO;
    while ((q & 1) == 0) {
        ss = ss + 1;
        q = q >> 1;
    }

    if (ss == BIGONE) {
        auto r1 = powMod(n, (p + 1) / 4);
        return Solution(r1, p - r1, true);
    }

    auto z = BIGTWO;
    while (ls(z) != p - 1) z = z + 1;
    auto c = powMod(z, q);
    auto r = powMod(n, (q + 1) / 2);
    auto t = powMod(n, q);
    auto m = ss;

    while (true) {
        if (t == 1) return Solution(r, p - r, true);
        auto i = BIGZERO;
        auto zz = t;
        while (zz != 1 && i < m - 1) {
            zz  = zz * zz % p;
            i = i + 1;
        }
        auto b = c;
        auto e = m - i - 1;
        while (e > 0) {
            b = b * b % p;
            e = e - 1;
        }
        r = r * b % p;
        c = b * b % p;
        t = t * c % p;
        m = i;
    }
}

void main() {
    auto pairs = [
        Pair(             10L,                13L),
        Pair(             56L,               101L),
        Pair(          1_030L,            10_009L),
        Pair(          1_032L,            10_009L),
        Pair(         44_402L,           100_049L),
        Pair(    665_820_697L,     1_000_000_009L),
        Pair(881_398_088_036L, 1_000_000_000_039L),
    ];

    foreach (pair; pairs) {
        auto sol = ts(pair.n, pair.p);

        writeln("n = ", pair.n);
        writeln("p = ", pair.p);
        if (sol.exists) {
            writeln("root1 = ", sol.root1);
            writeln("root2 = ", sol.root2);
        }
        else writeln("No solution exists");
        writeln();
    }

    auto bn = BigInt("41660815127637347468140745042827704103445750172002");
    auto bp = BIGTEN ^^ 50 + 577L;
    auto sol = ts(bn, bp);
    writeln("n = ", bn);
    writeln("p = ", bp);
    if (sol.exists) {
        writeln("root1 = ", sol.root1);
        writeln("root2 = ", sol.root2);
    }
    else writeln("No solution exists");
}

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

EchoLisp

(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)))))

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

' 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

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

' version 12-04-2017
' compile with: fbc -s console

#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

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.

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))
}

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.

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)
}

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.

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)
}

Haskell

Translation of: Python

import Data.List (genericTake, genericLength)
import Data.Bits (shiftR)

powMod :: Integer -> Integer -> Integer -> Integer
powMod m b e = go b e 1
  where
    go b e r
      | e == 0 = r
      | odd e  = go ((b*b) `mod` m) (e `div` 2) ((r*b) `mod` m)
      | even e = go ((b*b) `mod` m) (e `div` 2) r 

legendre :: Integer -> Integer -> Integer
legendre a p = powMod p a ((p - 1) `div` 2) 

tonelli :: Integer -> Integer -> Maybe (Integer, Integer)
tonelli n p | legendre n p /= 1 = Nothing
tonelli n p =
  let s = length $ takeWhile even $ iterate (`div` 2) (p-1)
      q = shiftR (p-1) s
  in if s == 1
    then let r = powMod p n ((p+1) `div` 4)
         in Just (r, p - r)
    else let z = (2 +) . genericLength
               $ takeWhile (\i -> p - 1 /= legendre i p)
               $ [2..p-1]
         in loop s
            ( powMod p z q )
            ( powMod p n $ (q+1) `div` 2 )
            ( powMod p n q )
  where
    loop m c r t
      | (t - 1) `mod` p == 0 = Just (r, p - r)
      | otherwise =
        let i = (1 +) . genericLength . genericTake (m - 2)
              $ takeWhile (\t2 -> (t2 - 1) `mod` p /= 0)
              $ iterate (\t2 -> (t2*t2) `mod` p)
              $ (t*t) `mod` p
            b = powMod p c (2^(m - i - 1))
            r' = (r*b)  `mod` p
            c' = (b*b)  `mod` p
            t' = (t*c') `mod` p
        in loop i c' r' t'

λ> tonelli 10 13
Just (7,6)
λ> tonelli 56 101
Just (37,64)
λ> tonelli 1030 10009
Just (1632,8377)
λ> tonelli 1032 10009
Nothing
λ> tonelli 44402 100049
Just (30468,69581)
λ> tonelli 665820697 1000000009
Just (378633312,621366697)
λ> tonelli 881398088036 1000000000039
Just (791399408049,208600591990)
λ> tonelli 41660815127637347468140745042827704103445750172002 $ (10^50)+577
Just (32102985369940620849741983987300038903725266634508,67897014630059379150258016012699961096274733366069)

J

Implementation:

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
)

Task examples:

   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

Java

Translation of: Kotlin

Works with: Java version 9

import java.math.BigInteger;
import java.util.List;
import java.util.Map;
import java.util.function.BiFunction;
import java.util.function.Function;

public class TonelliShanks {
    private static final BigInteger ZERO = BigInteger.ZERO;
    private static final BigInteger ONE = BigInteger.ONE;
    private static final BigInteger TEN = BigInteger.TEN;
    private static final BigInteger TWO = BigInteger.valueOf(2);
    private static final BigInteger FOUR = BigInteger.valueOf(4);

    private static class Solution {
        private BigInteger root1;
        private BigInteger root2;
        private boolean exists;

        Solution(BigInteger root1, BigInteger root2, boolean exists) {
            this.root1 = root1;
            this.root2 = root2;
            this.exists = exists;
        }
    }

    private static Solution ts(Long n, Long p) {
        return ts(BigInteger.valueOf(n), BigInteger.valueOf(p));
    }

    private static Solution ts(BigInteger n, BigInteger p) {
        BiFunction<BigInteger, BigInteger, BigInteger> powModP = (BigInteger a, BigInteger e) -> a.modPow(e, p);
        Function<BigInteger, BigInteger> ls = (BigInteger a) -> powModP.apply(a, p.subtract(ONE).divide(TWO));

        if (!ls.apply(n).equals(ONE)) return new Solution(ZERO, ZERO, false);

        BigInteger q = p.subtract(ONE);
        BigInteger ss = ZERO;
        while (q.and(ONE).equals(ZERO)) {
            ss = ss.add(ONE);
            q = q.shiftRight(1);
        }

        if (ss.equals(ONE)) {
            BigInteger r1 = powModP.apply(n, p.add(ONE).divide(FOUR));
            return new Solution(r1, p.subtract(r1), true);
        }

        BigInteger z = TWO;
        while (!ls.apply(z).equals(p.subtract(ONE))) z = z.add(ONE);
        BigInteger c = powModP.apply(z, q);
        BigInteger r = powModP.apply(n, q.add(ONE).divide(TWO));
        BigInteger t = powModP.apply(n, q);
        BigInteger m = ss;

        while (true) {
            if (t.equals(ONE)) return new Solution(r, p.subtract(r), true);
            BigInteger i = ZERO;
            BigInteger zz = t;
            while (!zz.equals(BigInteger.ONE) && i.compareTo(m.subtract(ONE)) < 0) {
                zz = zz.multiply(zz).mod(p);
                i = i.add(ONE);
            }
            BigInteger b = c;
            BigInteger e = m.subtract(i).subtract(ONE);
            while (e.compareTo(ZERO) > 0) {
                b = b.multiply(b).mod(p);
                e = e.subtract(ONE);
            }
            r = r.multiply(b).mod(p);
            c = b.multiply(b).mod(p);
            t = t.multiply(c).mod(p);
            m = i;
        }
    }

    public static void main(String[] args) {
        List<Map.Entry<Long, Long>> pairs = List.of(
            Map.entry(10L, 13L),
            Map.entry(56L, 101L),
            Map.entry(1030L, 10009L),
            Map.entry(1032L, 10009L),
            Map.entry(44402L, 100049L),
            Map.entry(665820697L, 1000000009L),
            Map.entry(881398088036L, 1000000000039L)
        );

        for (Map.Entry<Long, Long> pair : pairs) {
            Solution sol = ts(pair.getKey(), pair.getValue());
            System.out.printf("n = %s\n", pair.getKey());
            System.out.printf("p = %s\n", pair.getValue());
            if (sol.exists) {
                System.out.printf("root1 = %s\n", sol.root1);
                System.out.printf("root2 = %s\n", sol.root2);
            } else {
                System.out.println("No solution exists");
            }
            System.out.println();
        }

        BigInteger bn = new BigInteger("41660815127637347468140745042827704103445750172002");
        BigInteger bp = TEN.pow(50).add(BigInteger.valueOf(577));
        Solution sol = ts(bn, bp);
        System.out.printf("n = %s\n", bn);
        System.out.printf("p = %s\n", bp);
        if (sol.exists) {
            System.out.printf("root1 = %s\n", sol.root1);
            System.out.printf("root2 = %s\n", sol.root2);
        } else {
            System.out.println("No solution exists");
        }
    }
}

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

jq

Works with gojq and fq, two Go implementations of jq

Adapted from Wren

The Go implementations of jq provide indefinite-precision integer arithmetic.

See Modular_exponentiation for suitable jq definitions of `power/1` and `modPow/2` as used here.

include "rc-modular-exponentiation";  # see remark above

# If $j is 0, then an error condition is raised;
# otherwise, assuming infinite-precision integer arithmetic,
# if the input and $j are integers, then the result will be an integer.
def idivide($j):
  . as $i
  | ($i % $j) as $mod
  | ($i - $mod) / $j ;

def Solution(a;b;c):
  {"root1": a, "root2": b, "exists": c};

# pretty print a Solution
def pp:
  if .exists
  then "root1 = \(.root1)",
       "root2 = \(.root2)"
  else "No solution exists"
  end;

# Tonelli-Shanks
def ts($n; $p):
    def powModP($a; $e): $a | modPow($e; $p);

    def ls($a): powModP($a; ($p - 1) | idivide(2));

    if ls($n) != 1 then Solution(0; 0; false)
    else { q: ($p - 1), ss: 0}
    | until (.q % 2 != 0;
        .ss += 1
        | .q |= idivide(2) )
    | if .ss == 1
      then powModP(n; ($p+1) | idivide(4)) as $r1
      | Solution($r1; $p - $r1; true)
      else .z = 2
      | until ( ls(.z) == ($p - 1); .z += 1 )
      | .c = powModP(.z; .q)
      | .r = powModP($n; (.q+1) | idivide(2))
      | .t = powModP($n; .q)
      | .m = .ss
      | until (.emit;
          if .t == 1 then .emit = Solution(.r; $p - .r; true)
          else .i = 0
          | .zz = .t
          | until (.zz == 1 or .i >= (.m - 1);
              .zz = (.zz * .zz) % p
              | .i += 1 )
          | .b = .c
          | .e = .m - (1 + .i)
          | until (.e <= 0;
              .b = (.b * .b) % $p
              | .e += -1 )
          | .r = (.r * .b) % $p
          | .c = (.b * .b) % $p
          | .t = (.t * .c) % $p
          | .m = .i
	  end )
      | .emit
      end
    end;
    
def pairs: [
    [10, 13], [56, 101], [1030, 10009], [1032, 10009], [44402, 100049],
    [665820697, 1000000009], [881398088036, 1000000000039]
];

def task:
  pairs[] as [$n, $p]
  | ts($n; $p) as $sol
  | "n     = \($n)",
    "p     = \($p)",
    ($sol | pp),
    "";

def task2:
  def bn: 41660815127637347468140745042827704103445750172002;
  def bp: (10 | power(50)) + 577;
  ts(bn; bp) as $bsol
  | "n     = \(bn)",
    "p     = \(bp)",
    ( $bsol | pp );

task, task2

Output:

See Wren.

Julia

Works with: Julia version 0.6

Module:

module TonelliShanks

legendre(a, p) = powermod(a, (p - 1) ÷ 2, p)

function solve(n::T, p::T) where T <: Union{Int, Int128, BigInt}
    legendre(n, p) != 1 && throw(ArgumentError("$n not a square (mod $p)"))
    local q::T = p - one(p)
    local s::T = 0
    while iszero(q % 2)
        q ÷= 2
        s += one(s)
    end
    if s == one(s)
        r = powermod(n, (p + 1) >> 2, p)
        return r, p - r
    end
    local z::T
    for z in 2:(p - 1)
        p - 1 == legendre(z, p) && break
    end
    local c::T = powermod(z, q, p)
    local r::T = powermod(n, (q + 1) >> 1, p)
    local t::T = powermod(n, q, p)
    local m::T = s
    local t2::T = zero(p)
    while !iszero((t - 1) % p)
        t2 = (t * t) % p
        local i::T
        for i in Base.OneTo(m)
            iszero((t2 - 1) % p) && break
            t2 = (t2 * t2) % p
        end
        b = powermod(c, 1 << (m - i - 1), p)
        r = (r * b) % p
        c = (b * b) % p
        t = (t * c) % p
        m = i
    end
    return r, p - r
end

end  # module TonelliShanks

Main:

@show TonelliShanks.solve(10, 13)
@show TonelliShanks.solve(56, 101)
@show TonelliShanks.solve(1030, 10009)
@show TonelliShanks.solve(44402, 100049)
@show TonelliShanks.solve(665820697, 1000000009)
@show TonelliShanks.solve(881398088036, 1000000000039)
@show TonelliShanks.solve(41660815127637347468140745042827704103445750172002, big"10" ^ 50 + 577)

Output:

TonelliShanks.solve(10, 13) = (7, 6)
TonelliShanks.solve(56, 101) = (37, 64)
TonelliShanks.solve(1030, 10009) = (1632, 8377)
TonelliShanks.solve(44402, 100049) = (30468, 69581)
TonelliShanks.solve(665820697, 1000000009) = (378633312, 621366697)
TonelliShanks.solve(881398088036, 1000000000039) = (791399408049, 208600591990)
TonelliShanks.solve(@big_str("41660815127637347468140745042827704103445750172002"), @big_str("10") ^ 50 + 577) = (32102985369940620849741983987300038903725266634508, 67897014630059379150258016012699961096274733366069)

Kotlin

Translation of: Go

// 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")    
}

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

Nim

Based algorithm pseudo-code, referencing python 3.

proc pow*[T: SomeInteger](x, n, p: T): T = 
  var t = x mod p 
  var e = n 
  result = 1 
  while e > 0: 
    if (e and 1) == 1: 
      result = result * t mod p 
    t = t * t mod p 
    e = e shr 1 
 
proc legendre*[T: SomeInteger](a, p: T): T = pow(a, (p-1) shr 1, p) 
 
proc tonelliShanks*[T: SomeInteger](n, p: T): T =
  # Check that n is indeed a square.
  if legendre(n, p) != 1:
    raise newException(ValueError, "Not a square")
 
  # Factor out power of 2 from p-1.
  var q = p - 1
  var s = 0
  while (q and 1) == 0:
    s += 1
    q = q shr 1 
 
  if s == 1: 
    return pow(n, (p+1) shr 2, p)
 
  # Select a non-square z such as (z | p) = -1.
  var z = 2 
  while legendre(z, p) != p - 1: 
    z += 1
 
  var 
    c = pow(z, q, p)
    t = pow(n, q, p)
    m = s
  result = pow(n, (q+1) shr 1, p)
  while t != 1:
    var 
      i = 1
      z = t * t mod p 
    while z != 1 and i < m-1:
      i += 1
      z = z * z mod p 
 
    var b = pow(c, 1 shl (m-i-1), p)
    c = b * b mod p 
    t = t * c mod p 
    m = i 
    result = result * b mod p 
 
when isMainModule: 
  proc run(n, p: SomeInteger) = 
    try: 
      let r = tonelliShanks(n, p)
      echo r, " ", p-r
    except ValueError:
      echo getCurrentExceptionMsg()
 
  run(10, 13)
  run(56, 101)
  run(1030, 10009)
  run(1032, 10009)
  run(44402, 100049) 
  run(665820697, 1000000009)

output:

7 6
37 64
1632 8377
Not a square
30468 69581
378633312 621366697

OCaml

Translation of: Java

Library: zarith

An extra test case has been added for the `s = 1` branch.

let tonelli n p =
  let open Z in
  let two = ~$2 in
  let pp = pred p in
  let pph = pred p / two in
  let pow_mod_p a e = powm a e p in
  let legendre_p a = pow_mod_p a pph in

  if legendre_p n <> one then None
  else
    let s = trailing_zeros pp in
    if s = 1 then
      let r = pow_mod_p n (succ p / ~$4) in
      Some (r, p - r)
    else
      let q = pp asr s in
      let z =
        let rec find_non_square z =
          if legendre_p z = pp then z else find_non_square (succ z)
        in
        find_non_square two
      in
      let rec loop c r t m =
        if t = one then (r, p - r)
        else
          let mp = pred m in
          let rec find_i n i =
            if n = one || i >= mp then i else find_i (n * n mod p) (succ i)
          in
          let rec exp_pow2 b e =
            if e <= zero then b else exp_pow2 (b * b mod p) (pred e)
          in
          let i = find_i t zero in
          let b = exp_pow2 c (mp - i) in
          let c = b * b mod p in
          loop c (r * b mod p) (t * c mod p) i
      in
      Some
        (loop (pow_mod_p z q) (pow_mod_p n (succ q / two)) (pow_mod_p n q) ~$s)

let () =
  let open Z in
  [
    (~$9, ~$11);
    (~$10, ~$13);
    (~$56, ~$101);
    (~$1030, ~$10009);
    (~$1032, ~$10009);
    (~$44402, ~$100049);
    (~$665820697, ~$1000000009);
    (~$881398088036, ~$1000000000039);
    ( of_string "41660815127637347468140745042827704103445750172002",
      pow ~$10 50 + ~$577 );
  ]
  |> List.iter (fun (n, p) ->
         Printf.printf "n = %s\np = %s\n%!" (to_string n) (to_string p);
         match tonelli n p with
         | Some (r1, r2) ->
             Printf.printf "root1 = %s\nroot2 = %s\n\n%!" (to_string r1)
               (to_string r2)
         | None -> print_endline "No solution exists\n")

Output:

n = 9                
p = 11
root1 = 3
root2 = 8

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

Translation of: Raku

Library: ntheory

use bigint;
use ntheory qw(is_prime powmod kronecker);

sub tonelli_shanks {
    my($n,$p) = @_;
    return if kronecker($n,$p) <= 0;
    my $Q = $p - 1;
    my $S = 0;
    $Q >>= 1 and $S++ while 0 == $Q%2;
    return powmod($n,int(($p+1)/4), $p) if $S == 1;

    my $c;
    for $n (2..$p) {
        next if kronecker($n,$p) >= 0;
        $c = powmod($n, $Q, $p);
        last;
    }

    my $R = powmod($n, ($Q+1) >> 1, $p ); # ?
    my $t = powmod($n, $Q, $p );
    while (($t-1) % $p) {
        my $b;
        my $t2 = $t**2 % $p;
        for (1 .. $S) {
            if (0 == ($t2-1)%$p) {
                $b = powmod($c, 1 << ($S-1-$_), $p);
                $S = $_;
                last;
            }
            $t2 = $t2**2 % $p;
        }
        $R = ($R * $b) % $p;
        $c = $b**2 % $p;
        $t = ($t * $c) % $p;
    }
    $R;
}

my @tests = (
    (10, 13),
    (56, 101),
    (1030, 10009),
    (1032, 10009),
    (44402, 100049),
    (665820697, 1000000009),
    (881398088036, 1000000000039),
);

while (@tests) {
    $n = shift @tests;
    $p = shift @tests;
    my $t = tonelli_shanks($n, $p);
    if (!$t or ($t**2 - $n) % $p) {
        printf "No solution for (%d, %d)\n", $n, $p;
    } else {
        printf "Roots of %d are (%d, %d) mod %d\n", $n, $t, $p-$t, $p;
    }
}

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

Phix

Translation of: C#

Library: Phix/mpfr

with javascript_semantics 
include mpfr.e
 
function ts(string ns, ps)
    mpz n = mpz_init(ns),
        p = mpz_init(ps),
        t = mpz_init(),
        r = mpz_init(),
        pm1 = mpz_init(),
        pm2 = mpz_init()
    mpz_sub_ui(pm1,p,1)                 -- pm1 = p-1
    mpz_fdiv_q_2exp(pm2,pm1,1)          -- pm2 = pm1/2
    mpz_powm(t,n,pm2,p)                 -- t = mod(n^pm2,p)
    if mpz_cmp_si(t,1)!=0 then
        return "No solution exists"
    end if
    mpz q = mpz_init_set(pm1)
    integer ss = 0
    while mpz_even(q) do
        ss += 1
        mpz_fdiv_q_2exp(q,q,1)          -- q/=2
    end while
    if ss=1 then
        mpz_add_ui(t,p,1)
        mpz_fdiv_q_2exp(t,t,2)
        mpz_powm(r,n,t,p)               -- r = mod(n^((p+1)/4),p)
    else
        mpz z = mpz_init(2)
        while true do
            mpz_powm(t,z,pm2,p)         -- t = mod(z^pm2,p)
            if mpz_cmp(t,pm1)=0 then exit end if
            mpz_add_ui(z,z,1)           -- z+= 1
        end while
        mpz {b,c,zz} = mpz_inits(3)
        mpz_powm(c,z,q,p)               -- c = mod(z^q,p)
        mpz_add_ui(t,q,1)
        mpz_fdiv_q_2exp(t,t,1)
        mpz_powm(r,n,t,p)               -- r = mod(n^((q+1)/2),p)
        mpz_powm(t,n,q,p)               -- t = mod(n^q,p)
        integer m = ss
        while mpz_cmp_si(t,1) do        -- t!=1
            integer i = 0
            mpz_set(zz,t)
            while mpz_cmp_si(zz,1)!=0 and i<m-1 do
                mpz_powm_ui(zz,zz,2,p)  -- zz = mod(zz^2,p)
                i += 1
            end while
            mpz_set(b,c)
            integer e = m-i-1
            while e>0 do
                mpz_powm_ui(b,b,2,p)    -- b = mod(b^2,p)
                e -= 1
            end while
            mpz_mul(r,r,b)
            mpz_mod(r,r,p)              -- r = mod(r*b,p)
            mpz_powm_ui(c,b,2,p)        -- c = mod(b^2,p)
            mpz_mul(t,t,c)
            mpz_mod(t,t,p)              -- t = mod(t*c,p)
            m = i
        end while
    end if
    mpz_sub(p,p,r)
    return mpz_get_str(r)&" and "&mpz_get_str(p)
end function
 
constant tests = {{"10","13"},
                  {"56","101"},
                  {"1030","10009"},
                  {"1032","10009"},
                  {"44402","100049"},
                  {"665820697","1000000009"},
                  {"881398088036","1000000000039"},
                  {"41660815127637347468140745042827704103445750172002",
                   sprintf("1%s577",repeat('0',47))}} -- 10^50+577
 
for i=1 to length(tests) do
    string {p1,p2} = tests[i]   
    printf(1,"For n = %s and p = %s, %s\n",{p1,p2,ts(p1,p2)})
end for

Output:

For n = 10 and p = 13, 7 and 6
For n = 56 and p = 101, 37 and 64
For n = 1030 and p = 10009, 1632 and 8377
For n = 1032 and p = 10009, No solution exists
For n = 44402 and p = 100049, 30468 and 69581
For n = 665820697 and p = 1000000009, 378633312 and 621366697
For n = 881398088036 and p = 1000000000039, 791399408049 and 208600591990
For n = 41660815127637347468140745042827704103445750172002 and p = 100000000000000000000000000000000000000000000000577, 
        32102985369940620849741983987300038903725266634508 and 67897014630059379150258016012699961096274733366069

PicoLisp

Translation of: Go

# 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)))

Output:

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

Powershell

Translation of: Python

Works with: Powershell version 7

Function Invoke-ModuloExponentiation ([BigInt]$Base, [BigInt]$Exponent, $Modulo) {
    $Result = 1
    $Base = $Base % $Modulo
    If ($Base -eq 0) {return 0}
    
    While ($Exponent -gt 0) {
        If (($Exponent -band 1) -eq 1) {$Result = ($Result * $Base) % $Modulo}
        $Exponent = $Exponent -shr 1
        $Base = ($Base * $Base) % $Modulo
    }
    return ($Result % $Modulo)
}

Function Get-Legendre ([BigInt]$Integer, [BigInt]$Prime) {
    return (Invoke-ModuloExponentiation -Base $Integer -Exponent (($Prime - 1) / 2) -Modulo $Prime)
}

Function Invoke-TonelliShanks ([BigInt]$Integer, [BigInt]$Prime) {
    If ((Get-Legendre $Integer $Prime) -ne 1) {throw "$Integer not a square (mod $Prime)"}
    [bigint]$q = $Prime - 1
    $s = 0
    While (($q % 2) -eq 0) {
        $q = $q / 2
        $s++
    }
    If ($s -eq 1) {
        return (Invoke-ModuloExponentiation $Integer -Exponent (($Prime + 1) / 4) -Modulo $Prime)
    }
    For ($z = 2; [Bigint]::Compare($z, $Prime) -lt 0; $z++) {
        If ([BigInt]::Compare(($Prime - 1), (Get-Legendre $z $Prime)) -eq 0) {
            break
        }
    }
    $c = Invoke-ModuloExponentiation -Base $z -Exponent $q -Modulo $Prime
    $r = Invoke-ModuloExponentiation -Base $Integer -Exponent (($q + 1) / 2) -Modulo $Prime
    $t = Invoke-ModuloExponentiation -Base $Integer -Exponent $q -Modulo $Prime
    $m = $s
    $t2 = 0
    
    While ((($t - 1) % $Prime) -ne 0) {
        $t2 = $t * $t % $Prime
        Foreach ($i in (1..$m)) {
            If ((($t2 -1) % $Prime) -eq 0) {
                break
            }
            $t2 = Invoke-ModuloExponentiation -Base $t2 -Exponent 2 -Modulo $Prime
        }
        $b = Invoke-ModuloExponentiation -Base $c -Exponent ([Math]::Pow(2, ($m - $i - 1))) -Modulo $Prime
        $r = ($r * $b) % $Prime
        $c = ($b * $b) % $Prime
        $t = ($t * $c) % $Prime
        $m = $i
    }
    return $r
}

$TonelliTests = @(
    @{Integer = [BigInt]::Parse('10'); Prime = [BigInt]::Parse('13')},
    @{Integer = [BigInt]::Parse('56'); Prime = [BigInt]::Parse('101')},
    @{Integer = [BigInt]::Parse('1030'); Prime = [BigInt]::Parse('10009')},
    @{Integer = [BigInt]::Parse('44402'); Prime = [BigInt]::Parse('100049')},
    @{Integer = [BigInt]::Parse('665820697'); Prime = [BigInt]::Parse('1000000009')},
    @{Integer = [BigInt]::Parse('881398088036'); Prime = [BigInt]::Parse('1000000000039')},
    @{Integer = [BigInt]::Parse('41660815127637347468140745042827704103445750172002'); Prime = [BigInt]::Parse('100000000000000000000000000000000000000000000000577')}
)

$TonelliTests | Foreach-Object {
    $Result = Invoke-TonelliShanks @_
    [PSCustomObject]@{
        n = $_['Integer']
        p = $_['Prime']
        Roots = @($Result, ($_['Prime'] - $Result))
    }
} | Format-List

Output:

n     : 41660815127637347468140745042827704103445750172002
p     : 100000000000000000000000000000000000000000000000577
Roots : {32102985369940620849741983987300038903725266634508, 67897014630059379150258016012699961096274733366069}

n     : 41660815127637347468140745042827704103445750172002
p     : 100000000000000000000000000000000000000000000000577
Roots : {32102985369940620849741983987300038903725266634508, 67897014630059379150258016012699961096274733366069}

n     : 41660815127637347468140745042827704103445750172002
p     : 100000000000000000000000000000000000000000000000577
Roots : {32102985369940620849741983987300038903725266634508, 67897014630059379150258016012699961096274733366069}

n     : 41660815127637347468140745042827704103445750172002
p     : 100000000000000000000000000000000000000000000000577
Roots : {32102985369940620849741983987300038903725266634508, 67897014630059379150258016012699961096274733366069}

n     : 41660815127637347468140745042827704103445750172002
p     : 100000000000000000000000000000000000000000000000577
Roots : {32102985369940620849741983987300038903725266634508, 67897014630059379150258016012699961096274733366069}

n     : 41660815127637347468140745042827704103445750172002
p     : 100000000000000000000000000000000000000000000000577
Roots : {32102985369940620849741983987300038903725266634508, 67897014630059379150258016012699961096274733366069}

n     : 41660815127637347468140745042827704103445750172002
p     : 100000000000000000000000000000000000000000000000577
Roots : {32102985369940620849741983987300038903725266634508, 67897014630059379150258016012699961096274733366069}

Python

Translation of: EchoLisp

Works with: Python version 3

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

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

(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))))

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

Raku

(formerly Perl 6)

Works with: Rakudo version 2018.04

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

#  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";
}

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

REXX

Translation of: Python

The large numbers cannot reasonably be handled by the pow function shown here.

/* REXX (required by some interpreters) */
Numeric Digits 1000000
ttest ='[(10, 13), (56, 101), (1030, 10009), (44402, 100049)]'
Do While pos('(',ttest)>0
  Parse Var ttest '(' n ',' p ')' ttest
  r = tonelli(n, p)
  Say "n =" n "p =" p
  Say "          roots :" r (p - r)
  End
Exit

legendre: Procedure
  Parse Arg a, p
  return pow(a, (p - 1) % 2, p)

tonelli: Procedure
  Parse Arg n, p
  q = p - 1
  s = 0
  Do while q // 2 == 0
    q = q % 2
    s = s+1
    End
  if s == 1 Then
    return pow(n, (p + 1) % 4, p)
  Do z=2 To p
    if p - 1 == legendre(z, p) Then
      Leave
    End
  c = pow(z, q, p)
  r = pow(n, (q + 1) / 2, p)
  t = pow(n, q, p)
  m = s
  t2 = 0
  Do while (t - 1) // p <> 0
    t2 = (t * t) // p
    Do i=1 To m
      if (t2 - 1) // p == 0 Then
        Leave
      t2 = (t2 * t2) // p
      End
    y=2**(m - i - 1)
    b = pow(c, y, p)
    If b=10008 Then Trace ?R
    r = (r * b) // p
    c = (b * b) // p
    t = (t * c) // p
    m = i
    End
  return r
pow: Procedure
  Parse Arg x,y,z
  If y>0 Then
    p=x**y
  Else p=x
  If z>'' Then
    p=p//z
  Return p

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

Sidef

Translation of: Python

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}"
}

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

Visual Basic .NET

Translation of: C#

Imports System.Numerics

Module Module1

    Class Solution
        ReadOnly root1 As BigInteger
        ReadOnly root2 As BigInteger
        ReadOnly exists As Boolean

        Sub New(r1 As BigInteger, r2 As BigInteger, e As Boolean)
            root1 = r1
            root2 = r2
            exists = e
        End Sub

        Public Function GetRoot1() As BigInteger
            Return root1
        End Function

        Public Function GetRoot2() As BigInteger
            Return root2
        End Function

        Public Function GetExists() As Boolean
            Return exists
        End Function
    End Class

    Function Ts(n As BigInteger, p As BigInteger) As Solution
        If BigInteger.ModPow(n, (p - 1) / 2, p) <> 1 Then
            Return New Solution(0, 0, False)
        End If

        Dim q As BigInteger = p - 1
        Dim ss = BigInteger.Zero
        While (q Mod 2) = 0
            ss += 1
            q >>= 1
        End While

        If ss = 1 Then
            Dim r1 = BigInteger.ModPow(n, (p + 1) / 4, p)
            Return New Solution(r1, p - r1, True)
        End If

        Dim z As BigInteger = 2
        While BigInteger.ModPow(z, (p - 1) / 2, p) <> p - 1
            z += 1
        End While
        Dim c = BigInteger.ModPow(z, q, p)
        Dim r = BigInteger.ModPow(n, (q + 1) / 2, p)
        Dim t = BigInteger.ModPow(n, q, p)
        Dim m = ss

        Do
            If t = 1 Then
                Return New Solution(r, p - r, True)
            End If
            Dim i = BigInteger.Zero
            Dim zz = t
            While zz <> 1 AndAlso i < (m - 1)
                zz = zz * zz Mod p
                i += 1
            End While
            Dim b = c
            Dim e = m - i - 1
            While e > 0
                b = b * b Mod p
                e = e - 1
            End While
            r = r * b Mod p
            c = b * b Mod p
            t = t * c Mod p
            m = i
        Loop
    End Function

    Sub Main()
        Dim pairs = New List(Of Tuple(Of Long, Long)) From {
            New Tuple(Of Long, Long)(10, 13),
            New Tuple(Of Long, Long)(56, 101),
            New Tuple(Of Long, Long)(1030, 10009),
            New Tuple(Of Long, Long)(1032, 10009),
            New Tuple(Of Long, Long)(44402, 100049),
            New Tuple(Of Long, Long)(665820697, 1000000009),
            New Tuple(Of Long, Long)(881398088036, 1000000000039)
        }

        For Each pair In pairs
            Dim sol = Ts(pair.Item1, pair.Item2)
            Console.WriteLine("n = {0}", pair.Item1)
            Console.WriteLine("p = {0}", pair.Item2)
            If sol.GetExists() Then
                Console.WriteLine("root1 = {0}", sol.GetRoot1())
                Console.WriteLine("root2 = {0}", sol.GetRoot2())
            Else
                Console.WriteLine("No solution exists")
            End If
            Console.WriteLine()
        Next

        Dim bn = BigInteger.Parse("41660815127637347468140745042827704103445750172002")
        Dim bp = BigInteger.Pow(10, 50) + 577
        Dim bsol = Ts(bn, bp)
        Console.WriteLine("n = {0}", bn)
        Console.WriteLine("p = {0}", bp)
        If bsol.GetExists() Then
            Console.WriteLine("root1 = {0}", bsol.GetRoot1())
            Console.WriteLine("root2 = {0}", bsol.GetRoot2())
        Else
            Console.WriteLine("No solution exists")
        End If
    End Sub

End Module

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

Wren

Translation of: Kotlin

Library: Wren-dynamic

Library: Wren-big

import "./dynamic" for Tuple
import "./big" for BigInt

var Solution = Tuple.create("Solution", ["root1", "root2", "exists"])

var ts = Fn.new { |n, p|
    if (n is Num) n = BigInt.new(n)
    if (p is Num) p = BigInt.new(p)

    var powModP = Fn.new { |a, e| a.modPow(e, p) }

    var ls = Fn.new { |a| powModP.call(a, p.dec / BigInt.two) }

    if (ls.call(n) != BigInt.one) return Solution.new(BigInt.zero, BigInt.zero, false)
    var q = p.dec
    var ss = BigInt.zero
    while (q & BigInt.one == BigInt.zero) {
        ss = ss.inc
        q = q >> 1
    }
    if (ss == BigInt.one) {
        var r1 = powModP.call(n, p.inc / BigInt.four)
        return Solution.new(r1, p - r1, true)
    }
    var z = BigInt.two
    while (ls.call(z) != p.dec) z = z.inc
    var c = powModP.call(z, q)
    var r = powModP.call(n, q.inc/BigInt.two)
    var t = powModP.call(n, q)
    var m = ss
    while (true) {
        if (t == BigInt.one) return Solution.new(r, p - r, true)
        var i = BigInt.zero
        var zz = t
        while (zz != BigInt.one && i < m.dec) {
            zz = zz * zz % p
            i = i.inc
        }
        var b = c
        var e = m - i.inc
        while (e > BigInt.zero) {
            b = b * b % p
            e = e.dec
        }
        r = r * b % p
        c = b * b % p
        t = t * c % p
        m = i
    }
}

var pairs = [
    [10, 13], [56, 101], [1030, 10009], [1032, 10009], [44402, 100049],
    [665820697, 1000000009], [881398088036, 1000000000039]
]

for (pair in pairs) {
    var n = pair[0]
    var p = pair[1]
    var sol = ts.call(n, p)
    System.print("n     = %(n)")
    System.print("p     = %(p)")
    if (sol.exists) {
        System.print("root1 = %(sol.root1)")
        System.print("root2 = %(sol.root2)")
    } else {
        System.print("No solution exists")
    }
    System.print()
}

var bn = BigInt.new("41660815127637347468140745042827704103445750172002")
var bp = BigInt.ten.pow(50) + BigInt.new(577)
var bsol = ts.call(bn, bp)
System.print("n     = %(bn)")
System.print("p     = %(bp)")
if (bsol.exists) {
    System.print("root1 = %(bsol.root1)")
    System.print("root2 = %(bsol.root2)")
} else {
    System.print("No solution exists")
}

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

zkl

Translation of: EchoLisp

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
}

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));
}

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