Irrational base discrete weighted transform

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)^[1] 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.^[2][3]

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 ${\displaystyle M_{p}=2^{p}-1}$), is based on four key elements developed by Crandall and Fagin:^[4]

• Balanced-radix representations: Allows digits to be signed (e.g., in the range ${\displaystyle -W/2\leq x_{j}<W/2}$), which reduces error bounds.^[4]
• Variable-base digit representations: Allows each digit ${\displaystyle x_{j}}$ to have its own base ${\displaystyle W_{j}}$.^[4]
• Weighted cyclic convolutions: The multiplication is performed using a Discrete Weighted Transform (DWT).^[4]
• Irrational numeric bases: The base used for the transform is irrational.^[4]

This approach avoids the need for zero-padding the arrays and performs the multiplication modulo ${\displaystyle M_{p}}$ directly.^[4] The algorithm to compute the product ${\displaystyle xy{\pmod {M_{p}}}}$ is as follows:^[4]

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

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.^[4]

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.^[5] 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.^[6]

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.^[3]

Notes

1. For examples of empirical formulae and tables, see gpuowl and Prime95 "gwnum.c".