Cipolla's algorithm

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

Cipolla's algorithm

Solve x² ≡ n (mod p)

In computational number theory, Cipolla's algorithm is a technique for solving an equation of the form x² ≡ n (mod p), where p is an odd prime and x ,n ∊ Fp = {0, 1, ... p-1}.

To apply the algorithm we need the Legendre symbol, and arithmetic in Fp².

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 is a ≡ 0

Arithmetic in Fp²

Let ω a symbol such as ω² is a member of Fp and not a square, x and y members of Fp. The set Fp² is defined as {x + ω y }. The subset { x + 0 ω} of Fp² is Fp. Fp² is somewhat equivalent to the field of complex number, with ω analoguous to i, and i² = -1 . Remembering that all operations are modulo p, addition, multiplication and exponentiation in Fp² are defined as :

• (x1 + ω y1) + (x2 + ω y2) := (x1 + x2 + ω (y1 + y2))
• (x1 + ω y1) * (x2 + ω y2) := (x1*x2 + y1*y2*ω²) + ω (x1*y2 + x2*y1)
  • (0 + ω) * (0 + ω) := (ω² + 0 ω) ≡ ω² in Fp
• (x1 + ω y1) ^ n := (x + ω y) * (x + ω y) * ... ( n times) (1)

Algorithm pseudo-code

• Input : p an odd prime, and n ≠ 0 in Fp
• Step 0. Check that n is indeed a square  : (n | p) must be ≡ 1
• Step 1. Find, by trial and error, an a > 0 such as (a² - n) is not a square : (a²-n | p) must be ≡ -1.
• Step 2. Let ω² = a² - n. Compute, in Fp2 : (a + ω) ^ ((p + 1)/2) (mod p)

To compute this step, use a pair of numbers, initially [a,1], and use repeated "multiplication" which is defined such that [c,d] times [e,f] is (mod p) [ c*c + ω²*f*f, d*e + c*f ].

• Step 3. Check that the result is ≡ x + 0 * ω in Fp2, that is x in Fp.
• Step 4. Output the two positive solutions, x and p - x (mod p).
• Step 5. Check that x * x ≡ n (mod p)

Example from Wikipedia

n := 10 , p := 13
Legendre(10,13) → 1         // 10 is indeed a square
a := 2                      // try
ω² := a*a - 10             // ≡ 7 ≡ -6
Legendre (ω² , 13) → -1    // ok - not square
(2 + ω) ^ 7 → 6 + 0 ω      // by modular exponentiation (1)
                            // 6 and (13 - 6) = 7 are solutions
(6 * 6) % 13 → 10           // = n . Checked.

Task

Implement the above.

Find solutions (if any) for

• n = 10 p = 13
• n = 56 p = 101
• n = 8218 p = 10007
• n = 8219 p = 10007
• n = 331575 p = 1000003

Extra credit

• n 665165880 p 1000000007
• n 881398088036 p 1000000000039
• n = 34035243914635549601583369544560650254325084643201 p = 10^50 + 151

See also:

• Modular exponentiation
• Tonelli-Shanks algorithm

EchoLisp

