Cipolla's algorithm

From Rosetta Code
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:



11l

Translation of: Python
F convertToBase(n, b)
   I (n < 2)
      R [n]
   V temp = n
   V ans = [Int]()
   L (temp != 0)
      ans = [temp % b] [+] ans
      temp I/= b
   R ans

F cipolla(=n, p)
   n %= p
   I (n == 0 | n == 1)
      R [n, (p - n) % p]
   V phi = p - 1
   I (pow(BigInt(n), phi I/ 2, p) != 1)
      R [Int]()
   I (p % 4 == 3)
      V ans = Int(pow(BigInt(n), (p + 1) I/ 4, p))
      R [ans, (p - ans) % p]
   V aa = 0
   L(i) 1 .< p
      V temp = pow(BigInt((i * i - n) % p), phi I/ 2, p)
      I (temp == phi)
         aa = i
         L.break
   V exponent = convertToBase((p + 1) I/ 2, 2)
   F cipollaMult(ab, cd, w, p)
      V (a, b) = ab
      V (c, d) = cd
      R ((a * c + b * d * w) % p, (a * d + b * c) % p)
   V x1 = (aa, 1)
   V x2 = cipollaMult(x1, x1, aa * aa - n, p)
   L(i) 1 .< exponent.len
      I (exponent[i] == 0)
         x2 = cipollaMult(x2, x1, aa * aa - n, p)
         x1 = cipollaMult(x1, x1, aa * aa - n, p)
      E
         x1 = cipollaMult(x1, x2, aa * aa - n, p)
         x2 = cipollaMult(x2, x2, aa * aa - n, p)
   R [x1[0], (p - x1[0]) % p]

print(‘Roots of 2 mod 7: ’cipolla(2, 7))
print(‘Roots of 8218 mod 10007: ’cipolla(8218, 10007))
print(‘Roots of 56 mod 101: ’cipolla(56, 101))
print(‘Roots of 1 mod 11: ’cipolla(1, 11))
print(‘Roots of 8219 mod 10007: ’cipolla(8219, 10007))
Output:
Roots of 2 mod 7: [4, 3]
Roots of 8218 mod 10007: [9872, 135]
Roots of 56 mod 101: [37, 64]
Roots of 1 mod 11: [1, 10]
Roots of 8219 mod 10007: []

C

Translation of: FreeBasic
#include <stdbool.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

struct fp2 {
    int64_t x, y;
};

uint64_t randULong(uint64_t min, uint64_t max) {
    uint64_t t = (uint64_t)rand();
    return min + t % (max - min);
}

// returns a * b mod modulus
uint64_t mul_mod(uint64_t a, uint64_t b, uint64_t modulus) {
    uint64_t x = 0, y = a % modulus;

    while (b > 0) {
        if ((b & 1) == 1) {
            x = (x + y) % modulus;
        }
        y = (y << 1) % modulus;
        b = b >> 1;
    }

    return x;
}

//returns b ^^ power mod modulus
uint64_t pow_mod(uint64_t b, uint64_t power, uint64_t modulus) {
    uint64_t x = 1;

    while (power > 0) {
        if ((power & 1) == 1) {
            x = mul_mod(x, b, modulus);
        }
        b = mul_mod(b, b, modulus);
        power = power >> 1;
    }

    return x;
}

// miller-rabin prime test
bool isPrime(uint64_t n, int64_t k) {
    uint64_t a, x, n_one = n - 1, d = n_one;
    uint32_t s = 0;
    uint32_t r;

    if (n < 2) {
        return false;
    }

    // limit 2^63, pow_mod/mul_mod can't handle bigger numbers
    if (n > 9223372036854775808ull) {
        printf("The number is too big, program will end.\n");
        exit(1);
    }

    if ((n % 2) == 0) {
        return n == 2;
    }

    while ((d & 1) == 0) {
        d = d >> 1;
        s = s + 1;
    }

    while (k > 0) {
        k = k - 1;
        a = randULong(2, n);
        x = pow_mod(a, d, n);
        if (x == 1 || x == n_one) {
            continue;
        }
        for (r = 1; r < s; r++) {
            x = pow_mod(x, 2, n);
            if (x == 1) return false;
            if (x == n_one) goto continue_while;
        }
        if (x != n_one) {
            return false;
        }

    continue_while: {}
    }

    return true;
}

int64_t legendre_symbol(int64_t a, int64_t p) {
    int64_t x = pow_mod(a, (p - 1) / 2, p);
    if ((p - 1) == x) {
        return x - p;
    } else {
        return x;
    }
}

struct fp2 fp2mul(struct fp2 a, struct fp2 b, int64_t p, int64_t w2) {
    struct fp2 answer;
    uint64_t tmp1, tmp2;

    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) % p;
    tmp1 = mul_mod(a.x, b.y, p);
    tmp2 = mul_mod(a.y, b.x, p);
    answer.y = (tmp1 + tmp2) % p;

    return answer;
}

struct fp2 fp2square(struct fp2 a, int64_t p, int64_t w2) {
    return fp2mul(a, a, p, w2);
}

struct fp2 fp2pow(struct fp2 a, int64_t n, int64_t p, int64_t w2) {
    struct fp2 ret;

    if (n == 0) {
        ret.x = 1;
        ret.y = 0;
        return ret;
    }
    if (n == 1) {
        return a;
    }
    if ((n & 1) == 0) {
        return fp2square(fp2pow(a, n / 2, p, w2), p, w2);
    } else {
        return fp2mul(a, fp2pow(a, n - 1, p, w2), p, w2);
    }
}

void test(int64_t n, int64_t p) {
    int64_t a, w2;
    int64_t x1, x2;
    struct fp2 answer;

    printf("Find solution for n = %lld and p = %lld\n", n, p);
    if (p == 2 || !isPrime(p, 15)) {
        printf("No solution, p is not an odd prime.\n\n");
        return;
    }

    //p is checked and is a odd prime
    if (legendre_symbol(n, p) != 1) {
        printf(" %lld is not a square in F%lld\n\n", n, p);
        return;
    }

    while (true) {
        do {
            a = randULong(2, p);
            w2 = a * a - n;
        } while (legendre_symbol(w2, p) != -1);

        answer.x = a;
        answer.y = 1;
        answer = fp2pow(answer, (p + 1) / 2, p, w2);
        if (answer.y != 0) {
            continue;
        }

        x1 = answer.x;
        x2 = p - x1;
        if (mul_mod(x1, x1, p) == n && mul_mod(x2, x2, p) == n) {
            printf("Solution found: x1 = %lld, x2 = %lld\n\n", x1, x2);
            return;
        }
    }
}

int main() {
    srand((size_t)time(0));

    test(10, 13);
    test(56, 101);
    test(8218, 10007);
    test(8219, 10007);
    test(331575, 1000003);
    test(665165880, 1000000007);
    //test(881398088036, 1000000000039);

    return 0;
}
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 = 524868305, x2 = 475131702

C# 

using System;
using System.Numerics;

namespace CipollaAlgorithm {
    class Program {
        static readonly BigInteger BIG = BigInteger.Pow(10, 50) + 151;

