Module ch.randelshofer.fastdoubleparser
Package ch.randelshofer.fastdoubleparser

Class FftMultiplier

java.lang.Object
ch.randelshofer.fastdoubleparser.FftMultiplier
final class FftMultiplier extends Object
Provides methods for multiplying two BigIntegers using the FFT algorithm.

This code is based on bigint by Timothy Buktu.

References:

bigint, Copyright 2013 Timothy Buktu, 2-clause BSD License.
Note: We only use portions from this project, that have been marked with 2-clause BSD License in this file LICENSE.
github.com

  • Field Details

    • COS_0_25

      public static final double COS_0_25

    • SIN_0_25

      public static final double SIN_0_25

    • FFT_THRESHOLD

      private static final int FFT_THRESHOLD
      The threshold value for using floating point FFT multiplication. If the number of bits in each mag array is greater than the Toom-Cook threshold, and the number of bits in at least one of the mag arrays is greater than this threshold, then FFT multiplication will be used.
      See Also:

    • MAX_MAG_LENGTH

      private static final int MAX_MAG_LENGTH
      This constant limits mag.length of BigIntegers to the supported range.
      See Also:

    • ROOTS3_CACHE_SIZE

      private static final int ROOTS3_CACHE_SIZE
      for FFTs of length up to 3*2^19
      See Also:

    • ROOTS_CACHE2_SIZE

      private static final int ROOTS_CACHE2_SIZE
      for FFTs of length up to 2^19
      See Also:

    • TOOM_COOK_THRESHOLD

      private static final int TOOM_COOK_THRESHOLD
      The threshold value for using 3-way Toom-Cook multiplication.
      See Also:

    • ROOTS2_CACHE

      private static volatile FftMultiplier.ComplexVector[] ROOTS2_CACHE
      Sets of complex roots of unity. The set at index k contains 2^k elements representing all (2^(k+2))-th roots between 0 and pi/2. Used for FFT multiplication.

    • ROOTS3_CACHE

      private static volatile FftMultiplier.ComplexVector[] ROOTS3_CACHE
      Sets of complex roots of unity. The set at index k contains 3*2^k elements representing all (3*2^(k+2))-th roots between 0 and pi/2. Used for FFT multiplication.

  • Constructor Details

    • FftMultiplier

      FftMultiplier()

  • Method Details

    • bitsPerFftPoint

      static int bitsPerFftPoint(int bitLen)
      Returns the maximum number of bits that one double precision number can fit without causing the multiplication to be incorrect.
      Parameters:
      bitLen - length of this number in bits
      Returns:
      the maximum number of bits

    • calculateRootsOfUnity

      private static FftMultiplier.ComplexVector calculateRootsOfUnity(int n)
      Returns n-th complex roots of unity for the angles 0..pi/2, suitable for a transform of length n. They are used as twiddle factors and as weights for the right-angle transform. n must be 1 or an even number.

    • fft

      private static void fft(FftMultiplier.ComplexVector a, FftMultiplier.ComplexVector[] roots)
      Performs an FFT of length 2^n on the vector a. This is a decimation-in-frequency implementation.
      Parameters:
      a - input and output, must be a power of two in size
      roots - an array that contains one set of roots at indices log2(a.length), log2(a.length)-2, log2(a.length)-4, ... Each roots[s] must contain 2^s roots of unity such that roots[s][k] = e^(pi*k*i/(2*roots.length)), i.e., they must cover the first quadrant.

    • fft3

      private static void fft3(FftMultiplier.ComplexVector a0, FftMultiplier.ComplexVector a1, FftMultiplier.ComplexVector a2, int sign, double scale)
      Performs FFTs or IFFTs of size 3 on the vector (a0[i], a1[i], a2[i]) for each i. The output is placed back into a0, a1, and a2.
      Parameters:
      a0 - inputs / outputs for the first FFT coefficient
      a1 - inputs / outputs for the second FFT coefficient
      a2 - inputs / outputs for the third FFT coefficient
      sign - 1 for a forward FFT, -1 for an inverse FFT
      scale - 1 for a forward FFT, 1/3 for an inverse FFT

    • fftMixedRadix

      private static void fftMixedRadix(FftMultiplier.ComplexVector a, FftMultiplier.ComplexVector[] roots2, FftMultiplier.ComplexVector roots3)
      Performs an FFT of length 3*2^n on the vector a. Uses the 4-step algorithm to decompose the 3*2^n FFT into 2^n FFTs of length 3 and 3 FFTs of length 2^n. See https://www.nas.nasa.gov/assets/pdf/techreports/1989/rnr-89-004.pdf
      Parameters:
      a - input and output, must be 3*2^n in size for some n>=2
      roots2 - an array that contains one set of roots at indices log2(a.length/3), log2(a.length/3)-2, log2(a.length/3)-4, ... Each roots[s] must contain 2^s roots of unity such that roots[s][k] = e^(pi*k*i/(2*roots.length)), i.e., they must cover the first quadrant.
      roots3 - must be the same length as a and contain roots of unity such that roots[k] = e^(pi*k*i/(2*roots3.length)), i.e., they need to cover the first quadrant.

    • fromFftVector

      static BigInteger fromFftVector(FftMultiplier.ComplexVector fftVec, int signum, int bitsPerFftPoint)

    • getRootsOfUnity2

      private static FftMultiplier.ComplexVector[] getRootsOfUnity2(int logN)
      Returns sets of complex roots of unity. For k=logN, logN-2, logN-4, ..., the return value contains all k-th roots between 0 and pi/2.
      Parameters:
      logN - for a transform of length 2^logN

    • getRootsOfUnity3

      private static FftMultiplier.ComplexVector getRootsOfUnity3(int logN)
      Returns sets of complex roots of unity. For k=logN, logN-2, logN-4, ..., the return value contains all k-th roots between 0 and pi/2.
      Parameters:
      logN - for a transform of length 3*2^logN

    • ifft

      private static void ifft(FftMultiplier.ComplexVector a, FftMultiplier.ComplexVector[] roots)
      Performs an inverse FFT of length 2^n on the vector a. This is a decimation-in-time implementation.
      Parameters:
      a - input and output, must be a power of two in size
      roots - an array that contains one set of roots at indices log2(a.length), log2(a.length)-2, log2(a.length)-4, ... Each roots[s] must contain 2^s roots of unity such that roots[s][k] = e^(pi*k*i/(2*roots.length)), i.e., they must cover the first quadrant.

    • ifftMixedRadix

      private static void ifftMixedRadix(FftMultiplier.ComplexVector a, FftMultiplier.ComplexVector[] roots2, FftMultiplier.ComplexVector roots3)
      Performs an inverse FFT of length 3*2^n on the vector a. Uses the 4-step algorithm to decompose the 3*2^n FFT into 2^n FFTs of length 3 and 3 FFTs of length 2^n. See https://www.nas.nasa.gov/assets/pdf/techreports/1989/rnr-89-004.pdf
      Parameters:
      a - input and output, must be 3*2^n in size for some n>=2
      roots2 - an array that contains one set of roots at indices log2(a.length/3), log2(a.length/3)-2, log2(a.length/3)-4, ... Each roots[s] must contain 2^s roots of unity such that roots[s][k] = e^(pi*k*i/(2*roots.length)), i.e., they must cover the first quadrant.
      roots3 - must be the same length as a and contain roots of unity such that roots[k] = e^(pi*k*i/(2*roots3.length)), i.e., they need to cover the first quadrant.

    • multiply

      static BigInteger multiply(BigInteger a, BigInteger b)
      Returns a BigInteger whose value is (a * b).
      Parameters:
      a - value a
      b - value b
      Returns:
      this * val

    • multiplyFft

      static BigInteger multiplyFft(BigInteger a, BigInteger b)
      Multiplies two BigIntegers using a floating-point FFT.

      Floating-point math is inaccurate; to ensure the output of the FFT and IFFT rounds to the correct result for every input, the provably safe FFT error bounds from "Rapid Multiplication Modulo The Sum And Difference of Highly Composite Numbers" by Colin Percival, pg. 392 (fft.pdf) are used, the vector is "balanced" before the FFT, and accurate twiddle factors are used.

      This implementation incorporates several features compared to the standard FFT algorithm (Cooley Tukey FFT algorithm):

      • It uses a variant called right-angle convolution which weights the vector before the transform. The benefit of the right-angle convolution is that when multiplying two numbers of length n, an FFT of length n suffices whereas a regular FFT needs length 2n. This is because the right-angle convolution places half of the result in the real part and the other half in the imaginary part. See: Discrete Weighted Transforms And Large-Integer Arithmetic by Richard Crandall and Barry Fagin.
      • FFTs of length 3*2^n are supported in addition to 2^n.
      • Radix-4 butterflies; see https://www.nxp.com/docs/en/application-note/AN3666.pdf
      • Bernstein's conjugate twiddle trick for a small speed gain at the expense of (further) reordering the output of the FFT which is not a problem because it is reordered back in the IFFT.
      • Roots of unity are cached
      FFT vectors are stored as arrays of primitive doubles (two array elements are needed for representing one complex number). Storing them as arrays of primitive doubles instead of as MutableComplex objects is memory efficient, but in some cases below ~10^6 decimal digits, it hurts speed because it requires additional copying. Ideally this would be implemented using value types when they become available.
      Parameters:
      a - value a
      b - value b
      Returns:
      a*b

    • square

      static BigInteger square(BigInteger a)
      Returns a BigInteger whose value is (this<sup>2</sup>).
      Returns:
      this<sup>2</sup>

    • squareFft

      static BigInteger squareFft(BigInteger a)

    • toFftVector

      static FftMultiplier.ComplexVector toFftVector(byte[] mag, int fftLen, int bitsPerFftPoint)
      Converts this BigInteger into an array of complex numbers suitable for an FFT. Populates the real parts and sets the imaginary parts to zero.