Skip to content

Improve flint performance with multivariate polyratfun #838

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
jodavies wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Improve flint performance with multivariate polyratfun #838

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
a38f5ea
perf: flint: small improvements to flint interface + symbol cmp,norm
jodavies May 19, 2026
41798e1
perf: flint: use FORM lexicographic symbol order for polyratfun
jodavies May 26, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions check/fixes.frm
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1295,6 +1295,8 @@ assert succeeded?
assert result("F") =~ expr("1")
*--#] Issue126 :
*--#[ Issue128 :
* The result check here requires the old prf sign convention:
On oldprfsign;
* Rational arithmetic giving pi_
CF rat;
PolyRatFun rat;
Expand Down
9 changes: 9 additions & 0 deletions doc/manual/statements.tex
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3825,6 +3825,9 @@ \section{off}
statement~\ref{substafactarg}\index{factarg} to the new mode of version 4 or
later in which expressions in the argument of the mentioned function are
completely factored over the rationals. The default is off.}

\leftvitem{3.5cm}{oldprfsign\index{off!oldprfsign}}
\rightvitem{13cm}{\label{staoffoldprfsign} See \ref{staonoldprfsign}.}

\leftvitem{3.5cm}{parallel\index{off!parallel}}
\rightvitem{13cm}{Disallows the running of the program in parallel mode
Expand Down Expand Up @@ -4044,6 +4047,12 @@ \section{on}
version 4. This is a compatibility mode to allow oldprograms that rely on a
specific working of the FactArg statement to still run. The default is
off.}

\leftvitem{3.5cm}{oldprfsign\index{on!oldprfsign}}
\rightvitem{13cm}{\label{staonoldprfsign} Forces the numerator/denominator
sign convention of PolyRatFun~\ref{substapolyratfun} of \FORM\ versions prior
to 5.1, when using FLINT for polynomial arithmetic. Enabling this leads to a
loss of performance in multivariate cases.}

