user
**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
assistant
(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

```