diff --git a/dev/.documenter-siteinfo.json b/dev/.documenter-siteinfo.json index a5d8bd47..280b5acc 100644 --- a/dev/.documenter-siteinfo.json +++ b/dev/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.10.2","generation_timestamp":"2024-04-06T22:56:17","documenter_version":"1.3.0"}} \ No newline at end of file +{"documenter":{"julia_version":"1.10.2","generation_timestamp":"2024-04-14T22:26:36","documenter_version":"1.3.0"}} \ No newline at end of file diff --git a/dev/chebs.svg b/dev/chebs.svg index b4ddca0b..46998a4f 100644 --- a/dev/chebs.svg +++ b/dev/chebs.svg @@ -1,99 +1,99 @@ - + - - + - - + - - - - - - - - - - - - - - - - - - - - - - - - - -
Polynomials.AbstractUnivariatePolynomial — Type
AbstractUnivariatePolynomial{B,T,X} <: AbstractPolynomial{T,X}
 AbstractDenseUnivariatePolynomial{B,T,X} <: AbstractUnivariatePolynomial{B,T,X}
-AbstractLaurentUnivariatePolynomial{B,T,X} <: AbstractUnivariatePolynomial{B,T,X}

Abstract container types for polynomials with an explicit basis, B. AbstractDenseUnivariatePolynomial is for 0-based polynomials; AbstractLaurentUnivariatePolynomial is for polynomials with possibly negative powers of the indeterminate.

source
Polynomials.AbstractBasis — Type
AbstractBasis

Abstract type for specifying a polynomial basis.

source
+AbstractLaurentUnivariatePolynomial{B,T,X} <: AbstractUnivariatePolynomial{B,T,X}

Abstract container types for polynomials with an explicit basis, B. AbstractDenseUnivariatePolynomial is for 0-based polynomials; AbstractLaurentUnivariatePolynomial is for polynomials with possibly negative powers of the indeterminate.

source
Polynomials.AbstractBasis — Type
AbstractBasis

Abstract type for specifying a polynomial basis.

source
diff --git a/dev/extensions/index.html b/dev/extensions/index.html index c3c493cf..475c0813 100644 --- a/dev/extensions/index.html +++ b/dev/extensions/index.html @@ -10,4 +10,4 @@ c = [z(2d - 1) for i in 1:m]

then we can exploit its mutability with

MA.operate!(MA.add_mul, c, A, b)

to reduce the allocation down to 48 bytes due to 3 allocations.

These remaining allocations are due to the BigInt buffer used to store the result of intermediate multiplications. This buffer can be preallocated with:

buffer = MA.buffer_for(MA.add_mul, typeof(c), typeof(A), typeof(b))
 MA.buffered_operate!(buffer, MA.add_mul, c, A, b)

then the second line is allocation-free.

The MA.@rewrite macro rewrite an expression into an equivalent code that exploit the mutability of the intermediate results. For instance

MA.@rewrite(A1 * b1 + A2 * b2)

is rewritten into

c = MA.operate!(MA.add_mul, MA.Zero(), A1, b1)
 MA.operate!(MA.add_mul, c, A2, b2)

which is equivalent to

c = MA.operate(*, A1, b1)
-MA.mutable_operate!(MA.add_mul, c, A2, b2)
Note

Note that currently, only the Polynomial and Polynomials.PnPolynomial types implement the API and only

part of it is implemented

PolyCompat

While not an extension, the older Poly type that this package used prior to v0.7 is implemented as an alternate basis and provided on an opt-in bases by executing using Polynomials.PolyCompat. This is to provide support for older code bases.

+MA.mutable_operate!(MA.add_mul, c, A2, b2)
Note

Note that currently, only the Polynomial and Polynomials.PnPolynomial types implement the API and only

part of it is implemented

PolyCompat

While not an extension, the older Poly type that this package used prior to v0.7 is implemented as an alternate basis and provided on an opt-in bases by executing using Polynomials.PolyCompat. This is to provide support for older code bases.

