fft_small.h

Initialises R from the initial prime p and frees it. FLINT itself always uses p = 63 \cdot 2^{44} + 1 (0x0003f00000000001), the value get_default_mpn_ctx passes.

Returns the default context of the calling thread, creating it on the first call. flint_cleanup and flint_cleanup_master free it.

Returns at least n bytes of scratch memory from the context, aligned to 4096 bytes. The memory belongs to the context and must not be freed. A later call may move the buffer, so the pointer is valid only until the next call on R, and only one request can be outstanding.

Reserves the tail of the scratch buffer: the buffer is grown to hold head plus tail bytes (each rounded up to 4096), and the returned pointer to the last tail bytes stays valid until the release. Returns NULL if a reservation is already live. Reserving, like any growth, invalidates a pointer returned by mpn_ctx_fit_buffer.

Sets (z, a_n + b_n) to the product of (a, a_n) and (b, b_n), using R or the default context. Requires a_n, b_n \ge 1, in either order. If a == b and a_n = b_n, a square is computed.

Writes to (z, \mathit{hi} - \mathit{lo}) the limbs [\mathit{lo}, \mathit{hi}) of the product of (a, a_n) and (b, b_n). The result is exact for \mathit{lo} = 0; otherwise carries from below limb \mathit{lo} may be lost, making it smaller than the exact window by a few units. Limbs at or above a_n + b_n are zero.

Sets (z_1, a_n + x_n) to ax and (z_2, b_n + x_n) to bx, computing the transform of x once: three forward and two inverse transforms instead of four and two. Requires a_n, b_n, x_n \ge 1, in any order of size. z1 may overlap a and x but not b; z2 may overlap any input but not z1.

_fft_small_mpn_mul_precomp_size returns the number of doubles that the transform of an integer of x_n limbs needs for partners of up to yn_max limbs. _fft_small_mpn_mul_precomp_init stores the transform of (x, x_n) in storage, which must hold that many doubles and be aligned to 32 bytes, or in memory it allocates if storage is NULL. _fft_small_mpn_mul_precomp_clear frees that memory, but not storage provided by the caller. Require x_n, \mathtt{yn\_max} \ge 1.

Sets (z, x_n + y_n) to xy, where x is the integer whose transform P holds. Requires 1 \le y_n \le yn_max; z may overlap y.

Sets (z, z_h - z_l) to the coefficients [z_l, z_h) of the product of (a, a_n) and (b, b_n) modulo n. Requires a_n, b_n \ge 1, in either order, and reduced coefficients; any modulus n \ge 1 is allowed.

Transposed product: sets z_j = \sum_i a_i v_{i+j} \bmod n for 0 \le j < z_n, with (v, v_n) read as padded with zeros to length a_n + z_n - 1. This is the transpose of the multiplication by (a, a_n) of polynomials of length z_n; with a the reversal of f it is the window [a_n - 1, a_n - 1 + z_n) of the product fv. Requires a_n, v_n, z_n \ge 1, v_n \le a_n + z_n - 1 and reduced coefficients; the output must not overlap the inputs.

Sets (z, z_h - z_l) to the coefficients [z_l, z_h) of the product of (a, a_n) and (b, b_n) over \mathbb{Z} and returns 1, if eight primes suffice for the coefficient size; otherwise returns 0 and leaves z unchanged. Requires a_n, b_n \ge 1; z must hold initialised fmpz values.

Sets (q, a_n - b_n + 1) and (r, b_n - 1) to the quotient and remainder of (a, a_n) divided by (b, b_n) modulo n. Requires a_n \ge b_n \ge 2 and an invertible leading coefficient of b.

Reference implementations of the window product without transforms, used by the tests and profiles. _nmod_poly_mul_mid_classical evaluates each coefficient as a sum reduced after every term; _nmod_poly_mul_mid trims the operands to the window and splits a long operand into pieces of the length of the short one before calling it. Neither is used by FLINT's polynomial arithmetic, which has its own _nmod_poly_mulmid.

Initialises M with the transform of (b, b_n) for convolutions of length 2^d, d = depth \ge 4, in which at most len_bound products accumulate in one coefficient, and returns 1; returns 0 if more than eight primes would be needed. btrunc is the truncation of the transforms M will be used with, an upper bound for the output lengths.

