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

decouple basis from storage type #517

Merged
merged 65 commits into from
Aug 15, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
65 commits
Select commit Hold shift + click to select a range
6f38c91
Merge in initial code for explict basis
jverzani Jul 17, 2023
0015c2f
adding in explicit basis code
jverzani Jul 18, 2023
13d6fe3
WIP: tests pass
jverzani Jul 19, 2023
12ae7da
WIP: tests passing, cleaned up
jverzani Jul 20, 2023
0d90ff6
add polynomial_composition code path
jverzani Jul 20, 2023
dac8895
WIP
jverzani Jul 20, 2023
24a9959
WIP
jverzani Jul 21, 2023
2bbf302
WIP: tests passing
jverzani Jul 22, 2023
27b5c25
adjust to get tests passing
jverzani Jul 24, 2023
987b7a1
cleanup
jverzani Jul 24, 2023
89e8e9f
fix bug in scalar_add
jverzani Jul 24, 2023
16a6d03
WIP: cleanup
jverzani Jul 27, 2023
99937ae
WIP: coeffs
jverzani Jul 27, 2023
8b3bfb4
WIP: adjust hash, coeffs
jverzani Jul 27, 2023
266bf47
WIP: coeffs as before
jverzani Jul 27, 2023
db43bf3
WIP
jverzani Jul 28, 2023
ce64b61
WIP: coeffs tweak, mutable order
jverzani Jul 28, 2023
65f8e07
WIP: doc adjustments
jverzani Jul 28, 2023
b778db0
WIP: fiddle with promotion
jverzani Jul 28, 2023
555549e
WIP: identify DescriptorSystems issue
jverzani Jul 29, 2023
663cf93
WIP: More promotions to match old behaviours
jverzani Jul 29, 2023
e2c715f
WIP
jverzani Jul 29, 2023
b7a17c4
WIP: view
jverzani Aug 1, 2023
4c608f5
WIP: cleanup
jverzani Aug 2, 2023
eb986ab
WIP: cleanup
jverzani Aug 2, 2023
91bd228
WIP: cleanup
jverzani Aug 2, 2023
ae4ef09
WIP: more cleanup
jverzani Aug 2, 2023
f844a17
WIP: run doctests, adjust Chebyshev printing
jverzani Aug 2, 2023
0c13cda
WIP: convert
jverzani Aug 2, 2023
04017cd
WIP: more edits
jverzani Aug 3, 2023
0022fec
WIP: reshuffle
jverzani Aug 3, 2023
5123d41
WIP: change Polynomial type; migrate standard-basis.jl
jverzani Aug 3, 2023
c1bb6da
WIP: tests not working
jverzani Aug 3, 2023
a0958cf
derivative eltype
jverzani Aug 3, 2023
c5b15b7
get tests working
jverzani Aug 3, 2023
ee69d55
WIP: work to shift to new base type for most polys
jverzani Aug 4, 2023
95a3a88
doc fix=true
jverzani Aug 4, 2023
120bd1d
WIP: fixes for downstream
jverzani Aug 4, 2023
523c519
relax error type for nightly
jverzani Aug 4, 2023
31097a8
add skip, not broken for a test
jverzani Aug 4, 2023
e23774d
1.6 tests
jverzani Aug 4, 2023
a92bf2f
couldn't quite down JET; is this right version number?
jverzani Aug 5, 2023
77a8798
aqua test
jverzani Aug 5, 2023
2fe55dd
WIP, update docs, reduce redundancies
jverzani Aug 6, 2023
db21162
tweak test
jverzani Aug 6, 2023
ea6558a
WIP: more cleanup
jverzani Aug 7, 2023
a272fbe
don't test Aqua on 1.6
jverzani Aug 7, 2023
4bffe04
clean up
jverzani Aug 7, 2023
981ce90
test that copy is not alias
jverzani Aug 7, 2023
a93910e
merge master
jverzani Aug 7, 2023
5ddfc08
clean up identified invalidations
jverzani Aug 7, 2023
1b81927
remove tests of deprecated (now removed) functions
jverzani Aug 7, 2023
1193f15
Merge branch 'master' into bandaid
jverzani Aug 7, 2023
b79e354
reorg docs
jverzani Aug 8, 2023
c3c6047
fix sparse *
jverzani Aug 8, 2023
da20ed6
Merge branch 'bandaid' into bandaid-docs
jverzani Aug 8, 2023
4d87649
WIP: work on documentation
jverzani Aug 8, 2023
91406f6
stub version number
jverzani Aug 8, 2023
48e3865
reorg docs
jverzani Aug 8, 2023
993fa54
fix sparse *
jverzani Aug 8, 2023
c61c0ca
WIP: work on documentation
jverzani Aug 8, 2023
ce82631
Merge branch 'bandaid-bak' into bandaid
jverzani Aug 8, 2023
aec4d7f
SP test
jverzani Aug 8, 2023
521c86b
WIP: more cleanup
jverzani Aug 14, 2023
5f4c86f
version bump; signal breaking changes
jverzani Aug 15, 2023
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
11 changes: 8 additions & 3 deletions .github/workflows/downstream.yml
Original file line number Diff line number Diff line change
Expand Up @@ -14,12 +14,17 @@ jobs:
strategy:
fail-fast: false
matrix:
julia-version: [1,1.6]
julia-version: [1]
os: [ubuntu-latest]
package:
- {user: jverzani, repo: SpecialPolynomials.jl, group: All}
- {user: JuliaControl, repo: ControlSystems.jl, group: All}

