New package LaurentPolynomials

I want to announce two packages which go together, one for univariate and the other for multivariate polynomials. The first one is LaurentPolynomials

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 Tis 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.


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⁻¹
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)

julia> print(IOContext(stdout,:TeX=>true),p) # TeXable output (used in Pluto, IJulia)

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> 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}:

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

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

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> Pol(1)
Pol{Int64}: 1

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

julia> scalar(Pol(1))

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

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

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> 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

julia> p(0.5)

# 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> 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, Pols 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> 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).


Although you do mention the Polynomials package, let me state here explicitly that the Polynomials package does have LaurentPolynomial type 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> 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)

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> 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> 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)

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

1 Like

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.

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:

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.