newtonhensel.mws
Integer root extraction and perfectpower detection via
p
adic NewtonHensel lifting
Author: Carl Devore <devore@math.udel.edu>
8 January 2003
If you use this worksheet, please email me at
mailto://devore@math.udel.edu , to discuss it, or just to say "hi". If you find this worksheet useful, please consider supporting the work of a poor graduate student who spends thousands of dollars per year on books and computer equipment by sending me a small item from my Amazon.com wishlist
http://www.amazon.com/exec/obidos/registry/4SWHEIGN0FN9/ref=wl_s_3/10450624841726354 . You just have to select the items, and they are sent directly to me.
Keywords:
Newton's method, Hensel's lemma,
p
adic, perfectpower detection, root extraction
Purpose:
This worksheet can serve as an introduction to
p
adic methods, as a review of Newton's method, or as a discussion of the number theoretic problem of perfectpower detection.
The problems:
The computation problem (root extraction):
Given positive integers
n
and
r
we want to find if possible a positive integer
x
such that
. If no such
x
exists, then we do not care about getting an approximate solution.
The decision problem (perfectpower detection):
Given a positive integer
n
, we wish to detect whether there exists integers
x
and
r
such that
without necessarily finding a suitable
x
or
r.
However, in the code that follows, we will find the
x
and
r
.
Newton's method:
Recall that if
is a differentiable function, then the Newton iteration for finding an approximate solution to
is
.
As an example, consider the simple problem of computing a reciprocal. While it may seem trivial to use Newton's method for this
, this example is of fundamental importance and is used extensively in what follows, because it turns out that this Newton iteration can be expressed in a form that does not use division.
So that solving
is equivalent to computing
. The Newton iteration is
> 
G:= unapply(x  f(x)/diff(f(x),x), x);

But watch what happens when we expand that.
There is no division  just multiplication and subtraction. That idea is used in some computers to implement division, and it has many other uses. (See Knuth, 1998). In the general Newton's method, we need to divide by the derivative. We can avoid the division by performing a preliminary Newton step for reciprocating the derivative. While that is not very useful for floatingpoint computation, it is useful in modular arithmetic, where the concept of division is slightly different.
The
p
adic case:
It is a great surprise that the Newton iteration also works in modular arithmetic, even though the notion of derivative must be radically altered, and the usual graphical explanation of Newton's method in terms of finding the
x
intercept of the tangent line makes no sense. Indeed, in a sense, it works even better in modular arithmetic because it can give exact answers rather than approximations and because the number of iterations required to get an exact answer can be determined precisely.
Theorem 1 (NewtonHensel lifting of the modular reciprocal):
If
(
p
not necessarily prime) (thus
), then
.
Proof:
Since
p
divides
,
divides
. Thus
, and the conclusion follows.
Indeed, the idea can be extended to the Newton iteration for the roots of polynomials:
Theorem 2 (NewtonHensel lifting for modular polynomial roots):
Let
be a polynomial with integer coefficients,
, and
a unit modulo
p
. Then
.
Proof:
Consider the Taylor expansion of
f
:
.
When
x
and the coefficients of
f
are integers, the factorial denominators will always cancel with factors arising from differentiation, so the coefficients of the
y
's in the above expansion are integers. Substituting
and
and taking the result modulo
we get
.
Now substitute
. The division by
p
makes sense because we are assuming
, hence
is divisible by
p
. We get
.
So how does this help us solve problems over the integers? Once the modulus is larger than the possible root over the integers, we check whether the solution also works over the integers. Note that the modulus is squared in each step, so it grows very fast. This is analogous to the wellknown rule of thumb about Newton's method  that the number of digits of accuracy doubles with each iteration. In the
p
adic case, that rule of thumb is actually precise: the number of
p
adic digits of accuracy exactly doubles with each iteration.
In what follows, sometimes we need to factor a small prime out of an integer. The following procedure does that.
> 
FactorOut:= proc(n::posint, p::posint)
# return e and o such that n = p^e*o and o mod p <> 0.
local e,o,x;
o:= n;
for e from 0 do
if irem(o,p,'x') <> 0 then break end if;
o:= x
end do;
e,o
end proc:

To compute the square root of
n
, we use the Newton iteration
. The inverse cannot exist if
p
is even, thus we use
. That means we have to factor out 3's from
n,
and we use an initial solution
because 2 is not a square modulo 3. The inverse is updated via its own Newton iteration from Theorem 1.
> 
Isqrt:= proc(N::posint)
# Returns ``(sqrt(N))^2 if N is a perfect square, false otherwise,
# via 3adic Newton lifting.
local p, x, `/2/x`, n, e3;
# The 3adic algorithm requires that N not be a multiple of 3.
# So factor out all 3's.
e3, n:= FactorOut(N,3);
if irem(e3,2,'e3') <> 0 then return false end if;
# No square is congruent to 2 mod 3.
if irem(n,3)=2 then return false end if;
# Find the solution mod p=3^(2^k) until p > sqrt(n).
p:= 9; x:= 1; `/2/x`:= 2; # 1/2 mod 3 = 2
to ilog2(1+iquo(ilog2(n), 3)) do
# NewtonHensel lift for the solution of x^2 = n.
x:= mods(x  (x^2n)*`/2/x`, p);
# NewtonHensel lift for the next "denominator" `/2/x` in the line above.
`/2/x`:= mods(2*`/2/x`*(1x*`/2/x`), p);
p:= p^2
end do;
x:= mods(x  (x^2n)*`/2/x`, p);
if x^2 = n then ``(3^e3*abs(x))^2 else false end if
end proc:

