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.