Repository navigation
Faster F4 and Gröbner bases over prime fields up to 52 bits - #392
Merged
Merged
Conversation
With several threads, the serial parts of F4 dominated for systems like bayes148: sorting the columns of the matrix, inserting the new monomials of the reducer rows during symbolic preprocessing and moving the new columns into the basis hash table. The columns are now sorted blockwise in parallel and merged afterwards. For the insertions, the hash table is first searched in parallel; the missing monomials then get consecutive positions in the order in which a serial insertion would have added them and are written into the hash table in parallel, claiming their slots via compare-and-swap. Monomials missing in several rows of symbolic preprocessing are mapped to their first occurrence via a temporary concurrent hash table. Thus all hash tables have the same entries as before, which keeps the computation deterministic and the tracer usable across primes. For bayes148 over a 31 bit prime field this reduces the running time with 8 threads from 4.7 to 3.3 seconds, the other benchmarks are not affected. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Symbolic preprocessing searched the reducer of each column by checking the short divisor masks of all lead monomials of the basis. For bayes148 this means 11 billion mask checks of which only 6 million pass, so the search dominated symbolic preprocessing. At the start of symbolic preprocessing the lead monomials are now sorted by their masks and a binary trie on the sorted masks is built. Each node stores the bits all masks in its range have in common, so whole ranges whose masks are not contained in the mask of the column are skipped. The leaves store the lengths of the basis elements and the exponent vectors of their lead monomials contiguously, so the remaining exact divisibility checks do not need to access the basis. The reducer is chosen by the same rule as before (shortest, then newest), which does not depend on the order of the checks, so the chosen reducers do not change. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The reduction of the rows by the known pivots dominates the linear algebra. It is now done for groups of 8 rows at once: their dense rows are interleaved, so each term of a known pivot updates all rows of the group with one vector operation instead of a scattered update per row, and the pivot is read only once per group. This is used for the AB/CD reduction of the dense path and in the sparse path, where each remainder is afterwards reduced by the new pivots as before. Very sparse rows, which are reduced via a heap, the trace learning, normal forms and the final reduction keep the row wise reduction. The interreduction of the new pivots was serial. The new pivots are now interreduced as dense rows (in the sparse path if they fill at least a quarter of the triangle) in blocks of 64 pivots from the right to the left: the rows of a block are interreduced among themselves, then all rows left of the block are reduced by the block in parallel, two rows at once. The result is the reduced row echelon form, so it does not depend on the number of threads. All new kernels have AVX-512, AVX2, NEON and generic versions. The dense kernel for new pivots of the 32 bit linear algebra additionally gets AVX-512 and NEON versions. The AVX-512 and NEON versions are only compile tested. This roughly halves the running time over finite fields, e.g. kat13 over a 31 bit prime field takes 61.9 instead of 131.1 seconds on one thread and 13.6 instead of 25.0 seconds on eight threads. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Reducing the rows by the known pivots in groups of 8 rows slowed down the multi-modular computations over QQ with several threads: in the trace application many matrices have only 12 to 20 rows to be reduced, which gave only two or three groups for eight threads. Moreover, the interleaved kernel always works on all rows of a group, so for sparse rows, which rarely need the same pivots, most of its work is wasted. Groups are now only formed if each thread gets at least two full groups, otherwise the rows are reduced one by one as before; the AB/CD reduction gets its row wise loop back for this. In the sparse path only rows filling at least 1/32 of their remaining columns are put into groups, and the buffers for the groups are only allocated by the threads which form groups. With eight threads, eco13 over QQ now takes 10.7 instead of 10.8 seconds (median of five runs) and kat12 11.3 instead of 12.1 seconds, compared to before the grouping; the speedups over finite fields are kept. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Primes 2^32 <= p < 2^52 now use 64 bit coefficients (cf64_t). Products of coefficients do not fit into 64 bits anymore, so the dense rows are kept reduced to [0,p) and each update is reduced at once. The new kernels in la_ff_64.c use Montgomery reduction modulo 2^64 for scalar code, double precision with FMA for AVX2, AVX512 and NEON, and Montgomery reduction modulo 2^52 with vpmadd52luq/vpmadd52huq for AVX512 IFMA and AVX IFMA. The interleaved rows of the multi-row reduction stay in double precision during the whole reduction, which halves the number of instructions per update. They plug into la_ff_common.c via a few LA_WIDE hooks. export_f4() takes the field characteristic as uint64_t, for p >= 2^32 the input and output coefficients are int64_t. A test computes katsura-6 for the largest supported prime. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The learning phase of the multi-modular computation over QQ reduced each row on its own by the known pivots, since the multi-row kernel did not record which pivots were used for each row. The kernel now marks the pivots used for each row of a group in the reducer bit array of that row, so the learning phase can reduce rows in groups as well. This halves the time of the learning phase on dense systems like katsura or eco, e.g. 132s -> 63s for katsura-13 on one thread. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The final reduction of the basis now uses the usual matrix layout with the reducers on the left and the basis elements as rows to reduce, so the rows are reduced in parallel and in groups and then interreduced by the dense interreduction of new pivots. Since tail terms of a basis element may lie left of its lead term, the reduction starts at the smallest column of the row; rows without terms in known pivot columns become new pivots directly. The rows are built with the parallel batched hash insertion, which removes the serial bottleneck of this step. When learning the trace, dense matrices now also go through the dense linear algebra. The reducers used by each row are recorded in the sparse AB/CD step and the rows becoming pivots are tracked through the dense echelon form, so that the trace can be constructed as before. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…r algebra The probabilistic linear algebra (-l 44) reduces random linear combinations of the rows of each block until one of them reduces to zero. The first combination of a block is still reduced on its own, since most blocks do not give a new pivot. If it gives one and the rows of the block are dense enough, the next combinations are generated in groups of up to eight and reduced together by the known pivots with the multi-row kernel used by the exact linear algebra. Their remainders are then reduced one by one by the new pivots, and the block stops at the first remainder that becomes zero, as before. Sparse blocks keep the reduction of single combinations, there the passes over all columns of the interleaved combinations do not pay off. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
When a random linear form is used over QQ, msolve_trace_qq returned after the first modular step without freeing its data, so the trace, the meta data and the bases of that step leaked. It now frees them like the other early returns of the function. The caller only uses the degree of the minimal polynomial, which is stored before. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The F4 implementation handles primes up to 52 bits, but msolve rejected input files with a characteristic of 2^31 or more. The characteristic is now read as a 64 bit integer, and for characteristics between 2^32 and 2^52 the coefficients are stored as 64 bit integers in the new field cfs_64 of data_gens_ff_t. For such primes msolve computes Groebner bases and leading ideals (-g 1, -g 2) with the 64 bit F4 and prints 64 bit coefficients; solving, saturation, colon ideals, normal forms and signatures stop with an error message. Primes between 2^31 and 2^32 are rejected, neogb would handle them with its 32 bit arithmetic, which does not support them. msolve_julia still takes 32 bit coefficients and now rejects characteristics of 2^31 or more instead of misreading them. This changes data_gens_ff_t in the installed header msolve-data.h, field_char is now int64_t, thus the library version has to be bumped. The new test diff_kat6-52 checks the Groebner basis of katsura-6 over a 52-bit prime, its expected output agrees with the Groebner basis over QQ reduced modulo this prime. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Contributor
Author
|
Note that |
On Windows long is 32 bits, so atol() and strtol() overflow for the 52-bit characteristics and coefficients msolve now accepts. kat6-52 failed there because the characteristic was read with atol(). Use strtoll() for the characteristic and for all coefficients after the first term, which already used it. On Linux the behaviour is unchanged. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
clang requires the rounding mode of _mm256_round_pd() and _mm512_roundscale_pd() to be an integer constant expression; a local const int variable is accepted by gcc but not by clang, which broke the build on macOS for Intel. Use a macro instead. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Changing the field characteristic of export_f4() to uint64_t broke the interface used by AlgebraicSolving.jl, which passes it as a 32 bit integer: the upper half of the 64 bit argument is undefined then. Restore export_f4() with its previous signature and add export_f4_64() for primes up to 2^52, whose coefficients are int64_t for primes >= 2^32. 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.
Faster F4 and Gröbner bases over prime fields up to 52 bits
This branch further speeds up F4 over prime fields and over QQ (multi-modular
tracer), adds F4 for primes up to 52 bits, and lets the
msolvebinarycompute Gröbner bases for characteristics between 2^32 and 2^52.
Changes
monomials when searching reducers.
and interreduce the new pivots in parallel, also when learning the trace.
the trace.
(
-l 44).IFMA and NEON kernels).
msolveaccepts characteristics up to 2^52 for Gröbner bases(
-g 1,-g 2), see below.msolve_trace_qqwhen a random linear form is usedover QQ.
Prime fields with characteristic >= 2^31 in msolve
pwithp < 2^31(as before)or
2^32 <= p < 2^52. Forp >= 2^32msolve only computes Gröbner bases(
-g 1,-g 2), with any linear algebra option and number of threads.Without
-g, or together with saturation, colon ideals, normal forms orsignatures, msolve stops with an error message.
2^31 <= p < 2^32are rejected: neogb handles them with its 32 bitcode, whose intermediate
p*pdoes not fit into 63 bits.diff_kat6-52(katsura-6 over a 52-bit prime).ABI change: library version needs to be bumped
The public headers change the layout of installed structures, so this
breaks the binary interface of both libraries:
libmsolve(msolve-data.h): indata_gens_ff_tthe fieldfield_charis nowint64_t(wasint32_t) and the new fieldint64_t *cfs_64holds the coefficients for characteristics >= 2^31.get_data_from_file()takesint64_t *field_char.libneogb(data.h):bs_t,mat_t,trace_tandmd_thave newfields (64 bit coefficients
cf_64, trace cache), andmd_t.fcis nowuint64_t.libneogbgets the new functionexport_f4_64(): likeexport_f4(),but with a
uint64_tcharacteristicp < 2^52; forp >= 2^32theinput and output coefficients are
int64_t.The functions called by AlgebraicSolving.jl (
export_f4,export_nf,free_f4_julia_result_data,msolve_julia,free_msolve_julia_result_data,export_groebner_qq) keep exactly theirsignatures from master. In particular
export_f4still takes auint32_tcharacteristic: changing it touint64_twould have broken theccallin AlgebraicSolving.jl, which passes a 32 bit integer, so the upperhalf of the argument would be undefined. To use 52-bit primes from Julia,
AlgebraicSolving.jl can call
export_f4_64with aUInt64characteristic and
Int64coefficients.Both libraries use
LT_VERSIONfromconfigure.ac, currently4:0:0.Since interfaces changed incompatibly, it has to become
5:0:0before thenext release. This is not done in this PR.
msolve_juliaThe Julia interface is unchanged: it still takes
const uint32_t field_charand 32 bit coefficients, so it supports characteristic 0 and primes
< 2^31as before. A characteristic>= 2^31is now rejected at the startwith "Field characteristic must be 0 or < 2^31." and
exit(1), like othererrors in this function. Before, such a value was stored in the
int32_tfield and became negative, which was not handled. Supporting 52-bit primes
from Julia needs a new entry point with 64 bit characteristic and
coefficients.
Testing
make check: 73/73.input_files/with-t 1,-t 4,-g 1,-g 2,-P 1and-l 44(372 runs).basis over QQ reduced mod p (6 systems, 4 primes, 1 and 8 threads,
-l 2and-l 44): all equal.findings.
This PR was prepared with the assistance of generative AI (Claude Code);
the changes were reviewed and tested as described above.
Co-Authored-By: Claude Opus 5.5 noreply@anthropic.com