Skip to content

Commit

Permalink
minor doc adjustments (#357)
Browse files Browse the repository at this point in the history
  • Loading branch information
jverzani authored Oct 1, 2021
1 parent 34bd564 commit 31814ca
Show file tree
Hide file tree
Showing 8 changed files with 51 additions and 58 deletions.
2 changes: 1 addition & 1 deletion Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@ name = "Polynomials"
uuid = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
license = "MIT"
author = "JuliaMath"
version = "2.0.14"
version = "2.0.15"

[deps]
Intervals = "d8418881-c3e1-53bb-8760-2df7ec849ed5"
Expand Down
5 changes: 5 additions & 0 deletions docs/make.jl
Original file line number Diff line number Diff line change
@@ -1,6 +1,11 @@
using Documenter
using Polynomials


ENV["PLOTS_TEST"] = "true"
ENV["GKSwstype"] = "100"


DocMeta.setdocmeta!(Polynomials, :DocTestSetup, :(using Polynomials); recursive = true)

makedocs(
Expand Down
10 changes: 5 additions & 5 deletions docs/src/extending.md
Original file line number Diff line number Diff line change
Expand Up @@ -30,11 +30,11 @@ As always, if the default implementation does not work or there are more efficie
| `one`| | Convenience to find constant in new basis |
| `variable`| | Convenience to find monomial `x` in new basis|

Check out both the [`Polynomial`](@ref) and [`ChebyshevT`](@ref) for examples of this interface being extended.
Check out both the [`Polynomial`](@ref) and [`ChebyshevT`](@ref) for examples of this interface being extended.

## Example

The following shows a minimal example where the polynomial aliases the vector defining the coefficients.
The following shows a minimal example where the polynomial aliases the vector defining the coefficients.
The constructor ensures that there are no trailing zeros. The `@register` call ensures a common interface. This example subtypes `StandardBasisPolynomial`, not `AbstractPolynomial`, and consequently inherits the methods above that otherwise would have been required. For other bases, more methods may be necessary to define (again, refer to [`ChebyshevT`](@ref) for an example).

```jldoctest AliasPolynomial
Expand Down Expand Up @@ -75,7 +75,7 @@ julia> p(3)
142
```

For the `Polynomial` type, the default on operations is to copy the array. For this type, it might seem reasonable -- to avoid allocations -- to update the coefficients in place for scalar addition and scalar multiplication.
For the `Polynomial` type, the default on operations is to copy the array. For this type, it might seem reasonable -- to avoid allocations -- to update the coefficients in place for scalar addition and scalar multiplication.

Scalar addition, `p+c`, defaults to `p + c*one(p)`, or polynomial addition, which is not inplace without addition work. As such, we create a new method and an infix operator

Expand Down Expand Up @@ -134,7 +134,7 @@ julia> p ./= 2
AliasPolynomial(3 + 2*x + 3*x^2 + 4*x^3)
```

Trying to divide again would throw an error, as the result would not fit with the integer type of `p`.
Trying to divide again would throw an error, as the result would not fit with the integer type of `p`.

Now `p` is treated as the vector `p.coeffs`, as regards broadcasting, so some things may be surprising, for example this expression returns a vector, not a polynomial:

Expand All @@ -147,4 +147,4 @@ julia> p .+ 2
6
```

The [`Polynomials.PnPolynomial`](@ref) type implements much of this.
The unexported `Polynomials.PnPolynomial` type implements much of this.
8 changes: 5 additions & 3 deletions docs/src/index.md
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
# Polynomials.jl

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

To install the package, run

Expand Down Expand Up @@ -57,6 +57,8 @@ julia> p(1)
```

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.

### Arithmetic

The usual arithmetic operators are overloaded to work on polynomials, and combinations of polynomials and scalars.
Expand Down Expand Up @@ -101,7 +103,7 @@ ERROR: ArgumentError: Polynomials have different indeterminates
[...]
```

Except for operations involving constant polynomials.
Except for operations involving constant polynomials.

```jldoctest
julia> p = Polynomial([1, 2, 3], :x)
Expand Down Expand Up @@ -358,7 +360,7 @@ Most of the root finding algorithms have issues when the roots have
multiplicities. For example, both `ANewDsc` and `Hecke.roots` assume a
square free polynomial. For non-square free polynomials:

* The `Polynomials.Multroot.multroot` function is available (version v"1.2" or greater) for finding the roots of a polynomial and their multiplicities. This is based on work of Zeng.
* The `Polynomials.Multroot.multroot` function is available (version `v"1.2"` or greater) for finding the roots of a polynomial and their multiplicities. This is based on work of Zeng.

Here we see `IntervalRootsFindings.roots` having trouble isolating the roots due to the multiplicites:

Expand Down
5 changes: 2 additions & 3 deletions src/common.jl
Original file line number Diff line number Diff line change
Expand Up @@ -170,10 +170,9 @@ _wlstsq(vand, y, W::AbstractMatrix) = qr(vand' * W * vand) \ (vand' * W * y)
Returns the roots, or zeros, of the given polynomial.
This is calculated via the eigenvalues of the companion matrix. The `kwargs` are passed to the `LinearAlgeebra.eigvals` call.
!!! Note
For non-factored, standard basis polynomials the roots are calculated via the eigenvalues of the companion matrix. The `kwargs` are passed to the `LinearAlgeebra.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
Expand Down
2 changes: 1 addition & 1 deletion src/polynomials/ChebyshevT.jl
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@ julia> Polynomials.evalpoly(5.0, p, false) # bypasses the domain check done in p
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:
!!! note
The Chebyshev polynomials are also implemented in `ApproxFun`, `ClassicalOrthogonalPolynomials.jl`, `FastTransforms.jl`, and `SpecialPolynomials.jl`.
"""
Expand Down
46 changes: 17 additions & 29 deletions src/rational-functions/common.jl
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@ Default methods for basic arithmetic operations are provided.
Numeric methods to cancel common factors, compute the poles, and return the residues are provided.
!!! Note:
!!! note
Requires `VERSION >= v"1.2.0"`
"""
abstract type AbstractRationalFunction{T,X,P} end
Expand Down Expand Up @@ -131,8 +131,8 @@ end

Base.://(p::AbstractPolynomial,q::Number) = p // (q*one(p))
Base.://(p::Number, q::AbstractPolynomial) = (p*one(q)) // q


function Base.copy(pq::PQ) where {PQ <: AbstractRationalFunction}
p,q = pqs(pq)
rational_function(PQ, p, q)
Expand Down Expand Up @@ -165,7 +165,7 @@ _eltype(pq::Type{<:AbstractRationalFunction{T}}) where {T} = Polynomial{T}
_eltype(pq::Type{<:AbstractRationalFunction}) = Polynomial

"""
pqs(pq)
pqs(pq)
Return `(p,q)`, where `pq=p/q`, as polynomials.
"""
Expand All @@ -178,7 +178,7 @@ function Base.isinf(pq::AbstractRationalFunction)
p,q=pqs(pq)
iszero(q) && !iszero(p)
end
function Base.isnan(pq::AbstractRationalFunction)
function Base.isnan(pq::AbstractRationalFunction)
p,q= pqs(pq)
iszero(p) && iszero(q)
end
Expand All @@ -191,7 +191,7 @@ function Base.isone(pq::AbstractRationalFunction)
isconstant(p) && isconstant(q) && p == q
end



## -----

Expand Down Expand Up @@ -236,7 +236,7 @@ function variable(::Type{PQ}) where {PQ <: AbstractRationalFunction}
rational_function(PQ, x, one(x))
end


# use degree as largest degree of p,q after reduction
function degree(pq::AbstractRationalFunction)
pq′ = lowest_terms(pq)
Expand Down Expand Up @@ -355,7 +355,7 @@ function Base.:/(p::R, q::Number) where {R <: AbstractRationalFunction}
rational_function(R, p0, (p1*q))
end
Base.:/(p::AbstractPolynomial, q::PQ) where {PQ <: AbstractRationalFunction} = rational_function(PQ, p,one(p)) / q
function Base.:/(p::PQ, q::AbstractPolynomial) where {PQ <: AbstractRationalFunction}
function Base.:/(p::PQ, q::AbstractPolynomial) where {PQ <: AbstractRationalFunction}
p0,p1 = pqs(p)
rational_function(PQ,p0, p1*q)
end
Expand Down Expand Up @@ -418,17 +418,17 @@ function Base.divrem(pq::PQ; method=:numerical, kwargs...) where {PQ <: Abstract

d,r = divrem(p,q)
(d, rational_function(PQ, r, q))

end




"""
lowest_terms(pq::AbstractRationalFunction, method=:numerical)
Find GCD of `(p,q)`, `u`, and return `(p÷u)//(q÷u)`. Commonly referred to as lowest terms.
* `method`: passed to `gcd(p,q)`
* `kwargs`: passed to `gcd(p,q)`
Expand Down Expand Up @@ -475,11 +475,11 @@ end
If `p/q =d + r/q`, returns `d` and the residues of a rational fraction `r/q`.
First expresses `p/q =d + r/q` with `r` of lower degree than `q` through `divrem`.
First expresses `p/q =d + r/q` with `r` of lower degree than `q` through `divrem`.
Then finds the poles of `r/q`.
For a pole, `λj` of multiplicity `k` there are `k` residues,
For a pole, `λj` of multiplicity `k` there are `k` residues,
`rⱼ[k]/(z-λⱼ)^k`, `rⱼ[k-1]/(z-λⱼ)^(k-1)`, `rⱼ[k-2]/(z-λⱼ)^(k-2)`, …, `rⱼ[1]/(z-λⱼ)`.
The residues are found using this formula:
The residues are found using this formula:
`1/j! * dʲ/dsʲ (F(s)(s - λⱼ)^k` evaluated at `λⱼ` ([5-28](https://stanford.edu/~boyd/ee102/rational.pdf)).
Expand Down Expand Up @@ -521,19 +521,19 @@ julia> p′, q′ = lowest_terms(d);
julia> q′ ≈ (s-1) * (s+1)^2 # works, as q is monic
true
julia> p′ ≈ (-s^2 + s + 1)
julia> p′ ≈ (-s^2 + s + 1)
true
```
!!! Note:
!!! note
There are several areas where numerical issues can arise. The `divrem`, the identification of multiple roots (`multroot`), the evaluation of the derivatives, ...
"""
function residues(pq::AbstractRationalFunction; method=:numerical, kwargs...)


d,r′ = divrem(pq)
r = lowest_terms(r′; method=method, kwargs...)
b,a = pqs(r)
Expand Down Expand Up @@ -590,15 +590,3 @@ function partial_fraction(::Val{:residue}, r, s::PQ) where {PQ}
terms

end












31 changes: 15 additions & 16 deletions src/rational-functions/fit.jl
Original file line number Diff line number Diff line change
Expand Up @@ -11,15 +11,14 @@ Fit a rational function of the form `pq = (a₀ + a₁x¹ + … + aₘxᵐ) / (1
!!! Note:
!!! 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ₙ)`.
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.
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
Expand Down Expand Up @@ -69,7 +68,7 @@ julia> u(variable(pq)) # to see which polynomial is used
"""
function Polynomials.fit(::Type{PQ}, xs::AbstractVector{S}, ys::AbstractVector{T}, m, n; var=:x) where {T,S, PQ<:RationalFunction}


β₁,β₂ = gauss_newton(collect(xs), convert(Vector{float(T)}, ys), m, n)
P = eltype(PQ)
T′ = Polynomials._eltype(P) == nothing ? eltype(β₁) : eltype(P)
Expand Down Expand Up @@ -132,10 +131,10 @@ function pade_fit(p::Polynomial{T}, m::Integer, n::Integer; var=:x) where {T}
d = degree(p)
@assert (0 <= m) && (1 <= n) && (m + n <= d)

# could be much more perfomant
# could be much more perfomant
c = convert(LaurentPolynomial, p) # for better indexing
cs = [c[m+j-i] for j 1:n, i 0:n]

qs′ = cs[:, 2:end] \ cs[:,1]
qs = vcat(1, -qs′)

Expand Down Expand Up @@ -198,15 +197,15 @@ function J!(Jᵣ, xs::Vector{T}, β, n) where {T}
end
nothing
end


function gauss_newton(xs, ys::Vector{T}, n, m, tol=sqrt(eps(T))) where {T}

β = initial_guess(xs, ys, n, m)
model = make_model(n)

Jᵣ = zeros(T, length(xs), 1 + n + m)

Δ = norm(ys, Inf) * tol

ϵₘ = norm(ys - model(xs, β), Inf)
Expand All @@ -216,7 +215,7 @@ function gauss_newton(xs, ys::Vector{T}, n, m, tol=sqrt(eps(T))) where {T}

while no_steps < 25
no_steps += 1

r = ys - model(xs, β)
ϵ = norm(r, Inf)
ϵ < Δ && return (β[1:n+1], β[n+2:end])
Expand All @@ -227,12 +226,12 @@ function gauss_newton(xs, ys::Vector{T}, n, m, tol=sqrt(eps(T))) where {T}
J!(Jᵣ, xs, β, n)
Δᵦ = pinv(Jᵣ' * Jᵣ) * (Jᵣ' * r)
β .-= Δᵦ

end

@warn "no convergence; returning best fit of many steps"
return (βₘ[1:n+1], βₘ[n+2:end])
end
end


end

2 comments on commit 31814ca

@jverzani
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/45917

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v2.0.15 -m "<description of version>" 31814caaf1b1e55d2720512d36b296c9ebed7896
git push origin v2.0.15

Please sign in to comment.