Irrational base discrete weighted transform

Short description: Variant of fast Fourier 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) in the early 1990s using Mathematica.

The IBDWT is used in the Great Internet Mersenne Prime Search's client Prime95 to perform FFT multiplication, as well as in other programs implementing Lucas–Lehmer test, such as CUDALucas and Glucas.

Algorithm

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

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

This approach avoids the need for zero-padding the arrays and performs the multiplication modulo $M_p$ directly. The algorithm to compute the product $xy(\mathrm{mod}{M}_p)$ is as follows:

1. Choose a run length (signal-length) $N<p$.
2. Establish a variable base representation for the numbers. For example, $x=\sum _{j=0}^{N-1}{x}_j{2}^{\lceil pj/N\rceil }$. Each term is usually between 16 and 20 bits if using double-precision terms.
3. Define a weight-signal $a$ where each component $a_j=2^{\lceil pj/N\rceil -pj/N}$, approximated by floats in the interval [1, 2).
4. Compute the forward DWT for both numbers: ${𝒳}\leftarrow DW{T}_{(N,a)}x$ and ${𝒴}\leftarrow DW{T}_{(N,a)}y$. This is practically computed using a standard DFT (like an FFT) as $DW{T}_{(N,a)}x\triangleq DF{T}_{(N)}(ax)$.
5. Perform a component-wise product of the transformed arrays: ${𝒵}\leftarrow {𝒳}{𝒴}$.
6. Compute the inverse DWT: $z\leftarrow DW{T}_{(N,a)}^{-1}{𝒵}$. This is computed as $z=a^{-1}DF{T}_{(N)}^{-1}({𝒵})$.
7. Round the resulting components to the nearest integer: $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).
8. Adjust the resulting digits $\{z_n\}$ to restore the variable-base radix representation. This step handles carries and borrows. Single-step partial carrying is sufficient.

References

• 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