Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Doctests and documentation patch #571

Merged
merged 9 commits into from
Jun 7, 2024
Merged
Show file tree
Hide file tree
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
Add more doctests and syntax highlighting tags
  • Loading branch information
abhro committed May 18, 2024
commit 5356352918c979ccb78b8d3a7c29d869297b41d2
1 change: 1 addition & 0 deletions docs/Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -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"

Expand Down
28 changes: 14 additions & 14 deletions docs/src/extending.md
Original file line number Diff line number Diff line change
Expand Up @@ -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

```

Expand Down
42 changes: 20 additions & 22 deletions docs/src/index.md
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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))

Expand Down Expand Up @@ -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)

Expand All @@ -209,7 +208,7 @@ julia> roots(p)
to univariate polynomials that is more performant than `roots`. It
is based on an algorithm of Skowron and Gould.

```
```julia-repl
julia> import PolynomialRoots # import as `roots` conflicts

julia> p = fromroots(Polynomial, [1,2,3])
Expand All @@ -234,7 +233,7 @@ The roots are always returned as complex numbers.
package implements the algorithm in Julia, allowing the use of other
number types.

```
```julia-repl
julia> using AMRVW

julia> AMRVW.roots(float.(coeffs(p)))
Expand All @@ -249,11 +248,11 @@ The roots are returned as complex numbers.
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-repl
julia> filter(isreal, AMRVW.roots(rand(1001) .- 1/2))
2-element Vector{ComplexF64}:
0.993739974989572 + 0.0im
Expand All @@ -262,7 +261,7 @@ julia> filter(isreal, AMRVW.roots(rand(1001) .- 1/2))

* 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-repl
julia> import Hecke # import as `roots` conflicts

julia> Qx, x = Hecke.PolynomialRing(Hecke.QQ)
Expand All @@ -280,7 +279,7 @@ julia> Hecke.roots(q)

This next polynomial has 3 real roots, 2 of which are in a cluster; `Hecke` quickly identifies them:

```
```julia-repl
julia> p = -1 + 254*x - 16129*x^2 + 1*x^17
x^17 - 16129*x^2 + 254*x - 1

Expand All @@ -300,7 +299,7 @@ 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-repl
julia> using Polynomials, IntervalArithmetic

julia> import IntervalRootFinding # its `roots` method conflicts with `roots`
Expand All @@ -321,7 +320,7 @@ The output is a set of intervals. Those flagged with `:unique` are guaranteed to

* The `RealPolynomialRoots` package provides a function `ANewDsc` to find isolating intervals for the roots of a square-free polynomial, specified through its coefficients:

```
```julia-repl
julia> using RealPolynomialRoots

julia> st = ANewDsc(coeffs(p))
Expand All @@ -333,7 +332,7 @@ There were 3 isolating intervals found:

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...
Expand All @@ -343,7 +342,7 @@ julia> refine_roots(st)

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)

Expand All @@ -356,7 +355,7 @@ There were 3 isolating intervals found:

`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)
Expand All @@ -380,7 +379,7 @@ square free polynomial. For non-square free polynomials:

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)

Expand All @@ -397,7 +396,7 @@ julia> IntervalRootFinding.roots(x -> p(x), 0..10)

The `roots` function identifies the roots, but the multiplicities would need identifying:

```
```julia-repl
julia> roots(p)
5-element Vector{Float64}:
1.000000000000011
Expand All @@ -410,14 +409,14 @@ julia> roots(p)

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)
Expand All @@ -431,7 +430,7 @@ julia> IntervalRootFinding.roots(x -> q(x), 0..10)

Similarly:

```
```julia-repl
julia> ANewDsc(coeffs(q))
There were 3 isolating intervals found:
[2.62…, 3.62…]₂₅₆
Expand Down Expand Up @@ -817,7 +816,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))
```
Expand Down
1 change: 0 additions & 1 deletion docs/src/polynomials/chebyshev.md
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down
29 changes: 19 additions & 10 deletions src/common.jl
Original file line number Diff line number Diff line change
Expand Up @@ -105,7 +105,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.


"""
Expand Down Expand Up @@ -180,7 +182,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
Expand Down Expand Up @@ -230,7 +234,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)
Expand All @@ -241,7 +246,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)
Expand All @@ -251,7 +257,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)

Expand All @@ -263,15 +270,17 @@ 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


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)
Expand Down Expand Up @@ -1168,9 +1177,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))
Expand Down
Loading