Polynomial manipulation

[BLI]

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.)

[Jean_Michel]

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}}:
 -u²-uw-3ux-3wx+w+x³-3x²
 tux²-2tux+tu+2tx³-4tx²+2tx-1
 twx³-3twx²+3twx-tw-tx⁵+3tx⁴-3tx³+tx²+u+w+x
 d-nu-3nx+n

[nsajko]

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}}:
 -u-w+x²-2x-d⁻¹nux²+2d⁻¹nux-d⁻¹nu-2d⁻¹nx³+4d⁻¹nx²-2d⁻¹nx
 -wx²+2wx-w+x⁴-2x³+x²-d⁻¹nu²x²+2d⁻¹nu²x-d⁻¹nu²-d⁻¹nux⁴-d⁻¹nux³+5d⁻¹nux²-3d⁻¹nux-2d⁻¹nx⁵+2d⁻¹nx⁴+2d⁻¹nx³-2d⁻¹nx²

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.

[thofma]

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
Stacktrace:

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

[nsajko]

@thofma

Is it supposed to handle overflow?

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

[BLI]

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?

[nsajko]

@BLI

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

Definitely a separate topic. Open a new post.

[BLI]

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]

Symbolics.groebner_basis(System)

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.

[nsajko]

@BLI

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

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.

[BLI]

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

or:

-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.

[nsajko]

@BLI

x^2 + x*u+w

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)?

[BLI]

Yeah. Sorry for whitespace. No parenthesis.

[nsajko]

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
end

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

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)...]
end

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))
  end

  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))
  end

  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))
  end

  (s_x, s_u, s_w)
end

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))

[Jean_Michel]

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).

[Jean_Michel]

@nsajko

grobner_basis([P,N-D*n*d^-1])

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}}:
 -u²-uw-3ux-3wx+w+x³-3x²
 -u-w+x²-2x-d⁻¹nux²+2d⁻¹nux-d⁻¹nu-2d⁻¹nx³+4d⁻¹nx²-2d⁻¹nx
 du²+duw+2dux+7du+2dwx+6dw-6dx²+14dx-nu³-nu²w-5nu²x+12nu²-5nuwx+13nuw+21nux+7nu-6nwx²+38nwx-12nw+6nx²+14nx

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…

[Jean_Michel]

@thofma

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

“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.

[sumiya11]

@BLI

So here is the other case:

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

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

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

[sumiya11]

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
true

obtained with the help of Oscar.groebner_basis_with_transformation_matrix.

[nsajko]

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

[BLI]

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