Repository navigation
Split the solution set along equal coordinates and find separating elements from the multiplication matrices - #388
Open
mohabsafey wants to merge 17 commits into
Open
mohabsafey wants to merge 17 commits into
mohabsafey wants to merge 17 commits into
Conversation
Rational reconstruction in the lifting of rational parametrizations over QQ (linear forms, multiplication matrix witnesses, trace and determinant, eliminating polynomial and coordinates) now goes through FLINT's half-gcd based fmpq_reconstruct_fmpz_2 instead of the quadratic Euclidean algorithm, falling back to ratrecon/ratreconwden when the bounds do not satisfy 2ND < m. The wrappers introduced for the lifting of Groebner bases are moved from lifting-gb.c to lifting.c so that both liftings share them. Outputs are unchanged; the time spent in CRT and rational reconstruction decreases by 25 to 40% (kat11: 1.22s to 0.70s). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
When lifting rational parametrizations over QQ, rational reconstruction was called every nbdoit iterations, nbdoit being doubled whenever the reconstruction took more than 20% of the modular computations; the number of primes used thus depended on timings and on the machine load. As for the lifting of Groebner bases, the period is now rr primes, rr being the largest power of 2 <= nprimes/128, i.e. max(1, rr/nthrds) iterations, and the reconstruction is attempted once per iteration, after the last good prime, instead of after each prime of the iteration. Outputs are unchanged and the number of primes is now reproducible. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
When lifting rational parametrizations over QQ, every coefficient of the eliminating polynomial and of the coordinates was lifted by incremental CRT at each prime, which is quadratic in the number of primes, although these CRT images are only used once the trace and determinant witnesses are reconstructed and verified. The modular images of the parametrization are now stored for each prime (in chunks of 64 primes which are never moved), and FLINT's multi-CRT is applied to them only when a rational reconstruction is attempted, one coefficient at a time, right before this coefficient is reconstructed; no CRT is thus applied beyond the first failed reconstruction, and no array of CRT images is kept. The coordinates are lifted one at a time, trying in turn the three bounds which were tried before on all of them. The witnesses are still lifted by incremental CRT. Outputs and numbers of primes are unchanged on 38 benchmark systems (eco, kat, noon, cyclic, henrion, chandra, random, minrank, cp_d, crit and others). The time spent in CRT and rational reconstruction decreases by up to 4.2x (random2_9_9: 20.4s to 4.8s; chandra-11: 50.7s to 14.5s), the total time by up to 17% with one thread and by up to 32% with four, and peak memory is unchanged or slightly lower. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
As for the lifting of Groebner bases, the rational reconstruction bounds used when lifting parametrizations over QQ (witnesses, eliminating polynomial and coordinates) are now computed from modulus / 2^RR_MARGIN instead of the modulus, so that 2ND <= modulus / 2^16. This makes spurious reconstructions unlikely: on cyclic7, two spurious reconstructions of the eliminating polynomial, each rejected by the next prime, are avoided. Outputs are unchanged on 35 benchmark systems; the number of primes increases by 0 to 4 on all of them but random2_9_9 (2627 to 2643, the reconstruction period being 16 primes at that point) and cyclic7 (194 to 193). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
When lifting parametrizations over QQ, each coefficient of the eliminating polynomial and of the coordinates is now reconstructed by trying in turn the fixed bounds used so far (D = M^(1/16) for the eliminating polynomial; D = M^(1/32), balanced bounds and N = M^(1/16) for the coordinates, with M = modulus / 2^RR_MARGIN), then the bounds set from the denominator accumulated so far, g: D = sqrt(M / g) and N = M / D, as for the lifting of Groebner bases. CRT is applied only once per coefficient. Previously, each fixed bound had to succeed on a whole coordinate. Using only the bounds set from the guessed denominator increases the number of primes by up to 1.9x on several families (eco, noon, henrion, cyclic), whose numerators are much larger than their denominators; the fixed bounds are thus tried first. Outputs are unchanged on 35 benchmark systems and the number of primes never increases. It decreases on KLYZ22 (279 to 225, 18% faster), whose denominators are large, and on the smallest Katsura systems, where some coordinates were lifted later than the eliminating polynomial because no single fixed bound suited all their coefficients (kat8: 22 to 16, kat9: 40 to 31, about 20% faster). From kat10 to kat12 all coordinates are lifted together with the eliminating polynomial, which is the bottleneck, so the number of primes and timings are unchanged (kat12: 323 primes, 325s vs 321s); the same holds for eco8 to eco12, where the fixed bounds always succeed and the fallback is never used. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
mulmat_reconstruct allocated ncols integers for the row denominators but cleared only nrows of them, leaking ncols - nrows integers at each call when multiplication matrices are lifted (-L 1). Also remove twelve static functions of lifting.c which are no longer called: the former lifting of multiplication matrices (crt_mpz_matfglm_*, mpz_matfglm_*, mpq_matfglm_*, crt_lift_mat, rat_recon_matfglm, check_matrix_and_linear_forms), replaced by mulmat_reconstruct, as well as old_rat_recon_trace_det and check_linear_forms. The corresponding types are kept in msolve-data.h, which is installed. Outputs are unchanged, with and without -L 1. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…nator The coefficients of the monic eliminating polynomial share a common denominator (the leading coefficient of the primitive integer polynomial). On dense random, MinRank and critical point systems, this denominator is about as large as the numerators, so reconstructing a single witness coefficient (trace or determinant) needs about twice the number of bits of these numerators, and the eliminating polynomial cannot be attempted before. Four witness coefficients (the trace and determinant ones, and the first nonzero ones from degrees nsols/2 and nsols/4) are now also reconstructed simultaneously: a short vector (d, n_1, .., n_4) with n_j = d u_j mod m is computed by LLL (FLINT's fmpz_lll), which needs about bits(n) + bits(d) / 4 bits of modulus. The candidate is accepted when its entries are below 2^(4 bits(m) / 5 - RR_MARGIN) and when it matches the images of the next prime; d is then used as a guess for the eliminating polynomial, which can then be reconstructed with its usual bounds. As LLL is more expensive than rational reconstruction, it is tried only while no single witness is known, at 16 primes and then each time the number of primes grows by 5%. Outputs are unchanged on 35 benchmark systems. On the 18 dense random, MinRank, critical point (cp_d, crit) systems, the number of primes decreases by 32 to 38%, and the total time by 20 to 39% on those taking more than 0.3s (geometric mean speedup 1.44; random2_9_9: 2643 to 1651 primes, 81.7s to 56.1s; random3_6_6: 763 to 495 primes, 33.5s to 21.7s). On the other families the number of primes is unchanged and timings are within noise (+-3%), the LLL attempts made before a single witness is known costing a few hundredths of a second. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The eliminating polynomial is reconstructed from its highest degree down, carrying the lcm of the denominators already found, so most coefficients only bring a few new denominator bits and are bound by the size of their numerators. The bound D = M^(1/64) is now tried before D = M^(1/16), which lowers the number of primes needed by 2 to 5% on katsura and eco systems (kat12: 323 -> 313, eco13: 238 -> 231) with identical outputs. Random and minrank systems are unaffected. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
AddressSanitizer reports a few constant-size leaks in the multi-modular setup, independent of the number of primes. Free invalid_gens once the input is imported, free the per-thread copies of the linear-form data and of the current leading monomials in free_msolve_trace_qq_initial_data (zero-initializing the pointer arrays), free the monomial basis when initial_modular_step fails, and free the FGLM data of the first call before the second call in the non-shape-position case. Also free the parametrization on the failure paths of nmod_fglm_compute_trace_data, drop a name buffer that was immediately overwritten when adding a linear form, and count the terms of the added linear form in gens->nterms so their coefficients are cleared. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
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>
mohabsafey
marked this pull request as draft
October 31, 2026 21:03
…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>
mohabsafey
marked this pull request as ready for review
October 31, 2026 21:15
…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>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
This PR builds on #384 (lifting of parametrizations); until #384 is merged, its commits also appear here. The last eight commits belong to this PR.
It targets systems where the parametrization over QQ is much larger than the solution set requires. Either the chosen linear form makes F4 much harder, or the eliminating polynomial is a product of many factors of small height. Outputs are unchanged, except that some coordinates are now printed exactly.
Changes
nonradical-radicalshape, for which only one of its two solutions was returned.MSOLVE_SEPELEM=0restores the former behaviour.MSOLVE_SEPDENSE=0disables this path,MSOLVE_SEPDENSE=2forces it.-P 1,-P 2) or returned to Julia.MSOLVE_SPLIT=0disables the splitting.[1 / 2, 1 / 2]instead of an interval).-P 0,-P 1). This leak predates this PR.Benchmarks
The runs use
-t 1. Timings are in seconds, before → after the corresponding commit.The coplanarity system comes from the computation of the equilibria of image-based visual servoing controllers [1]. The fekete-valdes systems come from the characterization of logarithmic Fekete critical configurations [2].
-P 1-P 1-P 1-P 2-P 2-P 2On coplanarity, the former flow tried 15 linear forms and settled on a dense one,
(1, 4, 9, …, 169); its time is the CPU time, since that run shared the machine with another one. The sparse form found by the search halves the number of primes, but it makes F4 three times slower per prime: the F4 matrices have 292 million terms with the form against 93 million without. The multiplication-matrix path removes this cost. On n6p0idealVarsReordered, the extraction of the real roots goes from 10211 s to 0.01 s. On fekete-valdes n5, the former flow tried three linear forms, each with a full F4 run, before finding one that separates the solutions; the search finds a sparse form in 0.01 s. On G1 and square1.pZ, the term counts differ by less than 2% and the usual path is kept. On noon7, the search tests 61 forms in 73 s; with the dense test of the first commit, it took 985 s, which the last commit fixes. On eco, katsura and the random systems, no splitting occurs and nothing changes.[1] 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
[2] D. Armentano, L. Bentancur, F. Carrasco, M. Fiori, M. Valdés, M. Velasco. Characterization of Logarithmic Fekete Critical Configurations of at Most Six Points in All Dimensions. ISSAC 2025, pp. 309–318. https://doi.org/10.1145/3747199.3747575
Testing
make checkpasses (72/72).fekete-valdes-n5-qqis added as a test. The references for cyclic5, issue-230 and nonradical-radicalshape change because the search picks other linear forms, and nonradical-radicalshape now returns both solutions.-P 0,-P 1and-P 2, 1 and 2 threads, and with the multiplication-matrix path forced.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