        private static Tuple<BigInteger, BigInteger, bool> C(string ns, string ps) {
            BigInteger n = BigInteger.Parse(ns);
            BigInteger p = ps.Length > 0 ? BigInteger.Parse(ps) : BIG;

            // Legendre symbol. Returns 1, 0, or p-1
            BigInteger ls(BigInteger a0) => BigInteger.ModPow(a0, (p - 1) / 2, p);

            // Step 0: validate arguments
            if (ls(n) != 1) {
                return new Tuple<BigInteger, BigInteger, bool>(0, 0, false);
            }

            // Step 1: Find a, omega2
            BigInteger a = 0;
            BigInteger omega2;
            while (true) {
                omega2 = (a * a + p - n) % p;
                if (ls(omega2) == p - 1) {
                    break;
                }
                a += 1;
            }

            // Multiplication in Fp2
            BigInteger finalOmega = omega2;
            Tuple<BigInteger, BigInteger> mul(Tuple<BigInteger, BigInteger> aa, Tuple<BigInteger, BigInteger> bb) {
                return new Tuple<BigInteger, BigInteger>(
                    (aa.Item1 * bb.Item1 + aa.Item2 * bb.Item2 * finalOmega) % p,
                    (aa.Item1 * bb.Item2 + bb.Item1 * aa.Item2) % p
                );
            }

            // Step 2: Compute power
            Tuple<BigInteger, BigInteger> r = new Tuple<BigInteger, BigInteger>(1, 0);
            Tuple<BigInteger, BigInteger> s = new Tuple<BigInteger, BigInteger>(a, 1);
            BigInteger nn = ((p + 1) >> 1) % p;
            while (nn > 0) {
                if ((nn & 1) == 1) {
                    r = mul(r, s);
                }
                s = mul(s, s);
                nn >>= 1;
            }

            // Step 3: Check x in Fp
            if (r.Item2 != 0) {
                return new Tuple<BigInteger, BigInteger, bool>(0, 0, false);
            }

            // Step 5: Check x * x = n
            if (r.Item1 * r.Item1 % p != n) {
                return new Tuple<BigInteger, BigInteger, bool>(0, 0, false);
            }

            // Step 4: Solutions
            return new Tuple<BigInteger, BigInteger, bool>(r.Item1, p - r.Item1, true);
        }

        static void Main(string[] args) {
            Console.WriteLine(C("10", "13"));
            Console.WriteLine(C("56", "101"));
            Console.WriteLine(C("8218", "10007"));
            Console.WriteLine(C("8219", "10007"));
            Console.WriteLine(C("331575", "1000003"));
            Console.WriteLine(C("665165880", "1000000007"));
            Console.WriteLine(C("881398088036", "1000000000039"));
            Console.WriteLine(C("34035243914635549601583369544560650254325084643201", ""));
        }
    }
}
Output:
(6, 7, True)
(37, 64, True)
(9872, 135, True)
(0, 0, False)
(855842, 144161, True)
(475131702, 524868305, True)
(791399408049, 208600591990, True)
(82563118828090362261378993957450213573687113690751, 17436881171909637738621006042549786426312886309400, True)

D

Translation of: Kotlin
import std.bigint;
import std.stdio;
import std.typecons;

enum BIGZERO = BigInt(0);

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

alias Point = Tuple!(BigInt, "x", BigInt, "y");
alias Triple = Tuple!(BigInt, "x", BigInt, "y", bool, "b");

Triple c(string ns, string ps) {
    auto n = BigInt(ns);
    BigInt p;
    if (ps.length > 0) {
        p = BigInt(ps);
    } else {
        p = BigInt(10)^^50 + 151;
    }

    // Legendre symbol, returns 1, 0 or p - 1
    auto ls = (BigInt a) => modPow(a, (p-1)/2, p);

    // Step 0, validate arguments
    if (ls(n) != 1) return Triple(BIGZERO, BIGZERO, false);

    // Step 1, find a, omega2
    auto a = BIGZERO;
    BigInt omega2;
    while (true) {
        omega2 = (a * a + p - n) % p;
        if (ls(omega2) == p-1) break;
        a++;
    }

    // multiplication in Fp2
    auto mul = (Point aa, Point bb) => Point(
        (aa.x * bb.x + aa.y * bb.y * omega2) % p,
        (aa.x * bb.y + bb.x * aa.y) % p
    );

    // Step 2, compute power
    auto r = Point(BigInt(1), BIGZERO);
    auto s = Point(a, BigInt(1));
    auto nn = ((p+1) >> 1) % p;
    while (nn > 0) {
        if ((nn & 1) == 1) r = mul(r, s);
        s = mul(s, s);
        nn >>= 1;
    }

    // Step 3, check x in Fp
    if (r.y != 0) return Triple(BIGZERO, BIGZERO, false);

    // Step 5, check x * x = n
    if (r.x*r.x%p!=n) return Triple(BIGZERO, BIGZERO, false);

    // Step 4, solutions
    return Triple(r.x, p-r.x, true);
}

void main() {
    writeln(c("10", "13"));
    writeln(c("56", "101"));
    writeln(c("8218", "10007"));
    writeln(c("8219", "10007"));
    writeln(c("331575", "1000003"));
    writeln(c("665165880", "1000000007"));
    writeln(c("881398088036", "1000000000039"));
    writeln(c("34035243914635549601583369544560650254325084643201", ""));
}
Output:
Tuple!(BigInt, "x", BigInt, "y", bool, "b")(6, 7, true)
Tuple!(BigInt, "x", BigInt, "y", bool, "b")(37, 64, true)
Tuple!(BigInt, "x", BigInt, "y", bool, "b")(9872, 135, true)
Tuple!(BigInt, "x", BigInt, "y", bool, "b")(0, 0, false)
Tuple!(BigInt, "x", BigInt, "y", bool, "b")(855842, 144161, true)
Tuple!(BigInt, "x", BigInt, "y", bool, "b")(475131702, 524868305, true)
Tuple!(BigInt, "x", BigInt, "y", bool, "b")(791399408049, 208600591990, true)
Tuple!(BigInt, "x", BigInt, "y", bool, "b")(82563118828090362261378993957450213573687113690751, 17436881171909637738621006042549786426312886309400, true)

EchoLisp 

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

F#

The function

