Skip to content

Commit

Permalink
Remove MutableArithemetics code to new package PolynomialsMutableArit…
Browse files Browse the repository at this point in the history
…hmetics (#445)

* WIP

* edits

* move MutableArithmetics to separate package

* rm tests for MA

* replace tab with space
  • Loading branch information
jverzani authored Aug 15, 2022
1 parent 84ec58f commit 0e1edb9
Show file tree
Hide file tree
Showing 8 changed files with 35 additions and 110 deletions.
1 change: 1 addition & 0 deletions .github/workflows/downstream.yml
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@ jobs:
os: [ubuntu-latest]
package:
- {user: jverzani, repo: SpecialPolynomials.jl, group: All}
- {user: jverzani, repo: PolynomialsMutableArithmetics.jl, group: All}

steps:
- uses: actions/checkout@v2
Expand Down
4 changes: 1 addition & 3 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -2,15 +2,13 @@ name = "Polynomials"
uuid = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
license = "MIT"
author = "JuliaMath"
version = "3.1.8"
version = "3.2.0"

[deps]
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"

[compat]
MutableArithmetics = "0.3,1"
RecipesBase = "0.7, 0.8, 1"
julia = "1.6"

Expand Down
40 changes: 31 additions & 9 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -119,6 +119,7 @@ may be noticeable in some use cases. This is amplified when the coefficients
are for instance `BigInt` or `BigFloat` which are mutable themself.
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
Expand All @@ -127,34 +128,55 @@ 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
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 allocate approximately 70 KiB due to 4 k
allocations. If the resulting vector is already allocated, e.g.,

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.mutable_operate!(MA.add_mul, c, A, b)
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

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.mutable_buffered_operate!(buffer, MA.add_mul, c, A, 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
Expand All @@ -174,7 +196,7 @@ c = MA.operate(*, A1, b1)
MA.mutable_operate!(MA.add_mul, c, A2, b2)
```

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

### Integrals and Derivatives
Expand Down
1 change: 0 additions & 1 deletion src/Polynomials.jl
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,6 @@ include("common.jl")

# Polynomials
include("polynomials/standard-basis.jl")
include("polynomials/mutable-arithmetics.jl")
include("polynomials/Polynomial.jl")
include("polynomials/ImmutablePolynomial.jl")
include("polynomials/SparsePolynomial.jl")
Expand Down
4 changes: 2 additions & 2 deletions src/common.jl
Original file line number Diff line number Diff line change
Expand Up @@ -818,8 +818,8 @@ Base.:-(p::P) where {P <: AbstractPolynomial} = _convert(p, -coeffs(p))

Base.:*(p::AbstractPolynomial, c::Number) = scalar_mult(p, c)
Base.:*(c::Number, p::AbstractPolynomial) = scalar_mult(c, p)
Base.:*(c::T, p::P) where {T, P <: AbstractPolynomial{T}} = scalar_mult(c, p)
Base.:*(p::P, c::T) where {T, P <: AbstractPolynomial{T}} = scalar_mult(p, c)
Base.:*(c::T, p::P) where {T, X, P <: AbstractPolynomial{T,X}} = scalar_mult(c, p)
Base.:*(p::P, c::T) where {T, X, P <: AbstractPolynomial{T,X}} = scalar_mult(p, c)

# implicitly identify c::Number with a constant polynomials
Base.:+(c::Number, p::AbstractPolynomial) = +(p, c)
Expand Down
1 change: 0 additions & 1 deletion src/polynomials/Polynomial.jl
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,6 @@ struct Polynomial{T, X} <: StandardBasisPolynomial{T, X}
end

@register Polynomial
@register_mutable_arithmetic Polynomial



Expand Down
71 changes: 0 additions & 71 deletions src/polynomials/mutable-arithmetics.jl

This file was deleted.

23 changes: 0 additions & 23 deletions test/StandardBasis.jl
Original file line number Diff line number Diff line change
Expand Up @@ -1525,29 +1525,6 @@ end
end
end

import MutableArithmetics
const MA = MutableArithmetics

function alloc_test(f, n)
f() # compile
@test n == @allocated f()
end

@testset "Mutable arithmetics" begin
d = m = n = 4
p(d) = Polynomial(big.(1:d))
z(d) = Polynomial([zero(BigInt) for i in 1:d])
A = [p(d) for i in 1:m, j in 1:n]
b = [p(d) for i in 1:n]
c = [z(2d - 1) for i in 1:m]
buffer = MA.buffer_for(MA.add_mul, typeof(c), typeof(A), typeof(b))
@test buffer isa BigInt
c = [z(2d - 1) for i in 1:m]
MA.buffered_operate!(buffer, MA.add_mul, c, A, b)
@test c == A * b
@test c == MA.operate(*, A, b)
@test 0 == @allocated MA.buffered_operate!(buffer, MA.add_mul, c, A, b)
end

@testset "SparsePolynomial" begin
@test Polynomials.minimumexponent(SparsePolynomial) == typemin(Int)
Expand Down

2 comments on commit 0e1edb9

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

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 v3.2.0 -m "<description of version>" 0e1edb942a0d31570699b7ca16b92144d1e04bff
git push origin v3.2.0

Please sign in to comment.