Multivariate polynomial with symbolic values in the coefficients


With Sympy in Python, it is possible to define a polynomial such as

x^2 + a*x^2 + 2/3*y + b*y + x*z + a*x*z

with symbolic variables x,y,z and also symbolic numbers a,b in the coefficients. In other words this is a polynomial on the field Q[a,b] (not sure it is a field but it doesn’t matter). See my blog if this is not clear enough.

Is there a Julia package to do that? I checked MultivariatePolynomials and its ecosystem but didn’t find something like that.

I believe Symbolics utilizes MultivariatePolynomials internally, which you could use. But you can use SymPy using PyCall or PythonCall. See SymPy.jl for a wrapper, similar it seems to your Caracas example.

Thanks. But I was looking for a way faster than SymPy/caracas. With an improved method it takes 20 minutes for my example. This is not important, I’m just curious to know whether there is a Julia way.

Ah maybe there’s a way:

julia> using AbstractAlgebra

julia> using MultivariatePolynomials

julia> using TypedPolynomials

julia> S, (a,b,c) = PolynomialRing(QQ,["a","b","c"])
(Multivariate Polynomial Ring in a, b, c over Rationals, AbstractAlgebra.Generic.MPoly{Rational{BigInt}}[a, b, c])

julia> @polyvar x y
(x, y)

julia> p1 = MultivariatePolynomials.polynomial([a, b, c], [x, x^2*y, x*y])
(b)x²y + (c)xy + (a)x

I’m not able to add two such polynomials:

julia> DynamicPolynomials.@polyvar x y
(x, y)

julia> p = DynamicPolynomials.polynomial([a, b, c], [x, x^2*y, x*y])
(b)x²y + (c)xy + (a)x

julia> 2*p
(2*b)x²y + (2*c)xy + (2*a)x

julia> p + p
ERROR: MethodError: no method matching zero(::Type{AbstractAlgebra.Generic.MPoly{Rational{BigInt}}})

Here is a way to perform an addition but it is not convenient:

julia> DynamicPolynomials.polynomial(vcat(p.a, p.a), vcat(p.x, p.x))
(2*b)x²y + (2*c)xy + (2*a)x
julia> using AbstractAlgebra

julia> S, (a,b,c) = PolynomialRing(QQ, ["a","b","c"]);

julia> T, (x,y) = PolynomialRing(S, ["x","y"]);

julia> p = a*x + b*x^2*y + c*x*y
b*x^2*y + c*x*y + a*x

julia> p+p
2*b*x^2*y + 2*c*x*y + 2*a*x

Ah ok, nice! Thank you.

I can get the monomials, the coefficients, the terms, but how to get the powers?

julia> collect(coeffs(p))
3-element Vector{AbstractAlgebra.Generic.MPoly{Rational{BigInt}}}:

julia> collect(terms(p))
3-element Vector{AbstractAlgebra.Generic.MPoly{AbstractAlgebra.Generic.MPoly{Rational{BigInt}}}}:

julia> collect(monomials(p))
3-element Vector{AbstractAlgebra.Generic.MPoly{AbstractAlgebra.Generic.MPoly{Rational{BigInt}}}}:

I get it!

julia> map(degrees, terms(p))
3-element Vector{Vector{Int64}}:
 [2, 1]
 [1, 1]
 [1, 0]

Wow! The task I do takes 20 minutes with SymPy and only ~1 minute with Julia!!!

Blog post: Fast expansion of a polynomial with Julia.


A bit more idiomatic would be:

julia> collect(exponent_vectors(p))
3-element Vector{Vector{Int64}}:
 [2, 1]
 [1, 1]
 [1, 0]

P.S.: If you replace AbstractAlgebra with Nemo, it should be even faster. You should not have to change anything else.


I’ve written a macro @polyring to automate it, so that you can write

S = @polyring QQ[a, b, c]
T = @polyring S[x, y]

instead of the more verbose

S, (a,b,c) = PolynomialRing(QQ, ["a","b","c"])
T, (x,y) = PolynomialRing(S, ["x","y"])

Macro definition:

macro polyring(r) # use with `using AbstractAlgebra` or `using Nemo`
    @assert r.head == :ref && length(r.args) > 1
    coeff_ring = r.args[1] # if r is QQ[a, b, c], this is :QQ
    vars = r.args[2:end] # this is [:a, :b, :c]
    vars_string = string.(vars) # this is ["a", "b", "c"]
    vars_tuple = Expr(:tuple, vars...) # this is :((a, b, c))
        R, $(esc(vars_tuple)) = PolynomialRing($(esc(coeff_ring)), $vars_string)

P.S. with this definition, you can also just write

@polyring (@polyring QQ[a,b,c])[x, y]

if you only need to use a, b, c, x, y, but not S or T in the code above.

1 Like

Very nice! We actually thought about having something like that, but there are only that many good names for functions/macros and then it becomes a bit challenging to remember which function does which. Also we had some reservation about this, because it (silently) rebinds the variables a, b, c, x, y (without explicit assignment).

I agree that nesting the polynomial ring constructor is at the moment a bit awkward, because one always has to do the [1] business. One the other hand people would complain if we would not return the indeterminant.

FWIW, we have the following macro:

julia> Qx, x, y = @PolynomialRing(QQ, x[1:10], y[1:10])
(Multivariate Polynomial Ring in 20 variables x[1], x[2], x[3], x[4], ..., y[10] over Rational Field, fmpq_mpoly[x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10]], fmpq_mpoly[y[1], y[2], y[3], y[4], y[5], y[6], y[7], y[8], y[9], y[10]])

julia> Qx, x = @PolynomialRing(QQ, x[1:2, 1:2])
(Multivariate Polynomial Ring in x[1,1], x[2,1], x[1,2], x[2,2] over Rational Field, fmpq_mpoly[x[1,1] x[1,2]; x[2,1] x[2,2]])

It is quite handy if one has an example with multidimensional indices for the variables.