A new algorithm for generic modular reduction

[Yves Gallot]

I think that the following algorithm is faster than the current generic modular reduction included in gwnum for fast exponentiation.
The method is to compute Montgomery modular multiplication using R = 2ⁿ - 1 (a Mersenne number).

Let z = x * y mod p, where p is an odd integer. R is a Mersenne number 2ⁿ - 1 >= p. q = 1/p (mod R).
REDC algorithm is

• t = x * y (a wide multiplication N x N => 2N).
• t_lo = t mod R, t_hi = t / R.
• r = (t_lo * q) mod R (a multiplication N x N => N, q is constant)
• r' = (r * p) / R (a wide multiplication N x N => 2N, p is constant)

The result is t_hi - r'.
The number of multiplications is two (N x N => 2N) and one (N x N => N), rather than three (N x N => 2N).

We can go further. We have to compute t_lo = t mod R, t_hi = t / R. t_lo is the cyclic convolution of x and y. t can be calculated using a cyclic convolution and a negacyclic convolution.
Let (a_1·x + a_0)·(b_1·x + b_0) mod x^2 - 1 = (a_0·b_1 + a_1·b_0)·(x - 1) + (a_0·b_1 + a_1·b_0) + (a_0·b_0 + a_1·b_1) = h·(x - 1) + l.
If C = (a_0 + a_1)·(b_0 + b_1) and N = (a_0 - a_1)·(b_0 - b_1) we have l = C and h = (C - N)/2.
If the 2 wide multiplications are evaluated using this method, we don't have to compute t mod R, t / R and (r * p) / R.
The number of operations is now three cyclic convolutions mod 2ⁿ - 1 and two negacyclic convolutions mod 2ⁿ + 1.
All operations are now N x N => N where N is the size of the tested number p.
The multiplication modulo any number is five times slower than the multiplication modulo a Mersenne number of the same size.

I implemented a 2-prp test using this method and 64-bit integers. Of course big numbers and fast transforms must be used in practice.

[pepi37]

Is this algo only for Mersenne number, or can be used on others numbers and bases?

[paulunderwood]

Originally Posted by pepi37

Is this algo only for Mersenne number, or can be used on others numbers and bases?

It is not for Mersenne or simple k*2^n+-1. It is for other numbers which do not use "special" mod. It is for "generic" mod.

George, is this any good for GWNUM? If so, how long before it is implemented?

[rogue]

This could be helpful for primorial and factorial searches, such as what PrimeGrid is doing.

[Citrix]

This might be faster:-
https://www.ams.org/journals/mcom/20...03-01543-6.pdf

Another faster approach for k*2^n+1 (especially larger k) would be to use the Chinese remainder FFT method using k, R=2^n-1 and coprime R'=2^(n+1)-1. The modulus step can be computed by addition.

(Slight modification of section 7 in https://www.ams.org/journals/mcom/19...-1185244-1.pdf)

[Yves Gallot]

Originally Posted by Citrix

This might be faster:-
https://www.ams.org/journals/mcom/20...03-01543-6.pdf

Yes, the last step of variation 2 is a single negacyclic convolution.
The number of operations is two cyclic convolutions mod 2ⁿ - 1 (a·b and (a·b)·N') and two negacyclic convolutions mod 2ⁿ + 1 (a·b and m·N).
Then it is four times slower than the multiplication modulo a Mersenne number.

[Prime95]

Originally Posted by paulunderwood

George, is this any good for GWNUM? If so, how long before it is implemented?

It's a long-shot.

GWNUM is designed to do lots of operations using one FFT. Here, you are mixing results from 2 FFTs (cyclic and negacyclic).

[Yves Gallot]

Originally Posted by Prime95

GWNUM is designed to do lots of operations using one FFT. Here, you are mixing results from 2 FFTs (cyclic and negacyclic).

McLaughlin variation 1 should be an easiest option and is as fast as variation 2.

If g = 1 we have R = 2ⁿ - 1, Q' = 2ⁿ⁺¹ - 1: two cyclic transforms of the same size with different weights.
Note that the last step of variation 1 is not correct: we have 0 <= t < 2Q' then t - N if t >= N is not enough for the modular reduction.
A common optimisation of REDC is to set N' = 1/N (mod R) (and not -1/N (mod R)). We have here:

1. m = a·b·N' mod R.
2. s = ((2m + 1)·N - (2a)·b) mod Q′.
3. If a·b + (m + 1)·N = s (mod 2) then t = s else t = s + Q′.
4. We have 0 <= t < 2Q', return t mod N.

[ixfd64]

Could this benefit programs like mfaktc and mfakto?

[Mysticial]

Originally Posted by Prime95

Here, you are mixing results from 2 FFTs (cyclic and negacyclic).

Off topic for this thread, but figured I'd mention that I do something similar in y-cruncher as a memory/storage reduction measure for the absolute largest of multiplications. It is slower due to needing to do additional passes over the dataset, but it almost halves the scratch memory/storage required since you halve the transform lengths.

CRT of cyclic (2^N - 1) and negacyclic (2^N + 1) is mathematically the same as the first radix 2 reduction of a normal FFT. The difference is where you perform the carryout.

[henryzz]

Originally Posted by Mysticial

Off topic for this thread, but figured I'd mention that I do something similar in y-cruncher as a memory/storage reduction measure for the absolute largest of multiplications. It is slower due to needing to do additional passes over the dataset, but it almost halves the scratch memory/storage required since you halve the transform lengths.

CRT of cyclic (2^N - 1) and negacyclic (2^N + 1) is mathematically the same as the first radix 2 reduction of a normal FFT. The difference is where you perform the carryout.

Presumably this ~doubles the FFT length that will fit in a L3 cache. This could be very useful for some cache size/FFT length combinations.