Irrational base discrete weighted transform

Template:Short description In mathematics, the irrational base discrete weighted transform (IBDWT) is a variant of the fast Fourier transform using an irrational base; it was developed by Richard Crandall (Reed College), Barry Fagin (Dartmouth College) and Joshua Doenias (NeXT Software)<ref>Template:Cite journal</ref> in the early 1990s using Mathematica. It implies a fast, practical implementation of large-number modular multiplication on modern computers, at asymptotically 2× faster than non-modular FFT multiplication.<ref>Template:Cite web</ref><ref name=GrangerScott15/>

It is most notably used in the Great Internet Mersenne Prime Search.

Algorithm

The IBDWT method, as applied to the Lucas-Lehmer test for Mersenne primes (which requires repeated squaring modulo a Mersenne number <math>M_p = 2^p - 1</math>), is based on four key elements developed by Crandall and Fagin:<ref name="thall" />

• Balanced-radix representations: Allows digits to be signed (e.g., in the range <math>-W/2 \le x_j < W/2</math>), which reduces error bounds.<ref name="thall" />
• Variable-base digit representations: Allows each digit <math>x_j</math> to have its own base <math>W_j</math>.<ref name="thall" />
• Weighted cyclic convolutions: The multiplication is performed using a Discrete Weighted Transform (DWT).<ref name="thall" />
• Irrational numeric bases: The base used for the transform is irrational.<ref name="thall" />

This approach avoids the need for zero-padding the arrays and performs the multiplication modulo <math>M_p</math> directly.<ref name="thall" /> The algorithm to compute the product <math>xy \pmod{M_p}</math> is as follows:<ref name="thall" />

1. Choose a run length (signal-length) <math>N < p</math>.<ref name="thall" />
2. Establish a variable base representation for the numbers. For example, <math display=inline>x = \sum_{j=0}^{N-1} x_j 2^{\lceil pj/N \rceil}</math>.<ref name="thall" /> Each term is usually between 16 and 20 bits if using double-precision terms. The maximum size of each term is determined empirically,Template:Efn at least before the contribution of Percival 2003.
3. Define a weight-signal <math>a</math> where each component <math>a_j = 2^{\lceil pj/N \rceil - pj/N}</math>, approximated by floats in the interval [1, 2).<ref name="thall" />
4. Compute the forward DWT for both numbers: <math>\mathcal{X} \leftarrow DWT_{(N,a)}x</math> and <math>\mathcal{Y} \leftarrow DWT_{(N,a)}y</math>. This is practically computed using a standard DFT (like an FFT) as <math>DWT_{(N,a)}x \triangleq DFT_{(N)}(ax)</math>.<ref name="thall" />
5. Perform a component-wise product of the transformed arrays: <math>\mathcal{Z} \leftarrow \mathcal{X}\mathcal{Y}</math>.<ref name="thall" />
6. Compute the inverse DWT: <math>z \leftarrow DWT_{(N,a)}^{-1}\mathcal{Z}</math>. This is computed as <math>z = a^{-1} DFT_{(N)}^{-1}(\mathcal{Z})</math>.<ref name="thall" />
7. Round the resulting components to the nearest integer: <math>z \leftarrow \text{round}(z)</math>, optionally checking the roundoff error is no greater than 0.4 (greater indicates too many integer bits stuffed into each term).<ref name="thall" />
8. Adjust the resulting digits <math>\{z_n\}</math> to restore the variable-base radix representation. This step handles carries and borrows. Single-step partial carrying is sufficient.<ref name="thall" />

Applications

Double-precision IBDWT is used in the Great Internet Mersenne Prime Search's x86 client Prime95 to perform modular multiplication in the Lucas–Lehmer test and Fermat primarily tests. The prime95 IBDWT library gwnum is also used in programs such as PrimeGrid's LLR2 and PRST. It is chosen because x86 CPUs since Pentium 4 have so much double-precision floating-point computing power that it is much faster to multiply numbers using IBDWT than to do the so using a more straightforward integer FFT (NTT).

Double-precision IBDWT has also been ported to other CPU architectures in the form of Glucas. It has also been ported to GPUs in the form of CUDALucas, GPUowl, and PRPLL.<ref name=thall>Template:Cite web</ref>

IBDWT can also be done using integer arithmetic modulo 2⁶⁴-2³²+1, a number theoretic transform. This approach was first demonstrated by Nick Craig-Wood in ARMPrime.<ref>Template:Cite web</ref> This too has been ported to GPUs, providing an alternative for consumer GPUs with weak double-precision computing power but acceptable 32-bit integer power, especially Nvidia models from the 2020s boasting "1:1" or "1:2" 32-bit integer multiplication speed but "1:64" double-precision speed relative to 32-bit floating-point.<ref>Template:Cite web</ref>

Derived methods

Granger and Scott demonstrated using IBDWT-inspired "GRP (generalized repunit prime) multiplication" to accelerate eliptic curve cryptography over F(2⁵²¹-1), the P-521. This is a Karatsuba-like technique featuring a cyclic convolution similar to IBDWT.<ref name=GrangerScott15>Template:Multiref2</ref>

Further reading

• Richard Crandall, Barry Fagin: Discrete weighted transforms and large-integer arithmetic, Mathematics of Computation 62, 205, 305-324, January 1994 (PDF file)
• Richard Crandall: Topics in Advanced Scientific Computation, TELOS/Springer-Verlag
• Template:Cite journal