Sets (z, z_h - z_l) to the coefficients [z_l, z_h) of the product of (a, a_n) with the polynomial stored in M and returns 1, or returns 0 if the window needs a transform longer than 2^d. Requires a_n \ge 1.

Initialises M for the transposed products with (v, v_n) computed by _nmod_poly_mul_tr_precomp, in which at most len_bound products accumulate in one coefficient, and returns 1; returns 0 if no plan is available. M is freed with _mul_precomp_clear.

Sets z_j = \sum_i r_i v_{i+j} \bmod n for 0 \le j < z_n, with (v, v_n) the polynomial of M made by _mul_precomp_init_tr, and returns 1; returns 0, writing nothing, if r_n + z_n - 1 > v_n. Requires r_n, z_n \ge 1, at most len_bound nonzero coefficients in r and reduced coefficients.

Prepares the divisor (b, b_n) for quotients of length at most B_n, that is, for dividends of length at most b_n + B_n - 1. Requires b_n \ge 2 and an invertible leading coefficient.

Sets (q, a_n - b_n + 1) and (r, b_n - 1) to the quotient and remainder of (a, a_n) by the divisor prepared in M, and returns 1. Requires a_n \ge b_n. Returns 0 without computing anything if the quotient length a_n - b_n + 1 exceeds the precision B_n of the stored inverse or needs a longer transform than the one prepared.

Initialises P for computing the coefficients [z_l, z_h) of convolutions of length z_n over \mathbb{Z}/n\mathbb{Z} in which at most len_bound products, each less than 2^{\mathtt{prod\_bits}}, accumulate in one coefficient. Returns 0 if more than eight primes would be needed. Requires z_l < z_h \le z_n.

Initialises P for cyclic convolutions modulo x^N - 1, N = 2^d with d = depth \ge 4, writing the coefficients [0, z_h), where z_h \le N. len_bound and prod_bits are as for fft_small_plan_init_nmod; a full cyclic product accumulates up to N products. Returns 0 if more than eight primes would be needed.

Initialises P for products of nonnegative integers of at most an_max and bn_max limbs, of which at most len_bound are added in transform space. If is_signed is nonzero, the plan also allows subtractions and the signed exports. Returns 1 on success and 0 if no prime count and chunk size satisfy the bound.

Like fft_small_plan_init_mpn, with the regime forced for testing and profiling: np = 0 chooses automatically, np = 2 forces the two-prime geometry (and fails if its bound cannot be met), and any other value forces the profile selection with four to eight primes.

Initialises P for products of nonnegative integers of at most n_n = *nn limbs modulo B^{n_n} - 1. The length is rounded up to the nearest one for which the chunks split n_n limbs exactly into 2^d chunks, and the new value is written to *nn. Returns 0 if no geometry exists.

fft_small_plan_init_mpn_window is fft_small_plan_init_mpn for products of which only the limbs [\mathit{lo}, \mathit{hi}) are needed. fft_small_plan_mpn_nslots returns the number of chunk coefficients (slots) that the inverse transform computes for operands of at most an_max and bn_max limbs, and fft_small_plan_is_cyclic whether slots alias, that is, whether the transform is shorter than the product.

Frees the private transform context of a plan that transforms modulo n itself, and does nothing otherwise.

Initialise an operand for the plan P, with storage allocated by the operand or provided by the caller, and free it. The borrowed storage must hold fft_small_op_sizeof_data(P) bytes, be aligned to 4096 bytes and outlive the operand; fft_small_op_clear does not free it.

Sets X to the forward transform of (a, a_n), a polynomial over \mathbb{Z}/n\mathbb{Z} with reduced coefficients, read as zero-padded to length itrunc. Requires itrunc to be a multiple of 8 and at most P->ztrunc; normally itrunc is a_n rounded up to a multiple of 256, or 2^d if the depth d is less than 8.

Sets X to the forward transform of the integer (a, a_n), cut into chunks of P->bits bits. Requires an integer plan and a_n \ge 1, at most the operand length the plan was made for.

Like fft_small_fft_mpn, with the transforms truncated to their first otrunc points, which suffices for products of at most otrunc slots. Requires otrunc \le P->ztrunc, a multiple of 256 or equal to P->ztrunc; fft_small_fft_mpn is the case otrunc = P->ztrunc.

