Skip to content

Commit

Permalink
Merge pull request #142 from kalmarek/enh/poly-improvements
Browse files Browse the repository at this point in the history 
Polynomial performance improvements
  • Loading branch information
@Joel-Dahne
Joel-Dahne authored Nov 25, 2021
2 parents bbe59ae + 23695ca commit 57ee7e2
Show file tree
Hide file tree
Showing 4 changed files with 104 additions and 24 deletions.
4 changes: 2 additions & 2 deletions Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "Arblib"
uuid = "fb37089c-8514-4489-9461-98f9c8763369"
authors = ["Marek Kaluba <[email protected]>", "Sascha Timme <Sascha Timme <[email protected]>", "Joel Dahne <[email protected]>"]
version = "0.5.2"
version = "0.6.0"

[deps]
Arb_jll = "d9960996-1013-53c9-9ba4-74a4155039c3"
Expand All @@ -14,5 +14,5 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
[compat]
Arb_jll = "~200.2000"
FLINT_jll = "~200.800"
SpecialFunctions = "1.0"
SpecialFunctions = "1.0, 2"
julia = "1.3"
66 changes: 60 additions & 6 deletions src/poly.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,7 @@ Base.getindex(p::Union{Poly,Series}, I::AbstractRange{<:Integer}) = [p[i] for i

Base.@propagate_inbounds function Base.setindex!(
p::Union{ArbPoly,ArbSeries},
x::ArbLike,
x::Union{ArbLike,_BitSigned},
i::Integer,
)
@boundscheck checkbounds(p, i)
Expand All @@ -76,6 +76,31 @@ Base.@propagate_inbounds function Base.setindex!(p::Union{Poly,Series}, x, i::In
return x
end

"""
ref(p::Union{ArbPoly,ArbSeries,AcbPoly,AcbSeries}, i)
Similar to `p[i]` but instead of an `Arb` or `Acb` returns an `ArbRef`
or `AcbRef` which still shares the memory with the `i`-th entry of
`p[i]`.
For `ArbPoly` and `AcbPoly` it only allows accessing coefficients up
to the degree of the polynomial since higher coefficients are not
guaranteed to be allocated.
For `ArbSeries` and `AcbSeries` it allows accessing coefficients up to
the degree of the series. Note that this degree might not be the same
as the degree of the underlying polynomial (in case higher order
coefficients are zero) but the way they are constructed ensures that
the coefficients will always be initialised.
!!! Note: If you use this to change the coefficient in a way so that the degree
of the polynomial might change you need to normalise the polynomial
afterwards to make sure that Arb recognises the possibly new degree of
the polynomial. If the new degree is the same or lower this can be
done using [`Arblib.normalise!`](@ref). If the new degree is higher
you need to manually set the correct value of
`Arblib.cstruct(p).length` to be one higher than the new degree.
"""
Base.@propagate_inbounds function ref(p::Union{ArbPoly,ArbSeries}, i::Integer)
@boundscheck 0 <= i <= degree(p) || throw(BoundsError(p, i))
ptr = cstruct(p).coeffs + i * sizeof(arb_struct)
Expand All @@ -98,18 +123,32 @@ for (TPoly, TSeries) in [(:ArbPoly, :ArbSeries), (:AcbPoly, :AcbSeries)]

@eval function $TPoly(coeff; prec::Integer = _precision(coeff))
p = $TPoly(prec = prec)
p[0] = coeff
@inbounds p[0] = coeff
return p
end

@eval function $TPoly(coeffs::AbstractVector; prec::Integer = _precision(first(coeffs)))
@eval function $TPoly(
coeffs::Union{Tuple,AbstractVector};
prec::Integer = _precision(first(coeffs)),
)
p = fit_length!($TPoly(prec = prec), length(coeffs))
@inbounds for i = 1:length(coeffs)
p[i-1] = coeffs[i]
end
return p
end

# Add a specialised constructors for the common case of a tuple
# with two elements. This would for example be used when
# constructing a polynomial with a constant plus x, e.g
# ArbPoly((x, 1))
@eval function $TPoly(coeffs::Tuple{Any,Any}; prec::Integer = _precision(first(coeffs)))
p = fit_length!($TPoly(prec = prec), length(coeffs))
@inbounds p[0] = coeffs[1]
@inbounds p[1] = coeffs[2]
return p
end

@eval $TPoly(p::Union{$TPoly,$TSeries}; prec::Integer = precision(p)) =
set!($TPoly(prec = prec), p)
end
Expand All @@ -130,22 +169,37 @@ for (TSeries, TPoly) in [(:ArbSeries, :ArbPoly), (:AcbSeries, :AcbPoly)]

@eval function $TSeries(coeff; degree::Integer = 0, prec::Integer = _precision(coeff))
p = $TSeries(degree = degree, prec = prec)
p[0] = coeff
@inbounds p[0] = coeff
return p
end

@eval function $TSeries(
coeffs::AbstractVector;
coeffs::Union{Tuple,AbstractVector};
degree::Integer = length(coeffs) - 1,
prec::Integer = _precision(first(coeffs)),
)
p = fit_length!($TSeries(degree = degree, prec = prec), degree + 1)
p = $TSeries(degree = degree, prec = prec)
@inbounds for i = 1:min(length(coeffs), degree + 1)
p[i-1] = coeffs[i]
end
return p
end

# Add a specialised constructors for the common case of a tuple
# with two elements. This would for example be used when
# constructing a series with a constant plus x, e.g ArbSeries((x,
# 1))
@eval function $TSeries(
coeffs::Tuple{Any,Any};
degree::Integer = length(coeffs) - 1,
prec::Integer = _precision(first(coeffs)),
)
p = $TSeries(degree = degree, prec = prec)
@inbounds p[0] = coeffs[1]
@inbounds p[1] = coeffs[2]
return p
end

@eval $TSeries(p::Union{$TPoly,$TSeries}; prec::Integer = precision(p)) =
set!($TSeries(degree = degree(p), prec = prec), p)
end
Expand Down
23 changes: 15 additions & 8 deletions test/poly.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,23 +56,30 @@
end

@testset "Constructors" begin
@test TPoly() == TPoly(T(0)) == TPoly(T[0]) == zero(TPoly) == zero(TPoly())
@test TPoly(T(1)) == TPoly(T[1]) == one(TPoly) == one(TPoly())
@test TPoly() ==
TPoly(T(0)) ==
TPoly((0,)) ==
TPoly(T[0]) ==
zero(TPoly) ==
zero(TPoly())
@test TPoly(T(1)) == TPoly((1,)) == TPoly(T[1]) == one(TPoly) == one(TPoly())
@test Arblib.isx(TPoly((0, 1)))
@test Arblib.isx(TPoly(T[0, 1]))
@test TPoly(T[1, 2, 0]) == TPoly(T[1, 2])
@test TPoly(5) == TPoly(5.0) == TPoly([5.0]) == TPoly(T(5))
@test TPoly(TPoly([1, 2])) == TPoly(ArbSeries([1, 2])) == TPoly([1, 2])
@test TPoly((1, 2, 0)) == TPoly(T[1, 2, 0]) == TPoly((1, 2)) == TPoly(T[1, 2])
@test TPoly(5) == TPoly(5.0) == TPoly((5.0,)) == TPoly([5.0]) == TPoly(T(5))
@test TPoly(TPoly((1, 2))) == TPoly(ArbSeries((1, 2))) == TPoly((1, 2))

@test precision(TPoly(prec = 64)) == 64
@test precision(TPoly(T(0), prec = 64)) == 64
@test precision(TPoly((T(0),), prec = 64)) == 64
@test precision(TPoly(T[0], prec = 64)) == 64
@test precision(zero(TPoly(prec = 64))) == 64
@test precision(one(TPoly(prec = 64))) == 64
@test precision(TPoly(TPoly([1, 2], prec = 64))) == 64
@test precision(TPoly(ArbSeries([1, 2], prec = 64))) == 64
@test precision(TPoly(TPoly((1, 2), prec = 64))) == 64
@test precision(TPoly(ArbSeries((1, 2), prec = 64))) == 64

if TPoly == AcbPoly
@test TPoly(AcbSeries([1, 2])) == TPoly(ArbPoly([1, 2])) == TPoly([1, 2])
@test TPoly(AcbSeries((1, 2))) == TPoly(ArbPoly((1, 2))) == TPoly((1, 2))

@test precision(TPoly(AcbSeries(prec = 64))) == 64
@test precision(TPoly(ArbPoly(prec = 64))) == 64
Expand Down
35 changes: 27 additions & 8 deletions test/series.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,27 +57,46 @@
end

@testset "Constructors" begin
@test TSeries() == TSeries(T[0]) == zero(TSeries) == zero(TSeries())
@test TSeries() ==
TSeries(T(0)) ==
TSeries((0,)) ==
TSeries(T[0]) ==
zero(TSeries) ==
zero(TSeries())
@test TSeries(degree = 3) ==
TSeries(T(0), degree = 3) ==
TSeries((0,), degree = 3) ==
TSeries(T[0], degree = 3) ==
zero(TSeries(degree = 3))
@test TSeries(T[1]) == one(TSeries) == one(TSeries())
@test TSeries(T(1)) ==
TSeries((1,)) ==
TSeries(T[1]) ==
one(TSeries) ==
one(TSeries())
@test Arblib.isx(TSeries((0, 1)))
@test Arblib.isx(TSeries(T[0, 1]))
@test TSeries(T[1, 2, 0]) != TSeries(T[1, 2])
@test TSeries([5.0]) == TSeries([5]) == TSeries(T[5])
@test TSeries(TSeries([1, 2])) == TSeries(ArbPoly([1, 2])) == TSeries([1, 2])
@test TSeries((1, 2, 0)) ==
TSeries(T[1, 2, 0]) !=
TSeries((1, 2)) ==
TSeries(T[1, 2])
@test TSeries([5.0]) ==
TSeries([5]) ==
TSeries((5.0,)) ==
TSeries(T[5]) ==
TSeries(T(5))
@test TSeries(TSeries((1, 2))) == TSeries(ArbPoly((1, 2))) == TSeries((1, 2))

@test precision(TSeries(degree = 1, prec = 64)) == 64
@test precision(TSeries(0, degree = 1, prec = 64)) == 64
@test precision(TSeries((T(0),), prec = 64)) == 64
@test precision(TSeries([0], prec = 64)) == 64
@test precision(zero(TSeries(degree = 1, prec = 64))) == 64
@test precision(one(TSeries(degree = 1, prec = 64))) == 64
@test precision(TSeries(TSeries([1, 2], prec = 64))) == 64
@test precision(TSeries(ArbPoly([1, 2], prec = 64))) == 64
@test precision(TSeries(TSeries((1, 2), prec = 64))) == 64
@test precision(TSeries(ArbPoly((1, 2), prec = 64))) == 64

if TSeries == AcbSeries
@test TSeries(AcbPoly([1, 2])) == TSeries(ArbSeries([1, 2])) == TSeries([1, 2])
@test TSeries(AcbPoly((1, 2))) == TSeries(ArbSeries((1, 2))) == TSeries((1, 2))
@test Arblib.degree(TSeries(ArbSeries(degree = 4))) == 4
@test precision(TSeries(AcbPoly(prec = 64))) == 64
@test precision(TSeries(ArbSeries(prec = 64))) == 64
Expand Down

2 comments on commit 57ee7e2

@Joel-Dahne
Copy link
Collaborator Author

Choose a reason for hiding this comment

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

@JuliaRegistrator register()

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

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 v0.6.0 -m "<description of version>" 57ee7e2fd534aa4750603de801c780c8f072b734
git push origin v0.6.0

Please sign in to comment.