# New package PuiseuxPolynomials

**URL:** <https://discourse.julialang.org/t/new-package-puiseuxpolynomials/70948>\
**Category:** Package Announcements\
**Tags:** math, polynomials\
**Created:** [November 4, 2021, 12:54pm UTC](https://discourse.julialang.org/t/new-package-puiseuxpolynomials/70948 "2021-11-04T12:54:51Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![Jean\_Michel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jean_michel/32/8282_2.png) [@Jean\_Michel](https://discourse.julialang.org/u/Jean_Michel)\
**Post date:** [November 4, 2021, 12:54pm UTC](https://discourse.julialang.org/t/new-package-puiseuxpolynomials/70948/1 "2021-11-04T12:54:51Z")

</div>

I want to announce two new packages which go together. Here is the second one  
[PuiseuxPolynomials](https://github.com/jmichel7/PuiseuxPolynomials.jl)

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-repl
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 `Mvp`s 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.

```plaintext
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 `Mvp`s 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-repl
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`:

```plaintext
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-repl
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-repl
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-repl
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-repl
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-repl
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-repl
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-repl
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-repl
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-repl
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, `Mvp`s 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.