Flint interface

[jodavies]

I have a work-in-progress implementation of an interface to flint for polynomial arithmetic: https://github.com/jodavies/form/tree/flint-prf

So far it is enough to use for PolyRatFun without special modes (expand or modulus).

The results are very good. The univariate and multivariate tests in user.frm are 2x and 10x faster respectively. The forcer test (after making the integrals harder) is also ~1.6x faster. The harder you make these tests the better the benefit (presumably until expressions start being written to disk etc, and then the polynomial code takes a smaller % of the runtime).

Given the excellent performance, I think at least some amount of flint support should go into FORM5.

Currently the "sign convention" of numerator and denominator is not guaranteed to be the same as the poly code due to how it fixes its signs before an additional sort into FORM's highfirst ordering. It might make things easier to break "compatibility" with the poly code here to allow some flexibility in the interface.

Next I will add support for functions like gcd_ and factorization, then @tueda can run polybench, for eg.

[jodavies]

The code can run the GCD part of polybench nicely and shows performance nearly at FLINT's native level. When I try univariate tests FORM is "stuck" at 10^-3 seconds: I think due to TIMER_ returning values in ms?

Ubuntu 20.04's packaged libflint-dev does not contain all the required functions, making the github actions test fail to build. Compiling FLINT as part of the actions would be very time consuming. (I am testing using 24.04 on my own system). The tests pass, I edited one to remove dependence on the sign of the denominator currently.

[jodavies]

Without actually testing, I think 2.6.3 has the required functions.

FUEL's (https://gitlab.com/feynmanintegrals/fuel) use of flint found a bug which was fixed in 3.1.0 related to multiple threads calling a low-level function simultaneously which (I think) can happen in FORM also. So perhaps we should target >= 3.1.0 anyway.

[tueda]

So perhaps we should target >= 3.1.0 anyway.

OK, this is a high hurdle. But, for example, Homebrew supports it.

FYI:

• Ubuntu 24.04LTS - 3.0.1-3.1build1
• vcpkg - 2.9.0 (See also: [flint] update to 3.1.3-p1 microsoft/vcpkg#39293)

[jodavies]

Yes, it would not be so nice for a feature that is "enabled by default" and should be similarly easy to build FORM as eg, GMP and zlib support. At some point I will test some intermediate versions and make sure they run OK (potentially modulo "rare threading issues"). The github runner could be updated to 24.04 without difficulty (or run under 24.04 in addition, and the 20.04 tests simply don't use flint).

[tueda]

For reference: the multithreading bug in flint 3.0.x: flintlib/flint#1652

[jodavies]

For information, my understanding after discussion with a flint author is that if FORM can initialize contexts for the structs that it uses at startup, under mutex, then there should not be any problems later. This doesn't only mean fmpz_mpoly_ctx_init, but also anything that is used internally by those routines.

[jodavies]

Status update: I have made various bugfixes after some more testing, with a test which generates random polynomials and compares the output with the current poly code for everything that is implemented so far.

It was a bit tricky to work out how to reproduce the sign of the result of gcd_, poly does things a bit strangely there. See the comments in jodavies@33e4743 . Currently poly could even be considered buggy, since gcd_(1-x,1-x) gives 1 - x rather than -1 + x (i.e., with positive leading coefficient). gcd_(1-x,(1-x)*(1+x)) indeed gives - 1 + x, agreeing with Mathematica.

What should we do here?

The remaining bits to do are expression (un)factorization, and inverse_. inverse_ is a bit strange: it is not nice that it just terminates form if an inverse is not found, and it completely ignores denominators in the input. Thus inverse_(x^4-1,x^2) gives -1 (agreeing with Mathematica) but inverse_(x^4-1/2,x^2) gives -1 also.

[jodavies]

Maybe we can leave the GCD sign as it has been all of these years.

I've fixed the inverse_ bug in poly, now I should devise a good set of test cases...

(That fix should somehow go into 4.3.2 also. Also I made a few fixes to EndSort where it does not properly return the result size (but not in a way which causes problems in the rest of form).

[tueda]

Can the bug fix be a separate PR?

[jodavies]

Yes, I also found more to fix. poly does not compute this one correctly either:

Local x1 = x^2 + 1;
Local x2 = x^3 + 2*x + 1;
Local x3 = inverse_(x1,x2);
.sort
Local div = div_(x3*x1,x2);
Local rem = rem_(x3*x1,x2);

[tueda]

poly does not compute this one correctly either:

I think the output in this case is correct:

   x3 =
      1/2 - 1/2*x + 1/2*x^2;

   div =
       - 1/2 + 1/2*x;

   rem =
      1;

and the same as:

x1 = x^2 + 1;
x2 = x^3 + 2*x + 1;
{g, {x3, r2}} = PolynomialExtendedGCD[x1, x2, x];
x3 // Expand
div = PolynomialQuotient[x1 x3, x2, x]
rem = PolynomialRemainder[x1 x3, x2, x]

[jodavies]

Ah, on the master branch you are right. Then I have broken something, I'll look into it.