Set Z to the pointwise product AB, respectively A^2, of operands in the primal domain. Z is then in the product domain and may alias A or B.

Add or subtract the pointwise product AB, respectively A^2, of primal operands to Z, which must be in the product domain. The sqr versions read A once and are preferable to addmul with B = A.

Set Z to A + B, A - B and -A, for operands in the same domain; Z gets that domain. Sums of primal operands represent sums of the inputs, whose larger coefficients the plan's bound must allow for.

Applies the inverse transform in place to an operand in the product domain, at the truncation of the plan. The normalisation was applied by the pointwise products. The primes are processed in parallel from P->ztrunc \ge 3072.

Inverse transform of the first itrunc slots in place, with the truncation rounded up as for the plan. If scale is not NULL, scale[i] is a normalisation factor modulo the i-th prime that was not applied in the pointwise stage, and is applied after the transform. The slots [\mathtt{itrunc}, \mathtt{clear\_to}) are set to zero, for exports that read whole chunks. X is marked as being in the product domain.

The transposes of fft_small_ifft and fft_small_fft_nmod, for transposed products. fft_small_ifft_tr_nmod converts (a, a_n), read as zero-padded to P->ztrunc coefficients (so a_n \le P->ztrunc), and applies the transposed inverse transform; X is then in the primal domain. fft_small_fft_tr applies the transposed forward transform with truncation itrunc \le P->ztrunc in place to an operand in the product domain, whose first itrunc slots can then be exported.

Write the coefficients of the plan's window [z_l, z_h), respectively of the subwindow [z_l, z_h) \subseteq [P->zl, P->zh), reduced modulo n, to z, which receives z_h - z_l coefficients. The exact integer coefficients must lie in [0, P) for the product P of the plan's primes, as guaranteed by its bound for sums of products.

Writes to (z, z_n) the integer \sum_k s_k 2^{k \cdot \mathtt{bits}} reduced modulo B^{z_n}, where s_k \in [0, P) are the reconstructed slots. For a product of integers this is the product, exact if z_n is at least its length.

Reconstruct a signed integer from the first nslots slots, each lifted to the symmetric range (-P/2, P/2]. The absolute value is written to (z, z_n), zero-padded, and its sign to *sign. The trunc version writes the limbs from limb lo_limbs on, within (-1.5, 0.5) units of its lowest limb; with lo_limbs = 0 both are exact.

Write to (z, z_n) the limbs [\mathtt{lo\_limbs}, \mathtt{lo\_limbs} + z_n) of the two's complement value of the signed reconstruction, without resolving the sign; _low is the case lo_limbs = 0, which is the residue modulo B^{z_n} and exact. The value need not fit in the window.

Return the number of limbs that the signed reconstruction of nslots slots, respectively the unsigned one with up to len_bound accumulated products per slot, writes. An export into a shorter window must go through scratch memory, as fft_small_stage_begin arranges.

Initialise an empty pool, give it count slabs of slab_size bytes and a staging region of stage_bytes bytes, and free it. fit returns 0, leaving the pool empty, if the block can be neither reserved nor allocated as flags allow.

Take a free slab, for fft_small_op_init_borrowed, or NULL if none is free; give a slab back; and test whether an address lies in the pool's block. Taking and giving are thread safe.

Return scratch space for limbs limbs and give it back. The space is the staging region if it is free and large enough, otherwise a free slab if use_slabs is nonzero and one is large enough, otherwise a fresh allocation; concurrent requests never share it.

Staged conversion of a value whose reconstruction writes need limbs into a window of z_n limbs. begin returns z if z_n \ge need and scratch space from the pool otherwise (pool may be NULL, meaning a plain allocation). After exporting there, and possibly adjusting the value, end copies the normalised value of length len to (z, z_n), zero pads it, releases the scratch space, and returns 0 if it does not fit (the window is then undefined).

Write the limbs [\mathit{lo}, \mathit{hi}) of the real and imaginary parts of the product, respectively the square, as two's complement windows of \mathit{hi} - \mathit{lo} limbs. Each is within one unit of \lfloor x / B^{\mathit{lo}} \rfloor for the exact part x, and exact (the residue modulo B^{\mathit{hi}}) for \mathit{lo} = 0. Return 0 if no plan exists.

