fft.h

Sets (r1, n1 + n2) to the product of (i1, n1) and (i2, n2), choosing the transform parameters from tuning tables and calling mul_truncate_sqrt2 or, for large products, mul_mfa_truncate_sqrt2. Requires n_1, n_2 \ge 1, in either order. The output must not overlap the inputs; i1 == i2 gives a squaring. Intended for operands of thousands of limbs and more; assert builds reject products shorter than about 57 limbs. Uses threads for large products.

Sets (r1, n1 + n2) to the product of (i1, n1) and (i2, n2) with a truncated sqrt2 convolution of length at most 4n, n = 2^{\mathit{depth}}, modulo 2^{nw} + 1. The operands are split into chunks of b = \lfloor (nw - \mathit{depth} - 1)/2 \rfloor bits; if this gives j_1 and j_2 chunks, requires j_1 + j_2 - 1 \le 4n. Requires nw to be a multiple of \mathtt{FLINT\_BITS}. The pointwise products use flint_mpn_mulmod_2expp1_basecase. Single-threaded except for the splitting.

As mul_truncate_sqrt2, with the same requirements, but using the cache-friendly matrix Fourier algorithm, threads, and fft_mulmod_2expp1 for the pointwise products. Any w is allowed, but if the limb count nw/\mathtt{FLINT\_BITS} exceeds FFT_MULMOD_2EXPP1_CUTOFF it must be a value returned by fft_adjust_limbs. Requires \mathit{depth} \ge 2.

Multiplies the polynomials with coefficients ii and jj modulo p = 2^N + 1, N = \mathit{limbs} \cdot \mathtt{FLINT\_BITS}, with a truncated sqrt2 transform of length 4n, n = 2^{\mathit{depth}}, and writes the first trunc coefficients of the product, reduced, to ii. Requires 2n < \mathit{trunc} \le 4n (for \mathit{depth} \le 6 and an even w = N/n, 0 < \mathit{trunc} \le 4n suffices; other values can crash or hang) and that the product has at most trunc coefficients; for \mathit{trunc} = 4n the product is otherwise computed modulo x^{4n} - 1. Requires n to divide N and, if limbs exceeds FFT_MULMOD_2EXPP1_CUTOFF, limbs to be a value returned by fft_adjust_limbs. Both arrays have 4n entries, and the entries beyond each input must be zero. jj is overwritten unless jj == ii (squaring). The scratch arrays need one entry per thread.

As fft_convolution, but if sdctx is not NULL, the pointwise products modulo 2^N + 1 are computed by sd_fft_mpn_mulmod_2expp1 of fft_small with this context, which must have been initialised for N = \mathit{limbs} \cdot \mathtt{FLINT\_BITS} and the digit count returned by fft_adjust_limbs_sd_fft. With NULL it is fft_convolution.

Plans a convolution of which only the coefficients with indices in [\mathit{nlo}, \mathit{nhi}) of a product of length len_out are needed, where 0 \le \mathit{nlo} < \mathit{nhi} \le \mathit{len\_out}. Returns the depth and sets trunc for fft_convolution: either the truncated transform with \mathit{trunc} = \mathit{len\_out}, or a cyclic transform of length \mathit{trunc} = 4 \cdot 2^{\mathit{depth}} whose wrapped terms all fall below nlo, whichever processes fewer points. The truncated choice violates \mathit{trunc} > 2n when \mathit{len\_out} \le 2, so callers handle that case separately.

Replaces the coefficients jj by their reduced truncated sqrt2 transform, for use as the fixed operand of fft_convolution_precache with the same depth, limbs and trunc. The requirements are those of fft_convolution. Single-threaded.

As fft_convolution, where jj has been transformed by fft_precache with the same parameters. jj is not modified, so it can serve for any number of products. Single-threaded.

As fft_convolution_precache, with the optional pointwise engine of fft_convolution_sd_fft, whose pointwise products use threads.

Sets r to the product of the reduced residues i1 and i2 modulo 2^{nw} + 1, reduced, where nw is a multiple of \mathtt{FLINT\_BITS}. Up to FFT_MULMOD_2EXPP1_CUTOFF limbs (128 on all targets) it uses flint_mpn_mulmod_2expp1_basecase. Above, it uses a negacyclic FFT convolution, which requires the limb count nw/\mathtt{FLINT\_BITS} to be a value returned by fft_adjust_limbs; for other limb counts the result is wrong. r may alias the inputs. tt must have room for 2nw/\mathtt{FLINT\_BITS} limbs.

Returns the least limb count \ge limbs that fft_mulmod_2expp1 can handle: limbs itself if it is at most FFT_MULMOD_2EXPP1_CUTOFF, and otherwise a value rounded up to a suitable multiple, at most the next power of two.

