Mersenne trial division implementation

[mathPuzzles]

I was reading the great wiki on the stages of Mersenne testing. The initial filtering step is trial division. After some sieving and modular rejection, you still end up with various primes which are worthwhile to explicitly test as factors. And the wiki says "The remaining potential factors are tested using the [traditional generic] modular powering algorithm above." Is this correct? Is the best Mersenne hunting software literally using the classic power-mod evaluation, or is are there any enhancements since you know you're taking the modulus of a number with the very simple form of \(2^n\)?

Put another way, Is there a especially good way to efficiently evaluate \(2^n \mod m\) as opposed to a generic \(a^n \mod m\)?

I can imagine that you could speed up the square-and-multiply loop with an especially efficient Brauer \(2^k\)-ary ladder since you can build any \(2^b\) for \(b < \log_2 m\) for free by simply left shifting a single bit by \(b\). That should save about 1/3 of the modular multiplies, but you still need to do the \(\log_2 n\) squarings.

If \(n\) is large enough it would eventually be useful to use Montgomery multiplication instead, but you'd lose the elegant "shift a bit" ladder and again it'd be no different for \(2^n\) than a generic \(a^n\).

Are there any other speedup tricks or is it really like the wiki says and it's still really just the classic square-and-multiply modular powering algorithm?

[science_man_88]

Originally Posted by mathPuzzles

I was reading the great wiki on the stages of Mersenne testing. The initial filtering step is trial division. After some sieving and modular rejection, you still end up with various primes which are worthwhile to explicitly test as factors. And the wiki says "The remaining potential factors are tested using the [traditional generic] modular powering algorithm above." Is this correct? Is the best Mersenne hunting software literally using the classic power-mod evaluation, or is are there any enhancements since you know you're taking the modulus of a number with the very simple form of \(2^n\)?

Put another way, Is there a especially good way to efficiently evaluate \(2^n \mod m\) as opposed to a generic \(a^n \mod m\)?

I can imagine that you could speed up the square-and-multiply loop with an especially efficient Brauer \(2^k\)-ary ladder since you can build any \(2^b\) for \(b < \log_2 m\) for free by simply left shifting a single bit by \(b\). That should save about 1/3 of the modular multiplies, but you still need to do the \(\log_2 n\) squarings.

If \(n\) is large enough it would eventually be useful to use Montgomery multiplication instead, but you'd lose the elegant "shift a bit" ladder and again it'd be no different for \(2^n\) than a generic \(a^n\).

Are there any other speedup tricks or is it really like the wiki says and it's still really just the classic square-and-multiply modular powering algorithm?

they could use factoring of p-1 and then use the factors of p-1 to calculate it for 2^(p-1) and then for 2^p.

http://primes.utm.edu/glossary/xpage...entiation.html also suggest they could go left to right binary exponentiation but that's a time versus space tradeoff.

[mathPuzzles]

An extremely similar question comes up when doing the precompute necessary for Montgomery multiplication of any number, since you almost always use a computational base like \(R=2^{64i}\) for some convienient \(i\). The one-time precompute needs to find the modular inverse of R in the base in question, usually using the extended Euclidian algorithm. So, is there any special method or optimization to computing \((2^n)^{-1}\mod m?\)

[mathPuzzles]

Originally Posted by science_man_88;456967[url

http://primes.utm.edu/glossary/xpage/BinaryExponentiation.html[/url] also suggest they could go left to right binary exponentiation but that's a time versus space tradeoff.

Yep, to get the full advantage of the (free) ladder precompute, you'd work left to right in chunks of \(b\) bits, since those are the sizes where you can get the (free) evaluation by shifting. And then of course the actual multiplication of that ladder is really just shifting the second multiplicand, so it's not even a mulmod, it's just a mod. The bulk of the evaluation cost will be all the squarings, which as far as I can see are unavoidable.

[ewmayer]

Originally Posted by mathPuzzles

I was reading the great wiki on the stages of Mersenne testing. The initial filtering step is trial division. After some sieving and modular rejection, you still end up with various primes

No - lots of composites, too.

Originally Posted by mathPuzzles

An extremely similar question comes up when doing the precompute necessary for Montgomery multiplication of any number, since you almost always use a computational base like \(R=2^{64i}\) for some convienient \(i\). The one-time precompute needs to find the modular inverse of R in the base in question, usually using the extended Euclidian algorithm.

No eGCD needed.

[mathPuzzles]

Originally Posted by ewmayer

No eGCD needed.

That's very interesting and exactly what I'm trying to learn! So what is this method to compute \((2^n)^{-1}\mod m\) without eGCD? Perhaps some variant in the flavor of the binary GCD?

[ewmayer]

Originally Posted by mathPuzzles

That's very interesting and exactly what I'm trying to learn! So what is this method to compute \((2^n)^{-1}\mod m\) without eGCD? Perhaps some variant in the flavor of the binary GCD?

The mod-inverse needed by the Montgomery mod is not the one you cite, but rather an inverse of m (mod R), where R is the arithmetic radix convenient for hardware arithmetic, e.g. 2^64 (or some suitable power thereof which is > m) on a 64-bit CPU. To get that mod-inverse one needs just the fast Newtonian-iterative mod-inverse described in Montgomery's paper. (See also here for some improved-initial-guess tweaks.) The ensuing binary modpow then computes 2^n (mod m) times some power of R which scaling can be 'unwound' in log log n time if the true mod is needed - which it isn't for e.g. Mersenne-mod TF work.

[mathPuzzles]

ewmayer, you're right, I swapped the R and modulus value in my question. But your answer is also part of what I'm trying to learn! Now I have to study your paper. It's looks like part of the solution to my efficiency puzzle. I appreciate it!

[ewmayer]

Originally Posted by mathPuzzles

ewmayer, you're right, I swapped the R and modulus value in my question. But your answer is also part of what I'm trying to learn! Now I have to study your paper. It's looks like part of the solution to my efficiency puzzle. I appreciate it!

Algorithm E in the paper has an efficient modpow based on Montgomery-mul, but I encourage you to try to slog though the preceding material, as well. :) I included several worked 64-bit-modulus examples; hopefully you will find those useful.