diff --git a/dev/index.html b/dev/index.html index cd39215e..ed6b888c 100644 --- a/dev/index.html +++ b/dev/index.html @@ -385,4 +385,4 @@ P = FactoredPolynomial p,q = fromroots(P, [1,2,3]), fromroots(P, [2,3,3,0]) plot(p//q)
┌ Warning: Increasing Δ, terminating search
-└ @ Polynomials.Multroot ~/work/Polynomials.jl/Polynomials.jl/src/polynomials/multroot.jl:283

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.

+└ @ Polynomials.Multroot ~/work/Polynomials.jl/Polynomials.jl/src/polynomials/multroot.jl:283

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.

diff --git a/dev/objects.inv b/dev/objects.inv index fe2381e6..9db3d9bc 100644 Binary files a/dev/objects.inv and b/dev/objects.inv differ diff --git a/dev/polyfit.svg b/dev/polyfit.svg index ed9b749c..7566a37f 100644 --- a/dev/polyfit.svg +++ b/dev/polyfit.svg @@ -1,103 +1,103 @@ - + - - + - - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + - - - - - + - - + diff --git a/dev/polynomials/chebyshev/index.html b/dev/polynomials/chebyshev/index.html index e499a924..42b2cd1e 100644 --- a/dev/polynomials/chebyshev/index.html +++ b/dev/polynomials/chebyshev/index.html @@ -14,7 +14,7 @@ -4.5 julia> Polynomials.evalpoly(5.0, p, false) # bypasses the domain check done in p(5.0) -2088.0

The latter shows how to evaluate a ChebyshevT polynomial outside of its domain, which is [-1,1]. (For newer versions of Julia, evalpoly is an exported function from Base with methods extended in this package, so the module qualification is unnecessary.

Note

The Chebyshev polynomials are also implemented in ApproxFun, ClassicalOrthogonalPolynomials.jl, FastTransforms.jl, and SpecialPolynomials.jl.

source

The ChebyshevT type holds coefficients representing the polynomial $a_0 T_0 + a_1 T_1 + ... + a_n T_n$.

For example, the basis polynomial $T_4$ can be represented with ChebyshevT([0,0,0,0,1]).

Conversion

ChebyshevT can be converted to Polynomial and vice-versa.

julia> c = ChebyshevT([1, 0, 3, 4])
+2088.0

The latter shows how to evaluate a ChebyshevT polynomial outside of its domain, which is [-1,1]. (For newer versions of Julia, evalpoly is an exported function from Base with methods extended in this package, so the module qualification is unnecessary.

Note

The Chebyshev polynomials are also implemented in ApproxFun, ClassicalOrthogonalPolynomials.jl, FastTransforms.jl, and SpecialPolynomials.jl.

source

The ChebyshevT type holds coefficients representing the polynomial $a_0 T_0 + a_1 T_1 + ... + a_n T_n$.

For example, the basis polynomial $T_4$ can be represented with ChebyshevT([0,0,0,0,1]).

Conversion

ChebyshevT can be converted to Polynomial and vice-versa.

julia> c = ChebyshevT([1, 0, 3, 4])
 ChebyshevT(1⋅T_0(x) + 3⋅T_2(x) + 4⋅T_3(x))
 
 
@@ -22,4 +22,4 @@
 Polynomial(-2.0 - 12.0*x + 6.0*x^2 + 16.0*x^3)
 
 julia> convert(ChebyshevT, p)
-ChebyshevT(1.0⋅T_0(x) + 3.0⋅T_2(x) + 4.0⋅T_3(x))
+ChebyshevT(1.0⋅T_0(x) + 3.0⋅T_2(x) + 4.0⋅T_3(x)) diff --git a/dev/polynomials/polynomial/index.html b/dev/polynomials/polynomial/index.html index c9221e24..ae20dcfe 100644 --- a/dev/polynomials/polynomial/index.html +++ b/dev/polynomials/polynomial/index.html @@ -8,7 +8,7 @@ Polynomial(1 + 2*s + 3*s^2) julia> one(Polynomial) -Polynomial(1.0)source

Immutable Polynomial

Polynomials.ImmutablePolynomial — Type
ImmutablePolynomial{T, X, N}(coeffs)

Construct an immutable (static) polynomial from its coefficients a₀, a₁, …, aₙ, lowest order first, optionally in terms of the given variable x where x can be a character, symbol, or string.

If $p = a_n x^n + \ldots + a_2 x^2 + a_1 x + a_0$, we construct this through ImmutablePolynomial((a_0, a_1, ..., a_n)) (assuming a_n ≠ 0). As well, a vector or number can be used for construction.

The usual arithmetic operators are overloaded to work with polynomials as well as with combinations of polynomials and scalars. However, operations involving two non-constant polynomials of different variables causes an error. Unlike other polynomials, setindex! is not defined for ImmutablePolynomials.

As the degree of the polynomial (+1) is a compile-time constant, several performance improvements are possible. For example, immutable polynomials can take advantage of faster polynomial evaluation provided by evalpoly from Julia 1.4; similar methods are also used for addition and multiplication.

However, as the degree is included in the type, promotion between immutable polynomials can not promote to a common type. As such, they are precluded from use in rational functions.

Note

ImmutablePolynomial is not axis-aware, and it treats coeffs simply as a list of coefficients with the first index always corresponding to the constant term.

Examples

julia> using Polynomials
+Polynomial(1.0)
source

Immutable Polynomial

Polynomials.ImmutablePolynomial — Type
ImmutablePolynomial{T, X, N}(coeffs)

Construct an immutable (static) polynomial from its coefficients a₀, a₁, …, aₙ, lowest order first, optionally in terms of the given variable x where x can be a character, symbol, or string.

If $p = a_n x^n + \ldots + a_2 x^2 + a_1 x + a_0$, we construct this through ImmutablePolynomial((a_0, a_1, ..., a_n)) (assuming a_n ≠ 0). As well, a vector or number can be used for construction.

The usual arithmetic operators are overloaded to work with polynomials as well as with combinations of polynomials and scalars. However, operations involving two non-constant polynomials of different variables causes an error. Unlike other polynomials, setindex! is not defined for ImmutablePolynomials.

As the degree of the polynomial (+1) is a compile-time constant, several performance improvements are possible. For example, immutable polynomials can take advantage of faster polynomial evaluation provided by evalpoly from Julia 1.4; similar methods are also used for addition and multiplication.

However, as the degree is included in the type, promotion between immutable polynomials can not promote to a common type. As such, they are precluded from use in rational functions.

Note

ImmutablePolynomial is not axis-aware, and it treats coeffs simply as a list of coefficients with the first index always corresponding to the constant term.

Examples

julia> using Polynomials
 
 julia> ImmutablePolynomial((1, 0, 3, 4))
 ImmutablePolynomial(1 + 3*x^2 + 4*x^3)
@@ -17,7 +17,7 @@
 ImmutablePolynomial(1 + 2*s + 3*s^2)
 
 julia> one(ImmutablePolynomial)
-ImmutablePolynomial(1.0)
Note

This was modeled after StaticUnivariatePolynomials by @tkoolen.

source

Sparse Polynomial

Polynomials.SparsePolynomial — Type
SparsePolynomial{T, X}(coeffs::Dict{Int,T})

Polynomials in the standard basis backed by a dictionary holding the non-zero coefficients. For polynomials of high degree, this might be advantageous.

Examples:

julia> using Polynomials
+ImmutablePolynomial(1.0)
Note

This was modeled after StaticUnivariatePolynomials by @tkoolen.

source

Sparse Polynomial

Polynomials.SparsePolynomial — Type
SparsePolynomial{T, X}(coeffs::Dict{Int,T})

Polynomials in the standard basis backed by a dictionary holding the non-zero coefficients. For polynomials of high degree, this might be advantageous.

Examples:

julia> using Polynomials
 
 julia> P  = SparsePolynomial;
 
@@ -40,7 +40,7 @@
 SparsePolynomial(-1.0 + 1.0*x^1000000000)
 
 julia> p(1)
-0.0
Note

SparsePolynomial is an alias for MutableSparsePolynomial{StandardBasis}.

source

Laurent Polynomial

Polynomials.LaurentPolynomial — Type
LaurentPolynomial{T,X}(coeffs::AbstractVector, [m::Integer = 0], [var = :x])

A Laurent polynomial is of the form a_{m}x^m + ... + a_{n}x^n where m,n are integers (not necessarily positive) with m <= n.

The coeffs specify a_{m}, a_{m-1}, ..., a_{n}. The argument m represents the lowest exponent of the variable in the series, and is taken to be zero by default.

Laurent polynomials and standard basis polynomials promote to Laurent polynomials. Laurent polynomials may be converted to a standard basis polynomial when m >= 0,

Integration will fail if there is a x⁻¹ term in the polynomial.

Note

LaurentPolynomial is axis-aware, unlike the other polynomial types in this package.

Examples:

julia> using Polynomials
+0.0
Note

SparsePolynomial is an alias for MutableSparsePolynomial{StandardBasis}.

source

Laurent Polynomial

Polynomials.LaurentPolynomial — Type
LaurentPolynomial{T,X}(coeffs::AbstractVector, [m::Integer = 0], [var = :x])

A Laurent polynomial is of the form a_{m}x^m + ... + a_{n}x^n where m,n are integers (not necessarily positive) with m <= n.

The coeffs specify a_{m}, a_{m-1}, ..., a_{n}. The argument m represents the lowest exponent of the variable in the series, and is taken to be zero by default.

Laurent polynomials and standard basis polynomials promote to Laurent polynomials. Laurent polynomials may be converted to a standard basis polynomial when m >= 0,

Integration will fail if there is a x⁻¹ term in the polynomial.

Note

LaurentPolynomial is axis-aware, unlike the other polynomial types in this package.

Examples:

julia> using Polynomials
 
 julia> P = LaurentPolynomial;
 
@@ -87,7 +87,7 @@
 Polynomial(x)
 
 julia> x^degree(p) * p(x⁻¹) # reverses  coefficients
-LaurentPolynomial(3.0 + 2.0*x + 1.0*x²)
source

Factored Polynomial

Polynomials.FactoredPolynomial — Type
FactoredPolynomial{T,X}

A polynomial type that stores its data in a dictionary whose keys are the roots and whose values are the respective multiplicities along with a leading coefficient.

The structure is utilized for scalar multiplication, polynomial multiplication and powers, the finding of roots, and the identification of a greatest common divisor. For other operations, say addition, the operation is done after converting to the Polynomial type then converting back. (This requires the identification of the roots, so is subject to numeric issues.)

Examples

julia> using Polynomials
+LaurentPolynomial(3.0 + 2.0*x + 1.0*x²)
source

Factored Polynomial

Polynomials.FactoredPolynomial — Type
FactoredPolynomial{T,X}

A polynomial type that stores its data in a dictionary whose keys are the roots and whose values are the respective multiplicities along with a leading coefficient.

The structure is utilized for scalar multiplication, polynomial multiplication and powers, the finding of roots, and the identification of a greatest common divisor. For other operations, say addition, the operation is done after converting to the Polynomial type then converting back. (This requires the identification of the roots, so is subject to numeric issues.)

Examples

julia> using Polynomials
 
 julia> p = FactoredPolynomial(Dict([0=>1, 1=>2, 3=>4]))
 FactoredPolynomial(x * (x - 3)⁴ * (x - 1)²)
@@ -111,7 +111,7 @@
 FactoredPolynomial((x - 4.0000000000000036) * (x - 2.9999999999999942) * (x - 1.0000000000000002) * (x - 2.0000000000000018))
 
 julia> map(x->round(x, digits=12), q) # map works over factors and leading coefficient -- not coefficients in the standard basis
-FactoredPolynomial((x - 4.0) * (x - 2.0) * (x - 3.0) * (x - 1.0))
source

Rational Function

Polynomials.RationalFunction — Type
RationalFunction(p::AbstractPolynomial, q::AbstractPolynomial)
+FactoredPolynomial((x - 4.0) * (x - 2.0) * (x - 3.0) * (x - 1.0))
source

Rational Function

Polynomials.RationalFunction — Type
RationalFunction(p::AbstractPolynomial, q::AbstractPolynomial)
 p // q

Create a rational expression (p//q) from the two polynomials.

Common factors are not cancelled by the constructor, as they are for the base Rational type. The lowest_terms function attempts that operation.

For purposes of iteration, a rational function is treated like a two-element container.

Examples

julia> using Polynomials
 
 julia> p,q = fromroots(Polynomial, [1,2,3]), fromroots(Polynomial, [2,3,4])
@@ -133,4 +133,4 @@
 (36 - 132*x + 193*x^2 - 144*x^3 + 58*x^4 - 12*x^5 + x^6) // (576 - 1248*x + 1108*x^2 - 516*x^3 + 133*x^4 - 18*x^5 + x^6)
 
 julia> derivative(pq)
-(-108 + 180*x - 111*x^2 + 30*x^3 - 3*x^4) // (576 - 1248*x + 1108*x^2 - 516*x^3 + 133*x^4 - 18*x^5 + x^6)
Note

The RationalFunctions.jl package was a helpful source of ideas.

Note

The ImmutablePolynomial type can not be used for rational functions, as the type requires the numerator and denominator to have the exact same type.

source
+(-108 + 180*x - 111*x^2 + 30*x^3 - 3*x^4) // (576 - 1248*x + 1108*x^2 - 516*x^3 + 133*x^4 - 18*x^5 + x^6)
Note

The RationalFunctions.jl package was a helpful source of ideas.

Note

The ImmutablePolynomial type can not be used for rational functions, as the type requires the numerator and denominator to have the exact same type.

source diff --git a/dev/rational_function.svg b/dev/rational_function.svg index 23f8bb52..79aecd21 100644 --- a/dev/rational_function.svg +++ b/dev/rational_function.svg @@ -1,93 +1,93 @@ - + - - + - - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + diff --git a/dev/reference/index.html b/dev/reference/index.html index 2a34cf9b..9819f772 100644 --- a/dev/reference/index.html +++ b/dev/reference/index.html @@ -18,7 +18,7 @@ Polynomial(1 + 2*x - x^2 - 2*x^3) julia> q / 2 -Polynomial(0.5 - 0.5*x^2)

Inspection

Polynomials.degree — Function
degree(::AbstractPolynomial)

Return the degree of the polynomial, i.e. the highest exponent in the polynomial that has a nonzero coefficient.

For standard basis polynomials the degree of the zero polynomial is defined to be $-1$. For Laurent type polynomials, this is 0, or lastindex(p). The firstindex method gives the smallest power of the indeterminate for the polynomial. The default method assumes the basis polynomials, βₖ, have degree k.

source
Base.length — Function
length(::AbstractPolynomial)

The length of the polynomial.

source
Base.size — Function
size(::AbstractPolynomial, [i])

Returns the size of the polynomials coefficients, along axis i if provided.

source
Polynomials.domain — Function
Polynomials.domain(::Type{<:AbstractPolynomial})

Returns the domain of the polynomial.

source
Polynomials.mapdomain — Function
mapdomain(::Type{<:AbstractPolynomial}, x::AbstractArray)
+Polynomial(0.5 - 0.5*x^2)

Inspection

Polynomials.degree — Function
degree(::AbstractPolynomial)

Return the degree of the polynomial, i.e. the highest exponent in the polynomial that has a nonzero coefficient.

For standard basis polynomials the degree of the zero polynomial is defined to be $-1$. For Laurent type polynomials, this is 0, or lastindex(p). The firstindex method gives the smallest power of the indeterminate for the polynomial. The default method assumes the basis polynomials, βₖ, have degree k.

source
Base.length — Function
length(::AbstractPolynomial)

The length of the polynomial.

source
Base.size — Function
size(::AbstractPolynomial, [i])

Returns the size of the polynomials coefficients, along axis i if provided.

source
Polynomials.domain — Function
Polynomials.domain(::Type{<:AbstractPolynomial})

Returns the domain of the polynomial.

source
Polynomials.mapdomain — Function
mapdomain(::Type{<:AbstractPolynomial}, x::AbstractArray)
 mapdomain(::AbstractPolynomial, x::AbstractArray)

Given values of x that are assumed to be unbounded (-∞, ∞), return values rescaled to the domain of the given polynomial.

Examples

julia> using Polynomials
 
 julia> x = -10:10
@@ -26,15 +26,15 @@
 
 julia> extrema(mapdomain(ChebyshevT, x))
 (-1.0, 1.0)
-
source
Base.chop — Function
chop(::AbstractPolynomial{T};
-    rtol::Real = Base.rtoldefault(real(T)), atol::Real = 0))

Removes any leading coefficients that are approximately 0 (using rtol and atol with norm(p)). Returns a polynomial whose degree will guaranteed to be equal to or less than the given polynomial's.

source
Polynomials.chop! — Function
chop!(::AbstractPolynomial{T};
-    rtol::Real = Base.rtoldefault(real(T)), atol::Real = 0))

In-place version of chop

source
Base.truncate — Function
truncate(::AbstractPolynomial{T};
-    rtol::Real = Base.rtoldefault(real(T)), atol::Real = 0)

Rounds off coefficients close to zero, as determined by rtol and atol, and then chops any leading zeros. Returns a new polynomial.

source
Base.iszero — Function
iszero(p::AbstractPolynomial)

Is this a $0$ polynomial.

For most types, the $0$ polynomial is one with no coefficients (coefficient vector T[]), though some types have the possibility of trailing zeros. The degree of a zero polynomial is conventionally $-1$, though this is not the convention for Laurent polynomials.

source
Base.isreal — Function
isreal(p::AbstractPolynomial)

Determine whether a polynomial is a real polynomial, i.e., having only real numbers as coefficients.

See also: real

source
Base.real — Function
real(p::AbstractPolynomial)

Construct a real polynomial from the real parts of the coefficients of p.

See also: isreal

Note

This could cause losing terms in p. This method is usually called on polynomials like p = Polynomial([1, 2 + 0im, 3.0, 4.0 + 0.0im]) where you want to chop the imaginary parts of the coefficients of p.

source
Polynomials.isintegral — Function
isintegral(p::AbstractPolynomial)

Determine whether a polynomial is an integer polynomial, i.e., having only integers as coefficients.

source
Polynomials.ismonic — Function
ismonic(p::AbstractPolynomial)

Determine whether a polynomial is a monic polynomial, i.e., its leading coefficient is one.

source

Iteration

For the Polynomial type, a natural mapping between the polynomial $a_0 + a_1 x + a_2 x^2 + \cdots + a_n x^n$ with the coefficients $(a_0, a_1, \dots, a_n))$ leads to the view point of a polynomial being a $0$-based vector. Similarly, when the basis terms are not the standard basis. The coeffs method returns these coefficients in an iterable (a vector or tuple). For Laurent type polynomials, the coefficients between firstindex(p) and lastindex(p) are returned.

More generally, pairs(p) returns values i => aᵢ where the polynomial has terms $a_i T_i$ for the basis $T_i$. (For sparse polynomials these need not be in order and only terms where $a_i \ne 0$ are given.) The keys and values methods iterate over i and aᵢ.

The firstindex method refers to the lowest stored basis index, which due to offsets need not be 0. It will be no smaller than Polynomials.minimumexponent, which is the smalled allowed index for the polynomial type. The lastindex method refers to the last basis index. If the type allows trailing zeros (like ImmutablePolynomial) this will differ from the value returned by degree.

The getindex(p,i) method returns p_i or zero when out of bounds (if the element type of the polynomial has zero(T) defined). For mutable polynomials, the setindex!(p, val, i) method sets p[i] to val. This may extend the underlying storage container for some polynomial types. For ImmutablePolynomial the @set! macro from Setfield can be used with the typical setindex! notation.

The map(fn, p) method maps fn over the coefficients and returns a polynomial with the same polynomial type as p.

Base.chop — Function
chop(::AbstractPolynomial{T};
+    rtol::Real = Base.rtoldefault(real(T)), atol::Real = 0))

