mpn_extras.h

Set the n limbs of xp to zero, respectively to y. These are inline loops, cheaper than calls to memset or mpn_zero for the short arrays typical of FLINT.

Copy the n limbs of yp to xp, in increasing respectively decreasing order of addresses. The arrays may overlap if xp \le yp (for flint_mpn_copyi) or xp \ge yp (for flint_mpn_copyd).

Returns 1 if all limbs of (x, \mathit{xsize}) are zero, and 0 otherwise.

Returns 1 if the arrays (x, \mathit{xsize}) and (y, \mathit{xsize}) are equal limb by limb, and 0 otherwise.

Returns the representation of (x, \mathit{xn}), or of its negation if negative is 1, in the given base, which must be one supported by GMP (2 \le base \le 62). If res is NULL, a new string is allocated with flint_malloc; otherwise res must have room for all the digits, one extra digit, a minus sign and the terminating null.

Returns the double nearest towards zero to \pm (\mathit{ptr}, \mathit{size}) \cdot 2^{\mathit{exp}}, where the sign is that of sign. Requires size \ge 0 and a nonzero top limb if size > 0. Overflow gives an infinity of the right sign and underflow (below the denormal range) gives zero. This is a copy of GMP's internal mpn_get_d.

Prints (x, \mathit{xsize}) in decimal to stdout, followed by the binary digits of each limb, least significant bit first.

Sets (\mathit{rp}, n) to r + a + b modulo B^n and returns the carry (0, 1 or 2). Requires n \ge 4. Only available on x86-64 with the ADX instructions, where it is an assembly routine running two carry chains in parallel; FLINT_HAVE_NATIVE_mpn_2add_n_inplace is then defined.

Sets (\mathit{res}, n) to |x - y| and returns 0 if x \ge y and 1 if x < y. The output may alias either input.

flint_mpn_signed_sub_n for a fixed length from 2 to 8, with the subtraction done by straight-line carry chains (the NN_SUB macros).

Sets (\mathit{res}, n) to x / 2 rounded towards -\infty, where x is read as an n-limb two's complement integer (an arithmetic shift right by one bit). The output may alias the input.

Tries to compute y = a_1 x_1 - a_2 x_2, where x_1 and x_2 have n \ge 1 limbs (not necessarily normalised). If the result is nonnegative and fits in n limbs, it is written to (y, n) and its normalised length is returned; otherwise -1 is returned and (y, n) holds garbage. y may alias x_1 but not x_2.

Set (\mathit{res}, n) to -x, x + y and x - y modulo (m, n), for x and y reduced modulo m. The output may alias the inputs.

Versions of flint_mpn_addmod_n and flint_mpn_submod_n where the second operand has only yn \le n limbs (and is reduced modulo m).

Versions for two-limb moduli, written with two-limb arithmetic macros. _flint_mpn_addmod_2 requires m < B^2/2 (the top bit of m[1] clear), so that the sum cannot overflow two limbs.

Sets (z, \mathit{xn} + \mathit{yn}) to the product of (x, \mathit{xn}) and (y, \mathit{yn}) and returns its top limb z_{xn + yn - 1}. Requires \mathit{xn} \ge \mathit{yn} \ge 1; z must not overlap the inputs.

Sets (z, 2n) to the product of (x, n) and (y, n). Requires n \ge 1; z must not overlap the inputs.

Sets (z, 2n) to the square of (x, n) and returns its top limb. Requires n \ge 1; z must not overlap x.

Sets (z_1, \mathit{an} + \mathit{xn}) to ax and (z_2, \mathit{bn} + \mathit{xn}) to bx. All lengths must be at least 1, in any order; z_1 and z_2 must overlap neither each other nor any input. When both products are long enough, the transform of x is computed only once (Integer products with a common operand).

Returns nonzero if flint_mpn_mul_shared uses a shared transform for these lengths, that is, if fft_small is available and _fft_small_mpn_mul_shared is expected to beat two separate products. Callers that need the looser aliasing rules of _fft_small_mpn_mul_shared can call it directly when this returns nonzero.

The non-inline part of flint_mpn_mul, flint_mpn_mul_n and flint_mpn_sqr, used for sizes beyond the hard-coded routines. Same requirements.

Karatsuba (Toom-22) multiplication, a port of GMP's mpn_toom22_mul that recurses into flint_mpn_mul, so that FLINT's assembly basecases are used at the leaves. Sets (\mathit{pp}, \mathit{an} + \mathit{bn}) to the product, for \mathit{an} \ge \mathit{bn} with 4 \mathit{an} < 5 \mathit{bn} (as required by GMP). scratch must have room for 2 \mathit{an} + 2 \mathtt{FLINT\_BITS} limbs, or be NULL for internal allocation.

Toom-32 multiplication of unbalanced operands, a port of GMP's mpn_toom32_mul (evaluation at 0, \pm 1, \infty) built on flint_mpn_mul. Requires \mathit{bn} + 2 \le \mathit{an} and \mathit{an} + 6 \le 3 \mathit{bn} (GMP's conditions). scratch may be NULL. It is not currently used by flint_mpn_mul.

