mersenneforum.org  

Go Back   mersenneforum.org > Great Internet Mersenne Prime Search > Math

Reply
 
Thread Tools
Old 2009-08-17, 00:34   #1
SPWorley
 
Jul 2009

3110 Posts
Default Montgomery method in Prime Numbers: ACP

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?

Last fiddled with by SPWorley on 2009-08-17 at 00:36
SPWorley is offline   Reply With Quote
Old 2009-08-17, 03:44   #2
Robert Holmes
 
Robert Holmes's Avatar
 
Oct 2007

11010102 Posts
Default

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.
Robert Holmes is offline   Reply With Quote
Old 2009-08-17, 04:43   #3
SPWorley
 
Jul 2009

378 Posts
Default

Quote:
Originally Posted by Robert Holmes View Post
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 is offline   Reply With Quote
Old 2009-08-17, 17:41   #4
SPWorley
 
Jul 2009

3110 Posts
Default

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?
SPWorley is offline   Reply With Quote
Old 2009-08-17, 19:10   #5
Robert Holmes
 
Robert Holmes's Avatar
 
Oct 2007

2·53 Posts
Default

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.

Last fiddled with by Robert Holmes on 2009-08-17 at 19:13
Robert Holmes is offline   Reply With Quote
Old 2009-08-18, 17:27   #6
Ken_g6
 
Ken_g6's Avatar
 
Jan 2005
Caught in a sieve

1100010112 Posts
Default

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.)

Last fiddled with by Ken_g6 on 2009-08-18 at 17:28
Ken_g6 is offline   Reply With Quote
Reply

Thread Tools


Similar Threads
Thread Thread Starter Forum Replies Last Post
Base method of montgomery polyselect implementation csg Homework Help 2 2017-12-13 07:04
Reverse process sum of odd numbers from x to y. What is the fastest method? Alberico Lepore Alberico Lepore 12 2017-09-05 14:57
New method of finding large prime numbers georgelouis@mac Math 41 2011-01-25 21:06
Brent-Montgomery-te Riele numbers FactorEyes Factoring 23 2008-02-22 00:36
Elliptic Curve Method factoring - Fermat numbers philmoore Math 131 2006-12-18 06:27

All times are UTC. The time now is 08:12.

Sun May 9 08:12:39 UTC 2021 up 31 days, 2:53, 0 users, load averages: 2.45, 2.60, 2.83

Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2021, Jelsoft Enterprises Ltd.

This forum has received and complied with 0 (zero) government requests for information.

Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.2 or any later version published by the Free Software Foundation.
A copy of the license is included in the FAQ.