Amazingly enough, the above procedure compares favorably timewise with Maple's builtin
isqrt.
But the builtin command also computes the nearest integer approximation if the square root is not an integer.
Test on a number of 10,000 digits.
Maple's
iroot
is faster though.
> 
time(iroot(n,2,'exact'));

> 
evalb(op(op(1,Isqrt(n)))^2=n);

To compute the
r
th root of
n
for odd
r
, we use the Newton iteration
. The
is computed separately to avoid repetition. A separate Newton iteration is used to update the inverse. We use
and factor out all 2's from
n
. If
r
is even, we just use Isqrt as many times as needed before using the algorithm for odd exponents.
> 
Iroot:= proc(N::posint, R::posint)
# Returns ``(N^(1/R))^R if N is an Rth power, false otherwise,
# via 2adic Newton lifting.
local n, r, p, next_p, x, `r1`, `x^(r1)`, `/r/x^(r1)`, e, e2, e3, ee2;
if R=1 then return ``(N) end if;
n:= N;
# First, do as much of the computation as possible with Isqrt.
ee2, r:= FactorOut(R,2);
if ee2 > 0 then
e3, n:= FactorOut(n,3);
if irem(e3, R, 'e3') <> 0 then return false end if;
to ee2 do
n:= Isqrt(n);
if n=false then
return false
else
n:= op(op(1,n))
end if
end do
else
e3:= 0
end if;
# Find a solution mod p = 2^(2^k) until p > n^(1/r).
e2, n:= FactorOut(n,2);
if irem(e2, r, 'e2') <> 0 then return false end if;
p:= 4; x:= 1; `/r/x^(r1)`:= 1; `r1`:= r1; `x^(r1)`:= 1;
to ilog2(iquo(ilog2(n),r)+1) do
next_p:= p^2;
# NewtonHensel lift for the solution of x^r = n.
x:= mods(x  (x*`x^(r1)`n)*`/r/x^(r1)`, p);
# NewtonHensel lift for the denominator.
`x^(r1)`:= mods(x^`r1`, next_p);
`/r/x^(r1)`:= mods(`/r/x^(r1)`*(2r*`x^(r1)`*`/r/x^(r1)`), p);
p:= next_p
end do;
x:= mods(x  (x*`x^(r1)`n)*`/r/x^(r1)`, p);
if x^r = n then ``(2^e2*3^e3*x)^(r*2^ee2) else false end if
end proc:

Unlike Isqrt, Iroot is not competitive timewise with Maple's
iroot
, which is not builtin but uses extendedprecision floatingpoint computation. Nonetheless, Iroot is still fast. Note that
p
is always a power of 2. All of the mod operations could be done by bit manipulations, but Maple does not represent large integers in binary, nor does Maple have an adequate bittruncation command.
> 
time(iroot(n,3, 'exact'));

> 
evalb(op(op(1,Iroot(n,3)))^3=n);

To detect whether
n
is a perfect power, we simply check every possible prime exponent up to a logarithmic bound.
> 
Ispower:= proc(n::posint)
# Returns two items.
# The first is a true/false whether n is perfect power.
# The second is a partial factorization n = ``(2)^e*``(3)^f*``(a)^g
# where gcd(a,6) = 1 and g is maximal.
local e2,e3,o,B,b,factor,save_factor;
if n=1 then return true, ``(1)^2 end if;
e2,o:= FactorOut(n,2);
e3,o:= FactorOut(o,3);
if o=1 then
return evalb(igcd(e2,e3)>1), ``(2)^e2*``(3)^e3
end if;
save_factor:= ``(o);
# Upper bound on the remaining power (the minimal base is
# now 5).
B:= iquo(ilog2(o), 2);
b:= 2;
while b <= B do
factor:= Iroot(o,b);
if factor<>false then
if igcd(e2,e3,b) = b then
# Use recursive call (at minimal extra cost) to determine
# if the factor can be expressed to a higher power.
return true, ``(2)^e2*``(3)^e3*Ispower(op(op(1,factor)))[2]^b
else
save_factor:= factor
end if
end if;
b:= nextprime(b)
end do;
false
,``(2)^e2*``(3)^e3
*`if`(
save_factor::`^`
,Ispower(op(op(1,save_factor)))[2]^op(2,save_factor)
,save_factor
)
end proc:

If one was purely interested in perfectpower detection, then there is no need to return the partial factorization. The recursive call is only made to refine this factorization; whether or not
n
is a perfect power has already been determined with certainty when the recusive call is made. Some small extra efficiencies could be made to the above by sieving for the possible prime exponents.
See the paper by Berstein.
References:
1. Eric Bach and Jeffrey Shallit
Algorithmic Number Theory: Volume 1: Efficient Algorithms
(MIT Press, 1996).
2. Dario Bini and Victor Pan
Polynomial and Matrix Computations: Volume 1: Fundamental Algorithms
(Birkhauser, 1994).
3. Daniel J. Berstein "Detecting perfect powers in essentially linear time"
Mathematics of Computation
(67:223, pp. 125383, July 1998).
http://cr.yp.to/papers/powersams.pdf
4. Donald Ervin Knuth
The Art of Computer Programming: Volume 2: Seminumerical Algorithms
(3rd ed., Addison Wesley Longman, 1998).
5. Joachim von zur Gathen and Jurgen Gerhard
Modern Computer Algebra
(Cambridge University Press, 1999).