Removes any leading coefficients that are approximately 0 (using rtol and atol with norm(p)). Returns a polynomial whose degree will guaranteed to be equal to or less than the given polynomial's.

source
Polynomials.chop! — Function
chop!(::AbstractPolynomial{T};
+    rtol::Real = Base.rtoldefault(real(T)), atol::Real = 0))

In-place version of chop

source
Base.truncate — Function
truncate(::AbstractPolynomial{T};
+    rtol::Real = Base.rtoldefault(real(T)), atol::Real = 0)

Rounds off coefficients close to zero, as determined by rtol and atol, and then chops any leading zeros. Returns a new polynomial.

source
Base.iszero — Function
iszero(p::AbstractPolynomial)

Is this a $0$ polynomial.

For most types, the $0$ polynomial is one with no coefficients (coefficient vector T[]), though some types have the possibility of trailing zeros. The degree of a zero polynomial is conventionally $-1$, though this is not the convention for Laurent polynomials.

source
Base.isreal — Function
isreal(p::AbstractPolynomial)

Determine whether a polynomial is a real polynomial, i.e., having only real numbers as coefficients.

See also: real

source
Base.real — Function
real(p::AbstractPolynomial)

Construct a real polynomial from the real parts of the coefficients of p.

See also: isreal

Note

