Lift the parametrization over QQ p-adically by Newton's iteration by mohabsafey · Pull Request #391 · algebraic-solving/msolve

mohabsafey · 2026-11-06T08:03:35Z

This PR builds on #388 (separating elements and splitting) and on #394 (work of F4 recorded in the trace). Until they are merged, their commits also appear here. The last twenty-two commits belong to this PR.

Over the rationals, msolve recovers the parametrization of the solution set by computing it modulo many primes and recombining the results by CRT. This costs one Gröbner basis and one FGLM run per prime. This PR adds a p-adic alternative. The parametrization modulo a few primes p is lifted to p^N by Newton's iteration on the input system, then recombined by CRT and reconstructed over Q as before. A cost model, fed by the first modular runs, chooses between the two liftings. Outputs are unchanged.

The lifting is the global Newton iterator of Giusti, Lecerf and Salvy [1]. In their Section 4.6, "Lifting the integers", a geometric resolution computed modulo p is lifted to p^k, and its rational coefficients are then recovered. Here, it is applied to the parametrization that msolve computes modulo p by F4 and FGLM. Newton iteration and Hensel lifting are covered in [2], Chapters 9 and 15. The corrections of the last steps follow Dixon's p-adic solution of linear systems [3].

Related software

Kronecker [1, 4]. A Magma package by G. Lecerf, implementing [1]. It computes a geometric resolution modulo a prime and lifts it to Q with the global Newton iterator. Its page was last updated in 2002, and it requires Magma 2.9-2 or later.
TriangularizeModular [5]. A Maple 9.5 implementation on top of the RegularChains library. It computes a triangular decomposition modulo a prime, then lifts the triangular sets by Hensel lifting. It handles square systems that generate radical ideals.
geomsolvex. A C++ library of Mathemagix, which implements the Kronecker solver and the relaxed Hensel lifting of triangular sets of Lebreton [6].
Dixon's lifting [3]. The p-adic solution of nonsingular linear systems over Z by Dixon's method is implemented in IML [8], in LinBox [9] (DixonSolver) and in FLINT [10] (fmpz_mat_solve_dixon, fmpq_mat_solve_dixon). These libraries lift the solution of one linear system over Z from its inverse modulo p. Here, the corrections of the last Newton steps apply the same idea to the linear systems of Newton's iteration over (Z/p^N)[t]/(w), with J⁻¹ kept at a lower precision.

As far as we know, the Gröbner basis solvers have, up to now, recovered the coefficients over Q by multi-modular CRT, as msolve did. This PR may be the first use of the p-adic lifting of the parametrization in an open-source Gröbner basis solver.

What is new

To our knowledge, the following techniques appear neither in the papers above nor in the related software.

