# Multivariate polynomial with symbolic values in the coefficients

Hello,

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
``````
4 Likes

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}}}:
b
c
a

julia> collect(terms(p))
3-element Vector{AbstractAlgebra.Generic.MPoly{AbstractAlgebra.Generic.MPoly{Rational{BigInt}}}}:
b*x^2*y
c*x*y
a*x

julia> collect(monomials(p))
3-element Vector{AbstractAlgebra.Generic.MPoly{AbstractAlgebra.Generic.MPoly{Rational{BigInt}}}}:
x^2*y
x*y
x
``````

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.

2 Likes

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.

3 Likes

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]
``````

``````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 # 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))
quote
R, \$(esc(vars_tuple)) = PolynomialRing(\$(esc(coeff_ring)), \$vars_string)
R
end
end
``````

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 `` 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, x, x, x, ..., y over Rational Field, fmpq_mpoly[x, x, x, x, x, x, x, x, x, x], fmpq_mpoly[y, y, y, y, y, y, y, y, y, y])

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.