Modular Square Roots

Date: 03/24/2006 at 08:14:14
From: Philip
Subject: modular square roots

I am trying to find an algorithm for finding modular square roots 
where the modulus is either a square or an odd composite.

I've heard about something called Hensel lifting but I don't 
understand.  The example I found in a paper by Victor Shoup doesn't
seem to work but its probably my failure to understand the modular
arithmetic involved.

I do know that in the case of a square modulus (p^2) there can only 
be a square root if there is a square root mod p but I can't get 
much further than that.  In the other case I know that given a
factorization of p1e1.p2e2..... one can somehow use the Chinese
Remainder Theorem.

Date: 03/24/2006 at 17:42:31
From: Doctor Vogler
Subject: Re: modular square roots

Hi Philip,

Thanks for writing to Dr. Math.  The first thing to do is factor your
modulus, because the problem is much easier with smaller moduluses. 
Then you divide your problem into three steps:

  1) For each prime factor p, find a modular square root mod p.

  2) For each prime power p^e (if e > 1), use the modular square root
mod p and Hensel lifting to get a modular square root mod p^e.

  3) Use the Chinese Remainder Theorem to put the roots together.

Note that there are two square roots mod p for any prime (additive
inverses of one another), so you have two choices in step (1) for each
prime.  If you have k distinct prime factors, then you have will 2^k
total choices, and there will be 2^k final square roots.

But actually, the prime 2 behaves a little differently (which is
related to the fact that squares have exponent 2); see below.

Step 1 is, in some sense, the hardest part, because about the only
thing you can do is a big search.  The simplest way is to try all of
the numbers from 1 to p/2, square each one, and stop when you get the
right square.  (Why can you stop at p/2 instead of going all the way
to p?)  Knowing whether you started with a square mod p is a much
easier problem due to the nice technique known as quadratic
reciprocity.  But that's another issue entirely.

Hensel lifting is fairly simple.  In one sense, the idea is to use
Newton's method to get a better result.  That is, if p is an odd
prime, and

  r^2 = n (mod p),

then you can find the root mod p^2 by changing your first
"approximation" r to

  r - (r^2 - n)/(2r) (mod p^2).

It turns out that this method will double the exponent at each step,
so next you can change the new approximation r to

  r - (r^2 - n)/(2r) (mod p^4),

and so on.  Why does this work?  Consider that you already have

  r^2 = n (mod p),

and you want to compute the root mod p^2.  Reducing mod p, you find
that the new root has to be r mod p.  (Can you explain why?)  So we
need to find some number s so that

  (r + ps)^2 = n (mod p^2).

Expanding, we get

  r^2 + 2ps + p^2*s^2 = n (mod p^2)

and therefore

  r^2 + 2ps = n (mod p^2)

or

  2ps = n - r^2 (mod p^2).

We already know that n - r^2 is divisible by p, so we divide by p, and
also by 2 mod p^2, and we get

  s = (n - r^2)/(2p) (mod p),

which is our answer!  If you do the substitution and simplify, you'll
find that

  r + ps = r - (r^2 - n)/(2r) (mod p^2).

That's the idea of Hensel lifting.

For p = 2, you have to start your Hensel lifting from p^4 = 16 instead
of p, since you have problems dividing by 2.  In any case, the only
odd numbers that have square roots mod 8 are 1.  At least mod 16,
there are two choices (1 and 9).  But then if you have

  r^2 = n (mod 2^k)

and you want to find s so that

  (r + s*2^k)^2 = n (mod 2^(2k))

You expand and get

  r^2 + s*2^(k+1) + s*2^(2k) = n (mod 2^(2k))

and the third term is zero mod 2^(2k), so you solve for s and get

  s = (n - r^2)/2^(k+1) mod (2^(k-1)),

which is the same formula as before (so you can still use Newton's
method), but you need to start with n - r^2 divisible by 2^(k+1) or
else you won't be able to divide by 2^(k+1), and then you only get s
mod 2^(k-1), so you don't quite double the exponent at each iteration
of Newton's method.  You get from a root mod 2^(k+1) to a root mod
2^(2k-1).  Or, said differently, you get from a root mod 2^k to a root
mod 2^(2k-3).  That's why it's not helpful to start with k <= 3.

If you have any questions about this or need more help, please write
back and show me what you have been able to do, and I will try to
offer further suggestions.

- Doctor Vogler, The Math Forum
  http://mathforum.org/dr.math/ 

Date: 03/24/2006 at 23:21:54
From: Philip
Subject: Thank you (modular square roots)

Excellent - thank you.  You've made it seem very simple.  It seems to
work fine except in the division stage at each iteration we seem to
need to be working with real numbers (floating point) or the result is
off by one.  The idea of mixing reals with modular arithmetic is a bit
alien to me.