This task uses Extensible Prime Generator (F#) 

// Cipolla's algorithm. Nigel Galloway: June 16th., 2019
let Cipolla n g =
  let rec fN i g e l=match e with n when n=0I->i |_ when e%2I=1I->fN ((i*g)%l) ((g*g)%l) (e/2I) l |_-> fN i ((g*g)%l) (e/2I) l
  let rec fG g=match (n/g+g)>>>1 with n when bigint.Abs(g-n)>>>1<2I->n+1I |g->fG g
  let a,b=let rec fI i=let q=i*i-n in if fN 1I q ((g-1I)/2I) g>1I then (i,q) else fI (i+1I) in fI(fG (bigint(sqrt(double n))))
  let fE=Seq.unfold(fun(n,i)->Some((n,i),((n*n+i*i*b)%g,(2I*n*i)%g)))(a,1I)|>Seq.cache
  let rec fL Πn Πi α β=match 2I**α with
                        Ω when Ω<β->fL Πn Πi (α+1) β
                       |Ω when Ω>β->let n,i=Seq.item (α-1) fE in fL ((Πn*n+Πi*i*b)%g) ((Πn*i+Πi*n)%g) 0 (β-Ω/2I)    
                       |_->let n,i=Seq.item α fE in ((Πn*n+Πi*i*b)%g)
  if fN 1I n ((g-1I)/2I) g<>1I then None else Some(fL 1I 0I 0 ((g+1I)/2I))

The Task

let test=[(10I,13I);(56I,101I);(8218I,10007I);(8219I,10007I);(331575I,1000003I);(665165880I,1000000007I);(881398088036I,1000000000039I);(34035243914635549601583369544560650254325084643201I,10I**50+151I)] 
test|>List.iter(fun(n,g)->match Cipolla n g with Some r->printfn "Cipolla %A %A -> %A (%A) check %A" n g r (g-r) ((r*r)%g) |_->printfn "Cipolla %A %A -> has no result" n g)
Output:
Cipolla 10 13 -> 7 (6) check 10
Cipolla 56 101 -> 64 (37) check 56
Cipolla 8218 10007 -> 135 (9872) check 8218
Cipolla 8219 10007 -> has no result
Cipolla 331575 1000003 -> 144161 (855842) check 331575
Cipolla 665165880 1000000007 -> 475131702 (524868305) check 665165880
Cipolla 881398088036 1000000000039 -> 208600591990 (791399408049) check 881398088036
Cipolla 34035243914635549601583369544560650254325084643201 100000000000000000000000000000000000000000000000151 -> 17436881171909637738621006042549786426312886309400 (82563118828090362261378993957450213573687113690751) check 34035243914635549601583369544560650254325084643201
Real: 00:00:00.089, CPU: 00:00:00.090, GC gen0: 2, gen1: 0

Factor

Translation of: Sidef
Works with: Factor version 0.99 2020-08-14
USING: accessors assocs interpolate io kernel literals locals
math math.extras math.functions ;

TUPLE: point x y ;
C: <point> point

:: (cipolla) ( n p -- m )
    0 0 :> ( a! ω2! )
    [ ω2 p legendere -1 = ]
    [ a sq n - p rem ω2! a 1 + a! ] do until a 1 - a!

    [| a b |
        a x>> b x>> * a y>> b y>> ω2 * * + p mod
        a x>> b y>> * b x>> a y>> * + p mod <point>
    ] :> [mul]

    1 0 <point> :> r!
    a 1 <point> :> s!
    p 1 + -1 shift :> n!

    [ n 0 > ] [
        n odd? [ r s [mul] call r! ] when
        s s [mul] call s!
        n -1 shift n!
    ] while

    r y>> zero? r x>> f ? ;
     

: cipolla ( n p -- m/f )
    2dup legendere 1 = [ (cipolla) ] [ 2drop f ] if ;

! Task
{
    { 10 13 }
    { 56 101 }
    { 8218 10007 }
    { 8219 10007 }
    { 331575 1000003 }
    { 665165880 1000000007 }
    { 881398088036 1000000000039 }
    ${
        34035243914635549601583369544560650254325084643201
        10 50 ^ 151 +
    }
}
[
    2dup cipolla
    [ 2dup - [I Roots of ${3} are (${1} ${0}) mod ${2}I] ]
    [ [I No solution for (${}, ${})I] ] if* nl
] assoc-each
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

FreeBASIC

LongInt version

Had a close look at the EchoLisp code for step 2. Used the FreeBASIC code from the Miller-Rabin task for prime testing. 

' 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.
#Ifndef TRUE
    #Define FALSE 0
    #Define TRUE Not FALSE
#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 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

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

GMP version

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

#Include Once "gmp.bi"

Type fp2
    x As Mpz_ptr
    y As Mpz_ptr
End Type

Data "10", "13"
Data "56", "101"
Data "8218", "10007"
Data "8219", "10007"
Data "331575", "1000003"
Data "665165880", "1000000007"
Data "881398088036", "1000000000039"
Data "34035243914635549601583369544560650254325084643201"  ', 10^50 + 151

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

    Dim As fp2 r
    r.x = Allocate(Len(__Mpz_struct)) : Mpz_init(r.x)
    r.y = Allocate(Len(__Mpz_struct)) : Mpz_init(r.y)

    Mpz_mul (r.x, a.y, b.y)
    Mpz_mul (r.x, r.x, w2)
    Mpz_addmul(r.x, a.x, b.x)
    Mpz_mod (r.x, r.x, p)
    Mpz_mul (r.y, a.x, b.y)
    Mpz_addmul(r.y, a.y, b.x)
    Mpz_mod (r.y, r.y, p)

    Return r

End Function

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

    Return fp2mul(a, a, p, w2)

End Function

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

    If Mpz_cmp_ui(n, 0) = 0 Then
        Mpz_set_ui(a.x, 1)
        Mpz_set_ui(a.y, 0)
        Return a
    End If
    If Mpz_cmp_ui(n, 1) = 0 Then Return a
    If Mpz_cmp_ui(n, 2) = 0 Then Return fp2square(a, p, w2)
    If Mpz_tstbit(n, 0) = 0 Then
        Mpz_fdiv_q_2exp(n, n, 1) ' even
        Return fp2square(fp2pow(a, n, p, w2), p, w2)
    Else
        Mpz_sub_ui(n, n, 1)      ' odd
        Return fp2mul(a, fp2pow(a, n, p, w2), p, w2)
    End If

End Function

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

Dim As Long i
Dim As ZString Ptr zstr
Dim As String n_str, p_str

Dim As Mpz_ptr a, n, p, p2, w2, x1, x2
a  = Allocate(Len(__Mpz_struct)) : Mpz_init(a)
n  = Allocate(Len(__Mpz_struct)) : Mpz_init(n)
p  = Allocate(Len(__Mpz_struct)) : Mpz_init(p)
p2 = Allocate(Len(__Mpz_struct)) : Mpz_init(p2)
w2 = Allocate(Len(__Mpz_struct)) : Mpz_init(w2)
x1 = Allocate(Len(__Mpz_struct)) : Mpz_init(x1)
x2 = Allocate(Len(__Mpz_struct)) : Mpz_init(x2)

Dim As fp2 answer
answer.x = Allocate(Len(__Mpz_struct)) : Mpz_init(answer.x)
answer.y = Allocate(Len(__Mpz_struct)) : Mpz_init(answer.y)

For i = 1 To 8
    Read n_str
    Mpz_set_str(n, n_str, 10)
    If i < 8 Then
        Read p_str
        Mpz_set_str(p, p_str, 10)
    Else
        p_str = "10^50 + 151" ' set up last n
        Mpz_set_str(p, "1" + String(50, "0"), 10)
        Mpz_add_ui(p, p, 151)
    End If

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

    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

    ' p is checked and is a odd prime
    ' legendre symbol needs to be 1
    If Mpz_legendre(n, p) <> 1 Then
        Print n_str; " is not a square in F"; p_str
        Print
        Continue For
    End If

    Mpz_set_ui(a, 1)
    Do
        Do
            Do
                Mpz_add_ui(a, a, 1)
                Mpz_mul(w2, a, a)
                Mpz_sub(w2, w2, n)
            Loop Until Mpz_legendre(w2, p) = -1

            Mpz_set(answer.x, a)
            Mpz_set_ui(answer.y, 1)
            Mpz_add_ui(p2, p, 1)       ' p2 = p + 1
            Mpz_fdiv_q_2exp(p2, p2, 1) ' p2 = p2 \ 2 (p2 shr 1)

            answer = fp2pow(answer, p2, p, w2)

        Loop Until Mpz_cmp_ui(answer.y, 0) = 0
        Mpz_set(x1, answer.x)
        Mpz_sub(x2, p, x1)
        Mpz_powm_ui(a, x1, 2, p)
        Mpz_powm_ui(p2, x2, 2, p)
        If Mpz_cmp(a, n) = 0 AndAlso Mpz_cmp(p2, n) = 0 Then Exit Do
    Loop

    zstr = Mpz_get_str(0, 10, x1)
    Print "Solution found: x1 = "; *zstr;
    zstr = Mpz_get_str(0, 10, x2)
    Print ", x2 = "; *zstr
    Print
Next

Mpz_clear(x1) : Mpz_clear(p2) : Mpz_clear(p) : Mpz_clear(a) : Mpz_clear(n)
Mpz_clear(x2) : Mpz_clear(w2) : Mpz_clear(answer.x) : Mpz_clear(answer.y)

' 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: x1 = 6, x2 = 7

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 = 855842, x2 = 144161

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

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

Find solution for n = 34035243914635549601583369544560650254325084643201 and p = 10^50 + 151
Solution found: x1 = 17436881171909637738621006042549786426312886309400, x2 = 82563118828090362261378993957450213573687113690751

Go

int

Implementation following the pseudocode in the task description. 

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))
}
Output:
6 7 true
37 64 true
9872 135 true
0 0 false
855842 144161 true