<lang scheme> (lib 'struct) (lib 'types) (lib 'bigint)

test equality mod p

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

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

Arithmetic in Fp²

(struct Fp² ( x y ))

a + b

(define (Fp²-add Fp²:a Fp²:b p ω2) (Fp² (% (+ a.x b.x) p) (% (+ a.y b.y) p)))

a * b

(define (Fp²-mul Fp²:a Fp²:b p ω2) (Fp² (% (+ (* a.x b.x) (* ω2 a.y b.y)) p) (% (+ (* a.x b.y) (* a.y b.x)) p)))

a * a

(define (Fp²-square Fp²:a p ω2) (Fp² (% (+ (* a.x a.x) (* ω2 a.y a.y)) p) (% (* 2 a.x a.y) p)))

a ^ n

(define (Fp²-pow Fp²:a n p ω2) (cond ((= 0 n) (Fp² 1 0)) ((= 1 n) (Fp² a.x a.y)) ((= 2 n) (Fp²-mul a a p ω2)) ((even? n) (Fp²-square (Fp²-pow a (/ n 2) p ω2) p ω2)) (else (Fp²-mul a (Fp²-pow a (1- n) p ω2) p ω2))))

x^2 ≡ n (mod p) ?

(define (Cipolla n p)

check n is a square

(unless (= 1 (Legendre n p)) (error "not a square (mod p)" (list n p)))

iterate until suitable 'a' found

(define a (for ((t (in-range 2 p))) ;; t = tentative a #:break (= (1- p) (Legendre (- (* t t) n) p)) => t )) (define ω2 (- (* a a) n)) ;; (writeln 'a-> a 'ω2-> ω2 'ω-> 'ω) ;; (Fp² a 1) = a + ω (define r (Fp²-pow (Fp² a 1) (/ (1+ p) 2) p ω2)) ;; (writeln 'r r) (define x (Fp²-x r)) (assert (zero? (Fp²-y r))) ;; hope that ω has vanished (assert (mod= n (* x x) p)) ;; checking the result (printf "Roots of %d are (%d,%d) (mod %d)" n x (% (- p x) p) p)) </lang>

Output:

(Cipolla 10 13)
Roots of 10 are (6,7) (mod 13)
(% (* 6 6) 13) → 10 ;; checking

(Cipolla 56 101)
Roots of 56 are (37,64) (mod 101)

(Cipolla 8218 10007)
Roots of 8218 are (9872,135) (mod 10007)

Cipolla 8219 10007)
❌ error: not a square (mod p) (8219 10007)

(Cipolla 331575 1000003)
Roots of 331575 are (855842,144161) (mod 1000003)
(% ( * 855842 855842) 1000003) → 331575

FreeBASIC

Had a close look at the EchoLisp code for step 2. Used the FreeBASIC code from the Miller-Rabin task for prime testing. <lang freebasic>' version 08-04-2017 ' compile with: fbc -s console ' maximum for p is 17 digits to be on the save side

' TRUE/FALSE are built-in constants since FreeBASIC 1.04 ' But we have to define them for older versions.

1. Ifndef TRUE

   #Define FALSE 0
   #Define TRUE Not FALSE

1. EndIf

Type fp2

   x As LongInt
   y As LongInt

End Type

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 ' a mod modulus, but a is already smaller then 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

Function fp2mul(a As fp2, b As fp2, p As LongInt, w2 As LongInt) As fp2

   Dim As fp2 answer
   Dim As ULongInt tmp1, tmp2
   ' needs to be broken down in smaller steps to avoid overflow
   ' answer.x = (a.x * b.x + a.y * b.y * w2) Mod p
   ' answer.y = (a.x * b.y + a.y * b.x) Mod p
   tmp1 = mul_mod(a.x, b.x, p)
   tmp2 = mul_mod(a.y, b.y, p)
   tmp2 = mul_mod(tmp2, w2, p)
   answer.x = (tmp1 + tmp2) Mod p
   tmp1 = mul_mod(a.x, b.y, p)
   tmp2 = mul_mod(a.y, b.x, p)
   answer.y = (tmp1 + tmp2) Mod p

   Return answer

End Function

Function fp2square(a As fp2, p As LongInt, w2 As LongInt) As fp2

   Return fp2mul(a, a, p, w2)

End Function

Function fp2pow(a As fp2, n As LongInt, p As LongInt, w2 As LongInt) As fp2

   If n = 0 Then Return Type (1, 0)
   If n = 1 Then Return a
   If n = 2 Then Return fp2square(a, p, w2)
   If (n And 1) = 0 Then
       Return fp2square(fp2pow(a, n\2, p, w2), p , w2)
   Else
       Return fp2mul(a, fp2pow(a, n -1, p, w2), p, w2)
   End If

End Function

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

Data 10, 13, 56, 101, 8218, 10007,8219, 10007 Data 331575, 1000003, 665165880, 1000000007 Data 881398088036, 1000000000039

Randomize Timer Dim As LongInt n, p, a, w2 Dim As LongInt i, x1, x2 Dim As fp2 answer

For i = 1 To 7

   Read n, p
   Print
   Print "Find solution for n =";n ; " and p =";p

   If p = 2 OrElse Isprime(p,15) = FALSE Then
       Print "No solution, p is not a odd prime"
       Continue For
   End If

   ' p is checked and is a odd prime
   If legendre_symbol(n, p) <> 1 Then
       Print n; " is not a square in F";Str(p)
       Continue For
   End If

   Do 
       Do
           a = Rnd * (p -2) +2
           w2 = a * a - n
       Loop Until legendre_symbol(w2, p) = -1

       answer = Type(a, 1)
       answer = fp2pow(answer, (p +1) \ 2, p, w2)
       If answer.y <> 0 Then Continue Do

       x1 = answer.x : x2 = p - x1
       If mul_mod(x1, x1, p) = n AndAlso mul_mod(x2, x2, p) = n Then
           Print "Solution found: x1 ="; x1; ", "; "x2 ="; x2
           Exit Do
       End If
   Loop            ' loop until solution is found

Next

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

Output:

Find solution for n = 10 and p = 13
Solution found: x1 = 7, x2 = 6

Find solution for n = 56 and p = 101
Solution found: x1 = 37, x2 = 64

Find solution for n = 8218 and p = 10007
Solution found: x1 = 9872, x2 = 135

Find solution for n = 8219 and p = 10007
 8219 is not a square in F10007

Find solution for n = 331575 and p = 1000003
Solution found: x1 = 144161, x2 = 855842

Find solution for n = 665165880 and p = 1000000007
Solution found: x1 = 475131702, x2 = 524868305

Find solution for n = 881398088036 and p = 1000000000039
Solution found: x1 = 791399408049, x2 = 208600591990

Go

int

Implementation following the pseudocode in the task description. <lang go>package main

import "fmt"

func c(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)
   }
   // Step 0, validate arguments
   if ls(n) != 1 {
       return
   }
   // Step 1, find a, ω2
   var a, ω2 int
   for a = 0; ; a++ {
       // integer % in Go uses T-division, add p to keep the result positive
       ω2 = (a*a + p - n) % p
       if ls(ω2) == p-1 {
           break
       }
   }
   // muliplication in fp2
   type point struct{ x, y int }
   mul := func(a, b point) point {
       return point{(a.x*b.x + a.y*b.y*ω2) % p, (a.x*b.y + b.x*a.y) % p}
   }
   // Step2, compute power
   r := point{1, 0}
   s := point{a, 1}
   for n := (p + 1) >> 1 % p; n > 0; n >>= 1 {
       if n&1 == 1 {
           r = mul(r, s)
       }
       s = mul(s, s)
   }
   // Step3, check x in Fp
   if r.y != 0 {
       return
   }
   // Step5, check x*x=n
   if r.x*r.x%p != n {
       return
   }
   // Step4, solutions
   return r.x, p - r.x, true

}

func main() {

   fmt.Println(c(10, 13))
   fmt.Println(c(56, 101))
   fmt.Println(c(8218, 10007))
   fmt.Println(c(8219, 10007))
   fmt.Println(c(331575, 1000003))

}</lang>

Output:

6 7 true
37 64 true
9872 135 true
0 0 false
855842 144161 true

big.Int

Extra credit: <lang go>package main

import (

   "fmt"
   "math/big"

)

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

   if big.Jacobi(&n, &p) != 1 {
       return
   }
   var one, a, ω2 big.Int
   one.SetInt64(1)
   for ; ; a.Add(&a, &one) {
       // big.Int Mod uses Euclidean division, result is always >= 0
       ω2.Mod(ω2.Sub(ω2.Mul(&a, &a), &n), &p)
       if big.Jacobi(&ω2, &p) == -1 {
           break
       }
   }
   type point struct{ x, y big.Int }
   mul := func(a, b point) (z point) {
       var w big.Int
       z.x.Mod(z.x.Add(z.x.Mul(&a.x, &b.x), w.Mul(w.Mul(&a.y, &a.y), &ω2)), &p)
       z.y.Mod(z.y.Add(z.y.Mul(&a.x, &b.y), w.Mul(&b.x, &a.y)), &p)
       return
   }
   var r, s point
   r.x.SetInt64(1)
   s.x.Set(&a)
   s.y.SetInt64(1)
   var e big.Int
   for e.Rsh(e.Add(&p, &one), 1); len(e.Bits()) > 0; e.Rsh(&e, 1) {
       if e.Bit(0) == 1 {
           r = mul(r, s)
       }
       s = mul(s, s)
   }
   R2.Sub(&p, &r.x)
   return r.x, R2, true

}