Chooses the coefficient ring of a Schönhage–Strassen convolution with transform parameter n = n_outer and coefficients of at least tight limbs, taking the negacyclic engine of fft_small into account. If a ring N' = mb suitable for that engine is found, sets *m_out to its digit count m and returns N'/\mathtt{FLINT\_BITS}; otherwise sets *m_out to zero and returns fft_adjust_limbs(tight). Without fft_small it always does the latter.

Sets r1 to the reduced product of i1 and i2 modulo 2^N + 1, N = \mathit{r\_limbs} \cdot \mathtt{FLINT\_BITS}, with a negacyclic transform of length 2n, n = 2^{\mathit{depth}}, modulo 2^{nw} + 1 combined with a one-limb convolution. The inputs must be less than 2^N (top limbs zero), and depth and w must be as chosen by fft_mulmod_2expp1: nw = N/n, with N/(2n) a multiple of \mathtt{FLINT\_BITS}.

Sets r_k = \sum_{i + j = k} \mathit{ii}_i \mathit{jj}_j - \sum_{i + j = k + m} \mathit{ii}_i \mathit{jj}_j \bmod B for 0 \le k < m: the negacyclic convolution of two vectors of m limbs, with wrapping word arithmetic. r must not alias the inputs.

Split the integer (limbs, total_limbs) into consecutive chunks of coeff_limbs limbs, respectively bits bits, starting with the least significant, and write them to poly[0], poly[1], …, each zero-padded to \mathit{output\_limbs} + 1 limbs. Return the number of chunks written, \lceil \mathit{total\_limbs}/\mathit{coeff\_limbs} \rceil, respectively \lceil \mathtt{FLINT\_BITS} \cdot \mathit{total\_limbs}/\mathit{bits} \rceil; later entries of poly are not touched. Requires output_limbs to be at least the chunk size in limbs. May use threads.

Evaluate the polynomial with the length coefficients poly at B^{\mathit{coeff\_limbs}}, respectively 2^{\mathit{bits}}, and set (res, total_limbs) to the value, discarding the bits beyond total_limbs limbs. Each coefficient must be less than B^{\mathit{output\_limbs}} and have \mathit{output\_limbs} + 1 readable limbs. res must be zero on entry: the coefficients are added into it, but carries out of the window of each coefficient are lost; see Combining.

Adds the signed limb c to the (\mathit{limbs} + 1)-limb two's complement integer r, in place, without reduction. Called with r + k and limbs - k, it adds c B^k. Inline; when no carry leaves the first limb it costs a few instructions.

Reduces the residue t, whose top limb is a small signed value, to [0, p).

Sets z to -a \bmod p for a reduced residue a. The result is reduced.

Set t to i_1 \cdot 2^d, respectively i_1 / 2^d, modulo p, for 0 \le d < \mathtt{FLINT\_BITS} and a residue i_1 whose top limb is small. The result is not necessarily reduced.

Sets m to the signed integer represented by the \mathit{limbs} + 1 limbs of i in two's complement. The value is not reduced modulo p: to obtain the residue in [0, p), reduce i with mpn_normmod_2expp1 first. Used in tests.

Sets r to 2^{iw} i_1 \bmod p, for 0 \le iw \le N.

Sets r to \sqrt 2^{\,iw} i_1 \bmod p, for odd i < 2n and odd w, where nw = N. temp is scratch space of \mathit{limbs} + 1 limbs.

Sets t = (i_1 + i_2) B^x and u = (i_1 - i_2) B^y modulo p, for 0 \le x, y < \mathit{limbs}. The inputs are not modified.

Sets t = i_1 / B^x + i_2 / B^y and u = i_1 / B^x - i_2 / B^y modulo p, for 0 \le x, y < \mathit{limbs}. The inputs are overwritten when x \ne y.

The forward butterfly: sets s = i_1 + i_2 and t = 2^{iw} (i_1 - i_2) modulo p, for 0 \le iw < N. The inputs are not modified.

The inverse butterfly: sets s = i_1 + 2^{-iw} i_2 and t = i_1 - 2^{-iw} i_2 modulo p, for 0 \le iw < N. i2 is overwritten.

Sets s = i_1 + i_2 and t = \sqrt 2^{\,iw} (i_1 - i_2) modulo p, for odd i < 2n and odd w, where nw = N. temp is scratch space of \mathit{limbs} + 1 limbs.

Sets s = i_1 + \sqrt 2^{\,-iw} i_2 and t = i_1 - \sqrt 2^{\,-iw} i_2 modulo p, for odd i < 2n and odd w, where nw = N. i2 is overwritten; temp is scratch space of \mathit{limbs} + 1 limbs.

fft_radix2 replaces a by \hat a, the transform of length 2n with root 2^w. ifft_radix2 is its inverse up to scaling: it replaces \hat a by 2n \cdot a.

