Mpn middle products

[fredrik-johansson]

This PR adds flint_mpn_mulmid with the same contract as radix_mulmid (4 free size parameters an, bn, zlo, zhi; computing only the convolution window [zlo, zhi) without carry-in from lower limbs and discarding carry-out).

This can be used to implement optimized Newton division and square root code for mpn (as already done for radix integers) which should be quite a bit faster than the GMP functions. It could also be used to optimize arb_mul and similar functions (it is hard to base arb_mul directly on flint_mpn_mulhigh_n because inputs can be unbalanced; flint_mpn_mulmid doesn't have this limitation).

The fft_small-based multiplication gains the ability to compute a middle product directly. In principle this means halving the transform length in lucky cases, but in practice this almost never happens due to the block-size padding. So in almost all cases the only speedup comes from skipping the output CRT for unneeded limbs, which is non-negligible.

We also use the new truncating fft_small interface in flint_mpn_mulhigh_n and flint_mpn_mullow_n which speeds up these functions 10-20% for large n.

For smaller inputs the following algorithms are used:

• Schoolbook multiplication (implemented using mpn_mul_1 / mpn_addmul_1 and some extra carry propagation).

• David Harvey's balanced integer Karatsuba mulmid_n which is available in GMP as the private mpn_mulmid_n (or __gmpn_mulmid_n). Since this requires specific balanced lengths, it is used with zero-padding when sizes only roughly correspond. I first tried to port this to FLINT, but it depends on some special-purpose assembly routines and runs 10-20% slower with generic C fallbacks, so the easiest solution was to just call GMP. We could port the code to FLINT properly with all dependencies if the GMP developers decide to remove this internal function in the future.

• For sizes that almost look like a full mul, a mullow_n, or a mulhigh_n, we zero-pad or truncate and call flint_mpn_mul, flint_mpn_mullow_n or flint_mpn_mulhigh_n on temporary buffers.

The automatic dispatch between these algorithms is a little sketchy and wants more fine-tuning, especially for small sizes.

Generated mostly using Opus 4.8.

[albinahlback]

I don't think mulmid is necessary for fast reciprocal square root -- only mullo, sqrlo and mulhi. Essentially, you have something like
$x_{n + 1} = x_{n} + x_{n} \frac{1 - s x_{n}^2}{2}$
(This is normalized -- for a fast implementations you will use mixed precision, so you will have to denormalize with powers of two's)

For each iteration you have a bound for $1 - s x_{n}^2 \ll 1$ so you will only have to compute the low part of $s x_{n}^2$. Then $x_{n} (1 - s x_{n}^2)$ only requires a high multiplication (it is sloppy, but you don't need to be exact apart from the final iteration).

[fredrik-johansson]

The middle product is for the multiplication $s \cdot x_n^2$. For precision $n/2 \to n$ this is an $n \times n$ product where the least significant $n$ output limbs can be discarded and the most significant $n/2$ limbs are known and can be omitted. (So actually a bit different from the balanced middle product, and worth specifically tuning for.)

However, this reminds me that several of our square root functions don't actually use the optimal sequence...

[albinahlback]

I guess I was thinking of the lower cases where you want to avoid sign changes.