Module ch.randelshofer.fastdoubleparser

Package ch.randelshofer.fastdoubleparser

Class FftMultiplier

java.lang.Object

ch.randelshofer.fastdoubleparser.FftMultiplier

final class FftMultiplier
extends java.lang.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

Nested Class Summary

Nested Classes

Modifier and Type	Class	Description
(package private) static class	FftMultiplier.ComplexVector
(package private) static class	FftMultiplier.MutableComplex

Field Summary

Fields

Modifier and Type	Field	Description
static double	COS_0_25
private static int	FFT_THRESHOLD	The threshold value for using floating point FFT multiplication.
private static int	MAX_MAG_LENGTH	This constant limits mag.length of BigIntegers to the supported range.
private static int	ROOTS_CACHE2_SIZE	for FFTs of length up to 2^19
private static FftMultiplier.ComplexVector[]	ROOTS2_CACHE	Sets of complex roots of unity.
private static FftMultiplier.ComplexVector[]	ROOTS3_CACHE	Sets of complex roots of unity.
private static int	ROOTS3_CACHE_SIZE	for FFTs of length up to 3*2^19
static double	SIN_0_25
private static int	TOOM_COOK_THRESHOLD	The threshold value for using 3-way Toom-Cook multiplication.

Constructor Summary

Constructors

Constructor	Description
FftMultiplier()

Method Summary

Modifier and Type	Method	Description
(package private) 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.
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.
private static void	fft(FftMultiplier.ComplexVector a, FftMultiplier.ComplexVector[] roots)	Performs an FFT of length 2^n on the vector a.
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.
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.
(package private) static java.math.BigInteger	fromFftVector(FftMultiplier.ComplexVector fftVec, int signum, int bitsPerFftPoint)
private static FftMultiplier.ComplexVector[]	getRootsOfUnity2(int logN)	Returns sets of complex roots of unity.
private static FftMultiplier.ComplexVector	getRootsOfUnity3(int logN)	Returns sets of complex roots of unity.
private static void	ifft(FftMultiplier.ComplexVector a, FftMultiplier.ComplexVector[] roots)	Performs an inverse FFT of length 2^n on the vector a.
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.
(package private) static java.math.BigInteger	multiply(java.math.BigInteger a, java.math.BigInteger b)	Returns a BigInteger whose value is (a * b).
(package private) static java.math.BigInteger	multiplyFft(java.math.BigInteger a, java.math.BigInteger b)	Multiplies two BigIntegers using a floating-point FFT.
(package private) static java.math.BigInteger	square(java.math.BigInteger a)	Returns a BigInteger whose value is (this<sup>2</sup>).
(package private) static java.math.BigInteger	squareFft(java.math.BigInteger a)
(package private) static FftMultiplier.ComplexVector	toFftVector(byte[] mag, int fftLen, int bitsPerFftPoint)	Converts this BigInteger into an array of complex numbers suitable for an FFT.

Methods inherited from class java.lang.Object

clone, equals, finalize, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait

Field Detail

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:

Constant Field Values

MAX_MAG_LENGTH

private static final int MAX_MAG_LENGTH

This constant limits mag.length of BigIntegers to the supported range.

See Also:

Constant Field Values

ROOTS3_CACHE_SIZE

private static final int ROOTS3_CACHE_SIZE

for FFTs of length up to 3*2^19

See Also:

Constant Field Values

ROOTS_CACHE2_SIZE

private static final int ROOTS_CACHE2_SIZE

for FFTs of length up to 2^19

See Also:

Constant Field Values

TOOM_COOK_THRESHOLD

private static final int TOOM_COOK_THRESHOLD

The threshold value for using 3-way Toom-Cook multiplication.

See Also:

Constant Field Values

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 Detail

FftMultiplier

FftMultiplier()

Method Detail

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 java.math.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 java.math.BigInteger multiply(java.math.BigInteger a,
                                     java.math.BigInteger b)

Returns a BigInteger whose value is (a * b).

Parameters:

a - value a

b - value b

Returns:

this * val

multiplyFft

static java.math.BigInteger multiplyFft(java.math.BigInteger a,
                                        java.math.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 java.math.BigInteger square(java.math.BigInteger a)

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

Returns:

this<sup>2</sup>

squareFft

static java.math.BigInteger squareFft(java.math.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.