This could cause losing terms in p. This method is usually called on polynomials like p = Polynomial([1, 2 + 0im, 3.0, 4.0 + 0.0im]) where you want to chop the imaginary parts of the coefficients of p.

source
Polynomials.isintegral — Function
isintegral(p::AbstractPolynomial)

Determine whether a polynomial is an integer polynomial, i.e., having only integers as coefficients.

source
Polynomials.ismonic — Function
ismonic(p::AbstractPolynomial)

Determine whether a polynomial is a monic polynomial, i.e., its leading coefficient is one.

source

Iteration

For the Polynomial type, a natural mapping between the polynomial $a_0 + a_1 x + a_2 x^2 + \cdots + a_n x^n$ with the coefficients $(a_0, a_1, \dots, a_n))$ leads to the view point of a polynomial being a $0$-based vector. Similarly, when the basis terms are not the standard basis. The coeffs method returns these coefficients in an iterable (a vector or tuple). For Laurent type polynomials, the coefficients between firstindex(p) and lastindex(p) are returned.

More generally, pairs(p) returns values i => aᵢ where the polynomial has terms $a_i T_i$ for the basis $T_i$. (For sparse polynomials these need not be in order and only terms where $a_i \ne 0$ are given.) The keys and values methods iterate over i and aᵢ.

