Cipolla's algorithm
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:
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
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
<lang perl6># Legendre operator (𝑛│𝑝) sub infix:<│> (Int \𝑛, Int \𝑝 where (𝑝.is-prime and ?(𝑝 != 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, $a); for 0 .. * { $a = $_; $ω2 = ($a² - 𝑛) % 𝑝; last if ($ω2│𝑝) < 0; }
# 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 ...^ 0 { $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); if $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 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
Sage provides functions to construct a finite number field, modulo p*p, with polynomial x^2-ω², and this is precisely what we need. <lang sage> def Cipolla(n,p) :
if not is_prime(p) : print ("❌ %d is not a prime" % p) return False if legendre_symbol(n,p) != 1 : print ("❌ %d is not a square modulo %d" % (n,p)) return False for a in range (2,n) : if legendre_symbol(a*a-n,p) == -1 : break omega2 = a*a - n Q.<x> = PolynomialRing(GF(p)) R.<x> = GF(p*p,modulus = x^2-omega2) r = (a + x) ^ ((p+1)/2) return r, p-r
</lang>
- Output:
sage: Cipolla( 10,13) (6, 7) sage: Cipolla(56,101) (37, 64) sage: Cipolla(8218,10007) (9872, 135) sage: Cipolla (331575, 1000003) (855842, 144161) sage: Cipolla (8219,10007) ❌ 8219 is not a square modulo 10007 False
Sidef
<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
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)