Write the magnitudes of the limbs [\mathit{lo}, \mathit{hi}) of the real and imaginary parts of the product, respectively the square, to zr and zi (\mathit{hi} - \mathit{lo} limbs each, enough to hold the parts shifted right by \mathit{lo} limbs), and their lengths, negated for negative parts, to *zr_len and *zi_len. Exact for \mathit{lo} = 0 and otherwise within (-1.5, 0.5) units of the lowest limb. Return 0 if no plan exists.

Initialises C for multiplication modulo 2^N + 1 with m digits of b = N/m bits, taking the primes from R, normally get_default_mpn_ctx(), and frees it. N must be a positive multiple of FLINT_BITS, and m a power of two at least 16 dividing N, or 0 to let the function choose. Throws if N is not a positive multiple of FLINT_BITS or no configuration exists.

Initialise and free the scratch object for multiplications with C. The context is read-only during multiplication, so several threads can share it, each with its own scratch object.

Sets z = xy \bmod 2^N + 1. z may alias x or y; if x == y a square is computed with one forward transform.

Multiplication by a fixed operand. transform stores the transform of y in F, an array of transformed_size(C) doubles; mul_cached then sets z = xy \bmod 2^N + 1, like sd_fft_mpn_mulmod_2expp1 with one forward transform less.

Returns the number of primes the context uses for digits of b bits and m = 2^{\mathtt{depth}} digits, as described for the context, or -1 if there is none.

Initialise the context of the prime p = pp and free it. Throws if p fails fft_small_mulmod_satisfies_bounds, if p - 1 is odd or if 2^4 \nmid p - 1. Primality is assumed, not checked.

Make sure that the twiddle tables cover transforms of length 2^{\mathtt{depth}}, filling the missing tables under the context's mutex. Throw if 2^{\mathtt{depth}} \nmid p - 1. sd_fft_trunc and sd_ifft_trunc call this themselves.

Truncated forward transform of length 2^L in place: with f the polynomial with coefficients d[0], ..., d[itrunc - 1], sets d[sd_fft_ctx_trunc_index(L, i)] to a value congruent to f(w_i) for 0 \le i < \mathtt{otrunc}. The inputs must lie in (-3p, 3p), and so do the outputs.

Truncated inverse transform of length 2^L in place: if d[i] holds the value at position i of the forward transform of a polynomial g of length at most trunc, for i < \mathtt{trunc}, sets d[i] to a value congruent to 2^L g_i for i < \mathtt{trunc}. The inputs must lie in (-2p, 2p); the outputs lie in (-3p, 3p).

The transposes of sd_fft_trunc and sd_ifft_trunc with the same parameters, in place: for all vectors a and v, \langle \mathtt{sd\_fft\_trunc}(a), v \rangle = \langle a, \mathtt{sd\_fft\_trunc\_tr}(v) \rangle, and likewise for the inverse transform. sd_fft_trunc_tr reads the first otrunc entries, indexed like the output of sd_fft_trunc, and sets d_t to a value congruent to \sum_{i < \mathtt{otrunc}} d[\mathtt{sd\_fft\_ctx\_trunc\_index}(L, i)]\, w_i^t for 0 \le t < \mathtt{itrunc}. Each costs the same as the transform it transposes.

Set a_i \gets m a_i b_i, respectively m a_i^2, modulo p for i < 2^{\mathtt{depth}}, with 0 \le m \le p and depth at least 4. Inputs in (-3p, 3p) give outputs in (-3p/2, 3p/2), valid inputs for sd_ifft_trunc.

Return the root of unity w_j, w_{2j}, and -w_{2j}^{-1} modulo p, as integers in [-p/2, p/2]. Here w_0 = 1 and, for 2^{k-1} \le j < 2^k, w_j = \omega^{r} where \omega is a primitive 2^k-th root of unity and r is j with its k bits reversed; the roots for different k are chosen compatibly, \omega_{2^k} = \omega_{2^{k+1}}^2. The tables must cover the index.

The position of f(w_i) in the output of a forward transform of length 2^L: i itself for L < 4, and otherwise i with the two 2-bit fields of its low four bits exchanged, which transposes each group of 16 values as a 4 \times 4 matrix.

