Add sqrt(ZZ, ZZ) method (and sqrt for ZZ/p and GF p)

M2/Macaulay2/m2/galois.m2

@@ -228,6 +228,7 @@ GF(Ring) := GaloisField => opts -> (S) -> (
      F.cache = new CacheTable;
      F / F := (x,y) -> if y == 0 then error "division by zero" else x // y;
      F % F := (x,y) -> if y == 0 then x else 0_F;
+     sqrt F := x -> promote(tonelliShanks(lift(x, ZZ), F.order), F);
      F)
 
 random GaloisField := opts -> F -> (

M2/Macaulay2/m2/integers.m2

@@ -109,6 +109,36 @@ changeBase(String, ZZ)     := ZZ     => changeBase0
 changeBase(String, ZZ, ZZ) := String => (s, oldbase, newbase) -> (
     changeBase(changeBase(s, oldbase), newbase))
 
+-- Tonelli-Shanks algorithm (used by various sqrt methods)
+legendreSymbol = (n, p) -> powermod(n, (p - 1)//2, p)
+tonelliShanks = (n, p) -> (
+    if n % p == 0 then return 0;
+    if n % p == 1 then return 1;
+    if not isPrime p then error "expected a prime";
+    if legendreSymbol(n, p) != 1 then error "not a quadratic residue";
+    (q, s) := (p - 1, 0);
+    while even q do (
+	q //= 2;
+	s += 1);
+    z := 2;
+    while legendreSymbol(z, p) != p - 1 do z += 1;
+    m := s;
+    c := powermod(z, q, p);
+    t := powermod(n, q, p);
+    r := powermod(n, (q + 1)//2, p);
+    while true do (
+	if t == 0 then return 0;
+	if t == 1 then return (
+	    if 2*r > p then p - r
+	    else r);
+	i := 0;
+	while powermod(t, 2^i, p) != 1 do i += 1;
+	b := powermod(c, 2^(m - i - 1), p);
+	m = i;
+	c = powermod(b, 2, p);
+	t = (t * c) % p;
+	r = (r * b) % p))
+
 -----------------------------------------------------------------------------
 -- AtomicInt
 -----------------------------------------------------------------------------

M2/Macaulay2/m2/quotring.m2

@@ -103,6 +103,7 @@ ZZp Ideal := opts -> (I) -> (
 	  expression S := x -> expression rawToInteger raw x;
 	  fraction(S,S) := S / S := (x,y) -> if y === 0_S then error "division by zero" else x//y;
 	  S.frac = S;		  -- ZZ/n with n PRIME!
+	  sqrt S := x -> promote(tonelliShanks(lift(x, ZZ), n), S);
       savedQuotients#(typ, n) = S;
       lift(S,QQ) := opts -> liftZZmodQQ;
 	  S))

M2/Macaulay2/m2/typicalvalues.m2

@@ -319,6 +319,8 @@ taylor (log1p, (x,n) -> (
     s
     ))
 
+-- now that sqrt is a method function, we can finally install this
+sqrt(ZZ, ZZ) := tonelliShanks
 
 -- Local Variables:
 -- compile-command: "make -C $M2BUILDDIR/Macaulay2/m2 "

M2/Macaulay2/packages/Macaulay2Doc/operators.m2

@@ -107,6 +107,46 @@ RRi => { "an interval containing the square roots of the points of ", TT "I" }
      ///
      }
 
+doc ///
+  Key
+    (sqrt, ZZ, ZZ)
+  Headline
+    modular square root
+  Usage
+    sqrt(n, p)
+  Inputs
+    n:ZZ
+    p:ZZ -- a prime number
+  Outputs
+    :ZZ -- the square root of @VAR "n"@ modulo @VAR "p"@
+  Description
+    Text
+      This returns a solution of the equation $x^2\equiv n\pmod p$, provided
+      that @VAR "n"@ is a quadratic residue modulo @VAR "p"@.  If it is not,
+      then an error is raised.
+    Example
+      powermod(5, (41 - 1) // 2, 41) -- Euler's criterion
+      sqrt(5, 41)
+      powermod(13, 2, 41)
+    Text
+      Every quadratic residue modulo an odd prime has two square roots.  The
+      second square root may be obtained by subtracting the first from
+      @VAR "p"@.
+    Example
+      powermod(28, 2, 41)
+    Text
+      This method may also be used for finding square roots in finite prime
+      fields.
+    Example
+      sqrt(5_(ZZ/41))
+      sqrt(5_(GF 41))
+  References
+    Shanks, Daniel. "Five number-theoretic algorithms." @EM "Proceedings of
+    the Second Manitoba Conference on Numerical Mathematics (Winnipeg)"@, 1973.
+  SeeAlso
+    powermod
+///
+
 document {
      Key => {gcd,
 	 1:gcd,

M2/Macaulay2/packages/Macaulay2Doc/ov_examples.m2

@@ -10,6 +10,7 @@ Node
     times
     power
     powermod
+    (sqrt, ZZ, ZZ)
     lcm
     gcd
     gcdCoefficients

M2/Macaulay2/tests/normal/galois.m2

@@ -32,6 +32,9 @@ k = GF(2,3,Variable=>a)
 
 assert( length degree id_(k^1) == 0 )
 
+k = GF 41
+assert(sqrt(5_k)^2 == 5_k)
+
 end
 -- Local Variables:
 -- compile-command: "make -C $M2BUILDDIR/Macaulay2/packages/Macaulay2Doc/test galois.out"

M2/Macaulay2/tests/normal/integer.m2

@@ -79,6 +79,8 @@ assert try powermod(2, -1, 4) else true -- used to crash M2
 for b from 2 to 62 do assert(changeBase("10101", b, b) == "10101")
 assert(changeBase("0xdeadbeef", 0) == 0xdeadbeef)
 
+assert(powermod(sqrt(5, 41), 2, 41) == 5)
+
 end
 -- Local Variables:
 -- compile-command: "make -C $M2BUILDDIR/Macaulay2/packages/Macaulay2Doc/test integer.out"

M2/Macaulay2/tests/normal/quotientring.m2

@@ -14,6 +14,9 @@ assert( a_R^2 == 0 )
 R = QQ[a]/a^2[x][y]/y^3
 assert( a_R^2 == 0 )
 
+R = ZZ/41
+assert(sqrt(5_R)^2 == 5_R)
+
 end
 -- Local Variables:
 -- compile-command: "make -C $M2BUILDDIR/Macaulay2/packages/Macaulay2Doc/test quotientring.out"