func main() {

   var n, p big.Int
   n.SetInt64(665165880)
   p.SetInt64(1000000007)
   R1, R2, ok := c(n, p)
   fmt.Println(&R1, &R2, ok)

   n.SetInt64(881398088036)
   p.SetInt64(1000000000039)
   R1, R2, ok = c(n, p)
   fmt.Println(&R1, &R2, ok)

   n.SetString("34035243914635549601583369544560650254325084643201", 10)
   p.SetString("100000000000000000000000000000000000000000000000151", 10)
   R1, R2, ok = c(n, p)
   fmt.Println(&R1)
   fmt.Println(&R2)

}</lang>

Output:

475131702 524868305 true
791399408049 208600591990 true
82563118828090362261378993957450213573687113690751
17436881171909637738621006042549786426312886309400

J

Based on the echolisp implementation:

<lang J>leg=: dyad define

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

)

mul2=: conjunction define

 m| (*&{. + n**&{:), (+/ .* |.)

)

pow2=: conjunction define

 if. 0=y do. 1 0
 elseif. 1=y do. x
 elseif. 2=y do. x (m mul2 n) x
 elseif. 0=2|y do. (m mul2 n)~ x (m pow2 n) y%2
 elseif. do. x (m mul2 n) x (m pow2 n) y-1
 end.

)