The firstindex method refers to the lowest stored basis index, which due to offsets need not be 0. It will be no smaller than Polynomials.minimumexponent, which is the smalled allowed index for the polynomial type. The lastindex method refers to the last basis index. If the type allows trailing zeros (like ImmutablePolynomial) this will differ from the value returned by degree.

The getindex(p,i) method returns p_i or zero when out of bounds (if the element type of the polynomial has zero(T) defined). For mutable polynomials, the setindex!(p, val, i) method sets p[i] to val. This may extend the underlying storage container for some polynomial types. For ImmutablePolynomial the @set! macro from Setfield can be used with the typical setindex! notation.

The map(fn, p) method maps fn over the coefficients and returns a polynomial with the same polynomial type as p.

Polynomials.coeffs — Function
coeffs(::AbstractPolynomial)
 coeffs(::AbstractDenseUnivariatePolynomial)
-coeffs(::AbstractLaurentUnivariatePolynomial)

For a dense, univariate polynomial return the coefficients $(a_0, a_1, \dots, a_n)$ as an interable. This may be a vector or tuple, and may alias the polynomials coefficients.

For a Laurent type polynomial (e.g. LaurentPolynomial, SparsePolynomial) return the coefficients $(a_i, a_{i+1}, \dots, a_j)$ where $i$ is found from firstindex(p) and $j$ from lastindex(p).

For LaurentPolynomial and SparsePolynomial, the pairs iterator is more generically useful, as it iterates over $(i, p_i)$ possibly skipping the terms where $p_i = 0$.

Defaults to p.coeffs.

source
Base.pairs — Function
pairs(p::AbstractPolynomial)

Iterator over $(i, p_i)$ for each basis element, $\beta_i$, represented by the coefficients.

source
Base.values — Function
values(p::AbstractPolynomial)

Iterator over $p_i$s for each basis element, $\beta_i$, represented by the coefficients.

source
Base.keys — Function
keys(p::AbstractPolynomial)

Iterator over $i$s for each basis element, $\beta_i$, represented by the coefficients.

source
Base.firstindex — Function
firstindex(p::AbstractPolynomial)

The index of the smallest basis element, $eta_i$, represented by the coefficients. This is $0$ for a zero polynomial.

source
Base.lastindex — Function
lastindex(p::AbstractPolynomial)

The index of the largest basis element, $eta_i$, represented by the coefficients. May be $-1$ or $0$ for the zero polynomial, depending on the storage type.

source
Base.eachindex — Function
eachindex(p::AbstractPolynomial)

Iterator over all indices of the represented basis elements

source
Base.map — Function
map(fn, p::AbstractPolynomial, args...)

Transform coefficients of p by applying a function (or other callables) fn to each of them.

You can implement real, etc., to a Polynomial by using map. The type of p may narrow using this function.

source

Mathematical Functions

Base.zero — Function
zero(::Type{<:AbstractPolynomial})
-zero(::AbstractPolynomial)

Returns a representation of 0 as the given polynomial.

source
Base.one — Function
one(::Type{<:AbstractPolynomial})
-one(::AbstractPolynomial)

Returns a representation of 1 as the given polynomial.

source
Polynomials.variable — Function
variable(var=:x)
+coeffs(::AbstractLaurentUnivariatePolynomial)

For a dense, univariate polynomial return the coefficients $(a_0, a_1, \dots, a_n)$ as an interable. This may be a vector or tuple, and may alias the polynomials coefficients.

For a Laurent type polynomial (e.g. LaurentPolynomial, SparsePolynomial) return the coefficients $(a_i, a_{i+1}, \dots, a_j)$ where $i$ is found from firstindex(p) and $j$ from lastindex(p).

For LaurentPolynomial and SparsePolynomial, the pairs iterator is more generically useful, as it iterates over $(i, p_i)$ possibly skipping the terms where $p_i = 0$.

Defaults to p.coeffs.

source
Base.pairs — Function
pairs(p::AbstractPolynomial)

Iterator over $(i, p_i)$ for each basis element, $\beta_i$, represented by the coefficients.

