Fast Multiplication

About

Fast integer multiplication is an algorithm able to perform the multiplication of two n-digit integers using O(n log ⁡n) operations.

Modular multiplication of polynomials is the basis of this method.
An intricate method is often used: the product can be rewritten as a cyclic convolution and is evaluated using Fast Fourier Transforms. It is unnecessary and the purpose of this project is to demystify the fast multiplication for large integers.

Fast modular multiplication of polynomials is a divide-and-conquer algorithm, where factorization and the Chinese remainder theorem are applied recursively. This method is described here: starting from a simple implementation of the mathematical relation, a list of elementary programs illustrates how it can be extended and improved.

fastMul_rec.cpp

The recursion is applied to compute P(x)² mod x⁸ − 7, where P is a polynomial of degree 8−1. The field is the complex numbers.
P(x) mod x^m − r^m and P(x) mod x^m + r^m are calculated from P(x) mod x^2m − r^2m. The last step is Lagrange interpolation and we have P(x) mod x − r = P(r). P(r)² is a scalar. Finally, P(x)² mod x^2m − r^2m is calculated from P(x)² mod x^m − r^m and P(x)² mod x^m + r^m with the Chinese remainder theorem for polynomials.

fastMul_rec_235.cpp

The algorithm is extended to P(x)² mod x^{22 · 32 · 52} − 7. We have x^3m − r^3m = (x^m − r^m) · (x^m − jr^m) · (x^m − j²r^m) · , where j is a primitive root of x³ − 1. The equivalent relation is applied to x^5m − r^5m.

fastMul_rec_GF.cpp

Roots are calculated over a prime finite field of order p. There are two conditions:

• If P(x)² mod xⁿ − r is evaluated then a n^th root of r must be an element in Z/pZ.
• If the coefficients of P(x) are non-negative integers smaller than b (P(b) is the base-b representation of an integer) then P(x)² mod xⁿ − r over the natural numbers can be retrieved from the result over Z/pZ if r > 0 and p > n · (b − 1)². If r < 0 then outputs can be negative integers: the condition is p > 2 n · (b − 1)² to determine the sign of each coefficient.

Here n = 900, r = 7 and b = 1000. p = 913262401 is chosen: we have 243771734ⁿ = 7 (mod p) and p > n · (b − 1)² = 898200900. Rather than evaluating the square, cube or 5^th root of an element, the n^th root of unity and the n^th root of 7 are computed and roots needed for the algorithm are of the form 1^u/n · 7^v/n.

fastMul_rec_twisted.cpp

The recursion is based on the relation x^am − r^am = Prod_0≤k<a (x^m − α^kr^m), where a is prime and α is a primitive root of x^a − 1.
Let X = x^m and R = r^m. If X − α^kR is twisted into Y_k − R for each k then at each step all polynomial moduli are identical.
The transformation is a linear map. P(X) mod X − α^kR = P(X) mod α^−kX − R. Let Y_k = α^−kX, we have y_k = α^−k/mx.
If P(x) = Sum_0≤i<m c_i xⁱ then P(y_k) = Sum_0≤i<m c_i α^ik/m y_kⁱ. The coefficients of P(y_k) are c'_i = α^ik/m c_i. After the recursion, the coefficients of P(y_k) are untwisted into P(x) with c_i = α^−ik/m c'_i.

Twisting is the generalization of weighted transforms, where weighting is applied at each step.

With the previous version of the fast multiplication, n/2 multiplications by 1^u/n · r^v/n are computed at each step. With the twisted form, n/2 multiplications by r^v/n and n/2 multiplications by 1^u'/n are needed. The twisted transform is as efficient as the non-twisted form only if r = 1.

Twisted FFT was defined in Daniel J. Bernstein, Multidigit multiplication for mathematicians, 2001. It was restricted to r = 1 and n is a power of two.

fastMul_rec_FFT.cpp

The twisted fast multiplication is written for r = 1. Roots of r are equal to one and it can be seen that the twist factors α^ik/m are the twiddle factors of a Fast Fourier Transform.

The transform that splits P(x) mod xⁿ − 1 into P(x) mod x − ωⁱ, where ω is a primitive n^th root of unity and 0≤i<n is Gentleman-Sande FFT recursion. The inverse transform that merges P(x)² mod x − ωⁱ into P(x)² mod xⁿ − 1 is Cooley-Tukey FFT recursion. Because Gentleman-Sande and Cooley-Tukey diagrams are symmetric, data ordering of the output is irrelevant, even if an iterative FFT algorithm is implemented.

The non-twisted fast multiplication can also be used for the computation of P(x)² mod xⁿ − 1. It is a different algorithm.
The forward transform was found by Georg Bruun, "z-Transform DFT filters and FFTs, 1978, IEEE Transactions on Acoustics, Speech, and Signal Processing, 26 (1): 56-63. Here, the Discrete Fourier Transform is written as a Z-transform filter. A Discrete Fourier Transform is evaluated with n roots of unity along the z-domain's unit circle. But this restriction does not apply to Z-transforms. If the FFT based multiplication is limited to the polynomial xⁿ − 1 (cyclic convolution), Bruun's method can be applied to any polynomial.

The two algorithms are implemented: because P(x)² mod xⁿ − 1 is calculated, both are almost equivalent. But the method based on polynomial factorization can compute P(x)² mod xⁿ − r or P(x)² mod Φ_m(xⁿ), where Φ_m is the m^th cyclotomic polynomial, without an extra "weighted" step.

fastMul_GF.cpp

Fast multiplication based on a recursive polynomial factorization is an in-place algorithm. The roots needed at each step can be precomputed. An iterative (non-recursive) version of fastMul_rec_GF.cpp is implemented here. Note that roots are stored in 2-3-5-reversed order: it is a generalization of bit-reversal permutation based on decomposing n into its prime factors. Rather than a division by 2, 3 or 5 at each step of the reverse process, a single division by n is computed at the end of the algorithm.

The number of operations is equal to a FFT based multiplication, where the polynomial of the convolution is x⁹⁰⁰ − 1. But here the polynomial is x⁹⁰⁰ − 7.

fastMul_GFp2.cpp

A full implementation of the fast multiplication of numbers of the form b²ⁿ + 1. The probable primes of this form are generated.

The recursion is b²ⁿ − (−1) = (b²ⁿ⁻¹ − i) (b²ⁿ⁻¹ − (−i)) = ..., where i is a root of −1.
The finite field is of order p², where p is the Mersenne prime 2⁶¹ − 1. With the irreducible polynomial X² + 1, we have α² = −1: operations are similar to the complex numbers. Because i and its negative are the first roots of the recursion, we can use a half length transform over GF(p²).

Three multiplications are needed for a radix-4 butterfly. If the field is GF(p), 2×4 multiplications are necessary for a radix-4 step with 8 values. Four multiplications in GF(p) are required to compute a multiplication in GF(p²). The balance is GF(p): 8 multiplications, GF(p²): 12 multiplications.
But with a prime finite field, p must be of the form k 2ⁿ⁺¹ + 1 such that a (n+1)^th root of unity exists. With GF(p²), if p = 2^q − 1 the order of the multiplicative group is p² − 1 = 2^q+1 p: there is a primitive root of unity of order 2^q+1.
Using 64-bit integers, on the one hand two multiplications modulo 2⁶⁴ − 2³² + 1 on the other hand three multiplications modulo 2⁶¹ − 1. The fastest implementation depends on the architecture of the processor.