# New package LaurentPolynomials

**URL:** <https://discourse.julialang.org/t/new-package-laurentpolynomials/70947>\
**Category:** Package Announcements\
**Tags:** math, polynomials\
**Created:** [November 4, 2021, 12:51pm UTC](https://discourse.julialang.org/t/new-package-laurentpolynomials/70947 "2021-11-04T12:51:46Z")\
**Posts on this page:** 6\
**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:51pm UTC](https://discourse.julialang.org/t/new-package-laurentpolynomials/70947/1 "2021-11-04T12:51:46Z")

</div>

I want to announce two packages which go together, one for univariate and the other for multivariate polynomials. The first one is [LaurentPolynomials](https://github.com/jmichel7/LaurentPolynomials.jl)

I wanted to call this package `Pols`, 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 `LaurentPolynomials` following the suggestion  
of the moderator (reluctantly, since I did not want to reserve that name  
for my use).

This package, which depends on no other package, implements univariate Laurent polynomials (type `Pol{T}`), and univariate rational fractions (type `Frac{Pol{T}}`).

The initial motivation was to have a simple way to port GAP polynomials to Julia. The reasons for still having my own package are multiple:

- I need to have a simple and flexible interface, which I hope this provides.
- I need my polynomials to behave well when coefficients are in a ring, in which case I use pseudo-division and subresultant gcd.
- For my uses my polynomials are several times faster than those in the package `Polynomials`.
- I need my polynomials to work as well as possible with coefficients of types `T` where elements have a `zero` method but `T` itself does not have one (because `T` does not contain the necessary information). An example is modular arithmetic with a `BigInt` modulus which thus cannot be part of the type.
- a final justification is that `LaurentPolynomials` is used by `PuiseuxPolynomials`.

Laurent polynomials have the parametric type `Pol{T}`, where `T`is the type of the coefficients. They are constructed by giving a vector of coefficients of type `T`, and a valuation (an `Int`). We call true polynomials those whose valuation is `≥0`.

There is a current variable name (a `Symbol`) used to print polynomials nicely at the repl or in IJulia or Pluto. This name can be changed globally, or just changed for printing a given polynomial. But polynomials do not record individually with which symbol they should be printed.

**Examples**

```julia-repl
julia> Pol(:q) # define symbol used for printing and return Pol([1],1)
Pol{Int64}: q

julia> @Pol q # same as q=Pol(:q) useful to start session with polynomials
Pol{Int64}: q

julia> Pol([1,2]) # valuation is taken to be 0 if omitted
Pol{Int64}: 2q+1

julia> 2q+1 # same polynomial
Pol{Int64}: 2q+1

julia> Pol() # omitting all arguments gives Pol([1],1)
Pol{Int64}: q

julia> p=Pol([1,2,1],-1) # here the valuation is specified to be -1
Pol{Int64}: q+2+q⁻¹

julia> q+2+q^-1 # same polynomial
Pol{Int64}: q+2+q⁻¹

```

```plaintext
julia> print(p) # if not nice printing give an output which can be read back
Pol([1, 2, 1],-1)

# change the variable for printing just this time
julia> print(IOContext(stdout,:limit=>true,:varname=>"x"),p)
x+2+x⁻¹

julia> print(IOContext(stdout,:TeX=>true),p) # TeXable output (used in Pluto, IJulia)
q+2+q^{-1}

```

A polynomial can be taken apart with the functions `valuation`, `degree` and `getindex`. An index `p[i]` gives the coefficient of degree `i` of `p`.

```julia-repl
julia> valuation(p),degree(p)
(-1, 1)

julia> p[0], p[1], p[-1], p[10]
(2, 1, 1, 0)

julia> p[valuation(p):degree(p)]
3-element Vector{Int64}:
 1
 2
 1

julia> p[begin:end] # the same as the above line
3-element Vector{Int64}:
 1
 2
 1

julia> coefficients(p) # the same again
3-element Vector{Int64}:
 1
 2
 1

```

A polynomial is a _scalar_ if the valuation and degree are `0`. The function `scalar` returns the constant coefficient if the polynomial is a scalar, and `nothing` otherwise.

```julia-repl
julia> Pol(1)
Pol{Int64}: 1

julia> convert(Pol{Int},1) # the same thing
Pol{Int64}: 1

julia> scalar(Pol(1))
1

julia> convert(Int,Pol(1)) # the same thing
1

julia> Int(Pol(1)) # the same thing
1

julia> scalar(q+1) # nothing; convert would give an error

```

In arrays `Pol{T}` of different types `T` are promoted to the same type `T` (when the `T` involved have a promotion) and a number is promoted to a polynomial.

Usual arithmetic (`+`, `-`, `*`, `^`, `/`, `//`, `one`, `isone`, `zero`, `iszero`, `==`) works. Elements of type `<:Number` or of type `T` for a `Pol{T}` are considered as scalars for scalar operations on the coefficients.

```julia-repl
julia> derivative(p)
Pol{Int64}: 1-q⁻²

julia> p=(q+1)^2
Pol{Int64}: q²+2q+1

julia> p/2
Pol{Float64}: 0.5q²+1.0q+0.5

julia> p//2
Pol{Rational{Int64}}: (1//2)q²+(1//1)q+1//2

julia> p(1//2) # value of p at 1//2
9//4

julia> p(0.5)
2.25

# interpolation: find p taking values [2.0,1.0,3.0] at [1,2,3]
julia> Pol([1,2,3],[2.0,1.0,3.0])  
Pol{Float64}: 1.5q²-5.5q+6.0

```

Polynomials are scalars for broadcasting. They can be sorted (they have `cmp` and `isless` functions which compare the valuation and the coefficients), they can be keys in a `Dict` (they have a `hash` function).

The functions `divrem`, `div`, `%`, `gcd`, `gcdx`, `lcm`, `powermod` operate between true polynomials over a field, using the polynomial division. Over a ring it is better to use `pseudodiv` and `srgcd` instead of `divrem` and `gcd` (by default `gcd` between integer polynomials delegates to `srgcd`.

`exactdiv` does division (over a field or a ring) when it is exact, otherwise gives an error.

```julia-repl
julia> divrem(q^3+1,2q+1) # changes coefficients to field elements
(0.5q²-0.25q+0.125, 0.875)

julia> divrem(q^3+1,2q+1//1) # case of coeffcients already field elements
((1//2)q²+(-1//4)q+1//8, 7//8)

julia> pseudodiv(q^3+1,2q+1) # pseudo-division keeps the ring
(4q²-2q+1, 7)

julia> (4q^2-2q+1)*(2q+1)+7 # but multiplying back gives a multiple of the polynomial
Pol{Int64}: 8q³+8

julia> exactdiv(q+1,2.0) # exactdiv(q+1,2) would give an error
Pol{Float64}: 0.5q+0.5

```

Finally, `Pol`s have methods `conj`, `adjoint` which operate on coefficients, methods `positive_part`, `negative_part` and `bar` (useful for Kazhdan-Lusztig theory) and a method `randpol` to produce random polynomials.

Inverting polynomials is a way to get a rational fraction `Frac{Pol{T}}`s. Rational fractions are normalized so that the numerator and denominator are true polynomials prime to each other. They have the arithmetic operations `+`, `-` , `*`, `/`, `//`, `^`, `inv`, `one`, `isone`, `zero`, `iszero` (which can operate between a `Pol` or a `Number` and a `Frac{Pol{T}}`).

```julia-repl
julia> a=1/(q+1)
Frac{Pol{Int64}}: 1/(q+1)

julia> Pol(2/a) # convert back to `Pol`
Pol{Int64}: 2q+2

julia> numerator(a)
Pol{Int64}: 1

julia> denominator(a)
Pol{Int64}: q+1

julia> m=[q+1 q+2;q-2 q-3]
2×2 Matrix{Pol{Int64}}:
 q+1 q+2
 q-2 q-3

julia> n=inv(Frac.(m)) # convert to rational fractions to invert the matrix
2×2 Matrix{Frac{Pol{Int64}}}:
 (-q+3)/(2q-1) (-q-2)/(-2q+1)
 (q-2)/(2q-1) (q+1)/(-2q+1)

julia> map(x->x(1),n) # evaluate at 1 the inverse matrix
2×2 Matrix{Float64}:
  2.0 3.0
 -1.0 -2.0

julia> map(x->x(1;Rational=true),n) # evaluate at 1 using //
2×2 Matrix{Rational{Int64}}:
  2//1 3//1
 -1//1 -2//1

```

Rational fractions are also scalars for broadcasting and can be sorted (have `cmp` and `isless` methods).

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [November 4, 2021, 1:46pm UTC](https://discourse.julialang.org/t/new-package-laurentpolynomials/70947/2 "2021-11-04T13:46:29Z")

</div>

Although you do mention the [Polynomials](https://juliamath.github.io/Polynomials.jl) package, let me state here explicitly that the `Polynomials` package does have [LaurentPolynomial type](https://juliamath.github.io/Polynomials.jl/stable/polynomials/polynomial/#Laurent-Polynomial) too. See the code below. You also state that your package is much faster. Can you give some example where we can see it?

```julia
julia> using Polynomials

julia> p = LaurentPolynomial([1,1,1], -1, :q)
LaurentPolynomial(q⁻¹ + 1 + q)

julia> r = LaurentPolynomial([2,3,4], -2, :q)
LaurentPolynomial(2*q⁻² + 3*q⁻¹ + 4)

julia> p*r
LaurentPolynomial(4*q⁻⁴ + 12*q⁻³ + 25*q⁻² + 24*q⁻¹ + 16)

```

---

<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, 4:24pm UTC](https://discourse.julialang.org/t/new-package-laurentpolynomials/70947/3 "2021-11-04T16:24:28Z")

</div>

The motivation was not specially Laurent polynomials (I know it is available elsewhere) but  
correct behavior with coefficients in a ring. With respect to speed, everything seems faster. A random example:

```julia
julia> using Polynomials

julia> p=LaurentPolynomial([1,1,1],-1,:q)
LaurentPolynomial(q⁻¹ + 1 + q)

julia> @btime $p^10;
  917.176 ns (25 allocations: 3.16 KiB)

```

compared to

```julia
julia> using LaurentPolynomials

julia> p=Pol([1,1,1],-1)
Pol{Int64}: x+1+x⁻¹

julia> @btime $p^10;
  241.649 ns (6 allocations: 752 bytes)

```

---

<div class="post-metadata">

**Author:** ![jishnub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jishnub/32/33620_2.png) [@jishnub](https://discourse.julialang.org/u/jishnub)\
**Post date:** [November 9, 2021, 1:25pm UTC](https://discourse.julialang.org/t/new-package-laurentpolynomials/70947/4 "2021-11-09T13:25:32Z")

</div>

Would it be possible to include the performance improvements in `Polynomials`, or is there a fundamental difference in the approaches?

---

<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 9, 2021, 6:06pm UTC](https://discourse.julialang.org/t/new-package-laurentpolynomials/70947/5 "2021-11-09T18:06:46Z")

</div>

I have no idea I did not look at the code of `Polynomials`. In the other hand my code is I think easy to read and short.

---

<div class="post-metadata">

**Author:** ![j\_verzani](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/j_verzani/32/8551_2.png) [@j\_verzani](https://discourse.julialang.org/u/j_verzani)\
**Post date:** [November 9, 2021, 6:36pm UTC](https://discourse.julialang.org/t/new-package-laurentpolynomials/70947/6 "2021-11-09T18:36:08Z")

</div>

The `LaurentPolynomials` approach is more direct, and in initial testing speedier to construct polynomials as expected (the constructor does fewer checks, there is no copying involved, …) The `Polynomials` package is more generic, so different implementations are possible with different backends (a Dictionary for sparse polynomials, a tuple for immutable polynomials, …). I ran a very primitive benchmarking script I have for looking for regressions in `Polynomials` and got these numbers:

```julia
test	min P	IP	SP	LP	LP′	
p(x)	0.0ns >1000	1.0	>1000	>1000	>1000	
p+s	6.3ns 7.25	1.0	292.74	12.03	18.43	
p+q	11.8ns 8.02	1.0	204.25	36.25	8.02	
p+c	6.4ns 8.22	1.0	289.02	12.91	12.59	
p*s	55.9ns 1.0	119.3	41.08	5.74	2.61	
p*q	1.0μs 1.62	1.0	53.88	6.19	1.05	
p*c	66.3ns 1.0	100.32	34.35	4.83	1.92	
||q||	14.3ns 12.07	1.0	40.05	16.29	12.32	

```

They kind of surprised me. Each column is scaled by the speediest one, which in most cases are ImmutablePolynomials which take advantage of generated functions. The `Pol` constructor is the last column, labeled `LP'` I was expecting it to be faster than the `Polynomial` constructor here but they are about the same, save for the addition of a polynomial and a scalar (`s` is a scalar, `c` a constant polynomial). The `Pol` constructor is faster than the LP column, which is the `LaurentPolynomial` constructor for `Polynomials`. So at least until some cribbing from the LaurentPolynomial package is integrated into the Polynomials package, the `Pol` type has real performance benefits for Laurent Polynomials.