\leftvitem{3.5cm}{parallel\index{on!parallel}}
\rightvitem{13cm}{Allows the running of the program in parallel mode unless
Expand Down
1 change: 1 addition & 0 deletions sources/compcomm.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -135,6 +135,7 @@ static KEYWORDV onoffoptions[] = {
,{"oldfactarg", &(AC.OldFactArgFlag), 1, 0}
,{"memdebugflag", &(AC.MemDebugFlag), 1, 0}
,{"oldgcd", &(AC.OldGCDflag), 1, 0}
,{"oldprfsign", &(AC.OldPRFSignFlag), 1, 0}
,{"innertest", &(AC.InnerTest), 1, 0}
,{"wtimestats", &(AC.WTimeStatsFlag), 1, 0}
,{"sortreallocate", &(AC.SortReallocateFlag), 1, 0}
Expand Down
96 changes: 55 additions & 41 deletions sources/flintinterface.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -89,7 +89,7 @@ void flint::cleanup_master(void) {
#[ flint::divmod_mpoly :
*/
WORD* flint::divmod_mpoly(PHEAD const WORD *a, const WORD *b, const bool return_rem,
const WORD must_fit_term, const var_map_t &var_map) {
const WORD must_fit_term, const var_map_t &var_map, const bool sort_vars) {

flint::mpoly_ctx ctx(var_map.size());
flint::mpoly pa(ctx.d), pb(ctx.d), denpa(ctx.d), denpb(ctx.d);
Expand Down Expand Up @@ -133,10 +133,10 @@ WORD* flint::divmod_mpoly(PHEAD const WORD *a, const WORD *b, const bool return_
const uint64_t prev_size = 0;
const uint64_t out_size = return_rem ?
(uint64_t)flint::to_argument_mpoly(BHEAD NULL, with_arghead, must_fit_term, write, prev_size,
rem.d, var_map, ctx.d, scale.d)
rem.d, var_map, ctx.d, sort_vars, scale.d)
:
(uint64_t)flint::to_argument_mpoly(BHEAD NULL, with_arghead, must_fit_term, write, prev_size,
div.d, var_map, ctx.d, scale.d)
div.d, var_map, ctx.d, sort_vars, scale.d)
;
WORD* res = (WORD *)Malloc1(sizeof(WORD)*out_size, "flint::divrem_mpoly");

Expand All @@ -145,11 +145,11 @@ WORD* flint::divmod_mpoly(PHEAD const WORD *a, const WORD *b, const bool return_
write = true;
if ( return_rem ) {
(uint64_t)flint::to_argument_mpoly(BHEAD res, with_arghead, must_fit_term, write, prev_size,
rem.d, var_map, ctx.d, scale.d);
rem.d, var_map, ctx.d, sort_vars, scale.d);
}
else {
(uint64_t)flint::to_argument_mpoly(BHEAD res, with_arghead, must_fit_term, write, prev_size,
div.d, var_map, ctx.d, scale.d);
div.d, var_map, ctx.d, sort_vars, scale.d);
}

return res;
Expand Down Expand Up @@ -233,7 +233,7 @@ WORD* flint::divmod_poly(PHEAD const WORD *a, const WORD *b, const bool return_r
#[ flint::factorize_mpoly :
*/
WORD* flint::factorize_mpoly(PHEAD const WORD *argin, WORD *argout, const bool with_arghead,
const bool is_fun_arg, const var_map_t &var_map) {
const bool is_fun_arg, const var_map_t &var_map, const bool sort_vars) {

flint::mpoly_ctx ctx(var_map.size());
flint::mpoly arg(ctx.d), den(ctx.d), base(ctx.d);
Expand Down Expand Up @@ -276,8 +276,9 @@ WORD* flint::factorize_mpoly(PHEAD const WORD *argin, WORD *argout, const bool w
uint64_t output_size = 1;

// For finding the highest symbol, in FORM's lexicographic ordering
var_map_t var_map_inv;
for ( auto x: var_map ) {
vector<uint32_t> var_map_inv;
var_map_inv.resize(var_map.size());
for ( auto x : var_map ) {
var_map_inv[x.second] = x.first;
}

Expand All @@ -301,10 +302,10 @@ WORD* flint::factorize_mpoly(PHEAD const WORD *argin, WORD *argout, const bool w
fmpz_mpoly_get_term_exp_si((slong*)base_term_exponents.data(), base.d, (slong)j, ctx.d);

for ( size_t k = 0; k < var_map.size(); k++ ) {
if ( base_term_exponents[k] > 0 && ( var_map_inv.at(k) > max_var ||
( var_map_inv.at(k) == max_var && base_term_exponents[k] > max_pow ) ) ) {
if ( base_term_exponents[k] > 0 && ( var_map_inv[k] > max_var ||
( var_map_inv[k] == max_var && base_term_exponents[k] > max_pow ) ) ) {

max_var = var_map_inv.at(k);
max_var = var_map_inv[k];
max_pow = base_term_exponents[k];
base_sign[i] = fmpz_sgn(fmpz_mpoly_term_coeff_ref(base.d, j, ctx.d));
}
Expand All @@ -321,7 +322,7 @@ WORD* flint::factorize_mpoly(PHEAD const WORD *argin, WORD *argout, const bool w
const bool write = false;
for ( int64_t j = 0; j < exponent; j++ ) {
output_size += (uint64_t)flint::to_argument_mpoly(BHEAD NULL, with_arghead,
is_fun_arg, write, 0, base.d, var_map, ctx.d);
is_fun_arg, write, 0, base.d, var_map, ctx.d, sort_vars);
}
}
}
Expand Down Expand Up @@ -368,7 +369,7 @@ WORD* flint::factorize_mpoly(PHEAD const WORD *argin, WORD *argout, const bool w
const bool write = true;
for ( int64_t j = 0; j < exponent; j++ ) {
argout += flint::to_argument_mpoly(BHEAD argout, with_arghead, is_fun_arg, write,
argout-old_argout, base.d, var_map, ctx.d);
argout-old_argout, base.d, var_map, ctx.d, sort_vars);
}
}
// Final trailing zero to denote the end of the factors.
Expand Down Expand Up @@ -548,6 +549,7 @@ uint64_t flint::from_argument_mpoly(fmpz_mpoly_t poly, fmpz_mpoly_t denpoly, con
// Search for numerical or symbol denominators to create "denpoly".
flint::fmpz den_coeff, tmp;
fmpz_set_si(den_coeff.d, 1);
vector<uint64_t> exponents(var_map.size(), 0);
vector<uint64_t> neg_exponents(var_map.size(), 0);

for ( const WORD* term = args; term < arg_stop; term += term[0] ) {
Expand Down Expand Up @@ -599,7 +601,7 @@ uint64_t flint::from_argument_mpoly(fmpz_mpoly_t poly, fmpz_mpoly_t denpoly, con
const WORD* symbol_stop = term_stop - ABS(coeff_size);
const WORD* t = term;

vector<uint64_t> exponents(var_map.size(), 0);
fill(exponents.begin(), exponents.end(), 0);

t++; // skip over the total size entry
if ( t == symbol_stop ) {
Expand Down Expand Up @@ -913,7 +915,7 @@ void flint::fmpz_set_form(fmpz_t z, UWORD *a, WORD na) {
// If must_fit_term, this should be a TermMalloc buffer. Otherwise Malloc1 the buffer.
// For multi-variate cases.
WORD* flint::gcd_mpoly(PHEAD const WORD *a, const WORD *b, const WORD must_fit_term,
const var_map_t &var_map) {
const var_map_t &var_map, const bool sort_vars) {

flint::mpoly_ctx ctx(var_map.size());
flint::mpoly pa(ctx.d), pb(ctx.d), denpa(ctx.d), denpb(ctx.d), gcd(ctx.d);
Expand Down Expand Up @@ -988,7 +990,7 @@ WORD* flint::gcd_mpoly(PHEAD const WORD *a, const WORD *b, const WORD must_fit_t
const bool write = false;
const uint64_t prev_size = 0;
const uint64_t gcd_size = (uint64_t)flint::to_argument_mpoly(BHEAD NULL,
with_arghead, must_fit_term, write, prev_size, gcd.d, var_map, ctx.d);
with_arghead, must_fit_term, write, prev_size, gcd.d, var_map, ctx.d, sort_vars);

res = (WORD *)Malloc1(sizeof(WORD)*gcd_size, "flint::gcd_mpoly");
}
Expand All @@ -997,7 +999,7 @@ WORD* flint::gcd_mpoly(PHEAD const WORD *a, const WORD *b, const WORD must_fit_t
const bool write = true;
const uint64_t prev_size = 0;
flint::to_argument_mpoly(BHEAD res, with_arghead, must_fit_term, write, prev_size, gcd.d,
var_map, ctx.d);
var_map, ctx.d, sort_vars);

return res;
}
Expand Down Expand Up @@ -1125,20 +1127,25 @@ flint::var_map_t flint::get_variables(const vector <WORD *> &es, const bool with
else {
for ( WORD i = with_arghead ? ARGHEAD:0; with_arghead ? i < e[0]:e[i] != 0; i += e[i] ) {
num_terms++;
if ( i+1 < i+e[i]-ABS(e[i+e[i]-1]) && e[i+1] != SYMBOL ) {
const WORD coeff_size = e[i+e[i]-1];
const WORD symbols_size = e[i] - ABS(coeff_size);
if ( e[i+1] != SYMBOL && 1 < symbols_size ) {
MLOCK(ErrorMessageLock);
MesPrint("ERROR: polynomials and polyratfuns must contain symbols only");
MUNLOCK(ErrorMessageLock);
Terminate(1);
}

for ( WORD j = i+3; j<i+e[i]-ABS(e[i+e[i]-1]); j += 2 ) {
if ( !var_map.count(e[j]) ) {
var_map[e[j]] = num_vars++;
degrees.push_back(e[j+1]);
for ( WORD j = i+3; j < i+symbols_size; j += 2 ) {
const WORD symbol = e[j];
const WORD degree = e[j+1];
auto it = var_map.find(symbol);
if ( it == var_map.end() ) {
var_map[symbol] = num_vars++;
degrees.push_back(degree);
}
else {
degrees[var_map[e[j]]] = MaX(degrees[var_map[e[j]]], e[j+1]);
degrees[it->second] = MaX(degrees[it->second], degree);
}
}
}
Expand Down Expand Up @@ -1283,7 +1290,8 @@ WORD* flint::inverse_poly(PHEAD const WORD *a, const WORD *b, const var_map_t &v
*/
// Return a pointer to a buffer containing the product of the 0-terminated term lists at a and b.
// For multi-variate cases.
WORD* flint::mul_mpoly(PHEAD const WORD *a, const WORD *b, const var_map_t &var_map) {
WORD* flint::mul_mpoly(PHEAD const WORD *a, const WORD *b, const var_map_t &var_map,
const bool sort_vars) {

flint::mpoly_ctx ctx(var_map.size());
flint::mpoly pa(ctx.d), pb(ctx.d), denpa(ctx.d), denpb(ctx.d);
Expand Down Expand Up @@ -1320,12 +1328,12 @@ WORD* flint::mul_mpoly(PHEAD const WORD *a, const WORD *b, const var_map_t &var_
const bool must_fit_term = false;
const uint64_t prev_size = 0;
const uint64_t mul_size = (uint64_t)flint::to_argument_mpoly(BHEAD NULL,
with_arghead, must_fit_term, write, prev_size, pa.d, var_map, ctx.d, den.d);
with_arghead, must_fit_term, write, prev_size, pa.d, var_map, ctx.d, sort_vars, den.d);
res = (WORD*)Malloc1(sizeof(WORD)*mul_size, "flint::mul_mpoly");

write = true;
flint::to_argument_mpoly(BHEAD res, with_arghead, must_fit_term, write, prev_size, pa.d,
var_map, ctx.d, den.d);
var_map, ctx.d, sort_vars, den.d);

return res;
}
Expand Down Expand Up @@ -1387,7 +1395,7 @@ WORD* flint::mul_poly(PHEAD const WORD *a, const WORD *b, const var_map_t &var_m
*/
// Add the multi-variate FORM rational polynomials at t1 and t2. The result is written at out.
void flint::ratfun_add_mpoly(PHEAD const WORD *t1, const WORD *t2, WORD *out,
const var_map_t &var_map) {
const var_map_t &var_map, const bool sort_vars) {

flint::mpoly_ctx ctx(var_map.size());
flint::mpoly gcd(ctx.d), num1(ctx.d), den1(ctx.d), num2(ctx.d), den2(ctx.d);
Expand Down Expand Up @@ -1426,9 +1434,9 @@ void flint::ratfun_add_mpoly(PHEAD const WORD *t1, const WORD *t2, WORD *out,
const bool write = true;
// prev_size + 4, to account for final term size and coeff of "1/1"
out += flint::to_argument_mpoly(BHEAD out, with_arghead, must_fit_term, write, out-args_size+4,
num1.d, var_map, ctx.d);
num1.d, var_map, ctx.d, sort_vars);
out += flint::to_argument_mpoly(BHEAD out, with_arghead, must_fit_term, write, out-args_size+4,
den1.d, var_map, ctx.d);
den1.d, var_map, ctx.d, sort_vars);

*args_size = out - args_size + 1; // The +1 is to include the function ID
AT.WorkPointer = out;
Expand Down Expand Up @@ -1491,7 +1499,8 @@ void flint::ratfun_add_poly(PHEAD const WORD *t1, const WORD *t2, WORD *out,
*/
// Multiply and simplify occurrences of the multi-variate FORM rational polynomials found in term.
// The final term is written in place, with the rational polynomial at the end.
void flint::ratfun_normalize_mpoly(PHEAD WORD *term, const var_map_t &var_map) {
void flint::ratfun_normalize_mpoly(PHEAD WORD *term, const var_map_t &var_map,
const bool sort_vars) {

// The length of the coefficient
const WORD ncoeff = (term + *term)[-1];
Expand Down Expand Up @@ -1548,9 +1557,9 @@ void flint::ratfun_normalize_mpoly(PHEAD WORD *term, const var_map_t &var_map) {
const bool must_fit_term = true;
const bool write = true;
out += flint::to_argument_mpoly(BHEAD out, with_arghead, must_fit_term, write, out-term_size,
num1.d, var_map, ctx.d);
num1.d, var_map, ctx.d, sort_vars);
out += flint::to_argument_mpoly(BHEAD out, with_arghead, must_fit_term, write, out-term_size,
den1.d, var_map, ctx.d);
den1.d, var_map, ctx.d, sort_vars);

*args_size = out - args_size + 1; // The +1 is to include the function ID

Expand Down Expand Up @@ -1779,7 +1788,8 @@ void flint::ratfun_read_poly(const WORD *a, fmpz_poly_t num, fmpz_poly_t den) {
#define IFW(x) { if ( write ) {x;} }
uint64_t flint::to_argument_mpoly(PHEAD WORD *out, const bool with_arghead,
const bool must_fit_term, const bool write, const uint64_t prev_size, const fmpz_mpoly_t poly,
const var_map_t &var_map, const fmpz_mpoly_ctx_t ctx, const fmpz_t denscale) {
const var_map_t &var_map, const fmpz_mpoly_ctx_t ctx, const bool sort_vars,
const fmpz_t denscale) {

// out is modified later, keep the pointer at entry
const WORD* out_entry = out;
Expand All @@ -1797,8 +1807,9 @@ uint64_t flint::to_argument_mpoly(PHEAD WORD *out, const bool with_arghead,
}

// Create the inverse of var_map, so we don't have to search it for each symbol written
var_map_t var_map_inv;
for ( auto x: var_map ) {
vector<uint32_t> var_map_inv;
var_map_inv.resize(var_map.size());
for ( auto x : var_map ) {
var_map_inv[x.second] = x.first;
}

Expand Down Expand Up @@ -1956,8 +1967,10 @@ uint64_t flint::to_argument_mpoly(PHEAD WORD *out, const bool with_arghead,

if ( with_arghead ) {
IFW(*arg_size = out - arg_size);
if ( write ) {
// Sort into form highfirst ordering
if ( write && sort_vars ) {
// Sort into form highfirst ordering, if we have potentially re-ordered the variables by
// degree, in flint::get_variables. Otherwise, we can rely on FLINT's sorting with
// ORD_LEX and the variables ordered in FORM's lexicographic order.
flint::form_sort(BHEAD (WORD*)(out_entry));
}
}
Expand All @@ -1975,13 +1988,13 @@ uint64_t flint::to_argument_mpoly(PHEAD WORD *out, const bool with_arghead,
// If no denscale argument is supplied, just set it to 1 and call the usual function
uint64_t flint::to_argument_mpoly(PHEAD WORD *out, const bool with_arghead,
const bool must_fit_term, const bool write, const uint64_t prev_size, const fmpz_mpoly_t poly,
const var_map_t &var_map, const fmpz_mpoly_ctx_t ctx) {
const var_map_t &var_map, const fmpz_mpoly_ctx_t ctx, const bool sort_vars) {

flint::fmpz tmp;
fmpz_set_ui(tmp.d, 1);

uint64_t ret = flint::to_argument_mpoly(BHEAD out, with_arghead, must_fit_term, write,
prev_size, poly, var_map, ctx, tmp.d);
prev_size, poly, var_map, ctx, sort_vars, tmp.d);

return ret;
}
Expand Down Expand Up @@ -2012,8 +2025,9 @@ uint64_t flint::to_argument_poly(PHEAD WORD *out, const bool with_arghead,
}

// Create the inverse of var_map, so we don't have to search it for each symbol written
var_map_t var_map_inv;
for ( auto x: var_map ) {
vector<uint32_t> var_map_inv;
var_map_inv.resize(var_map.size());
for ( auto x : var_map ) {
var_map_inv[x.second] = x.first;
}

Expand Down
19 changes: 11 additions & 8 deletions sources/flintinterface.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,6 +65,7 @@ extern "C" {
// The bits of std that are needed:
using std::cout;
using std::endl;
using std::fill;
using std::map;
using std::swap;
using std::string;
Expand Down Expand Up @@ -132,10 +133,12 @@ namespace flint {
void cleanup(void);
void cleanup_master(void);

WORD* divmod_mpoly(PHEAD const WORD *, const WORD *, const bool, const WORD, const var_map_t &);
WORD* divmod_mpoly(PHEAD const WORD *, const WORD *, const bool, const WORD, const var_map_t &,
const bool);
WORD* divmod_poly(PHEAD const WORD *, const WORD *, const bool, const WORD, const var_map_t &);

WORD* factorize_mpoly(PHEAD const WORD *, WORD *, const bool, const bool, const var_map_t &);
WORD* factorize_mpoly(PHEAD const WORD *, WORD *, const bool, const bool, const var_map_t &,
const bool);
WORD* factorize_poly(PHEAD const WORD *, WORD *, const bool, const bool, const var_map_t &);

void form_sort(PHEAD WORD *);
Expand All @@ -147,30 +150,30 @@ namespace flint {
WORD fmpz_get_form(fmpz_t, WORD *);
void fmpz_set_form(fmpz_t, UWORD *, WORD);

WORD* gcd_mpoly(PHEAD const WORD *, const WORD *, const WORD, const var_map_t &);
WORD* gcd_mpoly(PHEAD const WORD *, const WORD *, const WORD, const var_map_t &, const bool);
WORD* gcd_poly(PHEAD const WORD *, const WORD *, const WORD, const var_map_t &);

var_map_t get_variables(const vector <WORD *> &, const bool, const bool);

WORD* inverse_poly(PHEAD const WORD *, const WORD *, const var_map_t &);

WORD* mul_mpoly(PHEAD const WORD *, const WORD *, const var_map_t &);
WORD* mul_mpoly(PHEAD const WORD *, const WORD *, const var_map_t &, const bool);
WORD* mul_poly(PHEAD const WORD *, const WORD *, const var_map_t &);

void ratfun_add_mpoly(PHEAD const WORD *, const WORD *, WORD *, const var_map_t &);
void ratfun_add_mpoly(PHEAD const WORD *, const WORD *, WORD *, const var_map_t &, const bool);
void ratfun_add_poly(PHEAD const WORD *, const WORD *, WORD *, const var_map_t &);

void ratfun_normalize_mpoly(PHEAD WORD *, const var_map_t &);
void ratfun_normalize_mpoly(PHEAD WORD *, const var_map_t &, const bool);
void ratfun_normalize_poly(PHEAD WORD *, const var_map_t &);

void ratfun_read_mpoly(const WORD *, fmpz_mpoly_t, fmpz_mpoly_t, const var_map_t &,
fmpz_mpoly_ctx_t);
void ratfun_read_poly(const WORD *, fmpz_poly_t, fmpz_poly_t);

uint64_t to_argument_mpoly(PHEAD WORD *, const bool, const bool, const bool, const uint64_t,
const fmpz_mpoly_t, const var_map_t &, const fmpz_mpoly_ctx_t);
const fmpz_mpoly_t, const var_map_t &, const fmpz_mpoly_ctx_t, const bool);
uint64_t to_argument_mpoly(PHEAD WORD *, const bool, const bool, const bool, const uint64_t,
const fmpz_mpoly_t, const var_map_t &, const fmpz_mpoly_ctx_t, const fmpz_t);
const fmpz_mpoly_t, const var_map_t &, const fmpz_mpoly_ctx_t, const bool, const fmpz_t);
uint64_t to_argument_poly(PHEAD WORD *, const bool, const bool, const bool, const uint64_t,
const fmpz_poly_t, const var_map_t &);
uint64_t to_argument_poly(PHEAD WORD *, const bool, const bool, const bool, const uint64_t,
Expand Down
Loading
Loading