big.Int

Extra credit: 

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)
}
Output:
475131702 524868305 true
791399408049 208600591990 true
82563118828090362261378993957450213573687113690751
17436881171909637738621006042549786426312886309400

J

Based on the echolisp implementation: 

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

Task examples: 

   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

Java

Translation of: Kotlin
Works with: Java version 8
import java.math.BigInteger;
import java.util.function.BiFunction;
import java.util.function.Function;

public class CipollasAlgorithm {
    private static final BigInteger BIG = BigInteger.TEN.pow(50).add(BigInteger.valueOf(151));
    private static final BigInteger BIG_TWO = BigInteger.valueOf(2);

    private static class Point {
        BigInteger x;
        BigInteger y;

        Point(BigInteger x, BigInteger y) {
            this.x = x;
            this.y = y;
        }

        @Override
        public String toString() {
            return String.format("(%s, %s)", this.x, this.y);
        }
    }

    private static class Triple {
        BigInteger x;
        BigInteger y;
        boolean b;

        Triple(BigInteger x, BigInteger y, boolean b) {
            this.x = x;
            this.y = y;
            this.b = b;
        }

        @Override
        public String toString() {
            return String.format("(%s, %s, %s)", this.x, this.y, this.b);
        }
    }

    private static Triple c(String ns, String ps) {
        BigInteger n = new BigInteger(ns);
        BigInteger p = !ps.isEmpty() ? new BigInteger(ps) : BIG;

        // Legendre symbol, returns 1, 0 or p - 1
        Function<BigInteger, BigInteger> ls = (BigInteger a)
            -> a.modPow(p.subtract(BigInteger.ONE).divide(BIG_TWO), p);

        // Step 0, validate arguments
        if (!ls.apply(n).equals(BigInteger.ONE)) {
            return new Triple(BigInteger.ZERO, BigInteger.ZERO, false);
        }

        // Step 1, find a, omega2
        BigInteger a = BigInteger.ZERO;
        BigInteger omega2;
        while (true) {
            omega2 = a.multiply(a).add(p).subtract(n).mod(p);
            if (ls.apply(omega2).equals(p.subtract(BigInteger.ONE))) {
                break;
            }
            a = a.add(BigInteger.ONE);
        }

        // multiplication in Fp2
        BigInteger finalOmega = omega2;
        BiFunction<Point, Point, Point> mul = (Point aa, Point bb) -> new Point(
            aa.x.multiply(bb.x).add(aa.y.multiply(bb.y).multiply(finalOmega)).mod(p),
            aa.x.multiply(bb.y).add(bb.x.multiply(aa.y)).mod(p)
        );

        // Step 2, compute power
        Point r = new Point(BigInteger.ONE, BigInteger.ZERO);
        Point s = new Point(a, BigInteger.ONE);
        BigInteger nn = p.add(BigInteger.ONE).shiftRight(1).mod(p);
        while (nn.compareTo(BigInteger.ZERO) > 0) {
            if (nn.and(BigInteger.ONE).equals(BigInteger.ONE)) {
                r = mul.apply(r, s);
            }
            s = mul.apply(s, s);
            nn = nn.shiftRight(1);
        }

        // Step 3, check x in Fp
        if (!r.y.equals(BigInteger.ZERO)) {
            return new Triple(BigInteger.ZERO, BigInteger.ZERO, false);
        }

        // Step 5, check x * x = n
        if (!r.x.multiply(r.x).mod(p).equals(n)) {
            return new Triple(BigInteger.ZERO, BigInteger.ZERO, false);
        }

        // Step 4, solutions
        return new Triple(r.x, p.subtract(r.x), true);
    }

    public static void main(String[] args) {
        System.out.println(c("10", "13"));
        System.out.println(c("56", "101"));
        System.out.println(c("8218", "10007"));
        System.out.println(c("8219", "10007"));
        System.out.println(c("331575", "1000003"));
        System.out.println(c("665165880", "1000000007"));
        System.out.println(c("881398088036", "1000000000039"));
        System.out.println(c("34035243914635549601583369544560650254325084643201", ""));
    }
}
Output:
(6, 7, true)
(37, 64, true)
(9872, 135, true)
(0, 0, false)
(855842, 144161, true)
(475131702, 524868305, true)
(791399408049, 208600591990, true)
(82563118828090362261378993957450213573687113690751, 17436881171909637738621006042549786426312886309400, true)

jq

Translation of: Wren

Works with gojq, the Go implementation of jq

A Point is represented by a numeric array of length two (i.e. [x,y]). 

# To take advantage of gojq's arbitrary-precision integer arithmetic:
def power($b): . as $in | reduce range(0;$b) as $i (1; . * $in);

# 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 bigBig: (10|power(50)) + 151;

# The remainder when $b raised to the power $e is divided by $m:
def modPow($b; $e; $m):
   if ($m == 1) then 0
   else {r: 1, b: ($b % $m), $e}
   | until (.e <= 0;
       if (.e % 2) == 1 then .r = (.r*.b) % $m else . end
       | .e |= idivide(2)
       | .b |= ((.*.) % $m) )
   | .r
   end;

def c($ns; $ps):
    $ns as $n
    | (if $ps != "" then $ps else bigBig end) as $p
 
    # Legendre symbol, returns 1, 0 or p - 1
    | def ls($a): modPow($a; ($p - 1) | idivide(2); $p);

    # multiplication in Fp2 where .omega2 comes from .
    def mul($aa; $bb):
      [($aa[0] * $bb[0] + $aa[1] * $bb[1] * .omega2) % $p,
       ($aa[0] * $bb[1] + $bb[0] * $aa[1]) % $p] ;

    # Step 0, validate arguments
    if (ls($n) != 1) then [0, 0, false]
    else
    # Step 1, find a, omega2
    { a: 0, stop: false }
    | until( .stop;
        .omega2 = ((.a * .a + $p - $n) % $p)
        | if ls(.omega2) == ($p - 1)
	  then .stop = true
          else .a += 1
          end )

    # Step 2, compute power
    | { r: [1, 0],
        s: [.a, 1],
        nn: (( ($p + 1) | idivide(2)) % $p),
	omega2  }
     | until (.nn <= 0;
        if (.nn % 2) == 1 then .r = mul(.r; .s) else . end
        | .s = mul(.s; .s)
        | .nn |= idivide(2)  )
 
    # Step 3, check x in Fp; or that x * x = n
    | if (.r[1] != 0) or ((.r[0] * .r[0]) % $p != $n) then [0, 0, false]
 
    # Step 4, solutions
      else  [.r[0], $p - .r[0], true]
      end
    end;

def exercise:
c(10; 13),
c(56; 101),
c(8218; 10007),
c(8219; 10007),
c(331575; 1000003),
c(665165880; 1000000007),
c(881398088036; 1000000000039),
c(34035243914635549601583369544560650254325084643201; "")
;

exercise
Output:
[6,7,true]
[37,64,true]
[9872,135,true]
[0,0,false]
[855842,144161,true]
[475131702,524868305,true]
[791399408049,208600591990,true]
[82563118828090362261378993957450213573687113690751,17436881171909637738621006042549786426312886309400,true]

Julia

Translation of: Perl
using Primes

function legendre(n, p)
    if p != 2 && isprime(p)
        x = powermod(BigInt(n), div(p - 1, 2), p)
        return x == 0 ? 0 : x == 1 ? 1 : -1
    end
    return -1
end

function cipolla(n, p)
    if legendre(n, p) != 1
        return NaN
    end
    a, w2 = BigInt(0), BigInt(0)
    while true
        w2 = (a^2 + p - n) % p
        if legendre(w2, p) < 0
            break
        end
        a += 1
    end
    r, s, i = (1, 0), (a, 1), p + 1
    while (i >>= 1) > 0
        if isodd(i)
            r = ((r[1] * s[1] + r[2] * s[2] * w2) % p, (r[1] * s[2] + s[1] * r[2]) % p)
        end
        s = ((s[1] * s[1] + s[2] * s[2] * w2) % p, (2 * s[1] * s[2]) % p)
    end
    return r[2] != 0 ? NaN : r[1]
