function newton(pol, root, derinv, p, m) assert Evaluate(pol, GF(p)!root) eq 0; k := 1; P := p; pold := Derivative(pol); assert Evaluate(pold, GF(p)!root)*derinv eq 1; vprintf User1: "root = %o\nderinv = %o\n\n", root, derinv; while 2*k lt m do k *:= 2; P := P^2; root := (root - Evaluate(pol, root)*derinv) mod P; derinv := (derinv*(2 - Evaluate(pold, root)*derinv)) mod P; vprintf User1: "root = %o\nderinv = %o\n\n", root, derinv; end while; root := (root - Evaluate(pol, root)*derinv) mod (P*p^(m-k)); vprintf User1: "root = %o\n\n", root; return root; end function; function nthroot(a, n) assert IsOdd(n); assert a gt 0; k := Valuation(a, 2); if k mod n ne 0 then return false, _; end if; a := ExactQuotient(a, 2^k); m := 1 + Floor(Log(2,a)/n); b := newton(P.1^n-a, 1, 1, 2, m) where P := PolynomialRing(Integers()); if b^n eq a then return true, b*2^(k div n); else return false, _; end if; end function; function sqrt(a) assert a gt 0; k := Valuation(a, 2); if IsOdd(k) then return false, _; end if; a := ExactQuotient(a, 2^k); m := Floor(Log(2,a)/2) - 1; if a mod 8 ne 1 then return false, _; end if; A := a div 8; if m le 0 then return true, (2*A+1)*2^(k div 2); end if; B := newton(Polynomial([-A,1,2]), A mod 2, 1, 2, m); b := 4*B + 1; if b^2 eq a then return true, b*2^(k div 2); end if; b := 2^(m+2) - b; if b^2 eq a then return true, b*2^(k div 2); end if; return false, _; end function; procedure checksqrt(len) a := 2*Random([2^(len-2), 2^(len-1)-1])+1; printf "len = %o:\n", len; t := Cputime(); a2 := a^2; t := Cputime(t); printf "a^2 ==> time = %6o\n", t; t := Cputime(); flag, b := sqrt(a2); assert flag and b eq a; t := Cputime(t); printf "sqrt ==> time = %6o\n", t; t := Cputime(); flag, b := IsSquare(a2); assert flag and b eq a; t := Cputime(t); printf "IsSquare ==> time = %6o\n", t; end procedure; // try checksqrt(10^7): // sqrt takes about 11 sec, // IsSquare sometimes 2.7 and sometimes 25 sec ?! procedure checknthroot(len, n) assert IsOdd(n); a := 2*Random([2^(len-2), 2^(len-1)-1])+1; printf "len = %o:\n", len; t := Cputime(); a2 := a^n; t := Cputime(t); printf "a^%o ==> time = %6o\n", n, t; t := Cputime(); flag, b := nthroot(a2, n); assert flag and b eq a; t := Cputime(t); printf "nthroot ==> time = %6o\n", t; t := Cputime(); flag, b := IsPower(a2, n); assert flag and b eq a; t := Cputime(t); printf "IsPower ==> time = %6o\n", t; end procedure;