Skip to content

Commit

Permalink
[WIP] Ring type (issue #349, #234) (#350)
Browse files Browse the repository at this point in the history
* drop T <: Number

* add `conv` alternative for `fastconv`

* restrict tests by VERSION

* doc updates
  • Loading branch information
jverzani authored Jul 7, 2021
1 parent 11457cf commit 60ba676
Show file tree
Hide file tree
Showing 15 changed files with 639 additions and 191 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.12"
version = "2.0.13"

[deps]
Intervals = "d8418881-c3e1-53bb-8760-2df7ec849ed5"
Expand Down
97 changes: 91 additions & 6 deletions docs/src/index.md
Original file line number Diff line number Diff line change
Expand Up @@ -161,7 +161,7 @@ By design, this is not type-stable; the returned roots may be real or complex.

The default `roots` function uses the eigenvalues of the
[companion](https://en.wikipedia.org/wiki/Companion_matrix) matrix for
a polynomial. This is an `𝑶(n^3)` operation.
a polynomial. This is an `𝑶(n^3)` operation.

For polynomials with `BigFloat` coefficients, the
`GenericLinearAlgebra` package can be seamlessly used:
Expand Down Expand Up @@ -204,7 +204,7 @@ julia> PolynomialRoots.roots(coeffs(p))
1.0000000000000002 + 0.0im
```

The roots are always returned as complex numbers.
The roots are always returned as complex numbers.


* The
Expand All @@ -226,7 +226,7 @@ julia> AMRVW.roots(float.(coeffs(p)))
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.
Expand Down Expand Up @@ -429,7 +429,7 @@ Fit a polynomial (of degree `deg`) to `x` and `y` using polynomial interpolation
using Plots, Polynomials
xs = range(0, 10, length=10)
ys = @. exp(-xs)
f = fit(xs, ys) # degree = length(xs) - 1
f = fit(xs, ys) # degree = length(xs) - 1
f2 = fit(xs, ys, 2) # degree = 2
scatter(xs, ys, markerstrokewidth=0, label="Data")
Expand Down Expand Up @@ -583,7 +583,7 @@ Polynomial(4 + 2*x + 3*x^2)

If `q` is non-constant, such as `variable(Polynomial, :y)`, then there would be an error due to the mismatched symbols. (The mathematical result would need a multivariate polynomial, not a univariate polynomial, as this package provides.)

The same conversion is done for polynomial multiplication: constant polynomials are treated as numbers; non-constant polynomials must have their symbols match.
The same conversion is done for polynomial multiplication: constant polynomials are treated as numbers; non-constant polynomials must have their symbols match.

There is an oddity -- though the following two computations look the same, they are technically different:

Expand All @@ -595,7 +595,7 @@ julia> one(Polynomial, :y) + one(Polynomial, :x)
Polynomial(2.0)
```

Both are constant polynomials over `Int`, but the first has the indeterminate `:y`, the second `:x`.
Both are constant polynomials over `Int`, but the first has the indeterminate `:y`, the second `:x`.

This technical difference causes no issues with polynomial addition or multiplication, as there constant polynomials are treated as numbers, but can be an issue when constant polynomials are used as array elements.

Expand Down Expand Up @@ -669,6 +669,91 @@ julia> [1 p; p 1] + [1 2one(q); 3 4] # array{P{T,:x}} + array{P{T,:y}}
Though were a non-constant polynomial with indeterminate `y` replacing
`2one(q)` above, that addition would throw an error.


## Non-number types for `T`

The coefficients of the polynomial may be non-number types, such as matrices or other polynomials, albeit not every operation is fully supported.

For example, a polynomial with matrix coefficients, might be constructed with:

```jldoctest non_number
julia> using Polynomials
julia> a,b,c = [1 0;2 1], [1 0; 3 1], [1 0; 4 1]
([1 0; 2 1], [1 0; 3 1], [1 0; 4 1])
julia> p = Polynomial([a,b,c])
Polynomial([1 0; 2 1] + [1 0; 3 1]*x + [1 0; 4 1]*x^2)
julia> q = derivative(p)
Polynomial([1 0; 3 1] + [2 0; 8 2]*x)
```

Various operations are available, `derivative` was shown above, here are the vector-space operations:

```jldoctest non_number
julia> 2p
Polynomial([2 0; 4 2] + [2 0; 6 2]*x + [2 0; 8 2]*x^2)
julia> p + q
Polynomial([2 0; 5 2] + [3 0; 11 3]*x + [1 0; 4 1]*x^2)
```

polynomial multiplication:

```jldoctest non_number
julia> p * q
Polynomial([1 0; 5 1] + [3 0; 18 3]*x + [3 0; 21 3]*x^2 + [2 0; 16 2]*x^3)
```

polynomial evaluation, here either with a scalar or a matrix:

```jldoctest non_number
julia> p(2)
2×2 Matrix{Int64}:
7 0
24 7
julia> p(b)
2×2 Matrix{Int64}:
3 0
18 3
```

But if the type `T` lacks support of some generic functions, such as `zero(T)` and `one(T)`, then there may be issues. For example, when `T <: AbstractMatrix` the output of `p-p` is an error, as the implementation assumes `zero(T)` is defined. For static arrays, this isn't an issue, as there is support for `zero(T)`. Other polynomial types, such as `SparsePolynomial` have less support, as some specialized methods assume more of the generic interface be implemented.

Similarly, using polynomials for `T` is a possibility:

```jldoctest non_number
julia> a,b,c = Polynomial([1],:y), Polynomial([0,1],:y), Polynomial([0,0,1],:y)
(Polynomial(1), Polynomial(y), Polynomial(y^2))
julia> p = Polynomial([a,b,c], :x)
Polynomial(Polynomial(1) + Polynomial(y)*x + Polynomial(y^2)*x^2)
julia> q = derivative(p)
Polynomial(Polynomial(y) + Polynomial(2*y^2)*x)
```

Again, much works:

```jldoctest non_number
julia> 2p
Polynomial(Polynomial(2) + Polynomial(2*y)*x + Polynomial(2*y^2)*x^2)
julia> p + q
Polynomial(Polynomial(1 + y) + Polynomial(y + 2*y^2)*x + Polynomial(y^2)*x^2)
julia> p(2)
Polynomial(1 + 2*y + 4*y^2)
julia> p(b)
Polynomial(1 + y^2 + y^4)
```

But much doesn't. For example, implicit promotion can fail. For example, the scalar multiplication `p * b` will fail, as the methods assume this is the fallback polynomial multiplication and not the intended scalar multiplication.


## Rational functions

The package provides support for rational functions -- fractions of polynomials (for most types). The construction of the basic type mirrors the construction of rational numbers.
Expand Down
46 changes: 39 additions & 7 deletions src/abstract.jl
Original file line number Diff line number Diff line change
@@ -1,17 +1,50 @@
export AbstractPolynomial



const SymbolLike = Union{AbstractString,Char,Symbol}

"""
AbstractPolynomial{T,X}
An abstract container for various polynomials.
An abstract type for various polynomials.
A polynomial type holds an indeterminate `X`; coefficients of type `T`, stored in some container type; and an implicit basis, such as the standard basis.
# Properties
- `coeffs` - The coefficients of the polynomial
# The type `T`
`T` need not be `T <: Number`, at the moment it is not restricted
Some `T`s will not be successful
* scalar mult: `c::Number * p::Polynomial` should be defined
* scalar mult: `c::T * p:Polynomial{T}` An ambiguity when `T <: AbstractPolynomial`
* scalar mult: `p:Polynomial{T} * c::T` need not commute
* scalar add/sub: `p::Polynomial{T} + q::Polynomial{T}` should be defined
* scalar sub: `p::Polynomial{T} - p::Polynomial{T}` generally needs `zeros(T,1)` defined for `zero(Polynomial{T})`
* poly mult: `p::Polynomial{T} * q::Polynomial{T}` Needs "`T * T`" defined (e.g. `Base.promote_op(*, Vector{Int}, Vector{Int}))` needs to be something.)
* poly powers: `p::Polynomial{T}^2` needs "`T^2`" defined
* implicit promotion: `p::Polynomial{T} + c::Number` needs `convert(T, c)` defined
* implicit promotion: `p::Polynomial{T} + c::T` ambiguity if `T <: AbstractPolynomial`
* evaluation: `p::Polynomial{T}(s::Number)`
* evaluation `p::Polynomial{T}(c::T)` needs `T*T` defined
* derivatives: `derivative(p::Polynomial{T})`
* integrals: `integrate(p::Polynomial{T})`
"""
abstract type AbstractPolynomial{T,X} end


## -----
# We want ⟒(P{α…,T}) = P{α…}; this default
# works for most cases
(P::Type{<:AbstractPolynomial}) = constructorof(P)
Expand All @@ -34,7 +67,7 @@ Polynomials.@register MyPolynomial
```
# Implementations
This will implement simple self-conversions like `convert(::Type{MyPoly}, p::MyPoly) = p` and creates two promote rules. The first allows promotion between two types (e.g. `promote(Polynomial, ChebyshevT)`) and the second allows promotion between parametrized types (e.g. `promote(Polynomial{T}, Polynomial{S})`).
This will implement simple self-conversions like `convert(::Type{MyPoly}, p::MyPoly) = p` and creates two promote rules. The first allows promotion between two types (e.g. `promote(Polynomial, ChebyshevT)`) and the second allows promotion between parametrized types (e.g. `promote(Polynomial{T}, Polynomial{S})`).
For constructors, it implements the shortcut for `MyPoly(...) = MyPoly{T}(...)`, singleton constructor `MyPoly(x::Number, ...)`, conversion constructor `MyPoly{T}(n::S, ...)`, and `variable` alternative `MyPoly(var=:x)`.
"""
Expand All @@ -44,15 +77,15 @@ macro register(name)
Base.convert(::Type{P}, p::P) where {P<:$poly} = p
function Base.convert(P::Type{<:$poly}, p::$poly{T,X}) where {T,X}
isconstant(p) && return constructorof(P){eltype(P),indeterminate(P)}(constantterm(p))
constructorof(P){eltype(P), indeterminate(P,p)}(coeffs(p))
constructorof(P){eltype(P), indeterminate(P,p)}(_coeffs(p))
end
Base.promote(p::P, q::Q) where {X, T, P <:$poly{T,X}, Q <: $poly{T,X}} = p,q
Base.promote_rule(::Type{<:$poly{T,X}}, ::Type{<:$poly{S,X}}) where {T,S,X} = $poly{promote_type(T, S),X}
Base.promote_rule(::Type{<:$poly{T,X}}, ::Type{S}) where {T,S<:Number,X} =
$poly{promote_type(T, S),X}
$poly(coeffs::AbstractVector{T}, var::SymbolLike = :x) where {T} =
$poly{T, Symbol(var)}(coeffs)
$poly{T}(x::AbstractVector{S}, var::SymbolLike = :x) where {T,S<:Number} =
$poly{T}(x::AbstractVector{S}, var::SymbolLike = :x) where {T,S} =
$poly{T,Symbol(var)}(T.(x))
function $poly(coeffs::G, var::SymbolLike=:x) where {G}
!Base.isiterable(G) && throw(ArgumentError("coeffs is not iterable"))
Expand All @@ -76,11 +109,11 @@ macro registerN(name, params...)
αs = tuple(esc.(params)...)
quote
Base.convert(::Type{P}, q::Q) where {$(αs...),T, P<:$poly{$(αs...),T}, Q <: $poly{$(αs...),T}} = q
Base.convert(::Type{$poly{$(αs...)}}, q::Q) where {$(αs...),T, Q <: $poly{$(αs...),T}} = q
Base.convert(::Type{$poly{$(αs...)}}, q::Q) where {$(αs...),T, Q <: $poly{$(αs...),T}} = q
Base.promote(p::P, q::Q) where {$(αs...),T, X, P <:$poly{$(αs...),T,X}, Q <: $poly{$(αs...),T,X}} = p,q
Base.promote_rule(::Type{<:$poly{$(αs...),T,X}}, ::Type{<:$poly{$(αs...),S,X}}) where {$(αs...),T,S,X} =
$poly{$(αs...),promote_type(T, S),X}
Base.promote_rule(::Type{<:$poly{$(αs...),T,X}}, ::Type{S}) where {$(αs...),T,X,S<:Number} =
Base.promote_rule(::Type{<:$poly{$(αs...),T,X}}, ::Type{S}) where {$(αs...),T,X,S<:Number} =
$poly{$(αs...),promote_type(T,S),X}

function $poly{$(αs...),T}(x::AbstractVector{S}, var::SymbolLike = :x) where {$(αs...),T,S}
Expand All @@ -98,4 +131,3 @@ macro registerN(name, params...)
(p::$poly)(x) = evalpoly(x, p)
end
end

Loading

3 comments on commit 60ba676

@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.

@JuliaRegister register

@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/40466

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.13 -m "<description of version>" 60ba67647ecbe9d9dced8733326f7360b03c724b
git push origin v2.0.13

Please sign in to comment.