end

const ctests = [(10, 13),
                (56, 101),
                (8218, 10007),
                (8219, 10007),
                (331575, 1000003),
                (665165880, 1000000007),
                (881398088036, 1000000000039),
                (big"34035243914635549601583369544560650254325084643201",
                    big"100000000000000000000000000000000000000000000000151")]

for (n, p) in ctests
   r = cipolla(n, p)
   println(r > 0 ? "Roots of $n are ($r, $(p - r)) mod $p." : "No solution for ($n, $p)")
end
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.

Kotlin

Translation of: Go
// version 1.2.0

import java.math.BigInteger

class Point(val x: BigInteger, val y: BigInteger)

val bigZero = BigInteger.ZERO
val bigOne  = BigInteger.ONE
val bigTwo  = BigInteger.valueOf(2L)
val bigBig  = BigInteger.TEN.pow(50) + BigInteger.valueOf(151L)

fun c(ns: String, ps: String): Triple<BigInteger, BigInteger, Boolean> {
    val n = BigInteger(ns)
    val p = if (!ps.isEmpty()) BigInteger(ps) else bigBig

    // Legendre symbol, returns 1, 0 or p - 1
    fun ls(a: BigInteger) = a.modPow((p - bigOne) / bigTwo, p)

    // Step 0, validate arguments
    if (ls(n) != bigOne) return Triple(bigZero, bigZero, false)

    // Step 1, find a, omega2
    var a = bigZero
    var omega2: BigInteger
    while (true) {
        omega2 = (a * a + p - n) % p
        if (ls(omega2) == p - bigOne) break
        a++
    }

    // multiplication in Fp2
    fun mul(aa: Point, bb: Point) =
        Point(
            (aa.x * bb.x + aa.y * bb.y * omega2) % p,
            (aa.x * bb.y + bb.x * aa.y) % p
        )

    // Step 2, compute power
    var r = Point(bigOne, bigZero)
    var s = Point(a, bigOne)
    var nn = ((p + bigOne) shr 1) % p
    while (nn > bigZero) {
        if ((nn and bigOne) == bigOne) r = mul(r, s)
        s = mul(s, s)
        nn = nn shr 1
    }

    // Step 3, check x in Fp
    if (r.y != bigZero) return Triple(bigZero, bigZero, false)

    // Step 5, check x * x = n
    if (r.x * r.x % p != n) return Triple(bigZero, bigZero, false)

    // Step 4, solutions
    return Triple(r.x, p - r.x, true)
}

fun main(args: Array<String>) {
    println(c("10", "13"))
    println(c("56", "101"))
    println(c("8218", "10007"))
    println(c("8219", "10007"))
    println(c("331575", "1000003"))
    println(c("665165880", "1000000007"))
    println(c("881398088036", "1000000000039"))
    println(c("34035243914635549601583369544560650254325084643201", ""))
}
Output:
(6, 7, true)
(37, 64, true)
(9872, 135, true)
(0, 0, false)
(855842, 144161, true)
(475131702, 524868305, true)
(791399408049, 208600591990, true)
(82563118828090362261378993957450213573687113690751, 17436881171909637738621006042549786426312886309400, true)

Mathematica /Wolfram Language 

ClearAll[Cipolla]
Cipolla[n_, p_] := Module[{ls, omega2, nn, a, r, s},
  ls = JacobiSymbol[n, p];
  If[ls != 1,
   {0, 0, False}
   ,
   a = 0;
   While[True,
    omega2 = Mod[a^2 + p - n, p];
    If[JacobiSymbol[omega2, p] == -1, Break[]];
    a++
    ];
   ClearAll[Mul];
   Mul[{ax_, ay_}, {bx_, by_}] := 
    Mod[{ax bx + ay by omega2, ax by + bx ay}, p];
   r = {1, 0};
   s = {a, 1};
   nn = Mod[BitShiftRight[p + 1, 1], p];
   While[nn > 0,
    If[BitAnd[nn, 1] == 1,
     r = Mul[r, s]
     ];
    s = Mul[s, s];
    nn = BitShiftRight[nn, 1];
    ];
   
   If[r[[2]] != 0, Return[{0, 0, False}]];
   If[Mod[r[[1]]^2, p] != n, Return[{0, 0, False}]];
   Return[{r[[1]], p - r[[1]], True}]
   ]
  ]
Cipolla[10, 13]
Cipolla[56, 101]
Cipolla[8218, 10007]
Cipolla[8219, 10007]
Cipolla[331575, 1000003]
Cipolla[665165880, 1000000007]
Cipolla[881398088036, 1000000000039]
Cipolla[34035243914635549601583369544560650254325084643201, 10^50 + 151]
Output:
{6,7,True}
{37,64,True}
{9872,135,True}
{0,0,False}
{855842,144161,True}
{475131702,524868305,True}
{791399408049,208600591990,True}
{82563118828090362261378993957450213573687113690751,17436881171909637738621006042549786426312886309400,True}

Nim

Translation of: Kotlin
Library: bignum
import options
import bignum

let
  Zero = newInt(0)
  One = newInt(1)
  BigBig = newInt(10)^50 + 15

type
  Point = tuple[x, y: Int]
  Solutions = (Int, Int)


proc exp(x, y, m: Int): Int =
  ## Missing function in "bignum" module.
  if m == 1: return Zero
  result = newInt(1)
  var x = x mod m
  var y = y.clone
  while not y.isZero:
    if y mod 2 == 1:
      result = result * x mod m
    y = y shr 1
    x = x * x mod m


proc c(ns, ps: string): Option[Solutions] =
  let n = newInt(ns)
  let p = if ps.len != 0: newInt(ps) else: BigBig

  # Legendre symbol: returns 1, 0 or p - 1.
  proc ls(a: Int): Int = a.exp((p - 1) div 2, p)

  # Step 0, validate arguments.
  if ls(n) != One: return none(Solutions)

  # Step 1, find a, omega2.
  var a = newInt(0)
  var omega2: Int
  while true:
    omega2 = (a * a + p - n) mod p
    if ls(omega2) == p - 1: break
    a += 1

  # Multiplication in Fp2.
  proc `*`(a, b: Point): Point =
    ((a.x * b.x + a.y * b.y * omega2) mod p,
     (a.x * b.y + b.x * a.y) mod p)

  # Step 2, compute power.
  var
    r: Point = (One, Zero)
    s: Point = (a, One)
    nn = ((p + 1) shr 1) mod p
  while not nn.isZero:
    if (nn and One) == One: r = r * s
    s = s * s
    nn = nn shr 1

  # Step 3, check x in Fp.
  if not r.y.isZero: return none(Solutions)

   # Step 5, check x * x = n.
  if r.x * r.x mod p != n: return none(Solutions)

   # Step 4, solutions.
  result = some((r.x, p - r.x))


when isMainModule:

  const Values = [("10", "13"), ("56", "101"), ("8218", "10007"),
                  ("8219", "10007"), ("331575", "1000003"),
                  ("665165880", "1000000007"), ("881398088036", "1000000000039"),
                  ("34035243914635549601583369544560650254325084643201", "")]

  for (n, p) in Values:
    let sols = c(n, p)
    if sols.isSome: echo sols.get()
    else: echo "No solutions."
Output:
(6, 7)
(37, 64)
(9872, 135)
No solutions.
(855842, 144161)
(475131702, 524868305)
(791399408049, 208600591990)
(82563118828090362261378993957450213573687113690751, 17436881171909637738621006042549786426312886309400)