Have I got this right??

Date: 03/25/2006 at 09:46:55
From: Doctor Vogler
Subject: Re: modular square roots

Hi Philip,

You're right to be concerned; in fact, reals and modular arithmetic
don't mix at all.

When I write

  a/b (mod n),

what I mean is a times the multiplicative inverse of b mod n.  For
example,

  1/2 = 4 (mod 7)

because 4 is the multiplicative inverse of 2 mod 7;

  2*4 = 1 (mod 7).

And 7 is its own multiplicative inverse mod 12 (check this), so

  3/7 = 3*7 = 21 = 9 = -3 (mod 12).

The usual rules of division apply, such as

  a/b + c/d = (ad + bc)/(bd)

and

  (ak)/(bk) = a/b

and so on.  But you have to be cautious about dividing by zero, since
there are more zeros mod n than there are in the reals.  In fact, you
have to be careful dividing by anything that isn't relatively prime to
the modulus.  You see, if the denominator isn't relatively prime to
the modulus, then it has no multiplicative inverse.  But we have
the rules:

  bx = ba (mod bm)

implies

  x = a (mod m),

which is a kind of dividing by b, but you have to divide the modulus
by b, too.  You'll see that happening when I divided by p.

Does that make sense?  Does it clear things up?

- Doctor Vogler, The Math Forum
  http://mathforum.org/dr.math/ 

Date: 04/21/2006 at 15:34:56
From: Doctor Vogler
Subject: Re: modular square roots

Hi Philip,

A little while ago, I said that to find a modular square root mod n,
you need to factor n and then do three steps:

  1) For each prime factor p, find a modular square root mod p.

  2) For each prime power p^e (if e > 1), use the modular square root
mod p and Hensel lifting to get a modular square root mod p^e.

  3) Use the Chinese Remainder Theorem to put the roots together.

Then I said that step 1 is, in some sense, the hardest part, because
about the only thing you can do is a big search.  I have since learned
that I was incorrect.  You can do better than a big search.

In fact, when p = 4k + 3 for some integer k, this is easy, since

  g^(k+1) mod p

is a square root of g mod p.

When p = 8k + 5 for some integer k, then this is only a little bit
harder, since you can take

  c = (2g)^k
  i = (2g)*(c^2) = (2g)^(2k+1)

and then

  g*c*(i - 1) mod p

is a square root of g mod p.

The case that is left is when p = 8k + 1 for some integer k, and this
case is harder, although there are still methods for doing this.  One
method is to do operations in the ring

  (Z/pZ)[x]/(x^2 - g)

or in the ring

  (Z/pZ)[x]/(x^2 + g).

Another method is a kind of partial discrete logarithm problem.  You
find a quadratic nonresidue mod p (half of them are, so they're not
hard to find), and raise it to the odd part of p-1 (that's p-1 after
you divide it by as many twos as you can).  You also raise your square
g to the odd part of p-1, and then you try to find the power of the
first that equals the second.

- Doctor Vogler, The Math Forum
  http://mathforum.org/dr.math/ 

Date: 04/22/2006 at 07:57:27
From: Philip
Subject: Thank you (modular square roots)

Many thanks and bless you.  This and your earlier advice have been
very helpful.

I needed to sort this out in connection with my project to implement
the multiple polynomial quadratic sieve (in python) which I have now
completed and which (so far as I can see) performs a great deal better
than the only other python implementation I know of in the public domain.

However, I'm still a little confused on moduli which are powers of 
two (I use precalculated roots for my current purposes).

Date: 04/23/2006 at 18:01:03
From: Doctor Vogler
Subject: Re: modular square roots

Hi Philip,

So you have a number n, and you want to find a square root mod 2^k. 
That is, you want to find a number r such that

  r^2 = n (mod 2^k).

In fact, you'd probably like to find all of them.  The tricky part of
this is that x^2 mod 2^k depends not on x mod 2^k but only on x mod
2^(k-1), when k > 1.  For example,

  x^2 = (x + 8)^2 mod 16.

That's because if

  x = y + z*2^(k-1)

for some integers y and z, then

  x^2 = y^2 + y*z*2^k + 2^(2k-2)
      = y^2 (mod 2^k).

So you'll recall that if n is a square mod p^k for some odd prime p,
then there are two square roots, r and -r, mod p^k.  But when p = 2,
then there are *four* square roots mod 2^k, namely

  r
  -r
  r + 2^(k-1)
  -r + 2^(k-1)

which is really just *two* square roots, r and -r, mod 2^(k-1).  Does
that make sense?

So we're actually only going to find the square roots mod 2^(k-1). 
There are two of those.