Sets (\mathit{rp}, n + 2) to the product of (\mathit{ap}, n) and the two-limb (\mathit{bp}, 2) and returns the top limb. Only available on ARM64 (FLINT_HAVE_NATIVE_mpn_mul_2), where it is an assembly routine used by flint_mpn_mul.

Sets (\mathit{res}, n) to an approximation of \lfloor uv / B^n \rfloor and returns an approximation of limb n - 1 of the product. Together they form a lower bound for \lfloor uv / B^{n-1} \rfloor that is wrong by a few units of the returned limb at most. Requires n \ge 1. The output may alias the inputs only for n at most FLINT_MPN_MULHIGH_FUNC_TAB_WIDTH.

The square u^2 with the contract of flint_mpn_mulhigh_n.

flint_mpn_mulhigh_n and flint_mpn_sqrhigh followed by normalisation: if the top bit of the high product is clear, the (n+1)-limb result (high limbs and returned limb) is shifted left by one bit. Returns the low limb in m1 and the shift (0 or 1) in m2. For normalised inputs the product is at least B^{2n}/4, so the result is normalised. This is the form needed for floating-point mantissas (mp_real, arf). The output must not alias the inputs.

Sets (\mathit{res}, n) to uv \bmod B^n and returns limb n of the product, that is, \lfloor uv / B^n \rfloor \bmod B (this extra limb may be removed from the interface in the future). Requires n \ge 1; res must not alias u.

Writes uv \bmod B^{n+1} to the low n + 1 limbs of res. The output must have room for 2n limbs: the high limbs may be used as scratch space, or the whole product may be written when this is fastest.

Writes a high product (with the error of flint_mpn_mulhigh_n) as n + 1 limbs to res + n - 1; the low n - 1 limbs of the 2n-limb output may be used as scratch space, or the full product may be written.

The algorithms behind flint_mpn_mulhigh_n, with its interface.

Rough high product by Mulders' recursive algorithm, as described by Harvey and Zimmermann [HZ2011]. Writes to \mathit{rp}_n, \ldots, \mathit{rp}_{2n-1} an approximation of the top n limbs of the product (respectively the square), which is a lower bound with an error of less than n units of \mathit{rp}_n. rp must have room for 2n limbs; the low n limbs are used as scratch space. The split points come from a tuned table flint_mpn_mulhigh_k_tab.

The algorithms behind flint_mpn_sqrhigh: dispatcher, quadratic basecase (on x86-64 with ADX an assembly routine requiring n \ge 8; elsewhere the multiplication basecase), Mulders' algorithm (from FLINT_MPN_SQRHIGH_MULDERS_CUTOFF, 90 limbs) and a full square (from FLINT_MPN_SQRHIGH_SQR_CUTOFF, 2000 limbs, or a middle product with fft_small from FLINT_MPN_SQRHIGH_FFT_SMALL_CUTOFF).

The non-inline part of flint_mpn_mulhigh_normalised and flint_mpn_sqrhigh_normalised, used beyond the hard-coded sizes.

The algorithms behind flint_mpn_mullow_n: dispatcher, quadratic basecase (up to FLINT_MPN_MULLOW_MULDERS_CUTOFF, 50 limbs; an assembly routine on x86-64 with ADX), Mulders' recursive short product (up to FLINT_MPN_MULHIGH_FFT_SMALL_CUTOFF or FLINT_MPN_MULHIGH_MUL_CUTOFF), and a full product (or an fft_small middle product). The low product is exact in every case, and so is the returned limb n. _flint_mpn_mullow_n_mulders_recursive writes the n low limbs only and uses rp beyond them as scratch space (it needs room for 2n limbs).

High and low products computed by flint_mpn_mulmid_fft_small. Only available with fft_small.

Sets (z, \mathit{zhi} - \mathit{zlo}) to the window [\mathit{zlo}, \mathit{zhi}) of ab with the contract above, choosing the fastest of the methods below for the shape. The operands may be given in either order; squaring is detected. z must not overlap the inputs.

The methods behind flint_mpn_mulmid, with its interface and contract: row-wise schoolbook multiplication in time O((\mathit{zhi} - \mathit{zlo}) \min(\mathit{an}, \mathit{bn})); a full product (flint_mpn_mul); a balanced low product (flint_mpn_mullow_n); a balanced high product (flint_mpn_mulhigh_n); GMP's balanced middle product (flint_mpn_mulmid_n); and the small-prime FFT, which only reconstructs the coefficients that reach the window. Each accepts any shape (padding internally) but is economical only in its own regime. flint_mpn_mulmid_via_n_padded exists only if FLINT_HAVE_NATIVE_mpn_mulmid_n is set, and flint_mpn_mulmid_fft_small only with fft_small.

Balanced middle product of (\mathit{ap}, 2n - 1) and (\mathit{bp}, n): writes the n + 2 limbs of \sum_{n-1 \le i + j \le 2n - 2} a_i b_j B^{i+j-n+1}, that is, the diagonals n - 1 to 2n - 2 of the product. Compared with the corresponding limbs of the full product only the carry from the lower diagonals is missing, so the low two limbs serve as guard limbs. This is GMP's internal mpn_mulmid_n, and is only declared if FLINT_HAVE_NATIVE_mpn_mulmid_n is set (when the build may call GMP internals and GMP exports the symbol).

