Getting Started With PolynomialRings.jl
Installation
Refer to the Julia website for details on installing Julia. As soon as you have, start it and type [
to get in package mode. Then, run
(v1.1) pkg> add PolynomialRings
to install PolynomialRings
and its dependencies. To test whether it worked, type
julia> using PolynomialRings
julia> @ring! Int[x,y]
@ring(Int64[x,y])
julia> (x + y) * (x - y)
x^2 + -y^2
If you see the same, you are all set!
Overview
Defining polynomial rings
The easiest way to define your polynomial rings is using the @ring!
macro:
julia> R = @ring! Int[x,y]
@ring(Int64[x,y])
This will create a type R
for your polynomials, and it will assign the polynomial x
to the variable x
and similarly for y
.
For a mathematically more pleasing look, try
julia> R = @ring! ℤ[x,y]
@ring(ℤ[x,y])
For entering ℤ, type \bbZ<tab>
in the Julia command line. (Juno and julia-vim
support this as well.) We support ℤ (arbitrary precision integers), ℚ (fractions of arbitrary precision integers), ℝ (arbitrary precision floating point) and ℂ (a + im*b
with both coefficients in ℝ).
If your variables have numbers instead of distinct names, you can use []
to represent that:
julia> R = @ring! ℤ[c[]]
@ring(ℤ[c[]])
This will make available the object c
, which you can use as follows:
julia> c1,c2,c3 = c[1], c[2], c[3]; c1,c2,c3 # or
(c[1], c[2], c[3])
julia> c1,c2,c3 = c[1:3]; c1,c2,c3 # or
(c[1], c[2], c[3])
julia> c1,c2,c3 = c[]; c1,c2,c3 # or
(c[1], c[2], c[3])
julia> # note that the following keeps state:
c1 = c(); c2 = c(); c3 = c(); c1,c2,c3
(c[1], c[2], c[3])
julia> c4 = c(); c5 = c(); c6 = c(); c4,c5,c6
(c[4], c[5], c[6])
Note that you cannot combine named and numbered variables in one ring definition. However, you can let one kind of ring be the base ring for another:
julia> R = @ring! ℤ[c[]][x,y]
@ring(ℤ[c[]][x,y])
A quick way to define a polynomial without first defining the ring is
julia> # will implicitly create the ring Int[x,y]
@polynomial x^2 - y^2
x^2 + -y^2
Arithmetic
The usual ring operations +
,-
,*
,^
work as you'd expect:
julia> @ring! ℤ[x,y]
@ring(ℤ[x,y])
julia> (x + y) * (x - y) == x^2 - y^2
true
We also support reduction operations between polynomials; for that, you can use the standard julia functions div
, rem
and divrem
. For example, divrem(f, g)
. You can also reduce with respect to a set of polynomials, e.g. divrem(f, [g1,g2])
.
For example, in the one-variable case, this is just the Euclidean algorithm:
julia> rem(x^2 - 1, x - 1)
0
If you prefer, you can also use the symbols ÷
for div
and %
for rem
:
julia> (x^2 - 1)÷(x - 1)
x + 1
Variables in your ring vs. variables in your script
It is common to use names such as f
,g
,h
for polynomials and names such as $x,y,z$ for the variables in your ring. For example, you might define
julia> R = @ring! ℤ[x,y,z]
@ring(ℤ[x,y,z])
julia> f = x^2 - x*y
x^2 + -x*y
In this situation, f
is a variable in your script, of type R
.
You might also define
julia> g(x,y) = x^2 - x*y
g (generic function with 1 method)
but this means something different. In this case, x
and y
are arguments to the function, and in its body, they take whatever value you pass to g
. For example:
julia> g(x,y)
x^2 + -x*y
julia> g(y,x)
-x*y + y^2
julia> g(y,y)
0
Maybe by now you wonder about x
and y
: are they Julia variables or just names? The answer is easiest to understand if we look at the @ring
macro. Note that this one does not have the !
in the end:
julia> using PolynomialRings
julia> S = @ring ℤ[x,y]
@ring(ℤ[x,y])
julia> x
ERROR: UndefVarError: x not defined
As you can see, we did define a type S
that contains polynomials with names $x$ and $y$ for the variables. However, in our script, the variable x
doesn't exist. The way to get the variable with name $x$ is to start with the symbol :x
, and convert it to S
. Here's how:
julia> S(:x)
x
The result is an object of type S
, much like how f
was an object of type R
. It represents a polynomial with just one term: $1x$.
Wouldn't it be practical if we would do x = S(:x)
and y = S(:y)
now? That way, we can use the Julia variable x
to refer to the polynomial $x$. Indeed, that's exactly what @ring!
(with !
) does!
In the next chapters, we will often pass variables as arguments. For example, we pass the variable in which we are doing an expansion, or the variable with respect to which we are taking a derivative. In those cases, we pass the variable as a symbol (e.g. :x
) to the function. For example, this works:
julia> diff(x^3, :x)
3*x^2
But this doesn't:
julia> diff(x^3, x)
3*x^2
In some cases, we offer a convenience macro. For example, the @deg
macro:
julia> deg(x^3*y^4, :x)
3
julia> deg(x^3*y^4, :y)
4
julia> @deg x^3*y^4 y
4
Expansions, coefficients, collecting monomials
The rings ℤ[a,b][c]
, ℤ[a,b,c]
and ℤ[b][a,c]
are canonically isomorphic, and we make it easy to switch perspective between them. For example, the different polynomials compare equal (using ==
) and can be easily convert
ed into each other. Type promotions also happen as you'd expect.
julia> R = @ring ℤ[a,b][c]
@ring(ℤ[a,b][c])
julia> T = @ring! ℤ[a,b,c]
@ring(ℤ[a,b,c])
julia> U = @ring ℤ[b][a,c]
@ring(ℤ[b][a,c])
julia> R(a*b + b*c + a*c)
(a + b)*c + a*b
julia> T(a*b + b*c + a*c)
a*b + a*c + b*c
julia> U(a*b + b*c + a*c)
a*c + b*a + b*c
Keep in mind that they don't have equal hash()
values, so don't rely on this for Set{Any}
and Dict{Any}
. Set{R}
and Dict{R}
should work, since type conversion will happen before hashing.
For seeing the constituent parts of a polynomial, use the @expand
macro. You need to specify in which variables you are expanding. After all, $(a+1)bc = abc + bc$, so the result from expanding in $a,b$ and $c$ is different from an expansion in just $b$ and $c$.
julia> @ring! ℤ[a,b,c]
@ring(ℤ[a,b,c])
julia> @expand (a*b*c + b*c) a b c
Base.Generator{Base.Generator{StepRange{Int64,Int64},PolynomialRings.Polynomials.var"#13#14"{@ring(ℤ[a,b,c])}},PolynomialRings.Expansions.var"#21#22"}(PolynomialRings.Expansions.var"#21#22"(), Base.Generator{StepRange{Int64,Int64},PolynomialRings.Polynomials.var"#13#14"{@ring(ℤ[a,b,c])}}(PolynomialRings.Polynomials.var"#13#14"{@ring(ℤ[a,b,c])}(a*b*c + b*c), 1:1:2))
julia> @expand (a*b*c + b*c) b c
Base.Generator{Array{(Term over @ring(ℤ[a]) in @degrevlex(b > c)),1},PolynomialRings.Expansions.var"#21#22"}(PolynomialRings.Expansions.var"#21#22"(), (Term over @ring(ℤ[a]) in @degrevlex(b > c))[(a + 1)*b*c])
For just obtaining a single coefficient, use @coefficient
.
julia> @coefficient (a*b*c + b*c) a^0*b^1*c^1
1
julia> @coefficient (a*b*c + b*c) b^1*c^1
a + 1
There is also corresponding functions expansion()
and coefficient()
. For those, you need to pass the variables as symbols. For example:
julia> coefficient(a*b*c + b*c, (0, 1, 1), :a, :b, :c)
1
julia> coefficient(a*b*c + b*c, (1, 1), :b, :c)
a + 1
Monomial orders
Functions such as leading_term
and divrem
have an implicit understanding of what the monomial order is. By default, if you create a ring, it will be ordered by degree, then reversely lexicographically. This is often called 'degrevlex'.
If you want to use a different order, you can specify that by creating the ring using the polynomial_ring
function:
julia> R,(x,y) = polynomial_ring(:x, :y, monomialorder=:lex)
(@ring(ℚ[x,y]), (x, y))
julia> PolynomialRings.monomialorder(R)
MonomialOrdering{:lex, typeof(@namingscheme((x,y)))}()
Built-in are :lex
, :deglex
and :degrevlex
. It is very easy to define your own order, though, and thanks to Julia's design, this doesn't incur any performance penalty. Read the documentation for PolynomialRings.MonomialOrderings.MonomialOrder
for details.
Gröbner basis computations
For computing a Gröbner basis for a set of polynomials, use the gröbner_basis
function. (For easier typing, groebner_basis
is an alias.)
You typically use this to figure out whether a polynomial is contained in the ideal generated by a set of other polynomials. For example, it is not obvious that $y^2$ is a member of the ideal $(x^5, x^2 + y, xy + y^2)$. Indeed, if one applies rem
, you will not find the expression of $y^2$ in terms of these polynomials:
julia> @ring! ℤ[x,y]
@ring(ℤ[x,y])
julia> I = [x^5, x^2 + y, x*y + y^2]
3-element Array{@ring(ℤ[x,y]),1}:
x^5
x^2 + y
x*y + y^2
julia> rem(y^2, I) # nonzero, even though y^2 ∈ (I)
y^2
However, if you compute a Gröbner basis first, you will:
julia> G = gröbner_basis(I)
3-element Array{@ring(ℚ[x,y]),1}:
x^2 + y
x*y + y^2
y^2
julia> rem(y^2, G) # now, it reduces to zero.
0
If you want to obtain the expression of $y^2$ in these elements, you can first use div
to obtain the (row) vector of factors expressing $y^2$ in G:
julia> div(y^2, G)
1×3 LinearAlgebra.Transpose{@ring(ℚ[x,y]),SparseArrays.SparseVector{@ring(ℚ[x,y]),Int64}}:
0 0 1//1
The gröbner_transformation
function gives a matrix tr
expressing G
in terms of I
:
julia> G, tr = gröbner_transformation(I); tr
3×3 SparseArrays.SparseMatrixCSC{@ring(ℚ[x,y]),Int64} with 7 stored entries:
[1, 1] = -1//1
[3, 1] = -1//1
[1, 2] = x^3 + y^2 + -y
[2, 2] = 1//1
[3, 2] = x^3 + y^2 + -y
[1, 3] = -x^2 + x + -y + 1//1
[3, 3] = -x^2 + x + -y
julia> div(y^2, G) * tr * I # gives back y^2
y^2
In other words, by looking at
julia> div(y^2, G) * tr
1×3 LinearAlgebra.Transpose{@ring(ℚ[x,y]),SparseArrays.SparseVector{@ring(ℚ[x,y]),Int64}}:
1//1 -x^3 + -y^2 + y x^2 + -x + y
we see that $y^2 = 1(x^5) + (y + xy - x^3)(x^2 + y) + -x(xy + y^2)$ which proves that $y^2 \in (I)$.
This is all summarized in the lift
function:
julia> lift(I, y^2)
1×3 LinearAlgebra.Transpose{@ring(ℚ[x,y]),SparseArrays.SparseVector{@ring(ℚ[x,y]),Int64}}:
1//1 -x^3 + -y^2 + y x^2 + -x + y
Using helper variables
(Be sure you understand Variables in your ring vs. variables in your script before reading this section.)
We now get to an important implementation detail. Imagine that you want to write a function that computes a derivative in the following way:
julia> function myderivative(f::RR, varsymbol) where RR <: Polynomial
varvalue = RR(varsymbol)
@ring! Int[ε]
return @coefficient f(; varsymbol => varvalue + ε) ε^1
end
myderivative (generic function with 1 method)
julia> myderivative(x^3 + x^2, :x)
3*x^2 + 2*x
(In fact, diff
is already built-in and has a more efficient implementation, but this example is for educational purposes.)
This works, but what about the following?
julia> @ring! ℚ[ε]
@ring(ℚ[ε])
julia> myderivative(ε^3 + ε^2, :ε)
0//1
This gives a wrong answer because of the naming clash inside myderivative
. You may be tempted to work around this as follows:
julia> function myderivative2(f::RR, varsymbol) where RR <: Polynomial
varvalue = RR(varsymbol)
ε = gensym()
_,(ε_val,) = polynomial_ring(ε)
return coefficient(f(;varsymbol => varvalue + ε_val), (1,), ε)
end
myderivative2 (generic function with 1 method)
julia> myderivative2(ε^3 + ε^2, :ε)
3//1*ε^2 + 2//1*ε
which gives the correct answer. Unfortunately, this is very inefficient:
julia> @time myderivative2(ε^3 + ε^2, :ε);
2.436041 seconds (5.06 M allocations: 265.733 MiB, 4.96% gc time)
julia> @time myderivative2(ε^3 + ε^2, :ε);
2.477522 seconds (5.06 M allocations: 265.800 MiB, 5.98% gc time)
The reason is that variable names are part of the type definition, and Julia specializes functions based on the type of its arguments. In this case, that means that for evaluating the @coefficient
call, and for the substitution call, all the code needs to be compiled every time you call myderivative
.
For this reason, we provide a function formal_coefficient(R)
which yields a variable that's guaranteed not to clash with the ring that you pass as an argument:
julia> function myderivative3(f::RR, varsymbol) where RR <: Polynomial
varvalue = RR(varsymbol)
ε_sym, ε_val = formal_coefficient(typeof(f))
return coefficient(f(;varsymbol => varvalue + ε_val), (1,), ε_sym)
end
myderivative3 (generic function with 1 method)
julia> @time myderivative3(ε^3 + ε^2, :ε); # first time is still slow (compiling)
2.382411 seconds (5.23 M allocations: 275.447 MiB, 4.48% gc time)
julia> @time myderivative3(ε^3 + ε^2, :ε); # but much faster the second time
0.000774 seconds (2.21 k allocations: 52.195 KiB)
julia> @time myderivative3(ε^3 + ε^2, :ε); # and the third
0.000502 seconds (2.21 k allocations: 52.195 KiB)
Free modules (arrays of polynomials)
For practical purposes, a free module (of finite rank) over a ring $R$ is just an array of polynomials in $R$. Many algorithms that work for polynomial rings also work for modules. For example, the function gröbner_basis
works just as well if we pass a vector of vectors of polynomials instead of a vector of polynomials:
julia> G = [[x^5-y,x^4],[x^3+y,y^3]];
julia> GG = gröbner_basis(G)
4-element Array{Array{@ring(ℚ[x,y]),1},1}:
[x^3 + y, y^3]
[x*y + -y^2, x^3*y^3 + -x^5 + -y^4]
[y^3 + y, -x^4*y^3 + -x^3*y^4 + x^6 + x^5*y + x^2*y^3 + x*y^4 + y^5 + -x^4]
[0, x^5*y^3 + -x^7 + -x^4*y + -y^4]
One can then use the functions rem
and div
to express a given vector as an $R$-linear combination of the others.
For these purposes, the leading term of a vector is defined to be the tuple $(i,t)$ where $i$ is the first nonzero index, and $t$ is the leading term of that nonzero element.
Base rings and base restriction / extension
Some operations need a field for a base ring. For example:
julia> R = @ring! ℤ[x]
@ring(ℤ[x])
julia> rem(2x^2, 3x + 1)
ERROR: InexactError: BigInt(2//3)
gives an error because we have to subtract $x^2 + \frac{2}{3}x$, which is not representable in $R$. We offer a convenience function base_extend
to extend to ℚ:
julia> rem(base_extend(2x^2), base_extend(3x + 1))
2//9
If you want, you can also extend to bigger base rings than the quotient field by passing that as an extra parameter. For example:
julia> f = base_extend(x^2 + 1, Complex{Rational{Int}})
x^2 + 1//1 + 0//1*im
julia> div(f, [x - im])
1×1 LinearAlgebra.Transpose{@ring(Complex{Rational{BigInt}}[x]),SparseArrays.SparseVector{@ring(Complex{Rational{BigInt}}[x]),Int64}}:
x + 0//1 + 1//1*im
By the way, if you are looking for an operation like rem
that stays in the integers, have a look at xrem
and friends. For example:
julia> xdiv(2x^2, 3x + 1)
(-9, -6*x + 2)
julia> xrem(2x^2, 3x + 1)
-2
This output signifies that -9
times the first argument is equal to -6x + 2
times the second argument plus -2
. The x
in the names xdiv
and xrem
is intended to represent cross-multiplying the leading coefficients.
Implementation of named and numbered variables
The difference between @ring Int[x1,x2,x3]
and @ring Int[x[]]
is not just the display name of the variables. In terms of implementation, the first version uses a Tuple
to represent the exponents, whereas the second version uses a SparseVector
. This means that for moderate number of variables, the former is often more efficient than the latter as tuples can often remain on the stack, saving allocation and garbage collection overhead. This stops being true when your exponents are very sparse, when the overhead of dealing with all the zeros in the tuple is worse than the overhead of garbage collection.
If you want to transform a set of polynomials from the latter representation to the former, use to_dense_monomials
. This is sometimes beneficial right before doing a computationally expensive operation, e.g. a Gröbner basis computation.
Frequently Asked Questions
Be the first to ask a question! Feel free to open an issue for it.