diff --git a/M2/Macaulay2/m2/galois.m2 b/M2/Macaulay2/m2/galois.m2 index bc28b70b98d..5486e972f4d 100644 --- a/M2/Macaulay2/m2/galois.m2 +++ b/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 -> ( diff --git a/M2/Macaulay2/m2/integers.m2 b/M2/Macaulay2/m2/integers.m2 index e014afbc673..8b777c81c33 100644 --- a/M2/Macaulay2/m2/integers.m2 +++ b/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 ----------------------------------------------------------------------------- diff --git a/M2/Macaulay2/m2/quotring.m2 b/M2/Macaulay2/m2/quotring.m2 index ab259a87c85..16675829f80 100644 --- a/M2/Macaulay2/m2/quotring.m2 +++ b/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)) diff --git a/M2/Macaulay2/m2/typicalvalues.m2 b/M2/Macaulay2/m2/typicalvalues.m2 index 82efde7c1d3..76256290af4 100644 --- a/M2/Macaulay2/m2/typicalvalues.m2 +++ b/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 " diff --git a/M2/Macaulay2/packages/Macaulay2Doc/operators.m2 b/M2/Macaulay2/packages/Macaulay2Doc/operators.m2 index dfe0df4a479..d2dac820895 100644 --- a/M2/Macaulay2/packages/Macaulay2Doc/operators.m2 +++ b/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, diff --git a/M2/Macaulay2/packages/Macaulay2Doc/ov_examples.m2 b/M2/Macaulay2/packages/Macaulay2Doc/ov_examples.m2 index c45c6e18119..de22a60dfe5 100644 --- a/M2/Macaulay2/packages/Macaulay2Doc/ov_examples.m2 +++ b/M2/Macaulay2/packages/Macaulay2Doc/ov_examples.m2 @@ -10,6 +10,7 @@ Node times power powermod + (sqrt, ZZ, ZZ) lcm gcd gcdCoefficients diff --git a/M2/Macaulay2/tests/normal/galois.m2 b/M2/Macaulay2/tests/normal/galois.m2 index f739628796e..ed2564834ec 100644 --- a/M2/Macaulay2/tests/normal/galois.m2 +++ b/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" diff --git a/M2/Macaulay2/tests/normal/integer.m2 b/M2/Macaulay2/tests/normal/integer.m2 index bcfa0c3f15e..cfe36f3487e 100644 --- a/M2/Macaulay2/tests/normal/integer.m2 +++ b/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" diff --git a/M2/Macaulay2/tests/normal/quotientring.m2 b/M2/Macaulay2/tests/normal/quotientring.m2 index 6d38d7d9d2f..30428ac62fe 100644 --- a/M2/Macaulay2/tests/normal/quotientring.m2 +++ b/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"