Polynomial manipulation

Thanks a lot for explanation. I will post another example.

I was trying to understand some thermodynamics which included cubic polynomials, and then some expressions which were fractions of polynomials. My expressions were more complex than those in the literature, so I was wondering what I did wrong. But then after some manual playing around, I found that I could reduce my polynomial fractions to those of the literature.

Essentially, I was wondering how some people could be so ingenious as to figure out those simplifications. And then I wondered: is there a Julia tool to do the same thing for me?

(I don’t know if I have time today… I really need to take an exercise walk… will hopefully be back in the evening.

I tried to install Nemo and Groebner via the Julia Pkg manager yesterday + update my packages. Don’t know what is wrong with my Windows 11 laptop… when I move away from the screen and the screen locker kicks in, it seems like the progress in the Julia REPL pauses… even though I have set the power settings to not pause with screen dimmed…so I need to finish the installation.)


My little PuiseuxPolynomials package (which has no dependencies) can do the computation:

julia> using PuiseuxPolynomials

julia> @Mvp t,n,d,x,u,w

julia> P = x^3 - 3x^2 - 3(u + w)*x - u^2 + (1 - u)*w
Mvp{Int64}: -u²-uw-3ux-3wx+w+x³-3x²

julia> N = x^2 - 2x - (u + w)
Mvp{Int64}: -u-w+x²-2x

julia> D = (2x + u) * (x - 1)^2
Mvp{Int64}: ux²-2ux+u+2x³-4x²+2x

julia> grobner_basis([P,d*N-D*n,t*D-1])
4-element Vector{Mvp{Int64, Int64}}:

Given that Puiseux polynomials include Laurent polynomials, which are able to represent some fractions, I hoped that PuiseuxPolynomials.jl would enable a more direct encoding of the system of equations. However, this doesn’t seem to work:

julia> grobner_basis([P,N-D*n*d^-1])
2-element Vector{Mvp{Int64, Int64}}:

I guess grobner_basis only supports the true polynomials, i.e., it doesn’t support Laurent polynomials?

BTW, @Jean_Michel, I made a PR to LaurentPolynomials.jl to fix the build on Nightly Julia.

That looks nice. Two quick question: Is it supposed to handle overflow? I tried it with some polynomials but got wrong results (arithmetic) and sometimes and error:

julia> grobner_basis([P,d*N-D*n,t*D-1])
ERROR: DivideError: integer division error

And regarding “no dependencies”, it also seems to depend on ModuleElts and LaurentPolynomials. Is this normal or did I maybe get an old version?

I believe it’s possible to avoid overflow simply by using BigInt coefficients, instead of Int.

I get an error with Nemo when trying to run it i a Jupyter Notebook (within VSCode):

julia> using Nemo

WARNING: could not import AbstractAlgebra.Solve into Nemo
ERROR: LoadError: UndefVarError: `@include_deprecated_bindings` not defined

It seems to work in the REPL, though.

  • Any clues?

Definitely a separate topic. Open a new post.

1 Like

Just for fun, I tried the code with Symbolics, too:

using Symbolics

@variables u, w, x, d, n, t# t, n, d, x, u, w

P = x^3 - 3*x^2 - 3*(u + w)*x - u^2 + (1 - u)*w
N = x^2 - 2*x - (u + w)
D = (2*x + u) * (x - 1)^2

System = [P, d*N - D*n, t*D - 1]


leading to:

9-element Vector{PolyForm{Real}}:
 -n - d + 3x*n + u*n
 w - 3(x^2) - 3w*x - 3u*x - u*w - (u^2) + (x^3)

Seems to work, too. I get the same result whether I change the order of defining variables.

OK – will check other polynomial now.

FTR, Symbolics.jl actually doesn’t do anything here on its own AFAIK. The polynomials are DynamicPolynomials.jl polynomials, which are supported by Groebner.jl.

1 Like

So here is the other case:

N_a = (x^2 + x*u+w)^2

D_a = (x^2 - 2*x - (u+w))

System_a = [P, d*N_a - D_a*n, t*D_a - 1]

when applying the Gröbner basis leads to:

-n - w*d + 3(x^2)*d + 3u*x*d + (u^2)*d


-n + (-w + 3x^2 + 3xu + u^2)d = 0 \\ \Downarrow \\ \frac{n}{d} = -w + 3x(x+u) + u^2

That is interesting, if correct… the expression I found in the literature is:

\frac{\left(x^{2}+u\cdot x+w\right)\left(2x+u\right)}{x-1}