p-adic lifting combined with Gröbner bases and CRT. The starting point is the parametrization that F4 and FGLM compute modulo p, with additional data on the per prime cost for these two steps, as well as the memory it uses. The parametrization is lifted at several primes, one per thread, each to a moderate precision p^N. The images are then recombined by CRT. Lifting several primes keeps the integers small and makes the method parallel, whereas a single prime would have to be lifted much further.
A cost model that chooses between the two liftings. It compares the cost of a Newton step with the cost of a modular run, both counted as operations. The work of F4 is read from the trace learned in the first modular run (Record the work of applying the trace when learning it #394), and that of FGLM, or of the dense route of Split the solution set along equal coordinates and find separating elements from the multiplication matrices #388, from its sizes. The same model chooses the precision N and the memory limit. It does not depend on the machine or on the number of threads, so the choice is reproducible. When the predicted gain is low, msolve keeps the multi-modular lifting, so it does not slow down.
Precision driven by the reconstruction. The number of bits of the coefficients over Q is not known in advance (even if bounds exist in the literature).
- The problem. The multi-modular lifting only has to add primes until the reconstruction succeeds. The p-adic lifting must also choose the precision of each prime. A Newton step costs about the same whatever the gain of precision, so steps to slightly higher precisions are expensive. A too high precision may waste time and memory.
- Attempts. The reconstruction is attempted whenever the modulus has grown by 10% since the last attempt. An attempt reconstructs the coefficients of w from the highest degree down and stops at the first failure, so a failed attempt is cheap. Each attempt records the number of bits b of the modulus and the fraction f(b) of the coefficients of w reconstructed.
- Estimate. The bits needed by the coefficients of w grow from the highest degree to a peak near the middle. So f(b) is extrapolated linearly from the last two attempts to f = 1/2, or a little further once that is reached. Once w is reconstructed, the coordinates B_v need a few percent more bits: 4% plus two digits are added.
- Target. The estimate gives the target precision of the primes being lifted. The remaining bits are split over as few batches of primes as possible, each prime being lifted at most to the most efficient precision of the cost model and within the memory limit. The precision is at least doubled at each step: on KLYZ12, the estimate grew from one attempt to the next (p^33, 77, 120, 175, 246, 271, 297), and doubling instead took the run from 841 s to 530 s. If the estimate is too low, the precision grows by a quarter.
- Memory limit. The estimate also enters the memory limit, which counts the residues the multi-modular lifting would hold. It is capped at twice the current number of bits, since the first estimates are often too high (katsura 13: 165000 bits estimated for 20600).
- Reused from Speed up the lifting of rational parametrizations over QQ #384. The reconstruction reuses, on the p-adic images, the techniques of the multi-modular lifting: a running common denominator, restarted for each coordinate when needed, and a denominator of w guessed by vector witnesses (LLL), with a safety margin in the reconstruction bounds.
A cheaper Newton step.
- In the last steps, J⁻¹ is not updated. dX is refined by Dixon-style corrections, whose operands are truncated to the precision of J⁻¹.
- F(X) and J(X) are evaluated through one table of monomials, and the last products are reduced once per entry instead of once per term.
- The products by J and J⁻¹ reuse transforms computed once per matrix.
Overdetermined systems. When there are more equations than variables, Newton's iteration uses a square subsystem of the sparsest equations whose Jacobian is invertible along the parametrization, if there is one. Otherwise it uses random combinations of all the equations.
Split eliminating polynomials. The inversion of J modulo (w, p) forms its pivots from random combinations of rows, which is needed when w splits.

Changes

p-adic lifting of the parametrization (src/msolve/hensel.c). The parametrization modulo p is written x_v = X_v(t) mod w(t), with the separating element t as the last variable. It is lifted by Newton's iteration in (Z/p^N)[t]/(w):
- dX = J⁻¹ F(X) and Δ = t − X_t + dX_t;
- w ← w + (w′ Δ mod w) and X_v ← X_v − dX_v + (X_v′ Δ mod w);
- J⁻¹ ← J⁻¹ (2I − J J⁻¹).
At p^N, each prime gives w and B_v = X_v w′ mod w, in the representation of splitting.c. They are recombined by CRT over the primes and reconstructed over Q with a running common denominator. The result is checked by evaluating the system at it modulo a prime that was not used. Until the check succeeds, more primes are lifted. Primes at which J is singular modulo (w, p) are dropped. After three such primes, msolve falls back to the multi-modular lifting.
Evaluation and solve.
- F(X) and J(X) are evaluated through a shared table of monomials, each obtained from a smaller one by one product modulo w. Products that no other monomial uses are summed without reduction and reduced once.
- Coefficients of F and J larger than p^N / 2 are reduced to the symmetric range modulo p^N before the evaluation, and their cost counts their limbs, so systems with large coefficients are neither slowed down nor mispriced.
- In the last steps J⁻¹ is not updated; instead, dX is refined by corrections at a decreasing precision, as in Dixon's lifting.
- The products by J and J⁻¹ use the transforms of fft_small, computed once per matrix (FLINT ≥ 3.2.0). Older FLINT releases keep fmpz_poly_mul.
Cost model and choice.
- The cost of a Newton step is counted in bit operations of the polynomial products, at the bit size of each step. The cost of a modular run counts the multiply-adds of F4 by the known pivots, an estimate of those among the new rows and the terms of its matrices, all recorded in the trace (Record the work of applying the trace when learning it #394), and the operations of FGLM or of the dense route. Both are independent of the machine and of the number of threads.
- The p-adic lifting is used when its predicted gain, at the most efficient precision up to p^1024, is at least 1.2. Otherwise, the multi-modular lifting is used as before. With -v 1 and above, msolve also prints the gain at the precision allowed by the memory limit at the start, which the run then exceeds as the limit grows.
- With the probabilistic linear algebra (-l 42, -l 44), no trace is learned and the multi-modular lifting is kept, unless --lifting 2 is given.
- Ideals that are not radical are not lifted (deg w < D).
Memory. The precision is bounded by a memory limit computed from the size of the modular computation and of the residues the multi-modular lifting would hold. The integers that FLINT caches per thread are freed after each Newton step, and the free memory is given back to the system (glibc). The reconstructed parametrization is stored as numerators over one common denominator per polynomial instead of one rational per coefficient, in the p-adic lifting and in the splitting along equal coordinates of Split the solution set along equal coordinates and find separating elements from the multiplication matrices #388. The rational values are freed as the parametrization is formed.
Overdetermined systems. Newton's iteration uses a square subsystem of the sparsest equations whose Jacobian is invertible modulo (w, p), if there is one. Otherwise it uses nv random combinations of all the equations. The check uses all the equations.
Inversion of J modulo (w, p). When w splits, J can be invertible while no entry of a column is invertible modulo w, each entry vanishing modulo a factor of w. The Gauss–Jordan elimination then forms its pivot from random combinations of rows. Without this, the lifting of jorge square1.sZ was abandoned and msolve fell back to the multi-modular lifting.
Threads. Each thread lifts its own primes. The first batch has one prime per thread.
Options.
- --lifting 0|1|2: cost model (default), multi-modular, or p-adic whenever applicable.
- --hensel-precision N, --hensel-cost-scale K and --hensel-memory MB override the parameters of the cost model. They are intended for experiments.
The options are described in the help and in a new subsection of the tutorial. With -v 2, the cost model and the Newton steps are printed as tables with a legend.
Tests. test/msolve/hensel_lifting.c covers the following (see Testing).
- The products by transforms, against fmpz_poly_mul.
- The inversion modulo (w, p), against the determinant, including split w.
- The lifting of systems with known solutions, with the variants listed under Testing. Each lifted parametrization is reconstructed over Q and compared with the parametrization modulo another prime.

Benchmarks

Settings:

Commit f94fd46 (the last one) for the systems whose p-adic peak memory was above the multi-modular one with 4674713, commit 4674713 for the others (marked *). The last three commits only change the memory of the reconstruction and speed it up; on the systems rerun, they reduced the peak memory by 12% to 40% and the time by 1% to 9%, with identical outputs.
Machine: Intel Core Ultra 9 285H (16 cores, one thread each, up to 5.4 GHz, 24 MB L3), 96 GB RAM, powersave governor.
Software: Ubuntu 26.04.1 (Linux 7.0), GCC 15.2 with -O2, FLINT 3.6.0, GMP 6.3.0, MPFR 4.2.2.
The machine was shared with other jobs, so timings vary by a few percent.
-P 2 -t 1 --random-seed 1, 3600 s limit (katsura 14 and eco 15 with the p-adic lifting: 2 hours).
"multi-modular" is --lifting 1 on the same build, which includes the faster F4 of Faster F4 and Gröbner bases over prime fields up to 52 bits #392. Times are wall-clock seconds, memory is the peak resident size.

All outputs that have a multi-modular reference are identical to it. Where both liftings finish, the p-adic lifting is 1.3 to 18 times faster than the multi-modular one, 3.9 times on average (geometric mean), and at least 3 times on 16 of the 22 systems; its peak memory is between 0.89 and 1.38 times the multi-modular one. On 16 other systems, the multi-modular lifting exceeds 3600 s while the p-adic one finishes in 5 to 44 minutes; the speed-ups there are at least 1.4 to 12.2, and about 33 and 41 by estimate on eco 15 and katsura 14.

Systems lifted, both liftings finished

system	p-adic	multi-modular	speed-up
KLYZ13	126.4 s, 419 MB *	2280.2 s, 441 MB	18.0
jorge rectangle1.sZ	240.4 s, 702 MB *	2865.2 s, 702 MB	11.9
jorge square2.sZ	153.3 s, 707 MB *	1618.9 s, 707 MB	10.6
jorge genericPlanar3_parallel.G1	134.6 s, 293 MB *	1148.1 s, 293 MB	8.5
coplanarity (visual servoing)	141.7 s, 706 MB *	958.7 s, 706 MB	6.8
jorge lozenge1.sZ	110.3 s, 601 MB *	637.0 s, 601 MB	5.8
katsura 12	45.7 s, 225 MB *	218.6 s, 232 MB	4.8
eco 13	39.6 s, 262 MB *	184.7 s, 293 MB	4.7
zimmermann	5.6 s, 49 MB *	24.9 s, 47 MB	4.5
jorge rectangle_parallel.G1	62.7 s, 353 MB *	276.4 s, 353 MB	4.4
crit_d_4_n7	846.3 s, 911 MB	2927.3 s, 939 MB	3.5
KLYZ22	6.9 s, 46 MB *	23.3 s, 45 MB	3.4
jorge square1.sZ	68.3 s, 691 MB *	223.2 s, 691 MB	3.3
random3_7_7	227.9 s, 390 MB	721.1 s, 327 MB	3.2
crit_deg_3_n_7_p_7	235.6 s, 415 MB	745.4 s, 327 MB	3.2
cp_d_3_n_7_p_7	239.8 s, 417 MB	719.5 s, 327 MB	3.0
random2_10_10	175.7 s, 375 MB	327.0 s, 334 MB	1.9
minrank_d_2_p_2_q_9_r_2	1527.3 s, 820 MB	2674.1 s, 596 MB	1.8
eco 12	6.5 s, 75 MB *	11.2 s, 82 MB	1.7
chandra 11	124.6 s, 226 MB	202.4 s, 249 MB	1.6
katsura 11	8.7 s, 63 MB	13.5 s, 71 MB	1.6
cp_d_3_n_7_p_2	391.1 s, 311 MB	492.1 s, 270 MB	1.3

The eco systems are run with their last two variables swapped.

The coplanarity and jorge systems come from the computation of the equilibria of image-based visual servoing controllers [7].

Systems lifted, multi-modular lifting over 3600 s

system	p-adic	speed-up
eco 14	296.0 s, 884 MB *	> 12.2
KLYZ12	299.0 s, 594 MB *	> 12.0
katsura 13	323.9 s, 663 MB	> 11.1
KLYZ23	398.6 s, 1103 MB *	> 9.0
jorge square_NOTparallel.G0	449.4 s, 344 MB	> 8.0
jorge genericPlanar3_NOTparallel.G1	546.1 s, 340 MB	> 6.6
henrion7	588.2 s, 805 MB	> 6.1
jorge rectangle_NOTparallel.G1	596.7 s, 371 MB	> 6.0
crit_deg_4_n_6_p_6	854.5 s, 885 MB	> 4.2
jorge square3.sZ	1078.9 s, 727 MB *	> 3.3
jorge rectangle2.sZ	1118.8 s, 743 MB *	> 3.2
jorge rectangle_NOTparallel.coplanarity	1129.7 s, 743 MB *	> 3.2
random2_11_11	1160.5 s, 1553 MB	> 3.1
gametwo7	1405.8 s, 1014 MB	> 2.6
eco 15	2413.7 s, 4769 MB *	> 1.5 (estimated ≈ 33)
katsura 14	2606.2 s, 2912 MB	> 1.4 (estimated ≈ 41)

For eco 15 and katsura 14, the estimates divide the time of the multi-modular lifting, extrapolated from the bits the p-adic run needed and the time of one modular run, by the p-adic time: about 660 primes of 121 s (22 hours) for eco 15 and 1440 primes (30 hours) for katsura 14.

Other systems

Multi-modular lifting kept. Seven systems of these benchmarks get a predicted gain below 1.2 and keep the multi-modular lifting: eco 11, vor1, crit_deg_3_n_6_p_6, random3_6_6, random2_9_9, crit_d_4_n6 and sot1. Their times are within 6% of --lifting 1. About fifty other systems of the benchmark set were triaged earlier and keep it too. The cost model only uses the first modular runs, which the multi-modular lifting needs anyway.
Both liftings over the limit. akt-2, steiner, minrank_d_2_p_3_q_5_r_2 and jorge generic1_NOTparallel.G0 exceed 3600 s with both liftings.

Testing

make check passes (76/76). hensel_lifting is a new test. It runs in about 2 s and covers:
- the products by transforms, against fmpz_poly_mul;
- the inversion modulo (w, p), against the determinant, including a split w (the case of square1.sZ) and singular matrices;
- the lifting of systems with known solutions: square, at two primes recombined by CRT, overdetermined with a square subsystem and with combinations, at a small prime with corrections, and in degree 70 with and without transforms (identical liftings).
Each lifted parametrization is reconstructed over Q and compared with the parametrization modulo another prime.
Injected bugs. The test fails when the combinations of rows in the inversion are disabled (84 checks fail), when the scaling of the transforms is dropped (38), or when the update of w in Newton's step is wrong (47). It also passes with HENSEL_NTT=0. Under valgrind, it shows no invalid access and no definite leak; the blocks reported as possibly lost are those of the integer cache of FLINT, which the test does not free at exit.
Options. -l 42 and -l 44 keep the multi-modular lifting by default and lift p-adically with --lifting 2. On cyclic 5, --split 0|2 with --lifting 1|2 give the same output.
Benchmarks. All outputs of the benchmarks above that have a reference are identical to it.

References

[1] M. Giusti, G. Lecerf, B. Salvy. A Gröbner free alternative for polynomial system solving. Journal of Complexity 17(1), pp. 154–211, 2001. https://doi.org/10.1006/jcom.2000.0571

[2] J. von zur Gathen, J. Gerhard. Modern Computer Algebra, 3rd edition. Cambridge University Press, 2013. Chapter 9 (Newton iteration) and Chapter 15 (Hensel lifting and factoring polynomials).

[3] J. D. Dixon. Exact solution of linear equations using p-adic expansions. Numerische Mathematik 40, pp. 137–141, 1982. https://doi.org/10.1007/BF01459082

[4] G. Lecerf. Kronecker, a Magma package for polynomial system solving. https://www.lix.polytechnique.fr/~lecerf/software/kronecker/

[5] X. Dahan, M. Moreno Maza, É. Schost, W. Wu, Y. Xie. Lifting techniques for triangular decompositions. ISSAC 2005, pp. 108–115.

[6] R. Lebreton. Relaxed Hensel lifting of triangular sets. Journal of Symbolic Computation, 2015. https://doi.org/10.1016/j.jsc.2014.09.012

[7] A. Colotti, J. García Fontán, A. Goldsztejn, S. Briot, F. Chaumette, O. Kermorgant, M. Safey El Din. Determination of All Stable and Unstable Equilibria for Image-Point-Based Visual Servoing. IEEE Transactions on Robotics 40, pp. 3406–3424, 2024. https://doi.org/10.1109/TRO.2024.3422050

[8] Z. Chen, A. Storjohann. A BLAS based C library for exact linear algebra on integer matrices. ISSAC 2005, pp. 92–99. https://cs.uwaterloo.ca/~astorjoh/iml.html

[9] J.-G. Dumas, T. Gautier, M. Giesbrecht, P. Giorgi, B. Hovinen, E. Kaltofen, B. D. Saunders, W. J. Turner, G. Villard. LinBox: A Generic Library for Exact Linear Algebra. ICMS 2002, pp. 40–50. https://linalg.org

[10] The FLINT team. FLINT: Fast Library for Number Theory. https://flintlib.org

Parts of this PR were written with the assistance of Claude Code (Claude Opus 5.5); all changes were reviewed and tested by the author.

🤖 Generated with Claude Code

mohabsafey · 2026-11-06T14:10:33Z

I swtiched it to draft because PRs on the F4 part impact the cost model. Let us wait for this one

When the parametrization over QQ fails because the default linear form does not separate the solutions, msolve used to try new linear forms one after the other, with a full F4 run for each attempt. It now builds the multiplication matrices modulo the prime of that run and searches for a separating element directly. The number of distinct solutions is the degree of the squarefree part of the minimal polynomial of a random linear form. The search first tries the single variables, then builds a sparse linear form with small coefficients greedily, adding one variable at a time until the minimal polynomial reaches the degree reached by a random form. The form found is added with a new variable as before. Its small coefficients lower the height of the parametrization, so fewer primes are needed. Setting MSOLVE_SEPELEM=0 restores the former behaviour. For a non-radical ideal, the eliminating polynomial of the form found has the degree of the squarefree part found by the search, and the parametrization is now accepted directly. This fixes nonradical-radicalshape, for which the former flow returned only one of the two solutions. The references for cyclic5 and issue-230 change because the search picks different linear forms. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

… slows F4 Adding the separating linear form with a new variable can make the Groebner basis much larger: on the coplanarity system (coming from image-based visual servoing systems), it grows from 358280 to 714997 terms, and the F4 learning time from 19.5 to 123 seconds. At the first prime, msolve now counts the nonzero terms of all the F4 matrices, with the form and for the original system, using a new counter in neogb. If there are more than 20% more terms with the form, the multi-modular steps run F4 on the original system and compute the parametrization of the linear form from its multiplication matrix: a Krylov sequence, Berlekamp-Massey for the eliminating polynomial, and the coordinates from the scalar sequences as in FGLM. These counts do not depend on the machine or on the number of threads, so neither does the choice. At the first prime, this parametrization is checked against the one from FGLM. It also handles non-radical ideals, returning the squarefree eliminating polynomial. MSOLVE_SEPDENSE=0 disables this path and MSOLVE_SEPDENSE=2 forces it. On coplanarity, the F4 matrices have 292 million terms with the form against 93 million without; the output is unchanged and the time goes from 9734 to 1668 seconds. On the fekete-valdes n6 systems (coming from the characterization of logarithmic Fekete critical configurations), the matrices have 27% more terms, and a prime costs 80 seconds instead of 128. On G1 and square1.pZ, the counts differ by less than 2% and the usual path is kept, which is also the faster one there. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

On some systems, the eliminating polynomial is a product of many factors of small height. On the fekete-valdes n6 system, it has 8924 bits, while its 392 factors over Q have at most 71 bits. At the first prime, msolve now splits the eliminating polynomial by taking gcds with the differences of the coordinates: two solutions are in the same component when the same pairs of their coordinates are equal. If the largest component has degree at most a quarter of the number of solutions, each component is lifted on its own, on either path (FGLM or the multiplication matrix). A component is identified across primes by its degree and by its pattern of equal coordinates, and it is done when its rational reconstruction is confirmed by one more prime. The parametrization of the whole solution set is then recombined over Q by CRT in Q[t]. If the components differ at three primes, msolve falls back to the usual lifting. MSOLVE_SPLIT=0 disables the splitting. On fekete-valdes n6, the output is unchanged: 8 primes are needed instead of 309, and the time goes from 29094 to 4093 seconds. The fekete-valdes n5 system is added as a test: its 32 components are lifted with 2 primes instead of 9. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

When the solution set splits along equalities of coordinates, each component now keeps its own parametrization over Q, with an integer eliminating polynomial of small degree and height. The real roots are isolated component by component, and the coordinates of each root are refined from the parametrization of its component instead of the one of the whole solution set. The points are then gathered and sorted by the value of the separating element, as before. The parametrization of the whole solution set is still recombined by CRT in Q[t] when it is printed (-P 1 or -P 2) and for the Julia interface, which always returns it; with -P 0, it is no longer built. On fekete-valdes n5, the 38 points are the same as with the global isolation, in the same order. Many coordinates are now exact: a linear component gives its rational root exactly, so the output prints [1 / 2, 1 / 2] instead of an interval around 1/2. The references of the test are updated accordingly. On fekete-valdes n6p0idealVarsReordered (392 components, 398 real roots), with -P 1 -t 1, the extraction of the real roots goes from 10211 s (5.8 s of isolation, the rest refining the coordinates) down to 0.01 s, and msolve from 14502 s down to 3996 s. On n6p0idealMaxRank3 (261 components, 267 real roots), it goes from 113 s down to 0.01 s, and msolve from 2706 s down to 2591 s. In both cases, the points are the same and come in the same order. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

When a linear form is added with a new variable, the last coordinate of each real point is hidden by decrementing its number of coordinates, so real_point_clear never freed it. AddressSanitizer reported these leaks with -P 0 and -P 1, for instance 664 bytes on cyclic5 and 1120 bytes on fekete-valdes n5. They are now cleared when the coordinate is hidden. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

The test added with the splitting of the solution set was listed in Makefile.am but not in AC_CONFIG_LINKS, so its script, input and reference outputs were missing from the tarball and make distcheck failed. They are now listed with the other tests. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

…jections The search for a separating element computed the degree of the minimal polynomial of each candidate form from its dense multiplication matrix, with a Wiedemann sequence: about 2 D^3 operations per form, D being the number of solutions. On noon6 (D = 717), the 32 forms tested took 9.7 of the 10.4 seconds of msolve, and on noon7 (D = 2187), the 61 forms took 985 seconds, while msolve took 383 seconds before the search was introduced. For a radical ideal, the random form computed first separates the solutions. Its parametrization x_v = X_v(t) mod Q is now computed once, and the degree for a form sum c_v x_v is the one of the minimal polynomial of sum c_v X_v in F_p[t]/(Q), from the sequence of the projections of its powers: 2 D products modulo Q. The degrees are the same up to an unlucky choice of the projections, so the same forms are found. Non-radical ideals keep the dense matrices. The search now takes 2.85 seconds on noon6 and 73 seconds on noon7, with the same forms and the same outputs. msolve takes 3.6 seconds on noon6 (5.6 before the search was introduced) and 128 seconds on noon7 (383). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

…he eliminating polynomial split_set_mpz_upoly copied each coordinate with its own length, so a coordinate whose leading coefficient vanishes lost it, while the usual lifting always writes deg(elim) coefficients. The parametrization was the same, but the output differed from the one of the usual lifting. This showed up on the coplanarity system of the visual servoing benchmarks, where four coordinates have degree 778 for an eliminating polynomial of degree 780. The coordinates are now padded with zeros to deg(elim) coefficients. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

…plitting by options The search for a separating element, the route of the multi-modular steps with a separating linear form and the splitting of the solutions along equalities of coordinates were controlled by the environment variables MSOLVE_SEPELEM, MSOLVE_SEPDENSE and MSOLVE_SPLIT. They are replaced by the options --separating-element (1: search, the default, 0: no search), --dense-route (0: chosen from the sizes of the F4 matrices, the default, 1: never through the multiplication matrix of the original system, 2: always when it applies) and --split (0: when the largest component has at most a quarter of the solutions, the default, 1: never, 2: whenever there are several components). Their values are kept in data_gens_ff_t, and they are described among the advanced options of the help, whose column of option names is widened for them. msolve no longer reads environment variables. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

Add a subsection on the options --separating-element, --dense-route and --split to the section on flags and options of the tutorial. It explains why the last variable must be separating, what each option controls, its default value, and whether the result depends on it. It also notes that the search for a separating element is not used when the staircase is rejected, and that forcing the splitting may slow down the computation, since the lifting of a large component costs about as much as that of all the solutions. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

The new test runs cyclic5-qq with the options --separating-element, --dense-route and --split, with one and two threads. On this system, the last variable is not separating, the search finds the linear form 0 4 3 2 1, and the multi-modular steps use the original system. The solutions split into 6 components, the largest of degree 20 out of 70, so the splitting is not used by default and --split 2 forces it. Without the search (--separating-element 0), the changes of variable order fail and the linear form 1 2 3 4 5 is added: the output is compared with a new reference. With the other options, it is compared with the reference of the existing test of cyclic5-qq. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

The F4 operations of the runs applying the trace are no longer counted: the linear algebra kernels do not increment the application counters, and the counts of a run were never added to the meta data of the caller, so the verbose output reported 0 Gops/sec for the learning and the application phases. Counting in the kernels would add work to their inner loops and needs care with several threads. The work of an application of the trace is the same for all primes, so it is now recorded in the trace when it is learned, outside the kernels. For each round, construct_trace adds to nmadd the lengths of the known pivots used by each row kept in the trace, that is the multiply-adds of the reductions of the rows by the known pivots, and to nterms the terms of the rows to be reduced and of the reducers used. add_lms_to_trace adds to nmadd_new the number of new elements of the round times their total length, an estimate of the multiply-adds of the reductions of the new rows by each other, which nmadd does not include. These quantities are meant for cost models choosing between strategies, such as the p-adic lifting of the parametrizations. In the verbose output of msolve_trace_qq and of the lifting of Gröbner bases, the rate of the application phase is computed from nmadd. The rate of the learning phase and the numbers of additions, multiplications and reductions printed at info level 3, which were always 0, are no longer printed. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

With probabilistic linear algebra (-l 42, 43 or 44), F4 learns no trace and the meta data has no trace to read the multiply-adds from, so msolve crashed when printing the rate of the application phase at info level 2. The rate is now printed only when a trace exists. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

For systems with at least as many equations as variables and simple solutions, the parametrization computed modulo a prime can be lifted p-adically by Newton's iteration instead of being recomputed at many primes. The approach is hybrid between a pure p-adic lifting and a multi-modular one: the parametrization is lifted modulo p^N at a few primes p, and the lifted images are combined by CRT. A single prime lifted to a high precision and many primes not lifted at all are its two extreme cases. The new file hensel.c lifts the eliminating polynomial and the coordinates together with the inverse of the Jacobian matrix, attempts rational reconstruction as the precision grows, using vector witnesses for the common denominator, and checks the result by evaluating the system modulo a prime that was not used. A cost model chooses the precision per prime, hence the balance between the number of primes and the lifting at each of them, and decides whether this lifting pays off compared with the multi-modular one. It is calibrated on measured costs: the operations and matrix terms of the F4 runs, which are now counted in md_t (max_mat_nterms, tracer_nr_*), and the matrix-vector products of FGLM. The precision is also bounded by a memory limit, computed from the sizes of the computation rather than from the memory of the machine, so that the computation does not depend on the machine it runs on. The limit is an estimate of the memory of the multi-modular computation: 8 bytes per term of the largest F4 matrix, plus the 28 D^2 bytes that the multi-modular lifting preallocates for D solutions, plus twice the size of the residues of the parametrization at the estimated number of bits. That estimate is bounded by twice the bits already obtained, so the limit grows as the lifting proceeds. The precision of a prime is then the largest one at which the data of Newton's iteration (values of the monomials of the system, the coordinates, the Jacobian matrix and its inverse) fit within the limit. The preallocations of the multi-modular lifting are skipped while this lifting runs and restored when it falls back. The new option --lifting selects the strategy: 0 (default) lets the cost model choose, 1 enforces the multi-modular lifting and 2 the p-adic lifting when it applies. The options --hensel-precision, --hensel-cost-scale and --hensel-memory fix the precision, scale the cost of the modular runs in the model and fix the memory limit, for experiments. They are described in a new section of the help. On the benchmarks, outputs are identical to the multi-modular lifting, with speedups from 1.3x to 18x where the p-adic lifting is chosen (katsura 13: 6853 s to 680 s, eco 13: 555 s to 82 s, coplanarity: 1561 s to 243 s), and peak memory within 1.65x of the multi-modular one. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

Add a subsection on the lifting of the parametrization over the rationals to the section on flags and options of the tutorial. It describes the multi-modular and p-adic liftings, and when the p-adic lifting applies, and explains that the latter is hybrid between a pure p-adic lifting and a multi-modular one. It details how its parameters are chosen: - the cost of a modular computation, from the operations of the F4 tracer and of FGLM; - the cost model of Newton's iteration; - the choice of the precision per prime and the predicted gain, with its threshold of 1.2; - the memory limit; - how the precision follows the estimate of the number of bits during the lifting. It also states that these choices do not depend on the machine and that both methods return the same parametrization. The flags --lifting, --hensel-precision, --hensel-cost-scale and --hensel-memory are documented with their default values, the parameters being set automatically by default. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

With several threads, the p-adic lifting runs in batches of primes, one prime per thread, each lifted by Newton's iteration. The precision reached by the primes of a batch is bounded by a memory limit. That limit was the estimated memory of the multi-modular computation with one thread, shared among the primes of the batch, so that with four threads each prime got about a quarter of the precision it would get alone, and two to four times more primes were needed. The F4 part of the memory reference is now counted once per thread for the limit of a batch, since the multi-modular computation also holds the F4 data of one prime per thread. The limit used to choose the lifting still counts it once, so that the choice does not depend on the number of threads. The first batch is also made of one prime per thread, the first prime and the images of the first multi-modular step, instead of the first prime alone. With one thread, the computation is unchanged. With four threads, eco 13 takes 32 s instead of 42 s (multi-modular with four threads: 167 s), katsura 12 32 s instead of 41 s (multi-modular with four threads: 112 s) and coplanarity 103 s (multi-modular with four threads: 509 s), with a peak memory within 1.4 times that of the multi-modular computation with four threads. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

In each step of Newton's iteration, J(X) dX = F(X) is solved with the inverse of the Jacobian matrix known only modulo p^jprec, by successive corrections. Each correction formed J^-1 s and J y with s and y at the full remaining precision, although only their jprec lowest digits are correct. The corrections now follow the p-adic solving of linear systems of Dixon (J. D. Dixon, Exact solution of linear equations using p-adic expansions, Numer. Math. 40 (1982), 137-141), in which each step multiplies the inverse modulo p by the residue modulo p only, and then updates the residue with the matrix. The same scheme is used, for instance, by IML, LinBox (DixonSolver) and FLINT (fmpz_mat_solve_dixon) for integer linear systems. Here, with the inverse modulo p^jprec, s and y are truncated to jprec digits: J^-1 s is formed modulo p^jprec, and only J y is formed at the remaining precision. This halves the time of the solve, the largest part of Newton's iteration (katsura 12: 34 s instead of 64 s for the whole iteration, 18 s instead of 34 s for the solve). The cost model counts the corrections in the same way, which raises the predicted gains by 5 to 20% and brings them closer to the measured ones. On the benchmarks where the p-adic lifting is chosen, outputs are identical and the computation time is reduced by up to 25% (katsura 13: 594 s instead of 680 s, katsura 12: 76 s instead of 102 s, eco 13: 70 s instead of 82 s), with the same peak memory; the choice of the lifting is unchanged on the other benchmarks. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

With the p-adic lifting, the time spent in the CRT and the rational reconstruction was not added to the reported total, which stayed at zero, and the timings of one modular run with the header of the steps were printed again at each batch of primes. The lifting now measures separately the elapsed times of Newton's iteration, of the CRT and rational reconstruction, and of the checks modulo an unused prime; each attempt of reconstruction reports its own times instead of the time elapsed since the start of the batch. The timings of one modular run and the header of the steps are printed once, the header naming the p-adic or the multi-modular steps. At the end, both liftings report the times of the learning run of F4 and of the rest of the first prime, the total times of F4 and FGLM in the modular runs, of Newton's iteration when it is used, of the CRT and rational reconstruction, of the checks and of the remaining work, which add up to the total time. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

When the solution set is split into components, or after the p-adic lifting when the parametrization is not requested (it then forms a single component), the real roots were isolated without any of the information printed at the verbosity levels, only their total number being reported. The information is now printed for each component, preceded by its number and degree when there are several, and the total is only reported in that case. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

Each attempt of reconstruction printed a long sentence, and the memory limit and the vector witnesses were reported on lines of their own. A legend is now printed once before the steps, and each attempt prints one line with the precision, the bits of the modulus and their estimate, the reconstructed coefficients, the bits of the denominator from the vector witnesses, the precision and memory allowed by the limit, and the times of the CRT and rational reconstruction and of the check. A line marks the start of each batch with the indices of its primes, and the timings of Newton's iteration at each prime are abbreviated as in the legend. In the final timings, those of F4 and FGLM in the modular runs are given as the number of primes times the elapsed time of one run, divided by the number of threads running them at once. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

The first batch of primes of the p-adic lifting is made of the first prime and of the images of the first modular run, which computes one image per thread: one of them was left out and discarded. It is now kept, and with one thread it makes the next batch without a modular run, which saves the time of one modular run (KLYZ23: 636 s instead of 685 s, the modular run taking 52 s). With several threads the image left out is computed alongside the others, and is still discarded. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

Newton's iteration needs the Jacobian matrix to be invertible at the solutions, which fails at the multiple ones: on an ideal which is not radical, the lifting was started, and given up after three primes at which the Jacobian matrix is singular. The parametrization is that of the radical, whose degree, the number of solutions, is smaller than the dimension of the quotient exactly when the ideal is not radical; the p-adic lifting is now skipped in that case before the first prime, with a message at the verbosity levels. On the benchmarks, this concerns bruno and phuoc1, where the usual lifting was already chosen by the cost model. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

With more equations than variables, Newton's iteration used nv combinations of all of them, each holding all the extra equations: on eco 12 times the ideal of a point, given by the 144 products of the equations and of the variables, they have 5600 terms instead of 89 for a square subsystem. The equations are now chosen from the one with the fewest terms at the first prime: an equation is kept if its row of the Jacobian matrix along the parametrization, reduced by those of the equations kept, has an entry invertible modulo the eliminating polynomial in a column without pivot. The solutions of the nv equations kept contain those of the system, at which their Jacobian matrix is then invertible, so that they are lifted as the combinations were; the combinations are used when the choice fails, as on coplanarity. On eco 12 times the ideal of a point, the p-adic lifting is now chosen, in 11.5 sec instead of 28.6 sec for the multi-modular one; on the cp_d benchmarks, whose equations have the same monomials, the cost is unchanged. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

The cost model and the system of Newton's iteration were printed on two long lines, holding about fifteen numbers. They are now printed as one block, with one quantity per line: the equations, the system used by Newton's iteration when there are more equations than variables, the terms and monomials of this system and of its Jacobian matrix, the operations of F4 and of FGLM in a modular run, the costs of a modular run and of the lifting of one prime, the gain at the most efficient precision and the one within the memory limit at the start, with the threshold above which the p-adic lifting is chosen. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

The values of the equations and of their Jacobian at the parametrization are computed from those of their monomials, each obtained from a predecessor by one product modulo w. The products without successor are only added to the equations, so that their reduction modulo w, about two thirds of their cost, can be avoided: they are kept unreduced and each polynomial they are added to is reduced once at the end. This is done for F or for J when they have more such products than polynomials to reduce (eco 12: 56 products for 10 polynomials in F, but 66 for 95 in J). The evaluation of F is 1.2 to 2.8 times faster (katsura 12: 18.7 s to 6.6 s, eco 13: 8.9 s to 3.8 s), and the whole computation up to 1.6 times faster with identical outputs. The cost model and the memory estimate account for it, so that the p-adic lifting is now chosen for chandra 11 (242 s to 148 s) and random3_6_6 (19.3 s to 12.9 s). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

In the last steps of Newton's iteration, where J^-1 is not updated, each step corrects its solution in rounds of products of J^-1 and of J by vectors, nv^2 products of polynomials each, whose operands were transformed again at every product although J^-1 stays the same over these steps and J over the rounds of a step. Their transforms modulo the 50-bit primes of fft_small are now computed once, so that a product by a vector costs nv transforms, nv^2 products of transforms and nv inverse transforms; the products of the updates of J^-1 are formed in the same way. The results are the same integers, reduced as before. The transforms take memory, which may lower the precision allowed by the memory limit (eco 12: p^7 instead of p^12), so that the cost model and the memory estimate account for them and the transforms are used only when they give the better gain within the memory limit. Since this limit grows with the bits of the residues, the choice is made again at each update of the limit, at the best precision up to it, and changed only for a gain of 5% (random2_10_10: products of polynomials at the start, by transforms from 54 MB on). With identical outputs, katsura 12 takes 54 s instead of 65 s, KLYZ13 161 s instead of 183 s, cp_d_3_n_7_p_7 240 s instead of 274 s, random2_10_10 198 s instead of 235 s, chandra 11 119 s instead of 148 s, and random2_9_9 is now lifted (33 s instead of 50 s). This needs the interface of fft_small of FLINT 3.2.0 or later; with older versions, or with HENSEL_NTT defined as 0, the products are those of fmpz_poly_mul. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

When w splits, every entry of a column of J can vanish modulo some factor of w while J is invertible, so the search for an invertible pivot failed and the p-adic lifting gave up at every prime (jorge square1.sZ: lifting abandoned, 342 s with the multi-modular lifting). When no entry of the column is invertible, random multiples of the rows below are added to the pivot row, which gives an invertible entry with high probability when J is invertible. square1.sZ is now lifted in 132 s with the same output. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

The test checks the main functions of hensel.c against independent computations: the products of polynomial matrices by transforms against fmpz_poly_mul, and the inversion of J modulo (w, p) against the determinant, including a split w where no entry of a column is invertible. It also lifts the parametrizations of systems with known solutions: square, at two primes, overdetermined with a square subsystem or with combinations, with the corrections of the solve, and in degree 70 with and without transforms. Each result is reconstructed over Q and compared with the parametrization modulo another prime. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

The buffers of the products by transforms of the p-adic lifting were allocated with aligned_alloc, which MinGW does not provide, so that msolve did not build on Windows. They are now allocated with flint_aligned_alloc and freed with flint_aligned_free, which FLINT implements on all platforms. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

…fting The cost model of hensel_choose weighed the F4 operations counted while applying the trace, which the linear algebra no longer counts, and the terms of the F4 matrices of the last run, which are no longer counted when the trace cache provides the matrices. It now uses the work recorded in the trace when it is learned (see construct_trace): the multiply-adds by the known pivots, the estimate of the multiply-adds among the new rows and the terms of the matrices, which do not depend on the prime nor on the number of threads. The counters of the runs applying the trace and their accumulation in finalize_f4 are removed, and the largest number of terms of an F4 matrix is updated only outside these runs, which may be made for several primes in parallel. On the dense route of the separating element, the parametrization is computed from the dense matrix of the multiplication by the separating form, and the model now counts its 2 D^3 multiply-adds, where it counted nothing. The constants of the model are fitted again on the modular runs of 48 systems, the F4 of master being faster: a multiply-add by a known pivot weighs 3.0 units of the model, one of the estimate among the new rows 5.2, a term of the matrices 44, a multiply-add of FGLM 1.5 and one of the dense route 7.6. With the probabilistic linear algebra (-l 42 or 44), F4 learns no trace, so the model has no measure of the work of F4: the multi-modular lifting is then used unless the p-adic one is forced with --lifting 2. The help of --lifting and the tutorial say so. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

…on's iteration The evaluation of F and of J multiplied the values of the monomials by the integer coefficients of the equations, reducing the sums modulo the current precision q only at the end. This is cheap for small coefficients, but for large ones each term costs a product by the whole coefficient. Coefficients larger than q / 2 in absolute value are now replaced by their symmetric residues modulo q, so that a term costs at most a product of coefficients of the size of q; the other ones, usually all of them, are kept, and the computation is unchanged for them. The cost model of the terms multiplies the part per bit by the number of limbs of each coefficient, at most that of q, which gives the former cost for coefficients of one limb. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

The gain deciding whether the p-adic lifting is used was the best precision per cost within the memory limit at the start of the lifting. This limit grows with the estimated number of bits, and the runs reach much higher precisions, so this gain underestimated the actual one: on chandra 11 it was 1.11, below the threshold 1.2, for a precision of p^6 at the start, while the run reaches about p^34 per prime and the p-adic lifting is 1.71 times faster than the multi-modular one. The decision now uses the gain without this limit (1.49 on chandra 11). The variant of the products, with or without the transforms, is still the one with the best gain within the memory limit, and both gains are printed. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

FLINT keeps the integers freed by Newton's iteration in a cache of each thread. The lifting emptied it only at the end of each batch of primes, after the CRT and the reconstruction attempts, which therefore ran with it. The cache is now emptied after each step of Newton's iteration, and on glibc the free pages are given back to the system with malloc_trim. On random2_10_10 the memory held during the lifting drops from 337 MB to 199 MB, about what the multi-modular lifting holds at that point, with the same run time. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

Each coefficient reconstructed by the p-adic lifting is a rational number holding its own denominator, which the polynomials of the parametrization share. Forming the parametrization kept all of them alongside their copies in the polynomials. The function that reads them last now frees them as it goes. The multi-modular lifting is unchanged. With both commits, the peak memory of the p-adic lifting drops by about 20% with identical outputs: random2_10_10 581 to 473 MB, chandra-11 372 to 299 MB, crit_d_4_n7 1515 to 1192 MB, minrank_d_2_p_2_q_9_r_2 1008 to 820 MB, against 333, 249, 939 and 596 MB for the multi-modular lifting. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

The reconstructed values of the parametrization were rational numbers, one per coefficient, each holding its own copy of the denominator that the eliminating polynomial and each B_v share. They are now numerators over one denominator per polynomial, in the p-adic lifting and in the splitting along equalities of coordinates. The rational reconstruction of the p-adic lifting writes the numerators over its running denominator, and the free pages left by its temporary integers are given back to the system between the polynomials. The outputs are identical; the CRT and rational reconstruction phase is 11% to 38% faster, and the peak memory of the p-adic lifting drops further: random2_10_10 473 to 375 MB, chandra-11 299 to 226 MB, crit_d_4_n7 1192 to 911 MB, against 333, 249 and 939 MB for the multi-modular lifting. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

mohabsafey marked this pull request as draft November 6, 2026 14:04

mohabsafey and others added 13 commits November 7, 2026 11:31

Adds missing co-author to separating.c

cde7e98

Adds missing co-author to splitting.c

197ab3f

mohabsafey mentioned this pull request Nov 7, 2026

Record the work of applying the trace when learning it #394

Merged

mohabsafey and others added 14 commits November 7, 2026 20:02

mohabsafey and others added 10 commits November 7, 2026 20:02

mohabsafey force-pushed the hensel-lifting branch from 8841c22 to f94fd46 Compare November 9, 2026 06:11

mohabsafey marked this pull request as ready for review November 9, 2026 06:29

Navigation Menu

Repository navigation

Lift the parametrization over QQ p-adically by Newton's iteration - #391

Lift the parametrization over QQ p-adically by Newton's iteration#391
mohabsafey wants to merge 37 commits into
algebraic-solving:masterfrom
mohabsafey:hensel-lifting

mohabsafey commented Nov 6, 2026 •

edited

Loading

Uh oh!

mohabsafey commented Nov 6, 2026

Uh oh!

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

1 participant

Conversation

mohabsafey commented Nov 6, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Related software

What is new

Changes

Benchmarks

Systems lifted, both liftings finished

Systems lifted, multi-modular lifting over 3600 s

Other systems

Testing

References

Uh oh!

mohabsafey commented Nov 6, 2026

Uh oh!

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

1 participant

mohabsafey commented Nov 6, 2026 •

edited

Loading