Bivariate resultants over nmod via geometric multipoint evaluation

[maelhos]

Following Issue #2605, this PR adds _gr_poly_resultant_multipoint / gr_poly_resultant_multipoint, a specialization for the resultant for bivariate polynomials, i.e. for the case where the coefficient ring is itself a polynomial ring. The coefficients of both inputs are evaluated at a geometric progression of points, the resultants of the resulting univariate polynomials are computed, and the result is interpolated back, using the geometric evaluation/interpolation added in #2449.

Currently only a base ring of nmod (word-size prime modulus) is supported; fmpz, fmpq and fmpz_mod can follow later (see discussion in the related issue). _gr_poly_resultant dispatches to it whenever it applies, replacing the Sylvester determinant that was previously used for these rings.

Follows the algorithm in PML's lzz_pXY, with the blockwise evaluation of @AntoineBak fork so that memory doesnt grow as n^3.

Two further refinements over PML: the evaluation points are scaled by a random constant (the FLINT geometric progression always starts at 1, so retrying alone would not avoid a leading coefficient vanishing at 1), and only the leading coefficient of the first argument is required not to vanish at the evaluation points.

The timing is unsurprisingly much better than before (on my small Zen 3 laptop)

deg_y	deg_x	sylvester	subres	new	syl/new	sub/new
4	4	0.00013	2.3e-05	1.59e-05	8.2	1.4
4	16	0.000744	0.000104	7.89e-05	9.4	1.3
4	64	0.00313	0.000475	0.000549	5.7	0.9
4	256	0.0107	0.00226	0.00221	4.8	1.0
8	4	0.00403	0.00016	5.33e-05	75.6	3.0
8	16	0.0203	0.000704	0.000495	41.0	1.4
8	64	0.078	0.00328	0.00191	40.8	1.7
8	256	0.285	0.0141	0.00791	36.0	1.8
16	4	0.106	0.00121	0.000431	245.9	2.8
16	16	0.599	0.00514	0.00188	318.6	2.7
16	64	2.1	0.0242	0.00738	284.6	3.3
16	256	-	0.099	0.0309	-	3.2
32	4	3.01	0.00959	0.00204	1477.5	4.7
32	16	17.3	0.0425	0.00799	2162.8	5.3
32	64	-	0.184	0.0313	-	5.9
32	256	-	0.782	0.127	-	6.2
64	4	89.6	0.0748	0.00871	10292.1	8.6
64	16	-	0.338	0.0344	-	9.8
64	64	-	1.44	0.14	-	10.3
64	256	-	6.25	0.58	-	10.8
128	4	-	0.638	0.0419	-	15.2
128	16	-	2.77	0.168	-	16.5
128	64	-	11.9	0.701	-	16.9
128	256	-	-	2.89	-	-
256	4	-	5.34	0.229	-	23.3
256	16	-	22.8	0.952	-	24.0
256	64	-	-	3.84	-	-
256	256	-	-	15.5	-	-

This PR is partially made using Claude Opus 5.

[vneiger]

Nice!

Concerning the evaluation at $c \cdot q^i$, could this maybe be done with several successive extrapolations?
https://flintlib.org/doc/nmod_poly.html#extrapolation
The first batch, done by evaluation, gives you a number batch of evaluations at the first powers of q (and batch is at least the polynomial length), so you could extrapolate with the right offset to get the evaluations at the points batch ... 2*batch-1, then extrapolate the latter to get the evaluations 2*batch ... 3*batch-1, etc. This would at least simplify the code a bit (no need to compute some specific powers, and no need to scale before evaluation), and "in theory" this should not impact the timings negatively but this needs to be checked.

[vneiger]

(note that this adds a dependency that is not wanted in case you plan to introduce multithreading in the evaluation step)

[maelhos]

Overall, I made a few tests, and it's about the same performance, though extrapolation is slightly worse in theory (it takes 2L + 2B scalar muls vs. 2L + B for the current version) for L, the total length, and B, the batch length. It also loses the opportunity of parallelization and seems to also be slightly worse for cache / branching. It's still not bad and makes the code easier to read for sure. I'd also argue that since we want to port this code to fmpz_mod and extrapolate is not available for that, we might want to keep it consistent. But overall I'm not opposed to it :)

[maelhos]

Guess crediting Vincent and Antoine now breaks the CI... windows is interesting

[vneiger]

