I want to announce two new packages which go together. Here is the second one
PuiseuxPolynomials
I wanted to call this package Mvps, since this is the name of the data
structure it provides, but could not register it under that name (too
short). I thus renamed it to PuiseuxPolynomials (reluctantly, since I did
not want to reserve that name for my use).
This package, which depends only on the packages LaurentPolynomials and ModuleElts, implements Puiseux polynomials, that is linear combinations of monomials of the type x₁^{a₁}… xₙ^{aₙ} where xᵢ are variables and aᵢ are exponents which can be arbitrary rational numbers (we use Puiseux polynomials with cyclotomic coefficients as splitting fields of cyclotomic Hecke algebras), and also implements multivariate rational fractions (type Frac{Mvp{T,Int}}). But I hope it is a good package for ordinary multivariate polynomials if you are only interested in that.
Some functions described below work only with polynomials where variables are raised to integral powers; we will refer to such objects as “Laurent polynomials”; some functions require further that variables are raised only to positive powers: we refer then to “true polynomials” (the numerator and denominator of Frac{Mvp{T,Int}} are true polynomials).
Puiseux polynomials have the parametric type Mvp{M,C} where M is the type of the monomials (Monomial{Int} for Laurent polynomials; Monomial{Rational{Int}} for more general Puisuex polynomials) and C is the type of the coefficients.
We first look at how to make Puiseux polynomials.
@Mvp x₁,…,xₙ
assigns to each Julia name xᵢ an Mvp representing an indeterminate suitable to build multivariate polynomials or rational fractions. Mvp(:x₁) creates the same Mvp without assigning it to variable x₁.
julia> @Mvp x,y
julia> (x+y^-1)^3
Mvp{Int64}: x³+3x²y⁻¹+3xy⁻²+y⁻³
julia> x+Mvp(:z)
Mvp{Int64}: x+z
julia> x^(1//2)
Mvp{Int64,Rational{Int64}}: x½
Mvp(x::Number) returns the multivariate polynomial with only a constant term, equal to x.
It is convenient to create Mvps using variables such as x,y above. To create them more directly, Monomial(:x=>1,:y=>-2) creates the monomial xy⁻², and then Mvp(Monomial(:x=>1,:y=>-2)=>3,Monomial()=>4) creates the Mvp 3xy⁻²+4. This is the way Mvp are printed in another context than the repl, IJulia or Pluto where they display nicely as show as above.
julia> print(3x*y^-2+4)
Mvp(Monomial(:x,:y => -2) => 3,Monomial() => 4)
julia> print(x^(1//2))
Mvp(Monomial(:x => 1//2) => 1)
Only monomials and one-term Mvps can be raised to a non-integral power; the Mvp with one term cm which is c times the monomial m can be raised to a fractional power of denominator d if and only if root(c,d) is defined (this is equivalent to c^{1//d} for floats);
julia> (4x)^(1//2)
Mvp{Int64,Rational{Int64}}: 2x½
julia> (2.0x)^(1//2)
Mvp{Float64,Rational{Int64}}: 1.4142135623730951x½
julia> root(2.0x)
Mvp{Float64,Rational{Int64}}: 1.4142135623730951x½
One may want to define root differently; for instance, in my other package Cycs I define square roots of rationals as cyclotomics; I also have implemented in Cycs arbitrary roots of roots of unity). Using Cycs:
julia> (2x)^(1//2)
Mvp{Cyc{Int64},Rational{Int64}}: √2x½
julia> (E(3)*x)^(2//3)
Mvp{Cyc{Int64},Rational{Int64}}: ζ₉²x⅔
There are various ways to take an Mvp apart. Below are the most direct; look also at the functions coefficient, coefficients, pairs, monomials, variables and powers.
julia> p=3x*y^-2+4
Mvp{Int64}: 3xy⁻²+4
julia> term(p,1)
xy⁻² => 3
julia> term(p,2)
=> 4
julia> length(p) # the number of terms
2
julia> m=first(term(p,1))
Monomial{Int64}:xy⁻²
julia> length(m) # how many variables in m
2
julia> m[:x] # power of x in m
1
julia> m[:y] # power of y in m
-2
The valuation and degree of an Mvp can be inspected globally or variable by variable.
julia> p
Mvp{Int64}: 3xy⁻²+4
julia> variables(p)
2-element Vector{Symbol}:
:x
:y
julia> degree(p),degree(p,:x),degree(p,:y)
(0, 1, 0)
julia> valuation(p),valuation(p,:x),valuation(p,:y)
(-1, 0, -2)
Terms are totally ordered in an Mvp by a monomial ordering (that is, an ordering on monomials so that x<y implies xz<yz for any monomials x,y,z). By default, the ordering is lex. The terms are in decreasing order, so that the first term is the highest. The orderings grlex and grevlex are also implemented.
An Mvp is a scalar if the valuation and degree are 0. The function scalar returns the constant coefficient if the Mvp is a scalar, and nothing otherwise.
Usual arithmetic (+, -, *, ^, /, //, one, isone, zero, iszero, ==) works. Elements of type <:Number are considered as scalars for scalar operations on the coefficients.
julia> p
Mvp{Int64}: 3xy⁻²+4
julia> p^2
Mvp{Int64}: 9x²y⁻⁴+24xy⁻²+16
julia> p/2
Mvp{Float64}: 1.5xy⁻²+2.0
julia> p//2
Mvp{Rational{Int64}}: (3//2)xy⁻²+2//1
One can evaluate an Mvp when setting the value of some variables by using the function call syntax (actually, the keyword syntax for the object used as a function)
julia> p=x+y
Mvp{Int64}: x+y
julia> p(x=2)
Mvp{Int64}: y+2
julia> p(x=2,y=x)
Mvp{Int64}: x+2
Note that an Mvp always evaluates to an Mvp, for consistency. You should use scalar on the result of evaluating all variables to get a number.
julia> p(x=1,y=2)
Mvp{Int64}: 3
julia> scalar(p(x=1,y=2))
3
julia> v=(x^(1//2))(x=2.0)
Mvp{Float64,Rational{Int64}}: 1.4142135623730951
julia> scalar(v)
1.4142135623730951
One can divide an Mvp by another when the division is exact, compute the gcd and lcm of two Mvp.
julia> exactdiv(x^2-y^2,x-y)
Mvp{Int64}: x+y
julia> (x+y)/(2x^2) # or by a monomial
Mvp{Float64}: 0.5x⁻¹+0.5x⁻²y
julia> (x+y)//(2x^2)
Mvp{Rational{Int64}}: (1//2)x⁻¹+(1//2)x⁻²y
julia> (x+y)/(x-y) # otherwise one gets a rational fraction
Frac{Mvp{Int64, Int64}}: (x+y)/(x-y)
Raising a non-monomial Laurent polynomial to a negative power returns a rational fraction. Rational fractions are normalized such that the numerator and denominators are true polynomials prime to each other. They have the arithmetic operations +, - , *, /, //, ^, inv, one, isone, zero, iszero (which can operate between an Mvp or a Number and a Frac{<:Mvp}).
julia> (x+1)^-2
Frac{Mvp{Int64, Int64}}: 1/(x²+2x+1)
julia> x+1/(y+1)
Frac{Mvp{Int64, Int64}}: (xy+x+1)/(y+1)
julia> 1-1/(y+1)
Frac{Mvp{Int64, Int64}}: y/(y+1)
One can evaluate a Frac when setting the value of some variables by using the function call syntax or the value function:
julia> ((x+y)/(x-y))(x=y+1)
Mvp{Float64}: 2.0y+1.0
julia> value((x+y)/(x-y),:x=>y+1;Rational=true) # use // for division
Mvp{Int64}: 2y+1
Frac are dissected using numerator and denominator. They are scalars for broadcasting and can be sorted (have cmp and isless methods).
julia> m=[x+y x-y;x+1 y+1]
2×2 Matrix{Mvp{Int64, Int64}}:
x+y x-y
x+1 y+1
julia> n=inv(Frac.(m))
2×2 Matrix{Frac{Mvp{Int64, Int64}}}:
(-y-1)/(x²-2xy-y²-2y) (x-y)/(x²-2xy-y²-2y)
(x+1)/(x²-2xy-y²-2y) (-x-y)/(x²-2xy-y²-2y)
julia> lcm(denominator.(n))
Mvp{Int64}: x²-2xy-y²-2y
Finally, Mvps have methods conj, adjoint which operate on coefficients, a derivative methods, and methods positive_part, negative_part and bar (useful for Kazhdan-Lusztig theory).
Despite the degree of generality of our polynomials, the speed is not too shabby. For the Fateman test f(f+1) where f=(1+x+y+z+t)^15, we take 4sec. According to the Nemo paper, Sagemath takes 10sec and Nemo takes 1.6sec.