source
Base.values — Function
values(p::AbstractPolynomial)

Iterator over $p_i$s for each basis element, $\beta_i$, represented by the coefficients.

source
Base.keys — Function
keys(p::AbstractPolynomial)

Iterator over $i$s for each basis element, $\beta_i$, represented by the coefficients.

source
Base.firstindex — Function
firstindex(p::AbstractPolynomial)

The index of the smallest basis element, $eta_i$, represented by the coefficients. This is $0$ for a zero polynomial.

source
Base.lastindex — Function
lastindex(p::AbstractPolynomial)

The index of the largest basis element, $eta_i$, represented by the coefficients. May be $-1$ or $0$ for the zero polynomial, depending on the storage type.

source
Base.eachindex — Function
eachindex(p::AbstractPolynomial)

Iterator over all indices of the represented basis elements

source
Base.map — Function
map(fn, p::AbstractPolynomial, args...)

Transform coefficients of p by applying a function (or other callables) fn to each of them.

You can implement real, etc., to a Polynomial by using map. The type of p may narrow using this function.

source

Mathematical Functions

Base.zero — Function
zero(::Type{<:AbstractPolynomial})
+zero(::AbstractPolynomial)

Returns a representation of 0 as the given polynomial.

source
Base.one — Function
one(::Type{<:AbstractPolynomial})
+one(::AbstractPolynomial)

Returns a representation of 1 as the given polynomial.

source
Polynomials.variable — Function
variable(var=:x)
 variable(::Type{<:AbstractPolynomial}, var=:x)
 variable(p::AbstractPolynomial, var=indeterminate(p))

Return the monomial x in the indicated polynomial basis. If no type is give, will default to Polynomial. Equivalent to P(var).

Examples

julia> using Polynomials
 
@@ -48,31 +48,31 @@
 2-element Vector{Float64}:
  -2.0
   3.0
-
source
Polynomials.basis — Function
basis(p::P, i::Int)
-basis(::Type{<:AbstractPolynomial}, i::Int, var=:x)

Return ith basis element for a given polynomial type, optionally with a specified variable.

source
Polynomials.basis — Function
basis(p::P, i::Int)
+basis(::Type{<:AbstractPolynomial}, i::Int, var=:x)

Return ith basis element for a given polynomial type, optionally with a specified variable.

source
Polynomials.fromroots — Function
fromroots(::AbstractVector{<:Number}; var=:x)
 fromroots(::Type{<:AbstractPolynomial}, ::AbstractVector{<:Number}; var=:x)

Construct a polynomial of the given type given the roots. If no type is given, defaults to Polynomial.

Examples

julia> using Polynomials
 
 julia> r = [3, 2]; # (x - 3)(x - 2)
 
 julia> fromroots(r)
-Polynomial(6 - 5*x + x^2)
source
fromroots(::AbstractMatrix{<:Number}; var=:x)
+Polynomial(6 - 5*x + x^2)
source
fromroots(::AbstractMatrix{<:Number}; var=:x)
 fromroots(::Type{<:AbstractPolynomial}, ::AbstractMatrix{<:Number}; var=:x)

Construct a polynomial of the given type using the eigenvalues of the given matrix as the roots. If no type is given, defaults to Polynomial.

Examples

julia> using Polynomials
 
 julia> A = [1 2; 3 4]; # (x - 5.37228)(x + 0.37228)
 
 julia> fromroots(A)
-Polynomial(-1.9999999999999998 - 5.0*x + 1.0*x^2)
source
Base.gcd — Function
gcd(a::AbstractPolynomial, b::AbstractPolynomial; atol::Real=0, rtol::Real=Base.rtoldefault)

Find the greatest common denominator of two polynomials recursively using Euclid's algorithm.

Examples

julia> using Polynomials
+Polynomial(-1.9999999999999998 - 5.0*x + 1.0*x^2)
source
Base.gcd — Function
gcd(a::AbstractPolynomial, b::AbstractPolynomial; atol::Real=0, rtol::Real=Base.rtoldefault)

Find the greatest common denominator of two polynomials recursively using Euclid's algorithm.

Examples

julia> using Polynomials
 
 julia> gcd(fromroots([1, 1, 2]), fromroots([1, 2, 3]))
 Polynomial(4.0 - 6.0*x + 2.0*x^2)
-
source
gcd(p1::StandardBasisPolynomial, p2::StandardBasisPolynomial; method=:euclidean, kwargs...)

Find the greatest common divisor.

By default, uses the Euclidean division algorithm (method=:euclidean), which is susceptible to floating point issues.

Passing method=:noda_sasaki uses scaling to circumvent some of these.

Passing method=:numerical will call the internal method NGCD.ngcd for the numerical gcd. See the docstring of NGCD.ngcd for details.

source
Polynomials.derivative — Function
derivative(::AbstractPolynomial, order::Int = 1)

Returns a polynomial that is the orderth derivative of the given polynomial. order must be non-negative.

source
Polynomials.integrate — Function
integrate(p::AbstractPolynomial)

Return an antiderivative for p

source
integrate(::AbstractPolynomial, C)

Returns the indefinite integral of the polynomial with constant C when expressed in the standard basis.

source
integrate(::AbstractPolynomial, a, b)

Compute the definite integral of the given polynomial from a to b. Will throw an error if either a or b are out of the polynomial's domain.

source
Polynomials.roots — Function
roots(::AbstractPolynomial; kwargs...)

Returns the roots, or zeros, of the given polynomial.

For non-factored, standard basis polynomials the roots are calculated via the eigenvalues of the companion matrix. The kwargs are passed to the LinearAlgebra.eigvals call.

Note

The default roots implementation is for polynomials in the standard basis. The companion matrix approach is reasonably fast and accurate for modest-size polynomials. However, other packages in the Julia ecosystem may be of interest and are mentioned in the documentation.

source
roots(pq::AbstractRationalFunction; kwargs...)