Overall, I made a few tests, and it's about the same performance, though extrapolation is slightly worse in theory (it takes 2L + 2B scalar muls vs. 2L + B for the current version) for L, the total length, and B, the batch length.

Ok, yes (these are negligible compared to the mulmid, except when the length is really small).

It also loses the opportunity of parallelization

Giving it a second thought, instead of doing the extrapolations with the same offset one after another, we may as well use different offsets B, 2*B, 3*B etc. always from the initial evaluations, and then we can parallelize efficiently.

and seems to also be slightly worse for cache / branching.

Agreed.

It's still not bad and makes the code easier to read for sure. I'd also argue that since we want to port this code to fmpz_mod and extrapolate is not available for that, we might want to keep it consistent. But overall I'm not opposed to it :)

Ok, thanks for having investigated the extrapolation version. I agree that consistency is a plus. There's no clear winner it seems, so, just pick the version you prefer!

For the rest, I'll read the code more carefully in the next few days.

[maelhos]

Ok, then I'll keep the current version and finish my prototype of fmpz and fmpq version. Though fmpz_mod will be for another PR as it requires geometric multipoint over fmpz_mod. So far I have very good stats (x486 lgtm) for my prototype over fmpz and fmpq...

$\mathbb Z$:

deg_y	deg_x	bits	subres	new	speedup	notes
4	4	8	4.42e-05	4.44e-05	1.0	below cutoff
4	4	1024	0.00675	0.00703	1.0	below cutoff
4	16	64	0.00257	0.00181	1.4
4	16	1024	0.0682	0.0268	2.5
4	64	8	0.00478	0.00174	2.7
4	64	1024	1.34	0.123	10.9
8	4	8	0.00113	0.000445	2.5
8	4	1024	0.184	0.0421	4.4
8	16	1024	1.22	0.216	5.6
8	64	1024	7.3	0.802	9.1
16	4	64	0.135	0.0195	6.9
16	4	1024	4.17	0.305	13.7
16	16	256	4.53	0.385	11.8
16	64	64	5.08	0.377	13.5
32	4	64	3.11	0.197	15.8
32	16	8	2.3	0.147	15.6
32	64	1024	-	50.2	-

$\mathbb Q$:

deg_y	deg_x	bits	subres	new	speedup	notes
4	4	64	0.0207	0.021	1.0	below cutoff
4	16	8	0.031	0.00277	11.2
4	16	256	11.1	0.425	26.2
4	64	8	9.82	0.0215	456.8
4	64	64	452	0.93	486.0
8	4	64	1.6	0.0383	41.9
8	4	256	21.9	0.304	72.0
8	16	8	0.535	0.0212	25.2
8	64	8	17.4	0.162	107.6

[maelhos]

I'm trying to use FFT-primes for multimodular fmpz to replace geometric multipoint by plain DFT which is sd_fft_trunc and the improvement is very mild... I'll also make the check for nmod version that if prime is FFT we should use sd_fft / fft_small (depending on range) instead of geometric multipoint.

@vneiger maybe we should make some kind of like "whatever multipoint" that depending on the primes either gives dft points (the right one) or geometric or in last resort subproduct tree. This could be usefull for some other things no ?

More generally we could implement this algorithm for generic dense $n$-variate polynomials, we can use Kronecker subsitution to get a complexity of $\mathcal O(d^{2n-1})$

[fredrik-johansson]

Noticed this is a comment:

(currently the Sylvester determinant, since polynomial rings do not report being a UFD)

Normally they do, but e.g. for nmod rings you may have to call gr_ctx_set_is_field to let them know that the modulus is prime. With this the generic dispatcher should already call the subresultant algorithm (or else we need to fix that).

[maelhos]

is there alrady a function that gives fft_primes ?

[maelhos]

A few things:

For the FFT-friendly idea: yes, but no. Basically, for nmod we have a x1.1-x1.2 speedup by using DFT if the prime is FFT.
BUT we are restricted to 50-bit primes, so in practice for multimodular fmpz doing a geometric multipoint with 62-bit primes actually wins; the FFT is ~0.9 times slower on average for input where the coeff in Z are ~100-200 bits.

Normally they do, but e.g. for nmod rings you may have to call gr_ctx_set_is_field to let them know that the modulus is prime. With this the generic dispatcher should already call the subresultant algorithm (or else we need to fix that).

It seems like it wasn't reported as far as I understand; see a97b017
The problem did not seem like it did not know that the modulus was prime but rather was the handler for the UFD function that was missing.