The offset 256 I of block I, the address d + 256 I, and the number 2^L of doubles of a transform of length 2^L.

Returns whether the floating-point modular multiplication modulo n maps products in (-2n^2, 2n^2) into (-n, n) and products in (-4n^2, 4n^2) into (-3n/2, 3n/2), by a bound that uses the actual error of the double 1/n. Every prime below 2^{49} passes, and about one 50-bit prime in 900 fails.

Returns whether a transform of length 2^{\mathtt{depth}} can be computed modulo n itself: 2 < n \le 2^{50}, n passes fft_small_mulmod_satisfies_bounds, depth is at most 40 and 2^{\max(\mathtt{depth}, 4)} \mid n - 1 (the context initialisation sd_fft_ctx_init_prime needs 2^4 \mid n - 1 for every depth). Primality is assumed.

Allocate storage for the data, which the caller fills in, and free it.

Pointers to the cofactor P / p_i (coeff_len limbs; the second form takes coeff_len as n for the compiler), to P (coeff_len limbs), and to the single limb (P / p_i) \bmod p_i.

Return the vectorised conversion of integers into chunks of bits bits modulo np primes, or NULL if none is instantiated, and the reconstruction for 3 \le np \le 8 primes. Used by the negacyclic engine.

For each prime l < np, reduce block I of the data d + l \cdot dstride to [0, p_l) and store it as integers at Xs + 256 l. The second form converts only the first prime and reduces the results further modulo the double p2 (with p2inv = 1/p_2).

Chooses the least number of primes k \le np_max whose product is at least c 2^e, starting at the first prime of P->R, and records the bound. Returns 0 if there is none.

Choose the depth and truncation for the window [z_l, z_h) of a convolution of length z_n with operands of truncated length at most xtrunc_max, and record the window. Require z_l < z_h \le z_n.

Like _fft_small_plan_set_window with the depth of P fixed. Returns 0 if the depth is below min_depth or the window needs a longer transform.

Computes the stride and the normalisation factors m[i] from the primes and the depth, which must be set.

Computes the window of the plan of the convolution described by A into z, using the scratch buffer of the plan's context.

The full product, its low z_h limbs, and its limbs [z_l, z_h), by the two-prime pipeline. Require a_n \ge b_n \ge 1 and z_l < z_h \le a_n + b_n. The window is exact for z_l = 0 and otherwise smaller than the exact window by 0 or 1 units of its lowest limb.

Chooses the least depth d \in [8, 30] and then the largest chunk size b \in [40, 45] with \lceil 64 a_n / b \rceil + \lceil 64 b_n / b \rceil - 1 \le 2^d and len_bound \cdot \lceil 64 b_n / b \rceil \cdot 2^{2b} < p_0 p_1. Returns 0 if there is none.

Writes the first len chunks of bits \le 45 bits of (x, x_n) as doubles to both d and d2, followed by zeros up to itrunc.

Writes to z the limbs [z_l, z_h) of \sum_j v_j 2^{j \cdot \mathtt{bits}} over nchunks chunk values v_j given by their residues modulo p_0 in d0 and modulo p_1 in d1, which it overwrites. Unless depth is UWORD_MAX, the residues are raw inverse transform outputs and are multiplied by 2^{-\mathtt{depth}}. navail8 is the number of entries available to the vectorised Garner step.

Return 2^k, and the minimum and maximum of two unsigned and of two signed words.

Return \lceil a/b \rceil, b \lceil a/b \rceil and b \lfloor a/b \rfloor. The first two compute a + b - 1 and are wrong if it overflows.

The bit length of x (0 for x = 0; _nz requires x \ne 0), the numbers of leading and trailing zero bits (FLINT_BITS for x = 0), and 2^{k} - 1 for a of bit length k.

Return \lceil \log_2 x \rceil and \lfloor \log_2 x \rfloor for x \ge 1; both are 0 for x = 0.

The bit length of (a, a_n), which need not be normalised, and the sign of a - b 2^e.

Adds (a, a_n) and the carry cf to (z, z_n) in place, propagating the carry through at most z_n limbs, and returns the carry out. Requires z_n \ge a_n \ge 1.