Perl

Translation of: Raku
Library: ntheory
use bigint;
use ntheory qw(is_prime);

sub Legendre {
    my($n,$p) = @_;
    return -1 unless $p != 2 && is_prime($p);
    my $x = ($n->as_int())->bmodpow(int(($p-1)/2), $p); # $n coerced to BigInt
    if    ($x==0) { return  0 }
    elsif ($x==1) { return  1 }
    else          { return -1 }
}

sub Cipolla {
    my($n, $p) = @_;
    return undef if Legendre($n,$p) != 1;

    my $w2;
    my $a = 0;
    $a++ until Legendre(($w2 = ($a**2 - $n) % $p), $p) < 0;

    my %r = ( x=> 1,  y=> 0 );
    my %s = ( x=> $a, y=> 1 );
    my $i = $p + 1;
    while (1 <= ($i >>= 1)) {
        %r = ( x => (($r{x} * $s{x} + $r{y} * $s{y} * $w2) % $p),
               y => (($r{x} * $s{y} + $s{x} * $r{y})       % $p)
             ) if $i % 2;
        %s = ( x => (($s{x} * $s{x} + $s{y} * $s{y} * $w2) % $p),
               y => (($s{x} * $s{y} + $s{x} * $s{y})       % $p)
             )
    }
    $r{y} ? undef : $r{x}
}

my @tests = (
    (10, 13),
    (56, 101),
    (8218, 10007),
    (8219, 10007),
    (331575, 1000003),
    (665165880, 1000000007),
    (881398088036, 1000000000039),
);

while (@tests) {
   $n = shift @tests;
   $p = shift @tests;
   my $r = Cipolla($n, $p);
   $r ? printf "Roots of %d are (%d, %d) mod %d\n", $n, $r, $p-$r, $p
      : print  "No solution for ($n, $p)\n"
}
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

Phix

Translation of: Kotlin
Library: Phix/mpfr
with javascript_semantics 
include mpfr.e 
 
procedure legendre(mpz r, a, p)
-- Legendre symbol, returns 1, 0 or p - 1 (in r)
    mpz_sub_ui(r,p,1)
    {} = mpz_fdiv_q_ui(r, r, 2)
    mpz_powm(r,a,r,p)
end procedure 
 
procedure mul_point(sequence a, b, mpz omega2, p)
-- (modifies a)
    mpz {xa,ya} = a,
        {xb,yb} = b,
        xaxb = mpz_init(),
        yayb = mpz_init(),
        xayb = mpz_init(),
        xbya = mpz_init()
    mpz_mul(xaxb,xa,xb)
    mpz_mul(yayb,ya,yb)
    mpz_mul(xayb,xa,yb)
    mpz_mul(xbya,xb,ya)
    mpz_mul(yayb,yayb,omega2)
    mpz_add(xaxb,xaxb,yayb)
    mpz_mod(xa,xaxb,p)      -- xa := mod(xaxb+yayb*omega2,p)
    mpz_add(xayb,xayb,xbya)
    mpz_mod(ya,xayb,p)      -- ya := mod(xayb+xbya,p)
    {xaxb,yayb,xayb,xbya} = mpz_free({xaxb,yayb,xayb,xbya})
end procedure
 
function cipolla(object no, po)
mpz n = mpz_init(no),
    p = mpz_init(po),
    t = mpz_init()
 
    -- Step 0, validate arguments
    legendre(t,n,p)
    if mpz_cmp_si(t,1)!=0 then return {"0","0","false"} end if
 
    -- Step 1, find a, omega2
    integer a = 0
    mpz omega2 = mpz_init(),
        pm1 = mpz_init()
    mpz_sub_ui(pm1,p,1)
    while true do
        mpz_sub(t,p,n)
        mpz_add_ui(t,t,a*a)
        mpz_mod(omega2,t,p)
        legendre(t,omega2,p)
        if mpz_cmp(t,pm1)=0 then exit end if
        a += 1
    end while
 
    -- Step 2, compute power
    sequence r = {mpz_init(1),mpz_init(0)},
             s = {mpz_init(a),mpz_init(1)}
    mpz nn = mpz_init()
    mpz_add_ui(nn,p,1)
    {} = mpz_fdiv_q_ui(nn, nn, 2)
    mpz_mod(nn,nn,p)
    while mpz_cmp_si(nn,0)>0 do
        if mpz_fdiv_ui(nn,2)=1 then
            mul_point(r,s,omega2,p)
        end if
        mul_point(s,s,omega2,p)
        {} = mpz_fdiv_q_ui(nn, nn, 2)
    end while
 
    -- Step 3, check x in Fp
    if mpz_cmp_si(r[2],0)!=0 then return {"0","0","false"} end if 
 
    -- Step 5, check x * x = n
    mpz_powm_ui(t,r[1],2,p)
    if mpz_cmp(t,n)!=0 then return {"0","0","false"} end if 
 
    -- Step 4, solutions
    mpz_sub(p,p,r[1])
    return {mpz_get_str(r[1]), mpz_get_str(p), "true"}
end function
 
constant tests = {{10,13},
                  {56,101},
                  {8218,10007},
                  {8219,10007},
                  {331575,1000003},
                  {665165880,1000000007},
                  {"881398088036","1000000000039"},
                  {"34035243914635549601583369544560650254325084643201",
                   "100000000000000000000000000000000000000000000000151"}}
 
for i=1 to length(tests) do
    object {n,p} = tests[i]
    ?{n,p,cipolla(n,p)}
end for

Obviously were you to use that in anger, you would probably rip out a few mpz_get_str() and return NULL/false rather than "0"/"false", etc.

Output:
{10,13,{"6","7","true"}}
{56,101,{"37","64","true"}}
{8218,10007,{"9872","135","true"}}
{8219,10007,{"0","0","false"}}
{331575,1000003,{"855842","144161","true"}}
{665165880,1000000007,{"475131702","524868305","true"}}
{"881398088036","1000000000039",{"791399408049","208600591990","true"}}
{"34035243914635549601583369544560650254325084643201","100000000000000000000000000000000000000000000000151",
 {"82563118828090362261378993957450213573687113690751","17436881171909637738621006042549786426312886309400","true"}}

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 mul ("A" B P W2)
   (let (A (copy "A")  B (copy B))
      (set
         "A"
         (%
            (+
               (* (car A) (car B))
               (* (cadr A) (cadr B) W2) )
            P )
         (cdr "A")
         (%
            (+
               (* (car A) (cadr B))
               (* (car B) (cadr A)) )
            P ) ) ) )