Sets (\mathit{out}, \mathit{khi} - \mathit{klo}) to the exact limbs [\mathit{klo}, \mathit{khi}) of the product xy, given that its low \mathit{klo} limbs are known and equal to (\mathit{kl}, \mathit{kl\_len}) (zero above kl_len). The window is computed from a windowed middle product starting three guard limbs below \mathit{klo}, whose deficit is less than B^2 units of its lowest limb, and the known low limbs tell which carry was lost. scratch must have room for \mathit{khi} limbs. This is used by Hensel division and the 2-adic square root, where the low part of a product is known in advance.

Sets z_r + z_i i = (a_r + a_i i)(b_r + b_i i), respectively (a_r + a_i i)^2. Each part has its own length, at least 1 limb, not necessarily normalised. A signed length is written for each output: the magnitude occupies |*zr_len| limbs and a negative value means that the part is negative. Limbs beyond the signed length are not significant. Each output must have room for \max(\mathit{arn}, \mathit{ain}) + \max(\mathit{brn}, \mathit{bin}) + 1 limbs (2 \max(\mathit{arn}, \mathit{ain}) + 1 for the square).

The individual algorithms behind the general functions, with the same interface; all accept any shape. The fft_small versions return 0, leaving the outputs untouched, when the method is unavailable or refuses the operands, and 1 on success.

High products: all parts have n limbs, and each output part receives exactly n + 1 limbs, zero padded, approximating the limbs [n, 2n] of the exact part, together with its sign.

Low products: zr and zi receive exactly n \ge 1 limbs, the real and imaginary parts of the product (respectively the square) modulo B^n as two's complement residues. No sign is returned: a negative part x is represented by B^n + x when |x| < B^n. The inputs are as for flint_mpn_mul_complex, and only their low n limbs are read. n may exceed the length of the product, the residues then being the sign-extended exact parts.

The individual algorithms behind the low products, with the conventions of the full products' individual algorithms.

Windowed products: zr and zi receive \mathit{hi} - \mathit{lo} limbs each, 0 \le \mathit{lo} < \mathit{hi}, the limbs [\mathit{lo}, \mathit{hi}) of the real and imaginary parts in two's complement: for each exact part x the window is \lfloor x / B^{\mathit{lo}} \rfloor + e \bmod B^{\mathit{hi} - \mathit{lo}} with a small error e. For \mathit{lo} = 0 these are the low products (e = 0); on the transformed path |e| \le 1; otherwise |e| < 4 (m + 1) B with m = \min(\max(\mathit{arn}, \mathit{ain}), \max(\mathit{brn}, \mathit{bin})). Inputs are as for the low products (only their low \mathit{hi} limbs are read).

The individual algorithms behind the windowed products, with the error bounds stated there.

Reversed dot products for any n_1 \ge n_2 \ge 1, \mathit{len} \ge 1 and s \in \{n_1 + n_2 - 1, n_1 + n_2, n_1 + n_2 + 1\}, using flint_mpn_mul and mpn_add_n. The scratch space must have n_1 + n_2 limbs. The sdiff version takes s - (n_1 + n_2 - 1), which should be a compile-time constant. The signed version handles the signs with branches and is competitive only when they are predictable. In the sign-magnitude version each coefficient is a sign limb (0 or 1) followed by the magnitude, so that the strides are n_1 + 1 and n_2 + 1; the output is in two's complement, and the scratch space must have n_1 + n_2 + 2s limbs.

Sets (\mathit{res}, (\mathit{nhi} - \mathit{nlo}) s) to the coefficients [\mathit{nlo}, \mathit{nhi}) of the product of (f, \mathit{flen}\, n_1) and (g, \mathit{glen}\, n_2), viewed as polynomials with n_1-limb respectively n_2-limb coefficients, computed modulo B^s by the classical algorithm. Requires n_1 \ge n_2 and s \in \{n_1 + n_2 - 1, n_1 + n_2, n_1 + n_2 + 1\}. The coefficients are unsigned for sgn = 0, two's complement for sgn = 1, and sign-magnitude for sgn = 2 (a sign limb followed by n_1 respectively n_2 magnitude limbs; the output is two's complement). The results are exact provided that the requested coefficients fit in s limbs. Squaring is detected when f == g and flen == glen.

