Modular arithmetic in floating point

Vector units multiply doubles much faster than they multiply 64-bit integers: AVX2 has no instruction for a full 64 \times 64-bit product, while it performs four double-precision fused multiply-adds per instruction. FLINT therefore does the modular arithmetic of its small-prime FFT (fft_small) on integers stored exactly in doubles. The functions are vec4d_mulmod, vec4d_reduce_to_pm1n and their relatives in machine_vectors.h, and the ranges stated there are proved below. The approach is of the kind studied by van der Hoeven, Lecerf and Quintin [HLQ2016]; exact linear algebra libraries such as FFLAS have long relied on integer computations in floating point [DGP2008].

Setting#

Doubles have a 53-bit significand. Let u = 2^{-53} be the unit roundoff and write \mathrm{fl}(x) for x rounded to the nearest double, with ties to even. We assume round to nearest throughout and ignore overflow and underflow, which cannot occur with the magnitudes below. The standard facts we need [MBD2018] are:

• Every integer x with |x| \le 2^{53} is a double. Consequently, an operation whose exact result is such an integer returns it exactly.
• For any real x in range, \mathrm{fl}(x) = x(1 + \delta) with |\delta| \le u.
• If h = \mathrm{fl}(ab) for doubles a and b, then ab - h is a double, so a fused multiply-add computes it exactly: \mathrm{fma}(a, b, -h) = ab - h.

Fix an integer modulus n \ge 1 and let \mathit{ninv} = \mathrm{fl}(1/n) = (1 + \delta_0)/n with |\delta_0| \le u. Residues are represented by integers in doubles, not necessarily reduced: a typical range is (-n, n), and intermediate sums may grow to a few multiples of n.

Rounding the quotient#

Every function computes an approximate quotient q \approx x/n and returns x - qn. The quotient is an integer obtained by rounding the product of a double h with \mathit{ninv}. Two ways are used. The mulmod and nmulmod functions, and all functions on NEON, always use native rounding; round_mul, and hence mulmod_fast and the reductions on AVX2, use the method selected per target by FLINT_FFT_SMALL_ROUND_USE_NATIVE:

• Native rounding: q = \mathrm{round}(\mathrm{fl}(h \cdot \mathit{ninv})) with a rounding instruction (vroundpd on x86, frintn on ARM). This rounds twice.
• Magic constant: with R = 3 \cdot 2^{51}, q = \mathrm{fl}(\mathrm{fl}(h \cdot \mathit{ninv} + R) - R), where the inner operation is one fused multiply-add. On Intel cores since Skylake the rounding instruction is slow, and this is faster.

For |x| > 2^{51} the lemma fails: the spacing of the doubles near x + R exceeds 1 when x > 2^{51}, and is below 1 when x < -2^{51}, in which case s - R need not even be an integer. Native rounding has no such limit. The difference between the two methods is in the error of the quotient:

Reduction#

This is vec4d_reduce_to_pm1n (and the identical vec4d_reduce_to_pm1no). Following it with "add n if negative" gives vec4d_reduce_to_0n.

The condition n \ge 4 is needed by the magic constant: for n = 2 or 3 and |a| close to 2^{53} the quotient exceeds 2^{51} and the result is wrong. It also gives |r| < n/2 + 2 \le n, so the result lies in (-n, n), inside the range [-n, n] promised by the name pm1n.

Multiplication#

vec4d_mulmod rounds natively in step 2 and vec4d_mulmod_fast uses round_mul, which is the magic constant on most Intel targets. The nmulmod variants compute -r the same way: they set \ell = \mathrm{fma}(-a, b, h) = h - ab and return \ell - \mathrm{fma}(-q, n, h).

The theorem holds for any n < 2^{50}, prime or not, odd or even. The ranges used by the transforms follow from it.

For 2^{49} \le n < 2^{50} the theorem gives only |r| < 1.25\,n for |ab| < 2n^2, and a finer analysis of the rounding in step 2, in the comments of fft_small/mulmod_satisfies_bounds.c, gives 9n/8. FLINT therefore checks each 50-bit prime individually with fft_small_mulmod_satisfies_bounds, which evaluates that finer bound with the actual error |n \cdot \mathit{ninv} - 1| in place of its maximum u. A prime passes if the bound gives |r| < n for |ab| < 2n^2 and |r| < 3n/2 for |ab| < 4n^2, and sd_fft_ctx_init_prime throws for a prime that fails. According to the same comments about one random 50-bit prime in 900 fails, and none of 49 bits or less, in agreement with Corollary 5.

Cheap reductions#

Between multiplications, values are brought back into a smaller range by comparisons and conditional additions, which are exact on integers below 2^{53}:

The vectorised versions use blendv on a sign bit or a comparison mask instead of branches. Testing the sign bit treats -0.0 as negative, so vec4d_reduce_pm1no_to_0n on AVX2 maps -0.0 to n, outside the stated output range; the result is still congruent to the input.

Use in FLINT#

• Small-prime FFT. fft_small computes number-theoretic transforms modulo primes p < 2^{50} with p \equiv 1 \pmod{2^k}, with all arithmetic in vec8d and vec4d. Twiddle factors are stored in (-p/2, p/2), the butterflies use mulmod_fast and reduce_to_pm1n, and the bound analysis in fft_small/mulmod_satisfies_bounds.c shows that the values stay within a few multiples of p, so that all sums are exact. That analysis needs the actual ranges of Theorem 3 and Theorem 4 in addition to the documented ones for the truncated inverse transform; it is carried out in Arithmetic and its ranges.
• Polynomial and integer multiplication. The conversions in and out of the transform use mulmod, reduce_to_0n and reduce_0n_to_pmhn. Through fft_small, this arithmetic underlies the multiplication of large integers (flint_mpn_mul) and of large polynomials in nmod_poly and fmpz_poly.
• Precomputed transforms. mpn_mod and the gr implementations of transformed polynomials multiply by fixed scalars with vec8d_mulmod.
• Matrices. nmod_mat_mul_u8 reduces its results with vec4d_reduce_to_0n on x86. The floating-point kernel nmod_mat_mul_fp50 uses its own versions of round_mul, mulmod_fast and reduce_to_pm1n with a fixed lane count per instruction set (nmod_mat/mul_fp_vec.h), with the magic constant except on NEON, and refers to the same analysis.