From 8b8055dcc9fea4ef2af0638fa574c6209f0a3716 Mon Sep 17 00:00:00 2001 From: abhro <5664668+abhro@users.noreply.github.com> Date: Fri, 7 Jun 2024 16:47:39 -0400 Subject: [PATCH] Doctests and documentation patch (#571) * Add lowest_terms doc for site Polynomials.RationalFunction makes a reference to lowest_terms * Add more doctests and syntax highlighting tags * Fix latex syntax in docstrings * Remove names from doctest blocks The named doctest blocks don't share scope and don't need to have the same environment * Fix spacing * Use markdown lists in docstrings * Fix more spacing * Use identity check for nothing-ness * Pull list text into list --- docs/Project.toml | 1 + docs/src/extending.md | 28 +- docs/src/index.md | 353 +++++++++--------- docs/src/polynomials/chebyshev.md | 1 - docs/src/polynomials/polynomial.md | 3 +- src/common.jl | 57 ++- src/legacy/Poly.jl | 4 +- src/legacy/pade.jl | 5 +- .../mutable-dense-laurent-polynomial.jl | 2 +- .../mutable-dense-polynomial.jl | 2 +- .../standard-basis/immutable-polynomial.jl | 13 +- .../standard-basis/laurent-polynomial.jl | 12 +- .../standard-basis/sparse-polynomial.jl | 4 +- .../standard-basis/standard-basis.jl | 32 +- src/rational-functions/fit.jl | 10 +- src/rational-functions/plot-recipes.jl | 3 +- src/rational-functions/rational-function.jl | 2 +- .../rational-transfer-function.jl | 24 +- src/show.jl | 7 - 19 files changed, 276 insertions(+), 287 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index 49ddb3c8..9e46871e 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -2,6 +2,7 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" DualNumbers = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74" GR_jll = "d2c73de3-f751-5644-a686-071e5b155ba9" +GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" diff --git a/docs/src/extending.md b/docs/src/extending.md index 748f38ad..4453ae92 100644 --- a/docs/src/extending.md +++ b/docs/src/extending.md @@ -283,24 +283,24 @@ We need to add in the vector-space operations: ```jldoctest new_container_type julia> function Base.:+(p::AliasPolynomialType{B,T,X}, q::AliasPolynomialType{B,S,X}) where {B,S,T,X} - R = promote_type(T,S) - n = maximum(degree, (p,q)) - cs = [p[i] + q[i] for i in 0:n] - AliasPolynomialType{B,R,X}(Val(false), cs) # save a copy - end + R = promote_type(T,S) + n = maximum(degree, (p,q)) + cs = [p[i] + q[i] for i in 0:n] + AliasPolynomialType{B,R,X}(Val(false), cs) # save a copy + end julia> function Base.:-(p::AliasPolynomialType{B,T,X}, q::AliasPolynomialType{B,S,X}) where {B,S,T,X} - R = promote_type(T,S) - n = maximum(degree, (p,q)) - cs = [p[i] - q[i] for i in 0:n] - AliasPolynomialType{B,R,X}(Val(false), cs) - end + R = promote_type(T,S) + n = maximum(degree, (p,q)) + cs = [p[i] - q[i] for i in 0:n] + AliasPolynomialType{B,R,X}(Val(false), cs) + end julia> function Base.map(fn, p::P) where {B,T,X,P<:AliasPolynomialType{B,T,X}} - cs = map(fn, p.coeffs) - R = eltype(cs) - AliasPolynomialType{B,R,X}(Val(false), cs) - end + cs = map(fn, p.coeffs) + R = eltype(cs) + AliasPolynomialType{B,R,X}(Val(false), cs) + end ``` diff --git a/docs/src/index.md b/docs/src/index.md index 005a8a43..ab47ede5 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -50,7 +50,6 @@ Polynomial(1 - x^2) julia> p(1) 0 - ``` The `Polynomial` constructor stores all coefficients using the standard basis with a vector. Other types (e.g. `ImmutablePolynomial`, `SparsePolynomial`, or `FactoredPolynomial`) use different back-end containers which may have advantage for some uses. @@ -117,7 +116,7 @@ Polynomial(2 + 2*x + 3*x^2) Arithmetic of different polynomial types is supported through promotion to a common type, which is typically the `Polynomial` type, but may be the `LaurentPolynomial` type when negative powers of the indeterminate are possible: -```jldoctext +```jldoctest julia> p, q = ImmutablePolynomial([1,2,3]), Polynomial([3,2,1]) (ImmutablePolynomial(1 + 2*x + 3*x^2), Polynomial(3 + 2*x + x^2)) @@ -184,7 +183,7 @@ a polynomial. This is an `𝑶(n^3)` operation. For polynomials with `BigFloat` coefficients, the `GenericLinearAlgebra` package can be seamlessly used: -``` +```jldoctest julia> p = fromroots(Polynomial{BigFloat}, [1,2,3]) Polynomial(-6.0 + 11.0*x - 6.0*x^2 + 1.0*x^3) @@ -209,20 +208,20 @@ julia> roots(p) to univariate polynomials that is more performant than `roots`. It is based on an algorithm of Skowron and Gould. -``` -julia> import PolynomialRoots # import as `roots` conflicts + ```julia-repl + julia> import PolynomialRoots # import as `roots` conflicts -julia> p = fromroots(Polynomial, [1,2,3]) -Polynomial(-6 + 11*x - 6*x^2 + x^3) + julia> p = fromroots(Polynomial, [1,2,3]) + Polynomial(-6 + 11*x - 6*x^2 + x^3) -julia> PolynomialRoots.roots(coeffs(p)) -3-element Vector{ComplexF64}: - 3.000000000000001 - 0.0im - 1.9999999999999993 + 0.0im - 1.0000000000000002 + 0.0im -``` + julia> PolynomialRoots.roots(coeffs(p)) + 3-element Vector{ComplexF64}: + 3.000000000000001 - 0.0im + 1.9999999999999993 + 0.0im + 1.0000000000000002 + 0.0im + ``` -The roots are always returned as complex numbers. + The roots are always returned as complex numbers. * The @@ -234,62 +233,62 @@ The roots are always returned as complex numbers. package implements the algorithm in Julia, allowing the use of other number types. -``` -julia> using AMRVW + ```julia-repl + julia> using AMRVW -julia> AMRVW.roots(float.(coeffs(p))) -3-element Vector{ComplexF64}: - 0.9999999999999997 + 0.0im - 2.0000000000000036 + 0.0im - 2.9999999999999964 + 0.0im -``` + julia> AMRVW.roots(float.(coeffs(p))) + 3-element Vector{ComplexF64}: + 0.9999999999999997 + 0.0im + 2.0000000000000036 + 0.0im + 2.9999999999999964 + 0.0im + ``` -The roots are returned as complex numbers. + The roots are returned as complex numbers. -Both `PolynomialRoots` and `AMRVW` are generic and work with -`BigFloat` coefficients, for example. + Both `PolynomialRoots` and `AMRVW` are generic and work with + `BigFloat` coefficients, for example. -The `AMRVW` package works with much larger polynomials than either -`roots` or `PolynomialRoots.roots`. For example, the roots of this 1000 -degree random polynomial are quickly and accurately solved for: + The `AMRVW` package works with much larger polynomials than either `roots` + or `PolynomialRoots.roots`. For example, the roots of this 1000 degree + random polynomial are quickly and accurately solved for: -``` -julia> filter(isreal, AMRVW.roots(rand(1001) .- 1/2)) -2-element Vector{ComplexF64}: - 0.993739974989572 + 0.0im - 1.0014677846996498 + 0.0im -``` + ```julia-repl + julia> filter(isreal, AMRVW.roots(rand(1001) .- 1/2)) + 2-element Vector{ComplexF64}: + 0.993739974989572 + 0.0im + 1.0014677846996498 + 0.0im + ``` * The [Hecke](https://github.com/thofma/Hecke.jl/tree/master/src) package has a `roots` function. The `Hecke` package utilizes the `Arb` library for performant, high-precision numbers: -``` -julia> import Hecke # import as `roots` conflicts + ```julia-repl + julia> import Hecke # import as `roots` conflicts -julia> Qx, x = Hecke.PolynomialRing(Hecke.QQ) -(Univariate Polynomial Ring in x over Rational Field, x) + julia> Qx, x = Hecke.PolynomialRing(Hecke.QQ) + (Univariate Polynomial Ring in x over Rational Field, x) -julia> q = (x-1)*(x-2)*(x-3) -x^3 - 6*x^2 + 11*x - 6 + julia> q = (x-1)*(x-2)*(x-3) + x^3 - 6*x^2 + 11*x - 6 -julia> Hecke.roots(q) -3-element Vector{Nemo.fmpq}: - 2 - 1 - 3 -``` + julia> Hecke.roots(q) + 3-element Vector{Nemo.fmpq}: + 2 + 1 + 3 + ``` -This next polynomial has 3 real roots, 2 of which are in a cluster; `Hecke` quickly identifies them: + This next polynomial has 3 real roots, 2 of which are in a cluster; `Hecke` quickly identifies them: -``` -julia> p = -1 + 254*x - 16129*x^2 + 1*x^17 -x^17 - 16129*x^2 + 254*x - 1 + ```julia-repl + julia> p = -1 + 254*x - 16129*x^2 + 1*x^17 + x^17 - 16129*x^2 + 254*x - 1 -julia> filter(isreal, Hecke._roots(p, 200)) # `_roots` not `roots` -3-element Vector{Nemo.acb}: - [0.007874015748031496052667730054749907629383970426203662570129818116411192289734968717460531379762086419 +/- 3.10e-103] - [0.0078740157480314960733165219137540296086246589982151627453855179522742093785877068332663198273096875302 +/- 9.31e-104] - [1.9066348541790688341521872066398429982632947292434604847312536201982593209326201234353468172497707769372732739429697289 +/- 7.14e-119] -``` + julia> filter(isreal, Hecke._roots(p, 200)) # `_roots` not `roots` + 3-element Vector{Nemo.acb}: + [0.007874015748031496052667730054749907629383970426203662570129818116411192289734968717460531379762086419 +/- 3.10e-103] + [0.0078740157480314960733165219137540296086246589982151627453855179522742093785877068332663198273096875302 +/- 9.31e-104] + [1.9066348541790688341521872066398429982632947292434604847312536201982593209326201234353468172497707769372732739429697289 +/- 7.14e-119] + ``` ---- @@ -300,75 +299,73 @@ To find just the real roots of a polynomial with real coefficients there are a f identifies real zeros of univariate functions and can be used to find isolating intervals for the real roots. For example, -``` -julia> using Polynomials, IntervalArithmetic - -julia> import IntervalRootFinding # its `roots` method conflicts with `roots` + ```julia-repl + julia> using Polynomials, IntervalArithmetic -julia> p = fromroots(Polynomial, [1,2,3]) -Polynomial(-6 + 11*x - 6*x^2 + x^3) - -julia> IntervalRootFinding.roots(x -> p(x), 0..10) -3-element Vector{IntervalRootFinding.Root{Interval{Float64}}}: - Root([0.999999, 1.00001], :unique) - Root([1.99999, 2.00001], :unique) - Root([2.99999, 3.00001], :unique) -``` + julia> import IntervalRootFinding # its `roots` method conflicts with `roots` + julia> p = fromroots(Polynomial, [1,2,3]) + Polynomial(-6 + 11*x - 6*x^2 + x^3) + julia> IntervalRootFinding.roots(x -> p(x), 0..10) + 3-element Vector{IntervalRootFinding.Root{Interval{Float64}}}: + Root([0.999999, 1.00001], :unique) + Root([1.99999, 2.00001], :unique) + Root([2.99999, 3.00001], :unique) + ``` -The output is a set of intervals. Those flagged with `:unique` are guaranteed to contain a unique root. + The output is a set of intervals. Those flagged with `:unique` are guaranteed to contain a unique root. * The `RealPolynomialRoots` package provides a function `ANewDsc` to find isolating intervals for the roots of a square-free polynomial, specified through its coefficients: -``` -julia> using RealPolynomialRoots - -julia> st = ANewDsc(coeffs(p)) -There were 3 isolating intervals found: -[2.62…, 3.62…]₂₅₆ -[1.5…, 2.62…]₂₅₆ -[-0.50…, 1.5…]₂₅₆ -``` - -These isolating intervals can be refined to find numeric estimates for the roots over `BigFloat` values. - -``` -julia> refine_roots(st) -3-element Vector{BigFloat}: - 2.99999999999999999999... - 2.00000000000000000000... - 1.00000000000000000000... -``` - -This specialized algorithm can identify very nearby roots. For example, returning to this Mignotte-type polynomial: - -``` -julia> p = SparsePolynomial(Dict(0=>-1, 1=>254, 2=>-16129, 17=>1)) -SparsePolynomial(-1 + 254*x - 16129*x^2 + x^17) - -julia> ANewDsc(coeffs(p)) -There were 3 isolating intervals found: -[1.5…, 3.5…]₅₃ -[0.0078740157480314960682066…, 0.0078740157480314960873178…]₁₃₉ -[0.0078740157480314960492543…, 0.0078740157480314960682066…]₁₃₉ -``` - -`IntervalRootFinding` has issues disambiguating the clustered roots of this example: - -``` -julia> IntervalRootFinding.roots(x -> p(x), 0..3.5) -7-element Vector{IntervalRootFinding.Root{Interval{Float64}}}: - Root([1.90663, 1.90664], :unique) - Root([0.00787464, 0.00787468], :unknown) - Root([0.00787377, 0.00787387], :unknown) - Root([0.00787405, 0.00787412], :unknown) - Root([0.00787396, 0.00787406], :unknown) - Root([0.00787425, 0.00787431], :unknown) - Root([0.00787394, 0.00787397], :unknown) -``` - -For this example, `filter(isreal, Hecke._roots(p))` also isolates the three real roots, but not quite as quickly. + ```julia-repl + julia> using RealPolynomialRoots + + julia> st = ANewDsc(coeffs(p)) + There were 3 isolating intervals found: + [2.62…, 3.62…]₂₅₆ + [1.5…, 2.62…]₂₅₆ + [-0.50…, 1.5…]₂₅₆ + ``` + + These isolating intervals can be refined to find numeric estimates for the roots over `BigFloat` values. + + ```julia-repl + julia> refine_roots(st) + 3-element Vector{BigFloat}: + 2.99999999999999999999... + 2.00000000000000000000... + 1.00000000000000000000... + ``` + + This specialized algorithm can identify very nearby roots. For example, returning to this Mignotte-type polynomial: + + ```julia-repl + julia> p = SparsePolynomial(Dict(0=>-1, 1=>254, 2=>-16129, 17=>1)) + SparsePolynomial(-1 + 254*x - 16129*x^2 + x^17) + + julia> ANewDsc(coeffs(p)) + There were 3 isolating intervals found: + [1.5…, 3.5…]₅₃ + [0.0078740157480314960682066…, 0.0078740157480314960873178…]₁₃₉ + [0.0078740157480314960492543…, 0.0078740157480314960682066…]₁₃₉ + ``` + + `IntervalRootFinding` has issues disambiguating the clustered roots of this example: + + ```julia-repl + julia> IntervalRootFinding.roots(x -> p(x), 0..3.5) + 7-element Vector{IntervalRootFinding.Root{Interval{Float64}}}: + Root([1.90663, 1.90664], :unique) + Root([0.00787464, 0.00787468], :unknown) + Root([0.00787377, 0.00787387], :unknown) + Root([0.00787405, 0.00787412], :unknown) + Root([0.00787396, 0.00787406], :unknown) + Root([0.00787425, 0.00787431], :unknown) + Root([0.00787394, 0.00787397], :unknown) + ``` + + For this example, `filter(isreal, Hecke._roots(p))` also isolates the three real roots, but not quite as quickly. ---- @@ -378,66 +375,65 @@ square free polynomial. For non-square free polynomials: * The `Polynomials.Multroot.multroot` function is available for finding the roots of a polynomial and their multiplicities. This is based on work of Zeng. -Here we see `IntervalRootFinding.roots` having trouble isolating the roots due to the multiplicities: - -``` -julia> p = fromroots(Polynomial, [1,2,2,3,3]) -Polynomial(-36 + 96*x - 97*x^2 + 47*x^3 - 11*x^4 + x^5) - -julia> IntervalRootFinding.roots(x -> p(x), 0..10) -335-element Vector{IntervalRootFinding.Root{Interval{Float64}}}: - Root([1.99999, 2], :unknown) - Root([1.99999, 2], :unknown) - Root([3, 3.00001], :unknown) - Root([2.99999, 3], :unknown) - ⋮ - Root([2.99999, 3], :unknown) - Root([2, 2.00001], :unknown) -``` - -The `roots` function identifies the roots, but the multiplicities would need identifying: - -``` -julia> roots(p) -5-element Vector{Float64}: - 1.000000000000011 - 1.9999995886034314 - 2.0000004113969276 - 2.9999995304339646 - 3.0000004695656672 -``` - - -Whereas, the roots along with the multiplicity structure are correctly identified with `multroot`: - -``` -julia> Polynomials.Multroot.multroot(p) -(values = [1.0000000000000004, 1.9999999999999984, 3.0000000000000018], multiplicities = [1, 2, 2], κ = 35.11176306900731, ϵ = 0.0) -``` - -The `square_free` function can help: - -``` -julia> q = Polynomials.square_free(p) -ANewDsc(q) -Polynomial(-0.20751433915978448 + 0.38044295512633425*x - 0.20751433915986722*x^2 + 0.03458572319332053*x^3) - -julia> IntervalRootFinding.roots(x -> q(x), 0..10) -3-element Vector{IntervalRootFinding.Root{Interval{Float64}}}: - Root([0.999999, 1.00001], :unique) - Root([1.99999, 2.00001], :unique) - Root([2.99999, 3.00001], :unique) -``` - -Similarly: - -``` -julia> ANewDsc(coeffs(q)) -There were 3 isolating intervals found: -[2.62…, 3.62…]₂₅₆ -[1.5…, 2.62…]₂₅₆ -[-0.50…, 1.5…]₂₅₆ -``` + Here we see `IntervalRootFinding.roots` having trouble isolating the roots due to the multiplicities: + + ```julia-repl + julia> p = fromroots(Polynomial, [1,2,2,3,3]) + Polynomial(-36 + 96*x - 97*x^2 + 47*x^3 - 11*x^4 + x^5) + + julia> IntervalRootFinding.roots(x -> p(x), 0..10) + 335-element Vector{IntervalRootFinding.Root{Interval{Float64}}}: + Root([1.99999, 2], :unknown) + Root([1.99999, 2], :unknown) + Root([3, 3.00001], :unknown) + Root([2.99999, 3], :unknown) + ⋮ + Root([2.99999, 3], :unknown) + Root([2, 2.00001], :unknown) + ``` + + The `roots` function identifies the roots, but the multiplicities would need identifying: + + ```julia-repl + julia> roots(p) + 5-element Vector{Float64}: + 1.000000000000011 + 1.9999995886034314 + 2.0000004113969276 + 2.9999995304339646 + 3.0000004695656672 + ``` + + Whereas, the roots along with the multiplicity structure are correctly identified with `multroot`: + + ```julia-repl + julia> Polynomials.Multroot.multroot(p) + (values = [1.0000000000000004, 1.9999999999999984, 3.0000000000000018], multiplicities = [1, 2, 2], κ = 35.11176306900731, ϵ = 0.0) + ``` + + The `square_free` function can help: + + ```julia-repl + julia> q = Polynomials.square_free(p) + ANewDsc(q) + Polynomial(-0.20751433915978448 + 0.38044295512633425*x - 0.20751433915986722*x^2 + 0.03458572319332053*x^3) + + julia> IntervalRootFinding.roots(x -> q(x), 0..10) + 3-element Vector{IntervalRootFinding.Root{Interval{Float64}}}: + Root([0.999999, 1.00001], :unique) + Root([1.99999, 2.00001], :unique) + Root([2.99999, 3.00001], :unique) + ``` + + Similarly: + + ```julia-repl + julia> ANewDsc(coeffs(q)) + There were 3 isolating intervals found: + [2.62…, 3.62…]₂₅₆ + [1.5…, 2.62…]₂₅₆ + [-0.50…, 1.5…]₂₅₆ + ``` ## Fitting a polynomial to arbitrary data @@ -817,7 +813,6 @@ julia> for (λ, rs) ∈ r # reconstruct p/q from output of `residues` end end - julia> d ((x - 4.0) * (x - 1.0000000000000002)) // ((x - 5.0) * (x - 2.0)) ``` diff --git a/docs/src/polynomials/chebyshev.md b/docs/src/polynomials/chebyshev.md index 3ef0187c..fc6f23cf 100644 --- a/docs/src/polynomials/chebyshev.md +++ b/docs/src/polynomials/chebyshev.md @@ -50,7 +50,6 @@ For example, the basis polynomial ``T_4`` can be represented with `ChebyshevT([0 julia> c = ChebyshevT([1, 0, 3, 4]) ChebyshevT(1⋅T_0(x) + 3⋅T_2(x) + 4⋅T_3(x)) - julia> p = convert(Polynomial, c) Polynomial(-2.0 - 12.0*x + 6.0*x^2 + 16.0*x^3) diff --git a/docs/src/polynomials/polynomial.md b/docs/src/polynomials/polynomial.md index 5442a5f7..d0dcfd47 100644 --- a/docs/src/polynomials/polynomial.md +++ b/docs/src/polynomials/polynomial.md @@ -41,6 +41,5 @@ FactoredPolynomial ## Rational Function ```@docs RationalFunction +lowest_terms ``` - - diff --git a/src/common.jl b/src/common.jl index 6d6468e4..b59990ad 100644 --- a/src/common.jl +++ b/src/common.jl @@ -24,9 +24,7 @@ export fromroots, Construct a polynomial of the given type given the roots. If no type is given, defaults to `Polynomial`. # Examples -```jldoctest common -julia> using Polynomials - +```jldoctest julia> r = [3, 2]; # (x - 3)(x - 2) julia> fromroots(r) @@ -49,9 +47,7 @@ fromroots(r::AbstractVector{<:Number}; var::SymbolLike = :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 -```jldoctest common -julia> using Polynomials - +```jldoctest julia> A = [1 2; 3 4]; # (x - 5.37228)(x + 0.37228) julia> fromroots(A) @@ -105,7 +101,9 @@ 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. +For fitting with a large degree, the Vandermonde matrix is exponentially +ill-conditioned. The `ArnoldiFit` type introduces an Arnoldi orthogonalization +that fixes this problem. """ @@ -180,7 +178,9 @@ _wlstsq(vand, y, W::AbstractMatrix) = (vand' * W * vand) \ (vand' * W * y) 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. +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 @@ -230,7 +230,8 @@ integrate(P::AbstractPolynomial) = throw(ArgumentError("`integrate` not implemen """ integrate(::AbstractPolynomial, C) -Returns the indefinite integral of the polynomial with constant `C` when expressed in the standard basis. +Returns the indefinite integral of the polynomial with constant `C` when +expressed in the standard basis. """ function integrate(p::P, C) where {P <: AbstractPolynomial} ∫p = integrate(p) @@ -241,7 +242,8 @@ end """ 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. +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. """ function integrate(p::AbstractPolynomial, a, b) P = integrate(p) @@ -251,7 +253,8 @@ end """ derivative(::AbstractPolynomial, order::Int = 1) -Returns a polynomial that is the `order`th derivative of the given polynomial. `order` must be non-negative. +Returns a polynomial that is the `order`th derivative of the given polynomial. +`order` must be non-negative. """ derivative(::AbstractPolynomial, ::Int) @@ -263,7 +266,9 @@ Return the critical points of `p` (real zeros of the derivative) within `I` in s * `p`: a polynomial -* `I`: a specification of a closed or infinite domain, defaulting to `Polynomials.domain(p)`. When specified, the values of `extrema(I)` are used with closed endpoints when finite. +* `I`: a specification of a closed or infinite domain, defaulting to + `Polynomials.domain(p)`. When specified, the values of `extrema(I)` are used + with closed endpoints when finite. * `endpoints::Bool`: if `true`, return the endpoints of `I` along with the critical points @@ -271,7 +276,7 @@ Return the critical points of `p` (real zeros of the derivative) within `I` in s Can be used in conjunction with `findmax`, `findmin`, `argmax`, `argmin`, `extrema`, etc. ## Example -``` +```julia x = variable() p = x^2 - 2 cps = Polynomials.critical_points(p) @@ -713,9 +718,7 @@ domain(::P) where {P <: AbstractPolynomial} = domain(P) Given values of x that are assumed to be unbounded (-∞, ∞), return values rescaled to the domain of the given polynomial. # Examples -```jldoctest common -julia> using Polynomials - +```jldoctest julia> x = -10:10 -10:10 @@ -746,14 +749,14 @@ minimumexponent(::Type{<:AbstractPolynomial}) = 0 """ firstindex(p::AbstractPolynomial) -The index of the smallest basis element, ``\beta_i``, represented by the coefficients. This is ``0`` for +The index of the smallest basis element, ``\\beta_i``, represented by the coefficients. This is ``0`` for a zero polynomial. """ Base.firstindex(p::AbstractPolynomial) = 0 # XXX() is a better default """ lastindex(p::AbstractPolynomial) -The index of the largest basis element, ``\beta_i``, represented by the coefficients. +The index of the largest basis element, ``\\beta_i``, represented by the coefficients. May be ``-1`` or ``0`` for the zero polynomial, depending on the storage type. """ Base.lastindex(p::AbstractPolynomial) = length(p) - 1 + firstindex(p) # not degree, which accounts for any trailing zeros @@ -934,12 +937,10 @@ Base.oneunit(p::P, args...) where {P <: AbstractPolynomial} = one(p, args...) 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`](@ref). Equivalent to `P(var)`. +Return the monomial `x` in the indicated polynomial basis. If no type is give, will default to [`Polynomial`](@ref). Equivalent to `P(var)`. # Examples -```jldoctest common -julia> using Polynomials - +```jldoctest julia> x = variable() Polynomial(x) @@ -950,7 +951,6 @@ julia> roots((x - 3) * (x + 2)) 2-element Vector{Float64}: -2.0 3.0 - ``` """ variable(::Type{P}) where {P <: AbstractPolynomial} = variable(⟒(P){eltype(P), indeterminate(P)}) @@ -1168,9 +1168,9 @@ function ⊕(P::Type{<:AbstractPolynomial}, p1::Dict{Int,T}, p2::Dict{Int,S}) wh # this allocates in the union -# for i in union(eachindex(p1), eachindex(p2)) -# p[i] = p1[i] + p2[i] -# end + #for i in union(eachindex(p1), eachindex(p2)) + # p[i] = p1[i] + p2[i] + #end for (i,pi) ∈ pairs(p1) @inbounds p[i] = pi + get(p2, i, zero(R)) @@ -1201,12 +1201,9 @@ Find the greatest common denominator of two polynomials recursively using # Examples -```jldoctest common -julia> using Polynomials - +```jldoctest julia> gcd(fromroots([1, 1, 2]), fromroots([1, 2, 3])) Polynomial(4.0 - 6.0*x + 2.0*x^2) - ``` """ function Base.gcd(p1::AbstractPolynomial{T}, p2::AbstractPolynomial{S}; kwargs...) where {T,S} diff --git a/src/legacy/Poly.jl b/src/legacy/Poly.jl index 999b4213..4e4d3be2 100644 --- a/src/legacy/Poly.jl +++ b/src/legacy/Poly.jl @@ -216,10 +216,10 @@ end polyint(p::Poly, C = 0) = integrate(p, C) polyint(p::Poly, a, b) = integrate(p, a, b) -polyint(p::AbstractPolynomial, args...) = error("`polyint` is a legacy name for use with `Poly` objects only. Use `integrate(p,...)`.") +polyint(p::AbstractPolynomial, args...) = error("`polyint` is a legacy name for use with `Poly` objects only. Use `integrate(p,...)`.") polyder(p::Poly, ord = 1) = derivative(p, ord) -polyder(p::AbstractPolynomial, args...) = error("`polyder` is a legacy name for use with `Poly` objects only. Use `derivative(p,[order=1])`.") +polyder(p::AbstractPolynomial, args...) = error("`polyder` is a legacy name for use with `Poly` objects only. Use `derivative(p,[order=1])`.") polyfit(x, y, n = length(x) - 1, sym=:x) = fit(Poly, x, y, n; var = sym) polyfit(x, y, sym::Symbol) = fit(Poly, x, y, var = sym) diff --git a/src/legacy/pade.jl b/src/legacy/pade.jl index b68e2b29..a13f4073 100644 --- a/src/legacy/pade.jl +++ b/src/legacy/pade.jl @@ -89,10 +89,9 @@ end Evaluate the Pade approximant at the given point. # Examples -```jldoctest pade +```jldoctest julia> using Polynomials, Polynomials.PolyCompat, SpecialFunctions - julia> p = Polynomial(@.(1 // BigInt(gamma(1:17)))); julia> pade = Pade(p, 8, 8); @@ -108,6 +107,4 @@ true padeval(PQ::Pade, x::Number) = PQ(x) padeval(PQ::Pade, x) = PQ.(x) - - end diff --git a/src/polynomial-container-types/mutable-dense-laurent-polynomial.jl b/src/polynomial-container-types/mutable-dense-laurent-polynomial.jl index 04dfc948..542b56a9 100644 --- a/src/polynomial-container-types/mutable-dense-laurent-polynomial.jl +++ b/src/polynomial-container-types/mutable-dense-laurent-polynomial.jl @@ -21,7 +21,7 @@ struct MutableDenseLaurentPolynomial{B,T,X} <: AbstractLaurentUnivariatePolynomi end i = findlast(!iszero, cs) - if i == nothing + if i === nothing xs = T[] else j = findfirst(!iszero, cs) diff --git a/src/polynomial-container-types/mutable-dense-polynomial.jl b/src/polynomial-container-types/mutable-dense-polynomial.jl index f3311ad1..e9658df1 100644 --- a/src/polynomial-container-types/mutable-dense-polynomial.jl +++ b/src/polynomial-container-types/mutable-dense-polynomial.jl @@ -2,7 +2,7 @@ MutableDensePolynomial{B,T,X} """ -struct MutableDensePolynomial{B,T,X} <: AbstractDenseUnivariatePolynomial{B,T, X} +struct MutableDensePolynomial{B,T,X} <: AbstractDenseUnivariatePolynomial{B,T,X} coeffs::Vector{T} function MutableDensePolynomial{B,T,X}(::Val{true}, cs::AbstractVector{S}, order::Int=0) where {B,T,X,S} if Base.has_offset_axes(cs) diff --git a/src/polynomials/standard-basis/immutable-polynomial.jl b/src/polynomials/standard-basis/immutable-polynomial.jl index 4dbf60b8..abb63f61 100644 --- a/src/polynomials/standard-basis/immutable-polynomial.jl +++ b/src/polynomials/standard-basis/immutable-polynomial.jl @@ -28,14 +28,13 @@ 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. + `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 -```jldoctest immutable_polynomials -julia> using Polynomials - +```jldoctest julia> ImmutablePolynomial((1, 0, 3, 4)) ImmutablePolynomial(1 + 3*x^2 + 4*x^3) @@ -47,7 +46,9 @@ ImmutablePolynomial(1.0) ``` !!! note - This was modeled after [StaticUnivariatePolynomials](https://github.com/tkoolen/StaticUnivariatePolynomials.jl) by `@tkoolen`. + This was modeled after + [StaticUnivariatePolynomials](https://github.com/tkoolen/StaticUnivariatePolynomials.jl) + by `@tkoolen`. """ const ImmutablePolynomial = ImmutableDensePolynomial{StandardBasis} diff --git a/src/polynomials/standard-basis/laurent-polynomial.jl b/src/polynomials/standard-basis/laurent-polynomial.jl index 504ad0fc..62ee0ef5 100644 --- a/src/polynomials/standard-basis/laurent-polynomial.jl +++ b/src/polynomials/standard-basis/laurent-polynomial.jl @@ -16,9 +16,7 @@ Integration will fail if there is a `x⁻¹` term in the polynomial. `LaurentPolynomial` is axis-aware, unlike the other polynomial types in this package. # Examples: -```jldoctest laurent -julia> using Polynomials - +```jldoctest julia> P = LaurentPolynomial; julia> p = P([1,1,1], -1) @@ -185,9 +183,7 @@ Call `p̂ = paraconj(p)` and `p̄` = conj(p)`, then this satisfies Examples: -```jldoctest laurent -julia> using Polynomials; - +```jldoctest julia> z = variable(LaurentPolynomial, :z) LaurentPolynomial(z) @@ -226,9 +222,7 @@ This satisfies for *imaginary* `s`: `conj(p(s)) = p̃(s) = (conj ∘ p)(s) = cco [reference](https://github.com/hurak/PolynomialEquations.jl#symmetrix-conjugate-equation-continuous-time-case) Examples: -```jldoctest laurent -julia> using Polynomials; - +```jldoctest julia> s = 2im 0 + 2im diff --git a/src/polynomials/standard-basis/sparse-polynomial.jl b/src/polynomials/standard-basis/sparse-polynomial.jl index ae2b7cfd..55d66dd3 100644 --- a/src/polynomials/standard-basis/sparse-polynomial.jl +++ b/src/polynomials/standard-basis/sparse-polynomial.jl @@ -10,9 +10,7 @@ advantageous. # Examples: ```jldoctest -julia> using Polynomials - -julia> P = SparsePolynomial; +julia> P = SparsePolynomial; julia> p, q = P([1,2,3]), P([4,3,2,1]) (SparsePolynomial(1 + 2*x + 3*x^2), SparsePolynomial(4 + 3*x + 2*x^2 + x^3)) diff --git a/src/polynomials/standard-basis/standard-basis.jl b/src/polynomials/standard-basis/standard-basis.jl index 0c172dd7..a5971506 100644 --- a/src/polynomials/standard-basis/standard-basis.jl +++ b/src/polynomials/standard-basis/standard-basis.jl @@ -577,7 +577,7 @@ Using constrained least squares, fit a polynomial of the type The degrees in `cs` and those in `J` should not intersect. Example -``` +```julia x = range(0, pi/2, 10) y = sin.(x) P = Polynomial @@ -647,7 +647,10 @@ This representation is useful for *evaluating* the polynomial, but does not lend # Returns -Returns an instance of `ArnoldiFit`. This object can be used to evaluate the polynomial. To manipulate the polynomial, the object can be `convert`ed to other polynomial types, though there may be some loss in accuracy when doing polynomial evaluations afterwards for higher-degree polynomials. +Returns an instance of `ArnoldiFit`. This object can be used to evaluate the +polynomial. To manipulate the polynomial, the object can be `convert`ed to other +polynomial types, though there may be some loss in accuracy when doing +polynomial evaluations afterwards for higher-degree polynomials. # Citations: @@ -665,7 +668,7 @@ Lei-Hong Zhang, Yangfeng Su, Ren-Cang Li. Accurate polynomial fitting and evalua # Examples: -``` +```julia f(x) = 1/(1 + 25x^2) N = 80; xs = [cos(j*pi/N) for j in N:-1:0]; p = fit(Polynomial, xs, f.(xs)); @@ -674,7 +677,7 @@ maximum(abs, p(x) - f(x) for x ∈ range(-1,stop=1,length=500)) # 3.304586010148 maximum(abs, q(x) - f(x) for x ∈ range(-1,stop=1,length=500)) # 1.1939520722092922e-7 ``` -``` +```julia N = 250; xs = [cos(j*pi/N) for j in N:-1:0]; p = fit(Polynomial, xs, f.(xs)); q = fit(ArnoldiFit, xs, f.(xs)); @@ -682,7 +685,7 @@ maximum(abs, p(x) - f(x) for x ∈ range(-1,stop=1,length=500)) # 3.553181862545 maximum(abs, q(x) - f(x) for x ∈ range(-1,stop=1,length=500)) # 8.881784197001252e-16 ``` -``` +```julia p = fit(Polynomial, xs, f.(xs), 10); # least-squares fit q = fit(ArnoldiFit, xs, f.(xs), 10); maximum(abs, q(x) - p(x) for x ∈ range(-1,stop=1,length=500)) # 4.6775083806238626e-14 @@ -736,9 +739,14 @@ end """ ArnoldiFit -A polynomial type produced through fitting a degree ``n`` or less polynomial to data ``(x_1,y_1),…,(x_N, y_N), N ≥ n+1``, This uses Arnoldi orthogonalization to avoid the exponentially ill-conditioned Vandermonde polynomial. See [`Polynomials.polyfitA`](@ref) for details. +A polynomial type produced through fitting a degree ``n`` or less polynomial to +data ``(x_1,y_1),…,(x_N, y_N), N ≥ n+1``, This uses Arnoldi orthogonalization to +avoid the exponentially ill-conditioned Vandermonde polynomial. See +[`Polynomials.polyfitA`](@ref) for details. -This is useful for polynomial evaluation, but other polynomial operations are not defined. Though these fitted polynomials may be converted to other types, for larger degrees this will prove unstable. +This is useful for polynomial evaluation, but other polynomial operations are +not defined. Though these fitted polynomials may be converted to other types, +for larger degrees this will prove unstable. """ struct ArnoldiFit{T, M<:AbstractArray{T,2}, X} <: AbstractPolynomial{T,X} @@ -766,7 +774,9 @@ Base.convert(::Type{P}, p::ArnoldiFit{T,M,X}) where {P <: AbstractPolynomial,T,M compensated_horner(p::P, x) compensated_horner(ps, x) -Evaluate `p(x)` using a compensation scheme of S. Graillat, Ph. Langlois, N. Louve [Compensated Horner Scheme](https://cadxfem.org/cao/Compensation-horner.pdf). Either a `Polynomial` `p` or its coefficients may be passed in. +Evaluate `p(x)` using a compensation scheme of +S. Graillat, Ph. Langlois, N. Louve [Compensated Horner Scheme](https://cadxfem.org/cao/Compensation-horner.pdf). +Either a `Polynomial` `p` or its coefficients may be passed in. The Horner scheme has relative error given by @@ -776,7 +786,11 @@ The compensated Horner scheme has relative error bounded by `|(p(x) - p̂(x))/p(x)| ≤ u + β(n) ⋅ u² ⋅ cond(p, x)`. -As noted, this reflects the accuracy of a backward stable computation performed in doubled working precision `u²`. (E.g., polynomial evaluation of a `Polynomial{Float64}` object through `compensated_horner` is as accurate as evaluation of a `Polynomial{Double64}` object (using the `DoubleFloat` package), but significantly faster. +As noted, this reflects the accuracy of a backward stable computation performed +in doubled working precision `u²`. (E.g., polynomial evaluation of a +`Polynomial{Float64}` object through `compensated_horner` is as accurate as +evaluation of a `Polynomial{Double64}` object (using the `DoubleFloat` package), +but significantly faster. Pointed out in https://discourse.julialang.org/t/more-accurate-evalpoly/45932/5. """ diff --git a/src/rational-functions/fit.jl b/src/rational-functions/fit.jl index 4c0d9799..aa6f13cc 100644 --- a/src/rational-functions/fit.jl +++ b/src/rational-functions/fit.jl @@ -31,12 +31,14 @@ Fit a rational function of the form `pq = (a₀ + a₁x¹ + … + aₘxᵐ) / (1 real poles and arbitrarily high approximation orders on any real interval, regardless of the distribution of the points") - The [RationalApproximations](https://github.com/billmclean/RationalApproximations) package also has implementations of the AAA algorithm. + The [RationalApproximations](https://github.com/billmclean/RationalApproximations) + package also has implementations of the AAA algorithm. - A python library, [polyrat](https://github.com/jeffrey-hokanson/polyrat), has implementations of other algorithms. + A python library, [polyrat](https://github.com/jeffrey-hokanson/polyrat), + has implementations of other algorithms. ## Example -``` +```julia-repl julia> x = variable(Polynomial{Float64}) Polynomial(1.0*x) @@ -87,7 +89,7 @@ Fit a Pade approximant (cf docstring for `Polynomials.pade_fit`) to `r`. Examples: -```jldoctext +```julia-repl julia> using Polynomials, PolynomialRatios julia> x = variable() diff --git a/src/rational-functions/plot-recipes.jl b/src/rational-functions/plot-recipes.jl index 8a27ffa4..04a3f7a0 100644 --- a/src/rational-functions/plot-recipes.jl +++ b/src/rational-functions/plot-recipes.jl @@ -36,7 +36,7 @@ function rational_function_trim(pq, a, b, xlims, ylims) a = isnothing(a) ? xlims[1] : a b = isnothing(b) ? xlims[2] : b - + if isnothing(a) && isnothing(b) u= poly_interval(p) v= poly_interval(q) @@ -64,4 +64,3 @@ function rational_function_trim(pq, a, b, xlims, ylims) xs, ys′ end - diff --git a/src/rational-functions/rational-function.jl b/src/rational-functions/rational-function.jl index 6c039f92..00c48f08 100644 --- a/src/rational-functions/rational-function.jl +++ b/src/rational-functions/rational-function.jl @@ -11,7 +11,7 @@ that operation. For purposes of iteration, a rational function is treated like a two-element container. ## Examples -``` +```julia-repl julia> using Polynomials julia> p,q = fromroots(Polynomial, [1,2,3]), fromroots(Polynomial, [2,3,4]) diff --git a/src/rational-functions/rational-transfer-function.jl b/src/rational-functions/rational-transfer-function.jl index 4c9cf35a..857fef3c 100644 --- a/src/rational-functions/rational-transfer-function.jl +++ b/src/rational-functions/rational-transfer-function.jl @@ -18,14 +18,14 @@ struct RationalTransferFunction{T,X,P<:AbstractPolynomial{T,X},Ts} <: AbstractRa new{T,X,P,Ts}(num,den) end function RationalTransferFunction{T,X,P,Ts}(num::P, den::P,ts::Union{Real,Nothing}) where{T,X,P<:AbstractPolynomial{T,X}, Ts} - check_den(den) - check_Ts(Ts,ts) + check_den(den) + check_Ts(Ts,ts) new{T,X,P,Ts}(num,den) end # can promote constants to polynomials too function RationalTransferFunction{T,X,P,Ts}(num::S, den::P, ts::Union{Real,Nothing}) where{S, T,X,P<:AbstractPolynomial{T,X}, Ts} check_den(den) - check_Ts(Ts,ts) + check_Ts(Ts,ts) new{T,X,P,Ts}(convert(P,num),den) end function RationalTransferFunction{T,X,P,Ts}(num::P,den::S, ts::Union{Real,Nothing}) where{S, T,X,P<:AbstractPolynomial{T,X}, Ts} @@ -69,7 +69,7 @@ end ## helpers for constructors # standardize Ts or throw error function standardize_Ts(Ts) - isnothing(Ts) || Ts >= 0 || Ts == -1 || + isnothing(Ts) || Ts >= 0 || Ts == -1 || throw(ArgumentError("Ts must be either a positive number, 0 (continuous system), or -1 (unspecified)")) Ts′ = isnothing(Ts) ? Ts : Float64(Ts) end @@ -103,7 +103,7 @@ end function promote_Ts(Ts1::Union{Float64,Nothing}, Ts2::Union{Float64,Nothing}) isnothing(Ts1) && (return Ts2) isnothing(Ts2) && (return Ts1) - Ts1 == Ts2 && (return Ts1) + Ts1 == Ts2 && (return Ts1) Ts1 == -1 && (Ts2 > 0 ? (return Ts2) : error("Sampling time mismatch")) Ts2 == -1 && (Ts1 > 0 ? (return Ts1) : error("Sampling time mismatch")) error("Sampling time mismatch") @@ -119,8 +119,8 @@ function Base.promote_rule(::Type{PQ}, ::Type{PQ′}) where {T,X,P,Ts,PQ <: Rati ts = promote_Ts(PQ, PQ′) RationalTransferFunction{S,Y,Q,Val(ts)} end - - + + @@ -135,15 +135,15 @@ function gain(pq::PQ) where {PQ <: RationalTransferFunction} p[end]/q[end] end - + # need to check here # """ - rt = adjoint(r) -Compute the adjoint `rt(λ)` of the rational transfer function `r(λ)` such that for + rt = adjoint(r) +Compute the adjoint `rt(λ)` of the rational transfer function `r(λ)` such that for `r(λ) = num(λ)/den(λ)` we have: - (1) `rt(λ) = conj(num(-λ))/conj(num(-λ))`, if `r.Ts = 0`; - (2) `rt(λ) = conj(num(1/λ))/conj(num(1/λ))`, if `r.Ts = -1` or `r.Ts > 0`. + 1. `rt(λ) = conj(num(-λ))/conj(num(-λ))`, if `r.Ts = 0`; + 2. `rt(λ) = conj(num(1/λ))/conj(num(1/λ))`, if `r.Ts = -1` or `r.Ts > 0`. """ function Base.adjoint(pq::PQ) where {PQ <: RationalTransferFunction} p,q = pqs(pq) diff --git a/src/show.jl b/src/show.jl index 1102ba2c..211232bf 100644 --- a/src/show.jl +++ b/src/show.jl @@ -222,9 +222,6 @@ parentheses. ```jldoctest julia> using Polynomials, DualNumbers - - - julia> Polynomial([Dual(1,2), Dual(3,4)]) Polynomial(1 + 2ɛ + 3 + 4ɛ*x) ``` @@ -232,8 +229,6 @@ Polynomial(1 + 2ɛ + 3 + 4ɛ*x) ```jldoctest julia> using DualNumbers, Polynomials - - julia> function Base.show_unquoted(io::IO, pj::Dual, indent::Int, prec::Int) if Base.operator_precedence(:+) <= prec print(io, "(") @@ -244,8 +239,6 @@ julia> function Base.show_unquoted(io::IO, pj::Dual, indent::Int, prec::Int) end end - - julia> Polynomial([Dual(1,2), Dual(3,4)]) Polynomial((1 + 2ɛ) + (3 + 4ɛ)*x) ```