Truncated FFTs: compute only \hat a_k for k < trunc. fft_truncate assumes a_j = 0 for j \ge trunc and uses those entries as scratch; fft_truncate1 reads all 2n inputs. Requires trunc even with 2 \le \mathit{trunc} \le 2n.

Inverse truncated FFTs. Given \hat a_k for k < trunc, they set entry j to 2n \cdot a_j for j < trunc. ifft_truncate assumes a_j = 0 for j \ge trunc and uses those entries as scratch; ifft_truncate1 instead requires entry j to hold 2n \cdot a_j for j \ge trunc on input, and overwrites it. Requires trunc even with 2 \le \mathit{trunc} \le 2n.

As fft_truncate and ifft_truncate for the transform of length 4n with root \sqrt 2^{\,w} over the same ring; the inverse returns 4n \cdot a_j. Requires trunc even with 2n < \mathit{trunc} \le 4n.

The weighted transform of length 2n for negacyclic convolution. fft_negacyclic multiplies a_j by \varphi^j, \varphi = \sqrt 2^{\,w}, and then applies fft_radix2, so that \hat a_k = a(\varphi^{2 \mathrm{rev}(k) + 1}) are the values at the roots of x^{2n} + 1. ifft_negacyclic applies ifft_radix2 and then multiplies entry j by \varphi^{-j}, which gives 2n \cdot a_j.

fft_butterfly_twiddle sets u = 2^{b_1}(s + t) and v = 2^{b_2}(s - t). ifft_butterfly_twiddle sets u = 2^{-b_1} s + 2^{-b_2} t and v = 2^{-b_1} s - 2^{-b_2} t, overwriting s and t. Both work modulo 2^N + 1, N = \mathit{limbs} \cdot \mathtt{FLINT\_BITS}, for 0 \le b_1, b_2 < 2N. They form the last layer of the forward and the first layer of the inverse column transforms, with the twiddle factors of the matrix Fourier algorithm merged into the shifts.

As fft_radix2 and ifft_radix2 on the entries ii[0], ii[is], …, ii[(2n - 1) is], with an extra twiddle: output k of the forward transform is multiplied by 2^{\mathit{ws} \cdot c \cdot (r + \mathit{rs} \cdot \mathrm{rev}(k))}, and input k of the inverse is divided by the same factor first.

The truncated versions of fft_radix2_twiddle and ifft_radix2_twiddle, in the sense of fft_truncate1 and ifft_truncate1.

As fft_truncate_sqrt2 and ifft_truncate_sqrt2 (the inverse returns 4n \cdot a_j), computed with the matrix Fourier algorithm on each half of length 2n. Requires 2n < \mathit{trunc} \le 4n with trunc a multiple of 2 n_1. The transformed values are stored in a different order than by fft_truncate_sqrt2, the order that ifft_mfa_truncate_sqrt2 expects. Single-threaded.

The three stages of a threaded convolution with the matrix Fourier algorithm, as used by fft_convolution and mul_mfa_truncate_sqrt2. fft_mfa_truncate_sqrt2_outer performs the sqrt2 layer and the column transforms of ii. fft_mfa_truncate_sqrt2_inner performs, one row at a time, the row transforms of ii and jj, the pointwise products with fft_mulmod_2expp1 and the inverse row transform. ifft_mfa_truncate_sqrt2_outer performs the inverse column transforms and the sqrt2 layer and divides by 4n, leaving reduced residues. The requirements are those of fft_mfa_truncate_sqrt2 and, for the pointwise products, those of fft_mulmod_2expp1. The scratch arrays need one entry per thread; temp is not used by the _inner function.

Sets (s, n) to x + y and (d, n) to x - y, each modulo B^n, and returns 2c + b, where c is the carry of the sum and b the borrow of the difference. The outputs may alias the inputs. Declared here for the butterflies; defined in mpn_extras.

As fft_convolution, but computing the transforms and the pointwise products in separate single-threaded passes instead of the fused threaded stages. Used to test fft_convolution.

As fft_mfa_truncate_sqrt2_inner, with the pointwise products computed by sd_fft_mpn_mulmod_2expp1 with the context sdctx, initialised for N = nw; each thread gets its own scratch for the context. Only defined in builds with fft_small.

For j_0 \le j < j_1, reduce ii[j] (and, in the first function, jj[j]) and set ii[j] to the product of ii[j] and jj[j] modulo 2^N + 1 with sd_fft_mpn_mulmod_2expp1, using threads. The _ro version does not modify jj, which must already be reduced. Without fft_small they do nothing; they are then never called.

As _fft_pointwise_sd_fft and _fft_pointwise_sd_fft_ro, for the indices j = \mathit{base} + \mathrm{rev}(s) \cdot \mathit{srt} + u with 0 \le s < \mathit{nrows} and 0 \le u < \mathit{srt}, where \mathrm{rev} reverses rvbits bits: the rows of the truncated half of a matrix Fourier transform.