cipolla=: dyad define

 assert. 1=1 p: y [ 'y must be prime'
 assert. 1= x leg y [ 'x must be square mod y'
 a=.1 
 whilst. (0 ~:{: r) do. a=. a+1
   while. 1>: leg&y@(x -~ *:) a do. a=.a+1 end.
   w2=. y|(*:a) - x
   r=. (a,1) (y pow2 w2) (y+1)%2
 end.
 if. x =&(y&|) *:{.r do.
   y|(,-){.r
 else.
   smoutput 'got ',":~.y|(,-){.r
   assert. 'not a valid square root'
 end.

)</lang>

Task examples:

<lang J> 10 cipolla 13 6 7

  56 cipolla 101

37 64

  8218 cipolla 10007

9872 135

  8219 cipolla 10007

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

  331575 cipolla 1000003

855842 144161

  665165880x cipolla 1000000007x

524868305 475131702

  881398088036x cipolla 1000000000039x

208600591990 791399408049

  34035243914635549601583369544560650254325084643201x cipolla (10^50x) + 151

17436881171909637738621006042549786426312886309400 82563118828090362261378993957450213573687113690751</lang>

Perl 6

Works with: Rakudo version 2016.10

Translation of: Sidef

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

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

}

1. a coordinate in a Field of p elements

class Fp {

   has Int $.x;
   has Int $.y;

}

sub cipolla ( Int \𝑛, Int \𝑝 ) {

   note "Invalid parameters ({𝑛}, {𝑝})"
     and return Nil if (𝑛│𝑝) != 1;
   my $ω2;
   my $a = 0;
   loop {
       last if ($ω2 = ($a² - 𝑛) % 𝑝)│𝑝 < 0;
       $a++;
   }

   # define a local multiply operator for Field coordinates
   multi sub infix:<*> ( Fp $a, Fp $b ){
       Fp.new: :x(($a.x * $b.x + $a.y * $b.y * $ω2) % 𝑝),
               :y(($a.x * $b.y + $b.x * $a.y)       % 𝑝)
   }

   my $r = Fp.new: :x(1),  :y(0);
   my $s = Fp.new: :x($a), :y(1);

   for (𝑝+1) +> 1, * +> 1 ... 1 {
       $r *= $s if $_ % 2;
       $s *= $s;
   }
   return Nil if $r.y;
   $r.x;

}

my @tests = (

   (10, 13),
   (56, 101),
   (8218, 10007),
   (8219, 10007),
   (331575, 1000003),
   (665165880, 1000000007),
   (881398088036, 1000000000039),
   (34035243914635549601583369544560650254325084643201,
     100000000000000000000000000000000000000000000000151)

);

for @tests -> ($n, $p) {

  my $r = cipolla($n, $p);
  say $r ?? "Roots of $n are ($r, {$p-$r}) mod $p"
         !! "No solution for ($n, $p)"

} </lang>

Output:

Roots of 10 are (6, 7) mod 13
Roots of 56 are (37, 64) mod 101
Roots of 8218 are (9872, 135) mod 10007
Invalid parameters (8219, 10007)
No solution for (8219, 10007)
Roots of 331575 are (855842, 144161) mod 1000003
Roots of 665165880 are (475131702, 524868305) mod 1000000007
Roots of 881398088036 are (791399408049, 208600591990) mod 1000000000039
Roots of 34035243914635549601583369544560650254325084643201 are (82563118828090362261378993957450213573687113690751, 17436881171909637738621006042549786426312886309400) mod 100000000000000000000000000000000000000000000000151

Sage

Works with: Sage version 7.6

<lang sage> def eulerCriterion(a, p):

   return -1 if pow(a, int((p-1)/2), p) == p-1 else 1

def cipollaMult(x1, y1, x2, y2, u, p):

   return ((x1*x2 + y1*y2*u) % p), ((x1*y2 + x2*y1) % p)

def cipollaAlgorithm(a, p, t):

   if not is_prime(p):
       print "❌ " + str(p) + " is not a prime."
       return False
   if eulerCriterion(a, p) == -1:
       print "❌ " + str(a) + " is not a quadratic residue modulo " + str(p)
       return False

   a = pow(a, 1, p)
   out = []

   for i in range(t):
       t = randrange(2, p)
       u = pow(t**2 - a, 1, p)
       while (eulerCriterion(u, p) == 1):
           t = randrange(2, p)
           u = pow(t**2 - a, 1, p)

       x0, y0 = t, 1
       x, y = t, 1
       for i in range(int((p + 1) / 2) - 1):
           x, y = cipollaMult(x, y, x0, y0, u, p)

       if x not in out:
           out.append(x)

   return sorted(out)

</lang>

Output:

sage: cipollaAlgorithm(10, 13, 30)
[6, 7]
sage: cipollaAlgorithm(56, 101, 30)
[37, 64]
sage: cipollaAlgorithm(8218, 10007, 30)
[135, 9872]
sage: cipollaAlgorithm(331575, 1000003, 30)
[144161, 855842]
sage: cipollaAlgorithm(8219, 10007, 30)
❌ 8219 is not a quadratic residue modulo 10007
False

Sidef

Translation of: Go

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

   legendre(n, p) == 1 || return nil

   var (a = 0, ω2 = 0)
   loop {
       ω2 = ((a*a - n) % p)
       if (legendre(ω2, p) == -1) {
           break
       }
       ++a
   }

   struct point { x, y }

   func mul(a, b) {
       point((a.x*b.x + a.y*b.y*ω2) % p, (a.x*b.y + b.x*a.y) % p)
   }

   var r = point(1, 0)
   var s = point(a, 1)

   for (var n = ((p+1) >> 1); n > 0; n >>= 1) {
       r = mul(r, s) if n.is_odd
       s = mul(s, s)
   }

   r.y == 0 ? r.x : nil

}

var tests = [

   [10, 13],
   [56, 101],
   [8218, 10007],
   [8219, 10007],
   [331575, 1000003],
   [665165880, 1000000007],
   [881398088036 1000000000039],
   [34035243914635549601583369544560650254325084643201, 10**50 + 151],

]

for n,p in tests {

   var r = cipolla(n, p)
   if (defined(r)) {
       say "Roots of #{n} are (#{r} #{p-r}) mod #{p}"
   } else {
       say "No solution for (#{n}, #{p})"
   }

}</lang>

Output:

Roots of 10 are (6 7) mod 13
Roots of 56 are (37 64) mod 101
Roots of 8218 are (9872 135) mod 10007
No solution for (8219, 10007)
Roots of 331575 are (855842 144161) mod 1000003
Roots of 665165880 are (475131702 524868305) mod 1000000007
Roots of 881398088036 are (791399408049 208600591990) mod 1000000000039
Roots of 34035243914635549601583369544560650254325084643201 are (82563118828090362261378993957450213573687113690751 17436881171909637738621006042549786426312886309400) mod 100000000000000000000000000000000000000000000000151

zkl

Translation of: EchoLisp

Uses lib GMP (GNU MP Bignum Library). <lang zkl>var [const] BN=Import("zklBigNum"); //libGMP fcn modEq(a,b,p) { (a-b)%p==0 } fcn Legendre(a,p){ a.powm((p - 1)/2,p) }

class Fp2{ // Arithmetic in Fp^2

  fcn init(_x,_y){ var [const] x=BN(_x), y=BN(_y) }	// two big ints
  //fcn add(b,p){ self((x + b.x)%p,(y + b.y)%p) }	// a + b
  fcn mul(b,p,w2){ self(( x*b.x + y*b.y*w2 )%p, (x*b.y + y*b.x) %p) } // a * b
  fcn square(p,w2){ mul(self,p,w2) }          	// a * a == self.mul(self,p,w2)
  fcn pow(n,p,w2){				// a ^ n
     if     (n==0)     self(1,0);
     else if(n==1)     self;
     else if(n==2)     square(p,w2);
     else if(n.isEven) pow(n/2,p,w2).square(p,w2);
     else 		mul(pow(n-1,p,w2),p,w2)
  }

}

fcn Cipolla(n,p){ n=BN(n); // x^2 == n (mod p) ?

  if(Legendre(n,p)!=1)   // check n is a square
     throw(Exception.AssertionError("not a square (mod p)"+vm.arglist));
  // iterate until suitable 'a' found (the first one found)
  a:=[BN(2)..p].filter1('wrap(a){ Legendre(a*a-n,p)==(p-1) });
  w2:=a*a - n;
  r:=Fp2(a,1).pow((p + 1)/2,p,w2);	    // (Fp2 a 1) = a + w2
  x:=r.x;
  _assert_(r.y==0,"r.y==0 : "+r.y);	    // hope that w has vanished
  _assert_(modEq(n,x*x,p),"modEq(n,x*x,p)"); // checking the result
  println("Roots of %d are (%d,%d)  (mod %d)".fmt(n,x,(p-x)%p,p));
  return(x,(p-x)%p);

}</lang> <lang zkl>foreach n,p in (T( T(10,13),T(56,101),T(8218,10007),T(8219,10007),T(331575,1000003), T(665165880,1000000007),T(881398088036,1000000000039), T("34035243914635549601583369544560650254325084643201", BN(10).pow(50) + 151) )){

  try{ Cipolla(n,p) }catch{ println(__exception) }

}</lang>

Output:

Roots of 10 are (6,7)  (mod 13)
Roots of 56 are (37,64)  (mod 101)
Roots of 8218 are (9872,135)  (mod 10007)
AssertionError(not a square (mod p)L(8219,10007))
Roots of 331575 are (855842,144161)  (mod 1000003)
Roots of 665165880 are (524868305,475131702)  (mod 1000000007)
Roots of 881398088036 are (208600591990,791399408049)  (mod 1000000000039)
Roots of 34035243914635549601583369544560650254325084643201 are (17436881171909637738621006042549786426312886309400,82563118828090362261378993957450213573687113690751)  (mod 100000000000000000000000000000000000000000000000151)