Repack to half-length fft_small convolution for tiny moduli in _nmod_poly_mul and mullow

[fredrik-johansson]

When the unreduced coefficients of an _nmod_poly product are small enough to fit in 50 bits, fft_small multiplies using a single FFT prime.

We can do even better when the moduli are really small by packing a linear polynomial into each coefficient. If the coefficients of the product (over $\mathbb{Z}$) are smaller than $M$, we can repack $a_0 + a_1 x + \ldots$ as $(a_0 + a_1 M) + (a_2 + a_3 M) x + \ldots$. The coefficients of the product of such polynomials will be quadratic polynomials in $M$, from which we can read off the coefficients of the original product.

It is clear that we can take $M = 2^{16}$ since $M^3$ is smaller than the 50-bit FFT primes used by fft_small. With a more careful analysis we can show that it is usually possible to work with $M = 2^{17}$.

Where this trick is applicable, currently for moduli $1 \le m \le 23$, it gives up to a 2x speedup by halving the convolution length. Plots demonstrating the speedup for _nmod_poly_mul and _nmod_poly_mullow are attached below.

Input lengths where this trick is applicable are roughly up to 250000 for $m = 2$, 77000 for $m = 3$, 21000 for $m = 5$, 9800 for $m = 7$, 3600 for $m = 11$, 1400 for $m = 17$ and 700 for $m = 23$.

This is definitely a hack: an optimized 32-bit FFT/NTT (or maybe something Schonhage-Strassen-like with many coefficients bit-packed in each word) should perform even better, and would allow much larger moduli and longer products. But this hack is easy to implement with the tools we have in FLINT right now, so we may as well use it.

Example improvement on a nontrivial benchmark problem: constructing GF($5^{3125}$) previously took 5.25 seconds, takes 3.97 seconds with this PR (1.32x speedup).

BTW, a new (to FLINT) trick used in the unpacking code is the 32-bit precomped remainder algorithm by Lemire, Kaser & Kurz which could be useful elsewhere in the nmod modules. This is even faster than the code generated by GCC for remainder by a compile-time constant.

[albinahlback]

Are these tabs architecture specific?

[fredrik-johansson]

Yes, and all the other tuning parameters in this file too. Note that these particular tabs were around before; I just moved them to a new file.

[albinahlback]

Yes, I did notice that. Do we have a tuning program somewhere that we can use at a later point?

[fredrik-johansson]

src/gr_poly/tune/cutoffs.c can generate this kind of table but it's not automatic.

[vneiger]

Nice!

BTW, a new (to FLINT) trick used in the unpacking code is the 32-bit precomped remainder algorithm by Lemire, Kaser & Kurz which could be useful elsewhere in the nmod modules. This is even faster than the code generated by GCC for remainder by a compile-time constant.

Is this different from using "Shoup precomputation" (as in #2061 ) but specialized to the 32-bit context?

[fredrik-johansson]

Nice!

BTW, a new (to FLINT) trick used in the unpacking code is the 32-bit precomped remainder algorithm by Lemire, Kaser & Kurz which could be useful elsewhere in the nmod modules. This is even faster than the code generated by GCC for remainder by a compile-time constant.

Is this different from using "Shoup precomputation" (as in #2061 ) but specialized to the 32-bit context?

Unless I missed something, the Shoup reduction still requires one conditional adjustment, but the Lemire et al. method doesn't: the high part of a product gives the exact remainder right away.