Return the zeros of the rational function (after cancelling commong factors, the zeros are the roots of the numerator.

source
Polynomials.fit — Function
fit(x, y, deg=length(x) - 1; [weights], var=:x)
-fit(::Type{<:AbstractPolynomial}, x, y, deg=length(x)-1; [weights], var=:x)

Fit the given data as a polynomial type with the given degree. Uses linear least squares to minimize the norm ||y - V⋅β||^2, where V is the Vandermonde matrix and β are the coefficients of the polynomial fit.

This will automatically scale your data to the domain of the polynomial type using mapdomain. The default polynomial type is Polynomial.

Weights

Weights may be assigned to the points by specifying a vector or matrix of weights.

When specified as a vector, [w₁,…,wₙ], the weights should be non-negative as the minimization problem is argmin_β Σᵢ wᵢ |yᵢ - Σⱼ Vᵢⱼ βⱼ|² = argmin_β || √(W)⋅(y - V(x)β)||², where, W the diagonal matrix formed from [w₁,…,wₙ], is used for the solution, V being the Vandermonde matrix of x corresponding to the specified degree. This parameterization of the weights is different from that of numpy.polyfit, where the weights would be specified through [ω₁,ω₂,…,ωₙ] = [√w₁, √w₂,…,√wₙ] with the answer solving argminᵦ | (ωᵢ⋅yᵢ- ΣⱼVᵢⱼ(ω⋅x) βⱼ) |^2.

When specified as a matrix, W, the solution is through the normal equations (VᵀWV)β = (Vᵀy), again V being the Vandermonde matrix of x corresponding to the specified degree.

(In statistics, the vector case corresponds to weighted least squares, where weights are typically given by wᵢ = 1/σᵢ², the σᵢ² being the variance of the measurement; the matrix specification follows that of the generalized least squares estimator with W = Σ⁻¹, the inverse of the variance-covariance matrix.)

large degree

For fitting with a large degree, the Vandermonde matrix is exponentially ill-conditioned. The ArnoldiFit type introduces an Arnoldi orthogonalization that fixes this problem.

source
fit(P::Type{<:StandardBasisPolynomial}, x, y, J, [cs::Dict{Int, T}]; weights, var)

Using constrained least squares, fit a polynomial of the type p = ∑_{i ∈ J} aᵢ xⁱ + ∑ cⱼxʲ where cⱼ are fixed non-zero constants

  • J: a collection of degrees to find coefficients for
  • cs: If given, a Dict of key/values, i => cᵢ, which indicate the degree and value of the fixed non-zero constants.

The degrees in cs and those in J should not intersect.

Example

x = range(0, pi/2, 10)
+
source
gcd(p1::StandardBasisPolynomial, p2::StandardBasisPolynomial; method=:euclidean, kwargs...)

Find the greatest common divisor.

By default, uses the Euclidean division algorithm (method=:euclidean), which is susceptible to floating point issues.

Passing method=:noda_sasaki uses scaling to circumvent some of these.

Passing method=:numerical will call the internal method NGCD.ngcd for the numerical gcd. See the docstring of NGCD.ngcd for details.

source
Polynomials.derivative — Function
derivative(::AbstractPolynomial, order::Int = 1)

Returns a polynomial that is the orderth derivative of the given polynomial. order must be non-negative.

source
Polynomials.integrate — Function
integrate(p::AbstractPolynomial)

Return an antiderivative for p

source
integrate(::AbstractPolynomial, C)

Returns the indefinite integral of the polynomial with constant C when expressed in the standard basis.

source
integrate(::AbstractPolynomial, a, b)

Compute the definite integral of the given polynomial from a to b. Will throw an error if either a or b are out of the polynomial's domain.

source
Polynomials.roots — Function
roots(::AbstractPolynomial; kwargs...)

Returns the roots, or zeros, of the given polynomial.

For non-factored, standard basis polynomials the roots are calculated via the eigenvalues of the companion matrix. The kwargs are passed to the LinearAlgebra.eigvals call.

Note

The default roots implementation is for polynomials in the standard basis. The companion matrix approach is reasonably fast and accurate for modest-size polynomials. However, other packages in the Julia ecosystem may be of interest and are mentioned in the documentation.

source
roots(pq::AbstractRationalFunction; kwargs...)

Return the zeros of the rational function (after cancelling commong factors, the zeros are the roots of the numerator.

source
Polynomials.fit — Function
fit(x, y, deg=length(x) - 1; [weights], var=:x)
+fit(::Type{<:AbstractPolynomial}, x, y, deg=length(x)-1; [weights], var=:x)

Fit the given data as a polynomial type with the given degree. Uses linear least squares to minimize the norm ||y - V⋅β||^2, where V is the Vandermonde matrix and β are the coefficients of the polynomial fit.

This will automatically scale your data to the domain of the polynomial type using mapdomain. The default polynomial type is Polynomial.

Weights

Weights may be assigned to the points by specifying a vector or matrix of weights.

When specified as a vector, [w₁,…,wₙ], the weights should be non-negative as the minimization problem is argmin_β Σᵢ wᵢ |yᵢ - Σⱼ Vᵢⱼ βⱼ|² = argmin_β || √(W)⋅(y - V(x)β)||², where, W the diagonal matrix formed from [w₁,…,wₙ], is used for the solution, V being the Vandermonde matrix of x corresponding to the specified degree. This parameterization of the weights is different from that of numpy.polyfit, where the weights would be specified through [ω₁,ω₂,…,ωₙ] = [√w₁, √w₂,…,√wₙ] with the answer solving argminᵦ | (ωᵢ⋅yᵢ- ΣⱼVᵢⱼ(ω⋅x) βⱼ) |^2.

When specified as a matrix, W, the solution is through the normal equations (VᵀWV)β = (Vᵀy), again V being the Vandermonde matrix of x corresponding to the specified degree.

(In statistics, the vector case corresponds to weighted least squares, where weights are typically given by wᵢ = 1/σᵢ², the σᵢ² being the variance of the measurement; the matrix specification follows that of the generalized least squares estimator with W = Σ⁻¹, the inverse of the variance-covariance matrix.)

large degree

For fitting with a large degree, the Vandermonde matrix is exponentially ill-conditioned. The ArnoldiFit type introduces an Arnoldi orthogonalization that fixes this problem.

source
fit(P::Type{<:StandardBasisPolynomial}, x, y, J, [cs::Dict{Int, T}]; weights, var)

Using constrained least squares, fit a polynomial of the type p = ∑_{i ∈ J} aᵢ xⁱ + ∑ cⱼxʲ where cⱼ are fixed non-zero constants

  • J: a collection of degrees to find coefficients for
  • cs: If given, a Dict of key/values, i => cᵢ, which indicate the degree and value of the fixed non-zero constants.

The degrees in cs and those in J should not intersect.

Example

x = range(0, pi/2, 10)
 y = sin.(x)
 P = Polynomial
 p0 = fit(P, x, y, 5)
 p1 = fit(P, x, y, 1:2:5)
 p2 = fit(P, x, y, 3:2:5, Dict(1 => 1))
-[norm(p.(x) - y) for p ∈ (p0, p1, p2)] # 1.7e-5, 0.00016, 0.000248
source
fit(::Type{RationalFunction}, xs::AbstractVector{S}, ys::AbstractVector{T}, m, n; var=:x)

Fit a rational function of the form pq = (a₀ + a₁x¹ + … + aₘxᵐ) / (1 + b₁x¹ + … + bₙxⁿ) to the data (x,y).

Note

This uses a simple implementation of the Gauss-Newton method to solve the non-linear least squares problem: minᵦ Σ(yᵢ - pq(xᵢ,β)², where β=(a₀,a₁,…,aₘ,b₁,…,bₙ).

A more rapidly convergent method is used in the LsqFit.jl package, and if performance is important, re-expressing the problem for use with that package is suggested.

Further, if an accurate rational function fit of adaptive degrees is of interest, the BaryRational.jl package provides an implementation of the AAA algorithm ("which offers speed, flexibility, and robustness we have not seen in other algorithms" Nakatsukasa, Sète, Trefethen) and one using Floater-Hormann weights Floater, Hormann ("that have no real poles and arbitrarily high approximation orders on any real interval, regardless of the distribution of the points")

The RationalApproximations package also has implementations of the AAA algorithm.

A python library, polyrat, has implementations of other algorithms.

Example

julia> x = variable(Polynomial{Float64})
+[norm(p.(x) - y) for p ∈ (p0, p1, p2)] # 1.7e-5, 0.00016, 0.000248
source
fit(::Type{RationalFunction}, xs::AbstractVector{S}, ys::AbstractVector{T}, m, n; var=:x)

Fit a rational function of the form pq = (a₀ + a₁x¹ + … + aₘxᵐ) / (1 + b₁x¹ + … + bₙxⁿ) to the data (x,y).

Note

This uses a simple implementation of the Gauss-Newton method to solve the non-linear least squares problem: minᵦ Σ(yᵢ - pq(xᵢ,β)², where β=(a₀,a₁,…,aₘ,b₁,…,bₙ).

A more rapidly convergent method is used in the LsqFit.jl package, and if performance is important, re-expressing the problem for use with that package is suggested.

Further, if an accurate rational function fit of adaptive degrees is of interest, the BaryRational.jl package provides an implementation of the AAA algorithm ("which offers speed, flexibility, and robustness we have not seen in other algorithms" Nakatsukasa, Sète, Trefethen) and one using Floater-Hormann weights Floater, Hormann ("that have no real poles and arbitrarily high approximation orders on any real interval, regardless of the distribution of the points")

The RationalApproximations package also has implementations of the AAA algorithm.

A python library, polyrat, has implementations of other algorithms.

Example

julia> x = variable(Polynomial{Float64})
 Polynomial(1.0*x)
 
 julia> pq = (1+x)//(1-x)
@@ -97,7 +97,7 @@
 4.440892098500626e-16
 
 julia> u(variable(pq)) # to see which polynomial is used
-(2.68328 + 0.447214*x - 1.78885*x^2 + 0.447214*x^3) // (2.68328 - 4.91935*x + 2.68328*x^2 - 0.447214*x^3)
source
fit(::Type{RationalFunction}, r::Polynomial, m, n; var=:x)

Fit a Pade approximant (cf docstring for Polynomials.pade_fit) to r.

Examples:

julia> using Polynomials, PolynomialRatios
+(2.68328 + 0.447214*x - 1.78885*x^2 + 0.447214*x^3) // (2.68328 - 4.91935*x + 2.68328*x^2 - 0.447214*x^3)
source
fit(::Type{RationalFunction}, r::Polynomial, m, n; var=:x)

Fit a Pade approximant (cf docstring for Polynomials.pade_fit) to r.

Examples:

julia> using Polynomials, PolynomialRatios
 
 julia> x = variable()
 Polynomial(x)
@@ -115,7 +115,7 @@
 0.001278729299871717
 
 julia> maximum(abs, exp(x) - fit(RationalFunction, ex, 2,2)(x) for x ∈ 0:.05:0.5)
-7.262205147950951e-5
source
Polynomials.vander — Function
vander(::Type{AbstractPolynomial}, x::AbstractVector, deg::Integer)

Calculate the pseudo-Vandermonde matrix of the given polynomial type with the given degree.

References

Vandermonde Matrix

source

Plotting

Polynomials can be plotted directly using Plots.jl or Makie.jl.

plot(::AbstractPolynomial; kwds...)

will automatically determine a range based on the critical points (roots, extrema and points of inflection).

plot(::AbstractPolynomial, a, b; kwds...)

will plot the polynomial within the range [a, b].

using Plots, Polynomials
+7.262205147950951e-5
source
Polynomials.vander — Function
vander(::Type{AbstractPolynomial}, x::AbstractVector, deg::Integer)

Calculate the pseudo-Vandermonde matrix of the given polynomial type with the given degree.

References

Vandermonde Matrix

source

Plotting

Polynomials can be plotted directly using Plots.jl or Makie.jl.

plot(::AbstractPolynomial; kwds...)

will automatically determine a range based on the critical points (roots, extrema and points of inflection).

plot(::AbstractPolynomial, a, b; kwds...)

will plot the polynomial within the range [a, b].

using Plots, Polynomials
 # T1, T2, T3, and T4:
 chebs = [
   ChebyshevT([0, 1]),
@@ -128,4 +128,4 @@
 p = plot(legend=false, label="")
 for (cheb, col) in zip(chebs, colors)
   plot!(cheb, c=col, lw=5)
-end

+end