I will check this.

The whitespace formatting is weird here? Is this supposed to be x^2 + x*u + w or is it supposed to be x^2 + x*(u+w)?

Yeah. Sorry for whitespace. No parenthesis.

The below code might be useful. It helps avoid repeating redundant work and gives multiple possible solutions by switching around the roles of the variables. It doesn’t find the solution you expect for your second problem. It finds these three solutions for your first problem:

n - 9//(u^3 + 6*u^2 + 12*u + 8)*d*x^2 + (3*u + 24)//(u^3 + 6*u^2 + 12*u + 8)*d*x + (-u^2 + 20*u + 27*w + 8)//(u^3 + 6*u^2 + 12*u + 8)*d
n - 1//(x^3 - 3*x^2 + 3*x - 1)*d*u + (-w - 1)//(x^3 - 3*x^2 + 3*x - 1)*d
n - 1//3//(x + 1//3*u - 1//3)*d

The code:

function denominator_constraint(t, d)
  t*d - 1

function numerator_denominator_constraints(t, n, d, num, den)
  nd = den*n - num*d
  (nd, denominator_constraint(t, d))

function problem_system(problem, t, n, d, vars)
  cs = problem.constraints(vars...)
  num = problem.numerator(vars...)
  den = problem.denominator(vars...)
  [cs..., numerator_denominator_constraints(t, n, d, num, den)...]

using Nemo, ParamPunPam

function all_systems_xuw(problem)
  s_x = let
    (R_param, (u, w)) = QQ["u", "w"]
    (_, (t, n, d, x)) = fraction_field(R_param)["t", "n", "d", "x"]
    problem_system(problem, t, n, d, (x, u, w))

  s_u = let
    (R_param, (x, w)) = QQ["x", "w"]
    (_, (t, n, d, u)) = fraction_field(R_param)["t", "n", "d", "u"]
    problem_system(problem, t, n, d, (x, u, w))

  s_w = let
    (R_param, (x, u)) = QQ["x", "u"]
    (_, (t, n, d, w)) = fraction_field(R_param)["t", "n", "d", "w"]
    problem_system(problem, t, n, d, (x, u, w))

  (s_x, s_u, s_w)

problem_1 = (
  constraints = ((x, u, w) -> (x^3 - 3*x^2 - 3*(u + w)*x - u^2 + (1 - u)*w,)),
  numerator = ((x, u, w) -> x^2 - 2*x - (u + w)),
  denominator = ((x, u, w) -> (2*x + u) * (x - 1)^2),

problem_1_solutions = map(paramgb, all_systems_xuw(problem_1))

As usual in Julia, there are two ways to detect overflows (which are not detected by default):

  • use Rational{Int} (fractions detect overflow)
  • use SafeInt (the nice package SaferIntegers is designed to detect overflows and has a rather low overhead compared to normal integers).

You can get slightly better results by changing the monomial order:

julia> grobner_basis([P,N-D*n*d^-1];lt=grlex)
3-element Vector{Mvp{Int64, Int64}}:

I guess grobner_basis only supports the true polynomials, i.e., it doesn’t support Laurent polynomials?

I do not know of a theoretical definition of Grobner bases for Laurent polynomials…

1 Like

“no dependencies” was a misstatement. I meant apart from these two additional packages, which themselves are small and have no other dependencies. The total amount of code in the 3 packages is rather small for the functionality.

Thanks, that’s interesting! Result seems correct:

using Nemo
R, (x,u,w) = QQ["x","u","w"]
N_a = (x^2 + x*u+w)^2
D_a = (x^2 - 2*x - (u+w))
P = x^3 - 3*x^2 - 3*(u + w)*x - u^2 + (1 - u)*w
Result = -w + 3x*(x+u) + u^2
# N_a / D_a = Result (mod P)
@assert mod(N_a, P) == mod(D_a * Result, P)

But I do not know which set of transformations one needs to apply to get this result form the initial expression. Maybe something like Gröbner/Standard Bases Over Fields · Oscar.jl (oscar-system.org) can help, will try

1 Like

Ah, here is the expression of the last result in terms of N_a, D_a, P:

julia> Result == (N_a + (2x + u)*P) // D_a

obtained with the help of Oscar.groebner_basis_with_transformation_matrix.


Apparently there’s yet another package implementing Gröbner bases: PolynomialRings.jl.

1 Like

Package description says “alpha status”, and the last work was done in April 2023??