- {user: andreasvarga, repo: DescriptorSystems.jl, group: All}
- {user: andreasvarga, repo: MatrixPencils.jl, group: All}
- {user: JuliaDSP, repo: DSP.jl, group: All}
- {user: tkluck, repo: GaloisFields.jl, group: All}
- {user: jverzani, repo: SpecialPolynomials.jl, group: All}
- {user: JuliaGNI, repo: QuadratureRules.jl, group: All}
- {user: JuliaGNI, repo: RungeKutta.jl, group: All}
steps:
- uses: actions/checkout@v3
- uses: julia-actions/setup-julia@v1
Expand Down
1 change: 1 addition & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -5,3 +5,4 @@ docs/site/
*.jl.mem

Manifest.toml
archive/
11 changes: 7 additions & 4 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -2,14 +2,15 @@ name = "Polynomials"
uuid = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
license = "MIT"
author = "JuliaMath"
version = "3.2.15"
version = "4.0.0"

[deps]
ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
MakieCore = "20f20a25-4f0e-4fdf-b5d1-57303727442b"
MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"

[weakdeps]
ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
Expand All @@ -26,19 +27,21 @@ ChainRulesCore = "1"
MakieCore = "0.6"
MutableArithmetics = "1"
RecipesBase = "0.7, 0.8, 1"
Setfield = "1"
julia = "1.6"

[extras]
Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
MakieCore = "20f20a25-4f0e-4fdf-b5d1-57303727442b"
ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a"
DualNumbers = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
MakieCore = "20f20a25-4f0e-4fdf-b5d1-57303727442b"
MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[targets]
test = ["ChainRulesCore", "DualNumbers", "LinearAlgebra", "SparseArrays", "OffsetArrays", "SpecialFunctions", "Test"]
test = ["Aqua", "ChainRulesCore", "DualNumbers", "LinearAlgebra", "SparseArrays", "OffsetArrays", "SpecialFunctions", "Test"]
322 changes: 2 additions & 320 deletions README.md
Original file line number Diff line number Diff line change
@@ -1,326 +1,8 @@
# Polynomials.jl

