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
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.
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))
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.
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.
Can you please tell me how you create this macro.
iit’s not working in my side
R, x, y = @PolynomialRing(QQ, x[1:3], y[1:3])
What is the error message? And what is the version you are using (copy the output of Pkg.status())?
I am using AbstractAlgebra and When I run R, x, y = @PolynomialRing(QQ, x[1:3], y[1:3]), I am getting
LoadError: UndefVarError: @PolynomialRing not defined.
Same with using Nemo.
The output of Pkg.status() is
Status `~/.julia/environments/v1.8/Project.toml`
⌃ [c3fe647b] AbstractAlgebra v0.30.8
[6e4b80f9] BenchmarkTools v1.3.2
[861a8166] Combinatorics v1.0.2
[86223c79] Graphs v1.8.0
[3e1990a7] Hecke v0.18.14
⌃ [7073ff75] IJulia v1.24.0
[093fc24a] LightGraphs v1.3.5
⌃ [2edaba10] Nemo v0.34.4
[f1435218] Oscar v0.12.1
[eebad327] PkgVersion v0.3.2
⌃ [91a5bcdd] Plots v1.38.14
[f27b6e38] Polynomials v3.2.13
⌃ [295af30f] Revise v3.5.2
⌃ [bcd08a7b] Singular v0.18.3
[2913bbd2] StatsBase v0.34.0
[0c5d862f] Symbolics v5.5.0
Info Packages marked with ⌃ have new versions available and may be upgradable.
Ah, I remember, we renamed it to @polynomial_ring. Sorry for the confusion.