Edit on GitHub

Polynomials.jl

Polynomials.jl is a Julia package that provides basic arithmetic, integration, differentiation, evaluation, and root finding over dense univariate polynomials.

To install the package, run

(v1.2) pkg> add Polynomials

The package can then be loaded into the current session using

using Polynomials

Quick Start

Construction and Evaluation

Construct a polynomial from its coefficients, lowest order first.

julia> Polynomial([1,0,3,4])
Polynomial(1 + 3*x^2 + 4*x^3)

An optional variable parameter can be added.

julia> Polynomial([1,2,3], :s)
Polynomial(1 + 2*s + 3*s^2)

Construct a polynomial from its roots.

julia> fromroots([1,2,3]) # (x-1)*(x-2)*(x-3)
Polynomial(-6 + 11*x - 6*x^2 + x^3)

Evaluate the polynomial p at x.

julia> p = Polynomial([1, 0, -1])
Polynomial(1 - x^2)

julia> p(1)
0

Arithmetic

The usual arithmetic operators are overloaded to work on polynomials, and combinations of polynomials and scalars.

julia> p = Polynomial([1,2])
Polynomial(1 + 2*x)

julia> q = Polynomial([1, 0, -1])
Polynomial(1 - x^2)

julia> 2p
Polynomial(2 + 4*x)

julia> 2 + p
Polynomial(3 + 2*x)

julia> p - q
Polynomial(2*x + x^2)

julia> p * q
Polynomial(1 + 2*x - x^2 - 2*x^3)

julia> q / 2
Polynomial(0.5 - 0.5*x^2)

julia> q ÷ p  # `div`, also `rem` and `divrem`
Polynomial(0.25 - 0.5*x)

Note that operations involving polynomials with different variables will error.

julia> p = Polynomial([1, 2, 3], :x)
Polynomial(1 + 2*x + 3*x^2)

julia> q = Polynomial([1, 2, 3], :s)
Polynomial(1 + 2*s + 3*s^2)

julia> p + q
ERROR: Polynomials must have same variable
[...]

Integrals and Derivatives

Integrate the polynomial p term by term, optionally adding constant term C. The degree of the resulting polynomial is one higher than the degree of p.

julia> integrate(Polynomial([1, 0, -1]))
Polynomial(1.0*x - 0.3333333333333333*x^3)

julia> integrate(Polynomial([1, 0, -1]), 2)
Polynomial(2.0 + 1.0*x - 0.3333333333333333*x^3)

Differentiate the polynomial p term by term. The degree of the resulting polynomial is one lower than the degree of p.

julia> derivative(Polynomial([1, 3, -1]))
Polynomial(3.0 - 2.0*x)

Root-finding

Return the roots (zeros) of p, with multiplicity. The number of roots returned is equal to the order of p. By design, this is not type-stable, the returned roots may be real or complex.

julia> roots(Polynomial([1, 0, -1]))
2-element Array{Float64,1}:
 -1.0
  1.0

julia> roots(Polynomial([1, 0, 1]))
2-element Array{Complex{Float64},1}:
 0.0 - 1.0im
 0.0 + 1.0im

julia> roots(Polynomial([0, 0, 1]))
2-element Array{Float64,1}:
 0.0
 0.0

Fitting arbitrary data

Fit a polynomial (of degree deg) to x and y using polynomial interpolation or a (weighted) least-squares approximation.

using Plots, Polynomials
xs = range(0, 10, length=10)
ys = exp.(-xs)
f = fit(xs, ys)  # fit(xs, ys, k)  for fitting a kth degreee polynomial

scatter(xs, ys, label="Data");
plot!(f, extrema(xs)..., label="Fit");
/home/travis/.julia/packages/GR/yMV3y/src/../deps/gr/bin/gksqt: error while loading shared libraries: libQt5Widgets.so.5: cannot open shared object file: No such file or directory
connect: Connection refused
GKS: can't connect to GKS socket application

GKS: Open failed in routine OPEN_WS
GKS: GKS not in proper state. GKS must be either in the state WSOP or WSAC in routine ACTIVATE_WS

Other bases

A polynomial, e.g. a_0 + a_1 x + a_2 x^2 + ... + a_n x^n, can be seen as a collection of coefficients, [a_0, a_1, ..., a_n], relative to some polynomial basis. The most familiar basis being the standard one: 1, x, x^2, ... Alternative bases are possible. The ChebyshevT polynomials are implemented, as an example. Instead of Polynomial or Polynomial{T}, ChebyshevT or ChebyshevT{T} constructors are used:

julia> p1 = ChebyshevT([1.0, 2.0, 3.0])
ChebyshevT(1.0⋅T_0(x) + 2.0⋅T_1(x) + 3.0⋅T_2(x))

julia> p2 = ChebyshevT{Float64}([0, 1, 2])
ChebyshevT(1.0⋅T_1(x) + 2.0⋅T_2(x))

julia> p1 + p2
ChebyshevT(1.0⋅T_0(x) + 3.0⋅T_1(x) + 5.0⋅T_2(x))

julia> p1 * p2
ChebyshevT(4.0⋅T_0(x) + 4.5⋅T_1(x) + 3.0⋅T_2(x) + 3.5⋅T_3(x) + 3.0⋅T_4(x))

julia> derivative(p1)
ChebyshevT(2.0⋅T_0(x) + 12.0⋅T_1(x))

julia> integrate(p2)
ChebyshevT(0.25⋅T_0(x) - 1.0⋅T_1(x) + 0.25⋅T_2(x) + 0.3333333333333333⋅T_3(x))

julia> convert(Polynomial, p1)
Polynomial(-2.0 + 2.0*x + 6.0*x^2)

julia> convert(ChebyshevT, Polynomial([1.0, 2,  3]))
ChebyshevT(2.5⋅T_0(x) + 2.0⋅T_1(x) + 1.5⋅T_2(x))

!!! Note The older Poly type that this package used prior to v0.7 is implemented as an alternate basis to provide support for older code bases. As of v1.0, this type will be only available by executing using Polynomials.PolyCompat.

Iteration

If its basis is implicit, then a polynomial may be seen as just a vector of coefficients. Vectors or 1-based, but, for convenience, polynomial types are 0-based, for purposes of indexing (e.g. getindex, setindex!, eachindex). Iteration over a polynomial steps through the basis vectors, e.g. a_0, a_1*x, ...

julia> as = [1,2,3,4,5]; p = Polynomial(as);

julia> as[3], p[2], collect(p)[3]
(3, 3, Polynomial(3*x^2))

Contributing

If you are interested in this project, feel free to open an issue or pull request! In general, any changes must be thoroughly tested, allow deprecation, and not deviate too far from the common interface. All PR's must have an updated project version, as well, to keep the continuous delivery cycle up-to-date.