Basic arithmetic, integration, differentiation, evaluation, and root finding over dense univariate [polynomials](https://en.wikipedia.org/wiki/Polynomial).
Basic arithmetic, integration, differentiation, evaluation, root finding, and fitting for
univariate [polynomials](https://en.wikipedia.org/wiki/Polynomial) in [Julia](https://julialang.org/).

[![](https://img.shields.io/badge/docs-stable-blue.svg)](https://JuliaMath.github.io/Polynomials.jl/stable)
[![CI](https://github.com/JuliaMath/Polynomials.jl/actions/workflows/ci.yml/badge.svg)](https://github.com/JuliaMath/Polynomials.jl/actions/workflows/ci.yml)
[![codecov](https://codecov.io/gh/JuliaMath/Polynomials.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/JuliaMath/Polynomials.jl)


## Installation

```julia
(v1.6) pkg> add Polynomials
```

This package supports Julia v1.6 and later.

## Available Types of Polynomials

* `Polynomial` –⁠ standard basis polynomials, $a(x) = a_0 + a_1 x + a_2 x^2 + … + a_n x^n$ for $n ≥ 0$.
* `ImmutablePolynomial` –⁠ standard basis polynomials backed by a [Tuple type](https://docs.julialang.org/en/v1/manual/functions/#Tuples-1) for faster evaluation of values
* `SparsePolynomial` –⁠ standard basis polynomial backed by a [dictionary](https://docs.julialang.org/en/v1/base/collections/#Dictionaries-1) to hold sparse high-degree polynomials
* `LaurentPolynomial` –⁠ [Laurent polynomials](https://docs.julialang.org/en/v1/base/collections/#Dictionaries-1), $a(x) = a_m x^m + … + a_n x^n$ for $m ≤ n$ and $m,n ∈ ℤ$. This is backed by an [offset array](https://github.com/JuliaArrays/OffsetArrays.jl); for example, if $m<0$ and $n>0$, we obtain $a(x) = a_m x^m + … + a_{-1} x^{-1} + a_0 + a_1 x + … + a_n x^n$
* `FactoredPolynomial` –⁠ standard basis polynomials, storing the roots, with multiplicity, and leading coefficient of a polynomial
* `ChebyshevT` –⁠ [Chebyshev polynomials](https://en.wikipedia.org/wiki/Chebyshev_polynomials) of the first kind
* `RationalFunction` - a type for ratios of polynomials.

## Usage

```julia
julia> using Polynomials
```

### Construction and Evaluation

Construct a polynomial from an array (a vector) of its coefficients, lowest order first.

```julia
julia> Polynomial([1,0,3,4])
Polynomial(1 + 3*x^2 + 4*x^3)
```

Optionally, the variable of the polynomial can be specified.

```julia
julia> Polynomial([1,2,3], :s)
Polynomial(1 + 2*s + 3*s^2)
```

Construct a polynomial from its roots.

```julia
julia> fromroots([1,2,3]) # (x-1)*(x-2)*(x-3)
Polynomial(-6 + 11*x - 6*x^2 + x^3)
```

Evaluate the polynomial `p` at `x`.

```julia
julia> p = Polynomial([1, 0, -1]);
julia> p(0.1)
0.99
```

### Arithmetic

Methods are added to the usual arithmetic operators so that they work on polynomials, and combinations of polynomials and scalars.

```julia
julia> p = Polynomial([1,2])
Polynomial(1 + 2*x)

julia> q = Polynomial([1, 0, -1])
Polynomial(1 - x^2)

julia> p - q
Polynomial(2*x + x^2)

julia> p = Polynomial([1,2])
Polynomial(1 + 2*x)

julia> q = Polynomial([1, 0, -1])
Polynomial(1 - x^2)

julia> 2p
Polynomial(2 + 4*x)

julia> 2+p
Polynomial(3 + 2*x)

julia> p - q
Polynomial(2*x + x^2)

julia> p * q
Polynomial(1 + 2*x - x^2 - 2*x^3)

julia> q / 2
Polynomial(0.5 - 0.5*x^2)

julia> q ÷ p # `div`, also `rem` and `divrem`
Polynomial(0.25 - 0.5*x)
```

Most operations involving polynomials with different variables will error.

```julia
julia> p = Polynomial([1, 2, 3], :x);
julia> q = Polynomial([1, 2, 3], :s);
julia> p + q
ERROR: ArgumentError: Polynomials have different indeterminates
```

#### Construction and Evaluation

While polynomials of type `Polynomial` are mutable objects, operations such as
`+`, `-`, `*`, always create new polynomials without modifying its arguments.
The time needed for these allocations and copies of the polynomial coefficients
may be noticeable in some use cases. This is amplified when the coefficients
are for instance `BigInt` or `BigFloat` which are mutable themselves.
This can be avoided by modifying existing polynomials to contain the result
of the operation using the [MutableArithmetics (MA) API](https://github.com/jump-dev/MutableArithmetics.jl).

Consider for instance the following arrays of polynomials
```julia
using Polynomials
d, m, n = 30, 20, 20
p(d) = Polynomial(big.(1:d))
A = [p(d) for i in 1:m, j in 1:n]
b = [p(d) for i in 1:n]
```

In this case, the arrays are mutable objects for which the elements are mutable
polynomials which have mutable coefficients (`BigInt`s).
These three nested levels of mutable objects communicate with the MA
API in order to reduce allocation.
Calling `A * b` requires approximately 40 MiB due to 2 M allocations
as it does not exploit any mutability.

Using

```julia
using PolynomialsMutableArithmetics
```

to register `Polynomials` with `MutableArithmetics`, then multiplying with:

```julia
using MutableArithmetics
const MA = MutableArithmetics
MA.operate(*, A, b)
```

exploits the mutability and hence only allocates approximately 70 KiB due to 4 k
allocations.

If the resulting vector is already allocated, e.g.,

```julia
z(d) = Polynomial([zero(BigInt) for i in 1:d])
c = [z(2d - 1) for i in 1:m]
```

then we can exploit its mutability with

```julia
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:

```julia
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
```julia
MA.@rewrite(A1 * b1 + A2 * b2)
```
is rewritten into
```julia
c = MA.operate!(MA.add_mul, MA.Zero(), A1, b1)
MA.operate!(MA.add_mul, c, A2, b2)
```
which is equivalent to
```julia
c = MA.operate(*, A1, b1)
MA.mutable_operate!(MA.add_mul, c, A2, b2)
```

*Note that currently, only the `Polynomial` type implements the API and it only
implements part of it.*

### Integrals and Derivatives

Integrate the polynomial `p` term by term, optionally adding a constant
term `k`. The degree of the resulting polynomial is one higher than the
degree of `p` (for a nonzero polynomial).

```julia
julia> integrate(Polynomial([1, 0, -1]))
Polynomial(1.0*x - 0.3333333333333333*x^3)

julia> integrate(Polynomial([1, 0, -1]), 2)
Polynomial(2.0 + 1.0*x - 0.3333333333333333*x^3)
```

Differentiate the polynomial `p` term by term. For non-zero
polynomials the degree of the resulting polynomial is one lower than
the degree of `p`.

```julia
julia> derivative(Polynomial([1, 3, -1]))
Polynomial(3 - 2*x)
```

### Root-finding


Return the roots (zeros) of `p`, with multiplicity. The number of
roots returned is equal to the degree of `p`. By design, this is not type-stable, the returned roots may be real or complex.

```julia
julia> roots(Polynomial([1, 0, -1]))
2-element Vector{Float64}:
-1.0
1.0

julia> roots(Polynomial([1, 0, 1]))
2-element Vector{ComplexF64}:
0.0 - 1.0im
0.0 + 1.0im

julia> roots(Polynomial([0, 0, 1]))
2-element Vector{Float64}:
0.0
0.0
```

### Fitting arbitrary data

Fit a polynomial (of degree `deg` or less) to `x` and `y` using a least-squares approximation.

```julia
julia> xs = 0:4; ys = @. exp(-xs) + sin(xs);

julia> fit(xs, ys) |> p -> round.(coeffs(p), digits=4) |> Polynomial
Polynomial(1.0 + 0.0593*x + 0.3959*x^2 - 0.2846*x^3 + 0.0387*x^4)

julia> fit(ChebyshevT, xs, ys, 2) |> p -> round.(coeffs(p), digits=4) |> ChebyshevT
ChebyshevT(0.5413⋅T_0(x) - 0.8991⋅T_1(x) - 0.4238⋅T_2(x))
```

Visual example:

![fit example](https://user-images.githubusercontent.com/14099459/70382587-9e055500-1902-11ea-8952-3f03ae08b7dc.png)

### Other methods

Polynomial objects also have other methods:

* For standard basis polynomials, 0-based indexing is used to extract
the coefficients of `[a0, a1, a2, ...]`; for mutable polynomials,
coefficients may be changed using indexing notation.

* `coeffs`: returns the coefficients

* `degree`: returns the polynomial degree, `length` is number of stored coefficients

* `variable`: returns the polynomial symbol as a polynomial in the underlying type

* `LinearAlgebra.norm`: find the `p`-norm of a polynomial

* `conj`: finds the conjugate of a polynomial over a complex field

* `truncate`: set to 0 all small terms in a polynomial;

* `chop` chops off any small leading values that may arise due to floating point operations.

* `gcd`: greatest common divisor of two polynomials.

* `Pade`: Return the
[Padé approximant](https://en.wikipedia.org/wiki/Pad%C3%A9_approximant) of order `m/n` for a polynomial as a `Pade` object.


## Related Packages

* [StaticUnivariatePolynomials.jl](https://github.com/tkoolen/StaticUnivariatePolynomials.jl) Fixed-size univariate polynomials backed by a Tuple

* [MultiPoly.jl](https://github.com/daviddelaat/MultiPoly.jl) for sparse multivariate polynomials

* [DynamicPolynomials.jl](https://github.com/JuliaAlgebra/DynamicPolynomials.jl) Multivariate polynomials implementation of commutative and non-commutative variables

* [MultivariatePolynomials.jl](https://github.com/JuliaAlgebra/MultivariatePolynomials.jl) for multivariate polynomials and moments of commutative or non-commutative variables

* [PolynomialRings.jl](https://github.com/tkluck/PolynomialRings.jl) A library for arithmetic and algebra with multi-variable polynomials.

* [AbstractAlgebra.jl](https://github.com/wbhart/AbstractAlgebra.jl), [Nemo.jl](https://github.com/wbhart/Nemo.jl) for generic polynomial rings, matrix spaces, fraction fields, residue rings, power series, [Hecke.jl](https://github.com/thofma/Hecke.jl) for algebraic number theory.

* [LaurentPolynomials.jl](https://github.com/jmichel7/LaurentPolynomials.jl) A package for Laurent polynomials.

* [CommutativeAlgebra.jl](https://github.com/KlausC/CommutativeRings.jl) the start of a computer algebra system specialized to discrete calculations with support for polynomials.

* [PolynomialRoots.jl](https://github.com/giordano/PolynomialRoots.jl) for a fast complex polynomial root finder. For larger degree problems, also [FastPolynomialRoots](https://github.com/andreasnoack/FastPolynomialRoots.jl) and [AMRVW](https://github.com/jverzani/AMRVW.jl). For real roots only [RealPolynomialRoots](https://github.com/jverzani/RealPolynomialRoots.jl).


* [SpecialPolynomials.jl](https://github.com/jverzani/SpecialPolynomials.jl) A package providing various polynomial types beyond the standard basis polynomials in `Polynomials.jl`. Includes interpolating polynomials, Bernstein polynomials, and classical orthogonal polynomials.

* [ClassicalOrthogonalPolynomials.jl](https://github.com/JuliaApproximation/ClassicalOrthogonalPolynomials.jl) A Julia package for classical orthogonal polynomials and expansions. Includes `chebyshevt`, `chebyshevu`, `legendrep`, `jacobip`, `ultrasphericalc`, `hermiteh`, and `laguerrel`. The same repository includes `FastGaussQuadrature.jl`, `FastTransforms.jl`, and the `ApproxFun` packages.


## Legacy code

As of v0.7, the internals of this package were greatly generalized and new types and method names were introduced. For compatibility purposes, legacy code can be run after issuing `using Polynomials.PolyCompat`.

## Contributing

If you are interested in contributing, feel free to open an issue or pull request to get started.
Loading