Set (\mathit{res}, (\mathit{flen} + \mathit{glen} - 1)\, \mathit{slimbs}) to the full product of (f, \mathit{flen}\, \mathit{nlimbs}) and (g, \mathit{glen}\, \mathit{nlimbs}), viewed as polynomials with nlimbs-limb coefficients (unsigned for sgn = 0, two's complement for sgn = 1), by the classical algorithm respectively Karatsuba's. The output coefficients are computed modulo B^{\mathit{slimbs}}, where slimbs is one of 2\,\mathit{nlimbs} - 1, 2\,\mathit{nlimbs}, 2\,\mathit{nlimbs} + 1. Squaring is detected when f == g and flen == glen.

Sets (\mathit{res}, (\mathit{nhi} - \mathit{nlo})\, \mathit{slimbs}) to the coefficients [\mathit{nlo}, \mathit{nhi}) of the product of (f, \mathit{flen}\, \mathit{nlimbs1}) and (g, \mathit{glen}\, \mathit{nlimbs2}), with the conventions of _flint_mpn_poly_mul_karatsuba except that the operands may have different numbers of limbs, with headrooms norm1 and norm2. slimbs must be at least \mathit{nlimbs1} + \mathit{nlimbs2} - 1, and at least \mathit{nlimbs1} + \mathit{nlimbs2} + 1 if either operand may get extended in the recursion (when its headroom is smaller than the recursion depth plus one). Squaring is detected when f == g, flen == glen and nlimbs1 == nlimbs2.

Chooses packing parameters that are correct and (as far as the tuning goes) optimal for computing the coefficients [\mathit{nlo}, \mathit{nhi}) of the product of polynomials of lengths len1 and len2 with _flint_mpn_poly_mulmid, given bounds in the format of _fmpz_vec_max_bits: the coefficients have at most |\mathit{bits1}| and |\mathit{bits2}| bits, and a negative bound indicates that negative coefficients may be present. squaring should be set if the operands are the same.

Returns the number of output limbs required by _flint_mpn_poly_mulmid for the given packed operands when the dot products have at most minlen terms and the Karatsuba cutoff is cutoff: the bound for the classical algorithm, or, if Karatsuba is used, the bound including a sign bit and the growth in the recursion, and \mathit{nlimbs1} + \mathit{nlimbs2} + 1 if an operand may get extended. Passing WORD_MAX as the cutoff gives the classical bound.

The Karatsuba cutoff (shorter length) used by _flint_mpn_poly_mulmid for coefficients of n_1 and n_2 limbs. _flint_mpn_poly_karatsuba_cutoff is the tuned value, from 32 for one-limb products down to 4 for large coefficients (64 down to 8 for squaring); _flint_mpn_poly_mulmid_cutoff returns it unless overridden with the global _flint_mpn_poly_mulmid_force_cutoff.

Whether the Karatsuba recursion may need to extend an operand with the given headroom by a limb, and whether _flint_mpn_poly_mulmid uses Karatsuba for dot products of minlen terms: never for sign-magnitude coefficients, and with tight headroom only from 1.5 times the cutoff.

Set (\mathit{res}, (\mathit{nhi} - \mathit{nlo})\, \mathit{slimbs}) to the coefficients [\mathit{nlo}, \mathit{nhi}) of the product (respectively the full product) of (f, \mathit{flen}) and (g, \mathit{glen}), packed as described by method (a biased operand is treated as unsigned), choosing between _flint_mpn_poly_mulmid_classical (always for sign-magnitude coefficients) and _flint_mpn_poly_mulmid_karatsuba, after discarding the input coefficients that cannot contribute to the window. slimbs must be at least _flint_mpn_poly_mulmid_slimbs (as computed by _flint_mpn_poly_mulmid_params). Squaring is detected when f == g, flen == glen and nlimbs1 == nlimbs2.

Euclidean division: sets (q, \mathit{an} - \mathit{bn} + 1) to \lfloor a/b \rfloor and (r, \mathit{bn}) to a - qb. The quotient must not overlap the inputs. As with GMP's mpn_tdiv_qr, r may be equal to a (the remainder then overwrites the low \mathit{bn} limbs of a) but must not overlap b. These replace GMP's mpn_tdiv_qr and are faster for all sizes.

Ceiling division: sets (q, \mathit{an} - \mathit{bn} + 2) to \lceil a/b \rceil and (r, \mathit{bn}) to qb - a \in [0, b). The quotient can overflow \mathit{an} - \mathit{bn} + 1 limbs (for example for a = B^2 - 1, b = B), so one extra limb, which is 0 or 1, is written. Aliasing as for flint_mpn_tdiv_qr.

Division rounding to nearest, ties to even: sets (q, \mathit{an} - \mathit{bn} + 2) to the integer nearest to a/b and (r, \mathit{bn}) to |a - qb| \le b/2, and returns the sign (-1, 0 or 1) of a - qb. As for ceiling division, one extra quotient limb (0 or 1) is written.

If b divides a, sets (q, \mathit{an} - \mathit{bn} + 1) to the quotient and returns 1; otherwise returns 0, the contents of q being undefined. This is a version of flint_mpn_divides based on flint_mpn_tdiv_qr.

If (\mathit{arrayg}, \mathit{limbsg}) divides (\mathit{array1}, \mathit{limbs1}), sets (q, \mathit{limbs1} - \mathit{limbsg} + 1) to the quotient and returns 1; otherwise returns 0. Requires \mathit{limbs1} \ge \mathit{limbsg} > 0 and a nonzero top limb of the divisor. temp must have room for limbsg limbs.

Approximate quotient: sets (q, \mathit{an} - \mathit{bn} + 1) to \lfloor a/b \rfloor or \lfloor a/b \rfloor + 1, like GMP's internal mpn_divappr_q. No aliasing is permitted. This follows the algorithm selection of flint_mpn_tdiv_q, so that it is never slower, but skips the final corrections.

Sets (q, \mathit{an} - \mathit{bn} + 1 + f) to \lfloor a B^f / b \rfloor or \lfloor a B^f / b \rfloor + 1, the approximate quotient with f fraction limbs. A negative f (with \mathit{an} - \mathit{bn} + 1 + f \ge 1) gives the top limbs of the quotient.

Correctly truncated reciprocal: sets (q, n - \mathit{xn} + 2) to \lfloor B^n / x \rfloor. Requires x_{xn-1} \ne 0 and n \ge \mathit{xn}. The top output limb is nonzero only when x is a power of B.

Approximate reciprocal with the interface of flint_mpn_inv: sets (q, n - \mathit{xn} + 2) to \lfloor B^n / x \rfloor or \lfloor B^n / x \rfloor + 1. It uses the algorithms of flint_mpn_inv (with the same cutoffs) without the correction steps, and above the Newton cutoffs _mp_real_inv_newton with one guard limb instead of three, so that it is never slower.

Exact division: sets (q, \mathit{an} - \mathit{bn} + 1) to a/b, assuming that b divides a (otherwise the result is meaningless).

Initialise pre for exact division by (b, \mathit{bn}) (with b_{bn-1} \ne 0), and free it. The init_prec version computes the inverse only to \min(\mathit{bn}' + 1, \mathit{qn}) limbs, which is much cheaper when all quotients are known to have at most \mathit{qn} limbs; longer quotients remain correct but are not accelerated.

Exact division by the fixed divisor of pre: sets (q, \mathit{an} - \mathit{bn} + 1) to a/b, assuming that b divides a. Quotients of at most \mathit{bn}' + 1 limbs cost a single flint_mpn_mullow_n; longer ones use block Hensel division with the stored inverse (one low and one high product per block), except for divisors b' of at most 7 limbs, where the register-based code of flint_mpn_divexact is faster. Only the low limbs of a that can influence the quotient are read.

Returns 1 if b divides a and 0 otherwise. Requires \mathit{bn} \ge 1 and b_{bn-1} \ne 0; a may have zero top limbs and \mathit{an} may be zero or less than \mathit{bn}. Used by fmpz_divisible.

Returns 1 if (x, \mathit{xsize}) is divisible by the odd limb d \ge 3 and 0 otherwise, using GMP's mpn_modexact_1_odd where available.

Returns a value below 2^{50} congruent to (a, n) modulo 2^{48} - 1, the quantity computed by GMP's internal mpn_mod_34lsub1 (which is called for inputs shorter than 48 limbs when available). Longer inputs use an AVX2 kernel when compiled with AVX2 support, or a portable carry-chain loop. Only available for 64-bit limbs.

Returns the 3/2 inverse \lfloor (B^3 - 1) / (d_1 B + d_0) \rfloor - B of the leading two limbs d_1, d_0 of a divisor, for use by the preinv1 functions. Requires d_1 \ge B/2. It is computed inline from the limb inverse of d_1, as GMP's invert_pi1.

Divides (a, m) by the normalised (b, n), returning the top limb of the quotient (0 or 1), writing the remaining m - n quotient limbs to q and the remainder in place to (a, n). Requires m \ge n \ge 2 and dinv = flint_mpn_preinv1(b[n - 1], b[n - 2]). This is _flint_mpn_divrem_preinv1 (schoolbook or divide and conquer division) with scratch space allocated as needed.

Divides (\mathit{up}, n) by the limb d, writing the quotient to (\mathit{qp}, n) and returning the remainder. Requires n \ge 1, d \ne 0, norm the number of leading zero bits of d and dinv = n_preinvert_limb(d). qp may alias up. This is equivalent to mpn_divrem_1(qp, 0, up, n, d) and faster for small n; for large n GMP's assembly is usually faster.

flint_mpn_divrem_1_preinv for dividends of exactly 2 and 3 limbs. The norm versions require a normalised divisor, the unnorm versions an unnormalised one (\mathit{norm} > 0).

Sets (\mathit{dinv}, n) to the precomputed inverse \lfloor B^{2n} / (d + 1) \rfloor - B^n of the normalised n-limb integer d (zero if d = B^n - 1), for use by the preinvn functions (multiword modular multiplication).

Sets (r, n) to a \bmod d, where (a, m) has m \ge n limbs and the normalised (d, n) has the inverse dinv from flint_mpn_preinvn. No aliasing of r with a is permitted. This is not always as fast as ordinary division.

Sets (q, m - n) to the quotient and (r, n) to the remainder of (a, m) by the normalised (d, n) with inverse dinv, returning the top quotient limb (0 or 1). Requires m \ge n; no aliasing is permitted. This is used by flint_mpn_tdiv_qr for short divisors with long dividends.

Sets (\mathit{res}, n) to x^{-1} \bmod B^n, for odd x. Newton iteration on middle products; up to FLINT_MPN_BINV_MULLOW_CUTOFF limbs full low products are used instead, and the cases n \le 3 are inlined.

Hensel division: for odd (b, \mathit{bn}), sets (q, n) to the 2-adic quotient with qb \equiv a \pmod{B^n}. If r is not NULL (it is NULL for bdiv_q), the remainder (a - qb)/B^n \bmod B^{\mathit{bn}} is written to (r, \mathit{bn}), so that a \equiv qb + B^n r \pmod{B^{n + \mathit{bn}}}; the division is exact with an n-limb quotient if and only if r = 0 and a has no nonzero limbs at positions \ge n + \mathit{bn}. q may alias a. Unlike GMP's mpn_bdiv_qr, the number n of quotient limbs is a parameter independent of \mathit{an} and \mathit{bn}.

The algorithms behind flint_mpn_bdiv_qr, with its interface: a one-limb divisor; the block algorithm running from the bottom up in blocks of \mathit{bn} limbs, each costing one low and one high \mathit{bn} \times \mathit{bn} product (used by the dispatcher when \mathit{bn} < n/2; the preinv version takes b^{-1} \bmod B^{\mathit{bn}}); and the Karp–Markstein form, which uses a half-precision inverse from flint_mpn_binv and one refinement step [KM1997].

If x \equiv 1 \pmod 8, set (\mathit{res}, n) to the 2-adic reciprocal square root y (with x y^2 \equiv 1 \pmod{B^n}), respectively the square root s (with s^2 \equiv x \pmod{B^n}), that is congruent to 1 modulo 4, and return 1. Otherwise return 0 without writing anything (an odd x is a square modulo B^n if and only if x \equiv 1 \pmod 8).

Sets (s, \lceil \mathit{an}/2 \rceil) to \lfloor \sqrt a \rfloor. If r is not NULL, the remainder a - s^2 \le 2s is written to (r, \lceil \mathit{an}/2 \rceil + 1), zero padded, and its number of limbs is returned; r must have room for \max(\mathit{an}, 2) limbs. If r is NULL, the return value is 0 if a is a perfect square and 1 otherwise (GMP's convention). Requires \mathit{an} \ge 1 and a_{an-1} \ne 0.

Checked exact square root, the analogue of flint_mpn_div: if a is a perfect square, sets (s, \lceil \mathit{an}/2 \rceil) to its square root and returns 1; otherwise returns 0, the contents of s being undefined. Requires a_{an-1} \ne 0.

Returns 1 if (a, \mathit{an}) is a perfect square and 0 otherwise, using the screening of flint_mpn_sqrt and allocating space for the root only when the screening passes. Requires a_{an-1} \ne 0. Used by fmpz_is_square.

The algorithms behind flint_mpn_sqrtrem, flint_mpn_sqrt and flint_mpn_is_square, with their interfaces: the dispatcher, the divide and conquer algorithm (requires \mathit{an} \ge 5 and needs no remainder space beyond \lceil \mathit{an}/2 \rceil + 1 limbs), the Newton algorithm, GMP's mpn_sqrtrem, and the checked root and square test after the test modulo 256.

Returns an upper bound for the number of limbs of x^e, including one limb of slack for the intermediate products of flint_mpn_pow. Requires x_{xn-1} \ne 0. For small results the bound e \cdot \mathrm{bits}(x) is used; otherwise e \log_2 x is evaluated in double precision from the top two limbs and rounded up, which is essentially tight also for small bases with large exponents.

Sets \mathit{res} to x^e by left-to-right binary exponentiation and returns the number of limbs of the result. res must have room for flint_mpn_pow_bound_limbs(x, xn, e) limbs and must not overlap x. Requires \mathit{xn} \ge 1 and x_{xn-1} \ne 0; e = 0 gives 1.

Drop-in replacements for the GMP functions of the same names, with identical semantics for signs, rounding and aliasing, and the limb-level work done by flint_mpn_tdiv_qr, flint_mpn_divexact and flint_mpn_sqrtrem. They are used by the fmpz division and square root functions for large operands.

Sets (r, n) to ab \bmod d, where d is the normalised modulus (d, n) with dinv = flint_mpn_preinv1(d[n - 1], d[n - 2]) (for n = 1, flint_mpn_preinv1(d[0], 0), which equals n_preinvert_limb(d[0])). If a, b and d have all been shifted left by norm bits so that d is normalised, the result carries the same shift. The product is reduced by _flint_mpn_divrem_preinv1.

Sets (r, n) to ab \bmod d for the normalised (d, n) with inverse dinv from flint_mpn_preinvn, with the shift convention of flint_mpn_mulmod_preinv1.

flint_mpn_mulmod_preinvn for n = 2, except that a, b and the result are not shifted: only d holds the shifted modulus.

Sets (\mathit{res}, n) to a^e \bmod d, where the exponent is the nonnegative (e, \mathit{en}) (not necessarily normalised), with the conventions of flint_mpn_mulmod_preinvn (a and the result are shifted, the exponent is not). Requires the unshifted modulus to be at least 2. res must not alias a or e.

Set (r, n) to ab + ef \bmod d with a single reduction, given d holding d 2^{\mathrm{norm}} and its inverse dinv; the operands and the result are unshifted. The second version is for n = 2.

Given a modulus of n limbs and shift norm, returns the best method for num multiplications by a fixed operand, accounting for the cost of the precomputation: MPN_MULMOD_PRECOND_NONE (use flint_mpn_mulmod_preinvn), MPN_MULMOD_PRECOND_SHOUP (use flint_mpn_mulmod_precond_shoup) or MPN_MULMOD_PRECOND_MATRIX (use flint_mpn_mulmod_precond_matrix).

Multiplication by a fixed operand 0 \le a < d with Shoup's method (Multiplication by a fixed operand (Shoup)): the precomputation writes the n-limb quotient \lfloor a B^n / d \rfloor to apre (the modulus given as dnormed = d 2^{\mathrm{norm}} with its inverse); the multiplication sets (\mathit{res}, n) to ab \bmod d with a high product and a low product, for b reduced modulo d. Here d is the unshifted modulus.

The matrix method of preconditioned multiplication (multiword modular multiplication): the precomputation writes a 2^{\mathrm{norm}} B^i \bmod d 2^{\mathrm{norm}} for 0 \le i < n (plus one junk limb) to apre, which must have room for flint_mpn_mulmod_precond_matrix_alloc(n) limbs; ab \bmod d is then a matrix-vector product with the limbs of b followed by one short reduction, in O(n^2) operations. The fmma version computes a_1 b_1 + a_2 b_2 \bmod d. Requires n \ge 2 and b reduced modulo d.

Multiplication modulo 2^b + 1 for the representation of the Schönhage–Strassen FFT (see fft.h): the inputs are y + c_1 2^b and z + c_0 2^b, where (y, n) and (z, n) with n = \lceil b/\mathtt{FLINT\_BITS} \rceil are reduced modulo 2^b and c holds the top bits (c_1 = 2 bit, c_0 = 1 bit), each input being at most 2^b. Sets (x, n) and returns the top bit of the product modulo 2^b + 1, fully reduced. tp must have room for 2n limbs.

Returns 1 if a is a square modulo the odd prime d, and 0 if it is not. This evaluates the Jacobi symbol (a/d) with GMP, in quasi-linear time; a return value of 0 means that the symbol is -1, so it also proves that a is not a square for an odd d that is not prime, while a return value of 1 proves nothing in that case.

Sets (\mathit{res}, n) to a square root of a modulo the odd prime d and returns 1. Returns 0, setting res to zero, if the Jacobi symbol rules a root out (which is conclusive whether or not d is prime), and -1, again setting res to zero, if the algorithm itself failed, which happens only for d not prime (or, with negligible probability, for a prime whose least quadratic nonresidue exceeds 2^{20}). The preinv version takes the inverse dinv of d 2^{\mathrm{norm}} from flint_mpn_preinvn and \mathrm{norm} = \mathrm{clz}(d_{n-1}), as an mpn_mod or fmpz_mod context holds them; dinv may be NULL.

Set (\mathit{arrayg}, \ell) to the greatest common divisor of (\mathit{array1}, \mathit{limbs1}) and (\mathit{array2}, \mathit{limbs2}) and return its length \ell. Unlike GMP's mpn_gcd, the inputs need not be odd or ordered by size, and they are not destroyed; they must be nonzero with nonzero top limbs. temp must provide \mathit{limbs1} + \mathit{limbs2} limbs of scratch space, or be NULL. The common power of 2 is removed first and GMP's mpn_gcd is called on the odd parts.

Divides (x, \mathit{xsize}) in place by 2^v, where v is the number of trailing zero bits, stores v in bits and returns the new length. x must be nonzero.

Divides (x, \mathit{xsize}) in place by the largest power p^v of (p, \mathit{psize}) that divides it, stores v in exp and returns the new length. x must be nonzero and p > 2. The divisibility is tested by ascending squares p, p^2, p^4, \ldots, which is efficient for large powers but has a high overhead, so it should not be used as the first stage of trial division.

Searches for a factor of (x, \mathit{xsize}) among the primes flint_primes[i] for start \le i < stop. Returns the first such i for which the prime divides x, and 0 if there is none. Requires start \ge 1, so that 0 is never a valid index.

Searches for factors of (x, \mathit{xsize}) among roughly the first num_primes primes, by GCDs with products of primes organised in a product tree; writes the indices in flint_primes of the primes found to factors and returns their number. Requires 0 \le num_primes \le 3512. If x fits in a small fmpz, it is factored completely instead. The product tree is computed once and cached globally (in a thread-safe way). This is the trial division used by fmpz_factor.

Initialise C for the given moduli. The version with flags (a bitwise or of FLINT_MPN_CRT_MOD and FLINT_MPN_CRT_CRT) only performs the precomputations needed for the selected operations, which saves time and memory when only one direction is needed; flint_mpn_crt_init selects both. Throws an exception if a modulus is zero, or if Chinese remaindering is selected and the moduli are not pairwise coprime. The precomputation costs O(M(N) \log N) for N total bits and uses O(N \log N) memory.

flint_mpn_crt_init2 with explicit tuning parameters: the approximate number of bits of a CRT basecase chunk, the number of bits at which the modular reduction switches to its basecase, and the number of limbs from which precomputed inverses are used for division. Intended for profiling.

Frees the memory used by C.

Writes the residues of the nonnegative (x, \mathit{xn}) modulo the moduli to out[0], ..., out[num_primes - 1]. The input may have any length, including zero, and need not be normalised. tmp must have room for C->tmp_limbs limbs, or be NULL for internal allocation. Requires C to be initialised for FLINT_MPN_CRT_MOD.

Reconstructs the integer x with the residues res[0], ..., res[num_primes - 1], each reduced modulo its modulus, and writes it as C->prod_len limbs, zero padded, to out. If sign is zero, x is the unique solution with 0 \le x < P and the return value is 0. Otherwise x is the unique solution with -P < 2x \le P, its absolute value is written, and the return value is 1 if x is negative. tmp is as for flint_mpn_multi_mod.

Vector versions. In the mod version, x is a packed array of len nonnegative integers of xn limbs each, and the residue of entry i modulo modulus l is written to out[l * out_stride + i] (so out_stride must be at least len). In the CRT version, the residue of entry i modulo modulus l is read from res[l * res_stride + i], entry i of the output is written to out + i * out_stride (with out_stride \ge C->prod_len), and if sign is nonzero the sign flags are written to negative[i] (negative may be NULL when sign is zero). These are faster than loops over the entries when the entries are small compared to P, since the precomputed tables are then traversed once per block of entries.

One-shot versions of flint_mpn_multi_mod and flint_mpn_multi_crt that need no precomputed structure. In the CRT version, out must have room for num_primes limbs; the length of the product P is written to outn, and if prod is not NULL, P itself is written there (with room for num_primes limbs). The return value is as for flint_mpn_multi_crt.

Returns the nonnegative (a, \mathit{an}) reduced modulo the j-th leaf modulus of C (the product of one or more consecutive batched moduli; leaf j is modulus j when all moduli exceed 32 bits). Requires an \le C->mod_pow_limbs. This exposes the basecase of flint_mpn_multi_mod; with a single modulus it is a faster replacement for mpn_mod_1 when many reductions with the same modulus are needed.

Sets (\mathit{rp}, \lceil n / \mathtt{FLINT\_BITS} \rceil) to a uniformly random integer in [0, 2^n).

Sets (\mathit{rp}, \mathit{xn}) to a uniformly random integer in [0, x). Requires x_{xn-1} \ne 0. All \mathit{xn} limbs are written even if the value is shorter.

Set (\mathit{rp}, n), respectively (\mathit{rp}, \lceil n/\mathtt{FLINT\_BITS} \rceil) with an n-bit value, to a random integer with long runs of zeros and ones in its binary representation. Such numbers trigger corner cases (carry propagation, near-overflow) more often than uniform ones and are used in tests.

The dispatcher behind flint_mpn_tdiv_qr and the individual algorithms for large operands, with its interface; in all of them r may be NULL. The Newton version views a and b as fixed-point numbers and computes m + 2 fraction limbs of a/b, where m = \mathit{an} - \mathit{bn} + 1; since the error is below 4 B^{-2} at the integer scale, the integer part is certified whenever the first fraction limb lies in [2, B-2], and the remainder then follows from a low product; otherwise the candidate quotient is corrected by O(1) steps using a full product. The unbalanced version (requires \mathit{an} > 2\mathit{bn} and \mathit{bn} \ge 3) divides blockwise, sharing one approximate inverse of b between the 2\mathit{bn} \times \mathit{bn} blocks. The preinv version requires (\mathit{binv}, \mathit{binvn} + 2) to be the output of _mp_real_inv_newton(binv, b, bn, binvn) with \mathit{binvn} \ge m + 2 and \mathit{an} \ge m + 2. The preinvn version uses flint_mpn_preinvn and flint_mpn_divrem_preinvn.

Euclidean division for 1 \le \mathit{bn} \le FLINT_MPN_DIV_SMALL_BN (7), with the conventions of flint_mpn_tdiv_qr (r may be NULL; no normalisation of b is needed). There is one code path per divisor length, in which the normalised divisor and the partial remainder are kept in registers and the dividend is normalised on the fly; each quotient limb takes a 3/2 division (2/1 for \mathit{bn} = 1) and, for \mathit{bn} \ge 3, the subtraction of a 1 \times (\mathit{bn} - 2) product, as in GMP's mpn_sbpi1_div_qr but without function calls or memory traffic. The cost is linear in \mathit{an}; beyond about 3\mathit{bn} limbs the schoolbook division is faster for \mathit{bn} \ge 4.

Division with a normalised divisor and the 3/2 inverse dinv = flint_mpn_preinv1 of its top two limbs, ported from GMP's internal mpn_sbpi1_*, mpn_dcpi1_*, mpn_tdiv_qr and mpn_div_q functions, with FLINT's multiplication in the divide and conquer steps [BurZie1998].

flint_mpn_inv without Newton inversion. It is used by _mp_real_inv_newton_basecase, so that the mutual calls between the Newton inversion and the division only ever go to lower precisions, whatever the cutoffs.

The part of flint_mpn_divexact for divisors of at least two limbs, and the version that always uses the Newton Hensel division flint_mpn_bdiv_q.

The parts of flint_mpn_divisible, for b = 2^v b' with b' odd, given that 2^v divides a. The small version (b' of at most FLINT_MPN_DIVEXACT_SMALL_BN limbs, \mathit{bn} \ge 2, v < 64, a_{an-1} \ne 0) streams the shifted a through a Hensel remainder computation with b' and the running remainder held in registers; b' divides a if and only if the final remainder is zero.