(de ci (N P)
   (and
      (=1 (legendre N P))
      (let
         (A 0
            W2 0
            R NIL
            S NIL )
         (loop
            (setq W2
               (% (- (+ (* A A) P) N) P) )
            (T (= (dec P) (legendre W2 P)))
            (inc 'A) )
         (setq R (list 1 0)  S (list A 1))
         (for
            (N
               (% (>> 1 (inc P)) P)
               (> N 0)
               (>> 1 N) )
            (and (bit? 1 N) (mul R S P W2))
            (mul S S P W2) )
         (=0 (cadr R))
         (=
            N
            (% (* (car R) (car R)) P) )
         (list (car R) (- P (car R))) ) ) )

(println (ci 10 13))
(println (ci 56 101))
(println (ci 8218 10007))
(println (ci 8219 10007))
(println (ci 331575 1000003))
(println (ci 665165880 1000000007))
(println (ci 881398088036 1000000000039))
(println (ci 34035243914635549601583369544560650254325084643201 (+ (** 10 50) 151)))
Output:
(6 7)
(37 64)
(9872 135)
NIL
(855842 144161)
(475131702 524868305)
(791399408049 208600591990)
(82563118828090362261378993957450213573687113690751 17436881171909637738621006042549786426312886309400)

Python 

#Converts n to base b as a list of integers between 0 and b-1
#Most-significant digit on the left
def convertToBase(n, b):
	if(n < 2):
		return [n];
	temp = n;
	ans = [];
	while(temp != 0):
		ans = [temp % b]+ ans;
		temp /= b;
	return ans;

#Takes integer n and odd prime p
#Returns both square roots of n modulo p as a pair (a,b)
#Returns () if no root
def cipolla(n,p):
	n %= p
	if(n == 0 or n == 1):
		return (n,-n%p)
	phi = p - 1
	if(pow(n, phi/2, p) != 1):
		return ()
	if(p%4 == 3):
		ans = pow(n,(p+1)/4,p)
		return (ans,-ans%p)
	aa = 0
	for i in xrange(1,p):
		temp = pow((i*i-n)%p,phi/2,p)
		if(temp == phi):
			aa = i
			break;
	exponent = convertToBase((p+1)/2,2)
	def cipollaMult((a,b),(c,d),w,p):
		return ((a*c+b*d*w)%p,(a*d+b*c)%p)
	x1 = (aa,1)
	x2 = cipollaMult(x1,x1,aa*aa-n,p)
	for i in xrange(1,len(exponent)):
		if(exponent[i] == 0):
			x2 = cipollaMult(x2,x1,aa*aa-n,p)
			x1 = cipollaMult(x1,x1,aa*aa-n,p)
		else:
			x1 = cipollaMult(x1,x2,aa*aa-n,p)
			x2 = cipollaMult(x2,x2,aa*aa-n,p)
	return (x1[0],-x1[0]%p)

print(f"Roots of 2 mod 7: {cipolla(2, 7)}")
print(f"Roots of 8218 mod 10007: {cipolla(8218, 10007)}")
print(f"Roots of 56 mod 101: {cipolla(56, 101)}")
print(f"Roots of 1 mod 11: {cipolla(1, 11)}")
print(f"Roots of 8219 mod 10007: {cipolla(8219, 10007)}")
Output:
Roots of 2 mod 7: (4, 3)
Roots of 8218 mod 10007: (9872, 135)
Roots of 56 mod 101: (37, 64)
Roots of 1 mod 11: (1, 10)
Roots of 8219 mod 10007: ()

Racket

Translation of: EchoLisp
#lang racket

(require math/number-theory)

;; math/number-theory allows us to parameterize a "current-modulus"
;; which obviates the need for p to be passed around constantly
(define (Cipolla n p) (with-modulus p (mod-Cipolla n)))

(define (mod-Legendre a)  
  (modexpt a (/ (sub1 (current-modulus)) 2)))
 
;; Arithmetic in Fp² 
(struct Fp² (x y))

(define-syntax-rule (Fp²-destruct* (a a.x a.y) ...)
  (begin (match-define (Fp² a.x a.y) a) ...)  )

;; a + b
(define (Fp²-add a b ω2)
  (Fp²-destruct* (a a.x a.y) (b b.x b.y))
  (Fp² (mod+ a.x b.x) (mod+ a.y b.y)))

;; a * b
(define (Fp²-mul a b ω2) 
  (Fp²-destruct* (a a.x a.y) (b b.x b.y))
  (Fp² (mod+ (* a.x b.x) (* ω2 a.y b.y)) (mod+ (* a.x b.y) (* a.y b.x))))
 
;; a * a	
(define (Fp²-square a ω2)
  (Fp²-destruct* (a a.x a.y))
  (Fp² (mod+ (sqr a.x) (* ω2 (sqr a.y))) (mod* 2 a.x a.y)))
 
;; a ^ n
(define (Fp²-pow a n ω2)
  (Fp²-destruct* (a a.x a.y))
  (cond 
    ((= 0 n) (Fp² 1 0))
    ((= 1 n) a)
    ((= 2 n) (Fp²-mul a a ω2))
    ((even? n) (Fp²-square (Fp²-pow a (/ n 2) ω2) ω2))
    (else (Fp²-mul a (Fp²-pow a (sub1 n) ω2) ω2))))
 
;; x^2 ≡ n (mod p) ?
(define (mod-Cipolla n) 
  ;; check n is a square
  (unless (= 1 (mod-Legendre n)) (error 'Cipolla "~a not a square (mod ~a)" n (current-modulus)))
  ;; iterate until suitable 'a' found
  (define a (for/first ((t (in-range 2 (current-modulus))) ;; t = tentative a
                        #:when (= (sub1 (current-modulus))
                                  (mod-Legendre (- (* t t) n))))
              t))
  (define ω2 (- (* a a) n))
  ;; (Fp² a 1) = a + ω
  (define r (Fp²-pow (Fp² a 1) (/ (add1 (current-modulus)) 2) ω2))
  (define x (Fp²-x r))
  (unless (zero? (Fp²-y r)) (error 'Cipolla "ω has not vanished")) ;; hope that ω has vanished
  (unless (mod= n (* x x)) (error 'Cipolla "result check failed")) ;; checking the result
  (values x (mod- (current-modulus) x)))

(define (report-Cipolla n p)
  (with-handlers ((exn:fail? (λ (x) (eprintf "Caught error: ~s~%" (exn-message x)))))
    (define-values (r1 r2) (Cipolla n p))
    (printf "Roots of ~a are (~a,~a)  (mod ~a)~%" n  r1 r2 p)))

(module+ test
  (report-Cipolla 10 13)
  (report-Cipolla 56 101)
  (report-Cipolla 8218 10007)
  (report-Cipolla 8219 10007)
  (report-Cipolla 331575 1000003)
  (report-Cipolla 665165880 1000000007)
  (report-Cipolla 881398088036 1000000000039)
  (report-Cipolla 34035243914635549601583369544560650254325084643201
                  100000000000000000000000000000000000000000000000151))
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)
Caught error: "Cipolla: 8219 not a square (mod 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)

Raku

(formerly Perl 6)

Works with: Rakudo version 2016.10
Translation of: Sidef
#  Legendre operator (𝑛│𝑝)
sub infix:<│> (Int \𝑛, Int \𝑝 where 𝑝.is-prime && (𝑝 != 2)) {
    given 𝑛.expmod( (𝑝-1) div 2, 𝑝 ) {
        when 0  {  0 }
        when 1  {  1 }
        default { -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)"
}
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
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(n, p):
    a = Mod(n, p)
    out = []

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

    if not is_prime(p):
        conglst = []                                    #congruence list
        crtlst = []
        factors = []

        for k in list(factor(p)):
            factors.append(int(k[0]))

        for f in factors:
            conglst.append(cipollaAlgorithm(a, f))

        for i in Permutations([0, 1] * len(factors), len(factors)).list():
            for j in range(len(factors)):
                crtlst.append(int(conglst[ j ][ i[j] ]))

            out.append(crt(crtlst, factors))
            crtlst = []

        return sorted(out)

    if pow(p, 1, 4) == 3:
        temp = pow(a, int((p+1)/4), p)
        return [temp, p - temp]


    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)

    out.extend([x, p - x])

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

Scala

Imperative solution

object CipollasAlgorithm extends App {
  private val BIG = BigInt(10).pow(50) + BigInt(151)

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

  private def c(ns: String, ps: String): Triple = {
    val (n, p) = (BigInt(ns), if (ps.isEmpty) BIG else BigInt(ps))

    // Legendre symbol, returns 1, 0 or p - 1
    def ls(a: BigInt) = a.modPow((p - BigInt(1)) / BigInt(2), p)

    // multiplication in Fp2
    def mul(aa: Point, bb: Point, omega2: BigInt) =
      new Point((aa.x * bb.x + aa.y * bb.y * omega2) % p, (aa.x * bb.y + (bb.x * aa.y)) % p)

    // Step 0, validate arguments
    if (ls(n) != BigInt(1)) new Triple(0, 0, false)
    else {
      // Step 1, find a, omega2
      var (a, flag, omega2) = (BigInt(0), true, BigInt(0))
      while (flag) {
        omega2 = (a * a + p - n) % p
        if (ls(omega2) == p - BigInt(1)) flag = false else a = a + BigInt(1)
      }

      // Step 2, compute power
      var (nn, r, s) = ((p + BigInt(1) >> 1) % p, new Point(BigInt(1), 0), new Point(a, BigInt(1)))
      while (nn > 0) {
        if ((nn & BigInt(1)) == BigInt(1)) r = mul(r, s, omega2)
        s = mul(s, s, omega2)
        nn = nn >> 1
      }
      // Step 3, check x in Fp
      if (r.y != 0) new Triple(0, 0, false)
      else // Step 5, check x * x = n
      if ((r.x * r.x) % p != n) new Triple(0, 0, false)
      else new Triple(r.x, p - r.x, true) // Step 4, solutions
    }
  }

  private class Point(val x: BigInt, val y: BigInt)

  private class Triple(val x: BigInt, val y: BigInt, val b: Boolean) {
    override def toString: String = f"($x%s, $y%s, $b%s)"
  }

}
Output:

See it running in your browser by ScalaFiddle (JavaScript, non JVM) or by Scastie (JVM).

Sidef

Translation of: Go
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})"
    }
}
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

Simpler implementation: 

func cipolla(n, p) {

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

    var (a, ω) = (
        0..Inf -> lazy.map {|a|
            [a, submod(a*a, n, p)]
        }.first_by {|t|
            kronecker(t[1], p) == -1
        }...
    )

    var r = lift(Mod(Quadratic(a, 1, ω), p)**((p+1)>>1))

    r.b == 0 ? r.a : nil
}

Visual Basic .NET

Translation of: C#
Imports System.Numerics

Module Module1

    ReadOnly BIG = BigInteger.Pow(10, 50) + 151

    Function C(ns As String, ps As String) As Tuple(Of BigInteger, BigInteger, Boolean)
        Dim n = BigInteger.Parse(ns)
        Dim p = If(ps.Length > 0, BigInteger.Parse(ps), BIG)

        ' Legendre symbol. Returns 1, 0, or p-1
        Dim ls = Function(a0 As BigInteger) BigInteger.ModPow(a0, (p - 1) / 2, p)

        ' Step 0: validate arguments
        If ls(n) <> 1 Then
            Return Tuple.Create(BigInteger.Zero, BigInteger.Zero, False)
        End If

        ' Step 1: Find a, omega2
        Dim a = BigInteger.Zero
        Dim omega2 As BigInteger
        Do
            omega2 = (a * a + p - n) Mod p
            If ls(omega2) = p - 1 Then
                Exit Do
            End If
            a += 1
        Loop

        ' Multiplication in Fp2
        Dim mul = Function(aa As Tuple(Of BigInteger, BigInteger), bb As Tuple(Of BigInteger, BigInteger))
                      Return Tuple.Create((aa.Item1 * bb.Item1 + aa.Item2 * bb.Item2 * omega2) Mod p, (aa.Item1 * bb.Item2 + bb.Item1 * aa.Item2) Mod p)
                  End Function

        ' Step 2: Compute power
        Dim r = Tuple.Create(BigInteger.One, BigInteger.Zero)
        Dim s = Tuple.Create(a, BigInteger.One)
        Dim nn = ((p + 1) >> 1) Mod p
        While nn > 0
            If nn Mod 2 = 1 Then
                r = mul(r, s)
            End If
            s = mul(s, s)
            nn >>= 1
        End While

        ' Step 3: Check x in Fp
        If r.Item2 <> 0 Then
            Return Tuple.Create(BigInteger.Zero, BigInteger.Zero, False)
        End If

        ' Step 5: Check x * x = n
        If r.Item1 * r.Item1 Mod p <> n Then
            Return Tuple.Create(BigInteger.Zero, BigInteger.Zero, False)
        End If

        ' Step 4: Solutions
        Return Tuple.Create(r.Item1, p - r.Item1, True)
    End Function

    Sub Main()
        Console.WriteLine(C("10", "13"))
        Console.WriteLine(C("56", "101"))
        Console.WriteLine(C("8218", "10007"))
        Console.WriteLine(C("8219", "10007"))
        Console.WriteLine(C("331575", "1000003"))
        Console.WriteLine(C("665165880", "1000000007"))
        Console.WriteLine(C("881398088036", "1000000000039"))
        Console.WriteLine(C("34035243914635549601583369544560650254325084643201", ""))
    End Sub

End Module
Output:
(6, 7, True)
(37, 64, True)
(9872, 135, True)
(0, 0, False)
(855842, 144161, True)
(475131702, 524868305, True)
(791399408049, 208600591990, True)
(82563118828090362261378993957450213573687113690751, 17436881171909637738621006042549786426312886309400, True)

Wren

Translation of: Kotlin
Library: Wren-big
Library: Wren-dynamic
import "./big" for BigInt
import "./dynamic" for Tuple

var Point = Tuple.create("Point", ["x", "y"])
var bigBig = BigInt.ten.pow(50) + BigInt.new(151)

var c = Fn.new { |ns, ps|
    var n = BigInt.new(ns)
    var p = (ps != "") ? BigInt.new(ps) : bigBig

    // Legendre symbol, returns 1, 0 or p - 1
    var ls = Fn.new { |a| a.modPow((p - BigInt.one) / BigInt.two, p) }

    // Step 0, validate arguments
    if (ls.call(n) != BigInt.one) return [BigInt.zero, BigInt.zero, false]

    // Step 1, find a, omega2
    var a = BigInt.zero
    var omega2
    while (true) {
        omega2 = (a * a + p - n) % p
        if (ls.call(omega2) == p - BigInt.one) break
        a = a.inc
    }

    // multiplication in Fp2
    var mul = Fn.new { |aa, bb|
        return Point.new((aa.x * bb.x + aa.y * bb.y * omega2) % p,
                         (aa.x * bb.y + bb.x * aa.y) % p)
    }

    // Step 2, compute power
    var r = Point.new(BigInt.one, BigInt.zero)
    var s = Point.new(a, BigInt.one)
    var nn = ((p + BigInt.one) >> 1) % p
    while (nn > BigInt.zero) {
        if ((nn & BigInt.one) == BigInt.one) r = mul.call(r, s)
        s = mul.call(s, s)
        nn = nn >> 1
    }

    // Step 3, check x in Fp
    if (r.y != BigInt.zero) return [BigInt.zero, BigInt.zero, false]

    // Step 5, check x * x = n
    if (r.x * r.x % p != n) return [BigInt.zero, BigInt.zero, false]

    // Step 4, solutions
    return [r.x, p - r.x, true]
}

System.print(c.call("10", "13"))
System.print(c.call("56", "101"))
System.print(c.call("8218", "10007"))
System.print(c.call("8219", "10007"))
System.print(c.call("331575", "1000003"))
System.print(c.call("665165880", "1000000007"))
System.print(c.call("881398088036", "1000000000039"))
System.print(c.call("34035243914635549601583369544560650254325084643201", ""))
Output:
[6, 7, true]
[37, 64, true]
[9872, 135, true]
[0, 0, false]
[855842, 144161, true]
[475131702, 524868305, true]
[791399408049, 208600591990, true]
[82563118828090362261378993957450213573687113690751, 17436881171909637738621006042549786426312886309400, true]

zkl

Translation of: EchoLisp

Uses lib GMP (GNU MP Bignum Library). 

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);
}
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) }
}
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)
Retrieved from "https://rosettacode.org/wiki/Cipolla%27s_algorithm?oldid=352952"