First of all, it would be nice to assume that n is odd, so let's
divide n by 2 until it is odd, so that

  n = n' * 2^m.

Then if m is odd (and less than k), then n has no square roots mod
2^k.  If m is even, then all we need to do is find the square root of
n' mod 2^(k-m) and multiply the root by 2^(m/2).  And n' is odd, so
what follows applies to n'.

So now we assume that n is odd.  Then if k = 1, we have exactly one
square root, namely r = 1 mod 2.  If k = 2, then we still have one
square root, namely r = 1 mod 2, but only if n = 1 mod 4.  If n = 3
mod 4, then it has no square roots.  The rest of the comments apply
for k >= 3.  In this case, 1 is the only odd square mod 8, which means
that n has no square roots mod 2^k unless n = 1 mod 8.  But if n = 1
mod 8, then there will be exactly two square roots mod 2^(k-1), which
will be additive inverses of one another, as I described above.

So how do we do Hensel lifting?  If we have an approximation x that
satisfies

  x^2 = n (mod 2^m)

for some m >= 3 (recall then that we only need to know x mod 2^(m-1)),
then we can change x to

  y = x - (x^2 - n)/(2x).

This is the formula from Newton's Method, possibly more familiar as

  y = (1/2)(x + n/x)

so that y is the average of x and n/x, a common way to approximate
square roots.  I'll find it easier to use the form

  y = (x^2 + n)/(2x)

in many cases.  In any case, it is important that

  (x^2 - n)/2

is an integer (because 2^m divides x^2 - n and m >= 3), so that we can
compute y mod 2^(2m-3).  To do this, we take the integer

  (x^2 - n)/2

and then multiply it by the multiplicative inverse of x mod 2^(2m-3),
and then subtract the product from x.  In a similar way, we could
multiply the integer

  (x^2 + n)/2

by the multiplicative inverse of x mod 2^(2m-3).  The result is

  y = (x^2 + n)/(2x) (mod 2^(2m-3)).

This determines

  y^2 = (x^2 + n)^2/(2x)^2 (mod 2^(2m-2))

and you can check that

  y^2 - n = (x^4 + 2nx^2 + n^2)/(4x^2) - n (mod 2^(2m-2))
          = (x^4 + 2nx^2 + n^2 - 4nx^2)/(4x^2) (mod 2^(2m-2))
          = (x^4 - 2nx^2 + n^2)/(4x^2) (mod 2^(2m-2))
          = (x^2 - n)^2/(4x^2) (mod 2^(2m-2)).

Since x^2 - n is divisible by 2^m, its square is divisible by 2^(2m).
After dividing by 4, it's still divisible by 2^(2m-2).  Since x is
odd, dividing by x^2 doesn't change that.  So that means that

  y^2 = n (mod 2^(2m-2)).

In a similar way, you can show that

  y = x (mod 2^(m-1)).

So this is how you use that:  You start with m = 3 and any odd number
for x (since x^2 = 1 = n mod 8 for all odd numbers x).  You'll get
different answers for x = 1 mod 4 and x = 3 mod 4, but they will be
negatives of one another, so you don't really need to do both.  You
compute y mod 8 and find that y^2 = n mod 16 and y = x mod 4.  The
second condition really only guarantees that you get two different
answers for x = 1 mod 4 and x = 3 mod 4, which, as I pointed out, will
be negatives of one another mod 2^(2m-3).  Then you change x to y and
m to 4 and recompute y mod 32 and get y^2 = n mod 64.  Then you change
x to y and m to 6 and recompute y mod 2^9 and get y^2 = n mod 2^10. 
And so on.  Each time you compute y mod 2^(m-1) and get y^2 = n mod
2^m, the m doubles minus 2 (changes to 2m-2).  (When p is odd, m
doubled each time, so this is marginally slower.)

For example, let's compute the square roots of 41 mod 2^k for various
numbers k.  So we have n = 41, and we start with

  x = 1 mod 4,

which gives

  1^2 = 41 (mod 8).

Then we change x to

  y = (x^2 + 41)/(2x)
    = (42)/2
    = 21
    = 5 (mod 8)

which has

  5^2 = 41 (mod 16).

Then we change x to

  y = (x^2 + 41)/(2x)
    = (25 + 41)/(2*5)
    = (66/2)/5
    = 33/5
    = 33*13 (mod 32)
    = 13 (mod 32)

which has

  13^2 = 41 (mod 64).

And so on.

- Doctor Vogler, The Math Forum
  http://mathforum.org/dr.math/ 

Date: 04/24/2006 at 01:59:44
From: Philip
Subject: Thank you (modular square roots)

Dear Dr. Vogler -

Thank you once again for a very clear and detailed explanation and for
all your help to date.

Philip