Montgomery method in Prime Numbers: ACP

[SPWorley]

I love the book Prime Numbers: a Computational Perspective.
But I'm going wild trying to get Montgomery reduction to work properly.

The problem is that this book uses a different style of explanation with wikipedia, and it's hard to sync them up, and I'm starting to doubt PN:ACP.

On page 449, there's algorithm 9.2.5.
The book's errata pdf says that the line "x=cd" should be "y=cd".
But I'm pretty sure that's wrong too. I think it should remain
"x=cd" and then the next line changed to:
"z=(x+N((xN') mod R))/R" as in theorem 9.2.1.

And on that same page, I also think the ladder algorithm (9.2.6) is wrong. Its notation makes it confusing if it's doing the digits in the reverse order or not.
But as a standard binary ladder, shouldn't the lines:

p=M(p,p)
if (y_i==1) p=M(p,x);

be swapped and changed to

if (y_i==1) p=M(p,x)
x=M(x,x)

?


I know I'm an idiot because I can't figure this out. But if you have a copy of PN:ACP on your shelf, could you look at page 449 and tell me what flavor of idiot I am.. one who can't understand the simple algorithm, or one who can't make a simple Montgomery reduction work.

(There's a third type of idiot, the one who uses wikipedia for math references, but let's not go there...)



A side question, does anyone have some example C or C++ code for using Montgomery reduction to compute x^y mod n which is designed for the simplest case of integer x,y,n all <2^32?

[Robert Holmes]

You appear to be correct about the first question --- it should be x = cd, then the transformation of 9.2.1 applied to x gives you y.

The powering algorithm is correct --- it is the same that is presented in 9.3.1. I believe you were thinking of the dual version of the square-and-multiply algorithm, also presented in 9.3.2.

As for sample code, check http://islab.oregonstate.edu/papers/j37acmon.pdf --- it presents various approaches for actual implementations of Montgomery multiplication.

[SPWorley]

Originally Posted by Robert Holmes

You appear to be correct about the first question --- it should be x = cd, then the transformation of 9.2.1 applied to x gives you y.

The powering algorithm is correct --- it is the same that is presented in 9.3.1. I believe you were thinking of the dual version of the square-and-multiply algorithm, also presented in 9.3.2.

Thanks for the verification!

But that laddering code is awkward. It doesn't match the 9.3.1 code, it unwraps one extra pass for the high bit, which makes it uselessly compute a full and potentially expensive Montgomery product to get the result p=1*1 (in Montgomery adjusted numbers) on that first pass. That made me doubt my understanding of it, especially as my code wasn't working.

Luckily, in the meantime, my code is now working... now to optimize.

Thanks again!

[SPWorley]

A followup question on Montgomery reduction.

Using a radix R=2^s is most common of course.
But the Montgomery method only works with N that are coprime to R..
so for the R=2^s case, N must be odd.

What's the solution for even N? You could break N into the product a power of 2 and an odd remainder, perform operations on both (the power of 2 mod version is especially easy), then use the Chinese Remainder theorem to combine the results.. but that's really inelegant and even with the efficiency of mod 2^32 computes, it's still computationally wasteful.

So what do you do when you need to perform many x^y mod N computes where N is even, and the Montgomery method fails for R=2^s?

[Robert Holmes]

Use the CRT, not unlike RSA decryption. If your modulus is even, it can be factored into N = p * 2^d. As p and 2^d are obviously coprime, one can perform the exponentiation mod p and 2^d separately and then paste it together using the CRT:

z0 = x^e mod p
z1 = x^e mod 2^d
z = z0 + p * ((z1-z0)* (p^-1) mod 2^d)

Needless to say, the arithmetic mod 2^d is trivial to perform in a computer. The inversion can be performed easily by Newton's iteration (see C&P Exercise 9.12)

Alternatively, you can use an alternative method, e.g. Barrett's reduction.

[Ken_g6]

SPWorley, you may also be interested in Alex Kruppa's Montgomery implementation, and what I did with it in this thread.

(This has nothing to do with even modulii.)