Mutable arithmetic
The high level interface can be combined with the low level wrapper to allow for efficient computations using mutable arithmetic. Making use of mutable arithmetic can in some cases have a significant impact on the performance of the code, in particular for multithreaded code (where avoiding to run the GC is generally more important) and code running many iterations of simple computations.
In the future it would be nice to have an interface to MutableArithmetics.jl, see #118.
The following methods are useful for mutating part of a value
Arblib.radref — Function
radref(x::ArbLike, prec = precision(x))Return a MagRef referencing the radius of x.
Arblib.midref — Function
midref(x::ArbLike, prec = precision(x))Return an ArfRef referencing the midpoint of x.
Arblib.realref — Function
realref(z::AcfLike, prec = precision(z))Return an ArfRef referencing the real part of z.
realref(z::AcbLike, prec = precision(z))Return an ArbRef referencing the real part of z.
Arblib.imagref — Function
imagref(z::AcfLike, prec = precision(z))Return an ArfRef referencing the imaginary part of z.
imagref(z::AcbLike, prec = precision(z))Return an ArbRef referencing the imaginary part of z.
Arblib.ref — Function
ref(v::Union{ArbVectorOrRef,AcbVectorOrRef}, i)Similar to v[i] but instead of an Arb or Acb returns an ArbRef or AcbRef which still shares the memory with the i-th entry of v.
ref(A::Union{ArbMatrixOrRef,AcbMatrixOrRef}, i, j)Similar to A[i,j] but instead of an Arb or Acb returns an ArbRef or AcbRef which still shares the memory with the (i,j)-th entry of A.
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 shares the memory with the i-th coefficient of p.
Using ref for reading coefficients is always safe, but if the coefficient is mutated then care has to be taken. See the comment further down for how to handle mutation.
It only allows accessing coefficients that are allocated. For ArbPoly and AcbPoly this is typically all coefficients up to the degree of the polynomial, but can be higher if for example Arblib.fit_length! is used. For ArbSeries and AcbSeries all coefficients up to the degree of the series are guaranteed to be allocated, even if the underlying polynomial has a lower degree.
If the coefficient is mutated in a way that the degree of the polynomial changes then one needs to also update the internally stored length of the polynomial.
- If the new degree is the same or lower this can be achieved with
Arblib.normalise!(p) - If the new degree is higher you need to manually set the length. This can be achieved with
whereArblib.set_length!(p, len) Arblib.normalise!(p)lenis one higher than (an upper bound of) the new degree.
See the extended help for more details.
Extended help
Here is an example were the leading coefficient is mutated so that the degree is lowered.
julia> p = ArbPoly([1, 2], prec = 64) # Polynomial of degree 1
1.0 + 2.0⋅x
julia> Arblib.zero!(Arblib.ref(p, 1)) # Set leading coefficient to 0
0
julia> Arblib.degree(p) # The degree is still reported as 1
1
julia> length(p) # And the length as 2
2
julia> p # Printing gives weird results
1.0 +
julia> Arblib.normalise!(p) # Normalising the polynomial updates the degree
1.0
julia> Arblib.degree(p) # This is now correct
0
julia> p # And so is printing
1.0Here is an example when a coefficient above the degree is mutated.
julia> q = ArbSeries([1, 2, 0], prec = 64) # Series of degree 3 with leading coefficient zero
1.0 + 2.0⋅x + 𝒪(x^3)
julia> Arblib.one!(Arblib.ref(q, 2)) # Set the leading coefficient to 1
1.0
julia> q # The leading coefficient is not picked up
1.0 + 2.0⋅x + 𝒪(x^3)
julia> Arblib.degree(q.poly) # The degree of the underlying polynomial is still 1
1
julia> Arblib.set_length!(q, 3) # After explicitly setting the length to 3 it is ok
1.0 + 2.0⋅x + 1.0⋅x^2 + 𝒪(x^3)Examples
Compare computing $\sqrt{x^2 + y^2}$ using mutable arithmetic with the default.
julia> using Arblib, BenchmarkToolsjulia> x = Arb(1 // 3)[0.33333333333333333333333333333333333333333333333333333333333333333333333333333 +/- 4.78e-78]julia> y = Arb(1 // 5)[0.20000000000000000000000000000000000000000000000000000000000000000000000000000 +/- 3.89e-78]julia> res = zero(x)0julia> f(x, y) = sqrt(x^2 + y^2)f (generic function with 1 method)julia> f!(res, x, y) = begin Arblib.sqr!(res, x) Arblib.fma!(res, res, y, y) return Arblib.sqrt!(res, res) endf! (generic function with 1 method)julia> @benchmark f($x, $y) samples=10000 evals=500BenchmarkTools.Trial: 7499 samples with 500 evaluations per sample. Range (min … max): 545.656 ns … 559.319 μs ┊ GC (min … max): 0.00% … 79.97% Time (median): 570.902 ns ┊ GC (median): 0.00% Time (mean ± σ): 1.342 μs ± 17.675 μs ┊ GC (mean ± σ): 39.58% ± 3.06% ▁▃▃▇█▇▆▄▄▆▆▄▄▂▁ ▁▁ ▁▃▃▂▁ ▁▂▁ ▂ █████████████████▇▇▇▆████▇▇▇█████▇████▇▅▅▄▅▅▄▅▅▅▄▃▅▄▄▄▅▅▆▄▅▅▅ █ 546 ns Histogram: log(frequency) by time 808 ns < Memory estimate: 416 bytes, allocs estimate: 8.julia> @benchmark f!($res, $x, $y) samples=10000 evals=500BenchmarkTools.Trial: 10000 samples with 500 evaluations per sample. Range (min … max): 331.496 ns … 622.138 ns ┊ GC (min … max): 0.00% … 0.00% Time (median): 352.916 ns ┊ GC (median): 0.00% Time (mean ± σ): 354.935 ns ± 11.749 ns ┊ GC (mean ± σ): 0.00% ± 0.00% ▁▁▆▇█▇▆▅▄▃▅▅▄▅▃ ▁ ▂▁▂▂▂▃▃▃▃▃▃▃▅▆▇██████████████████▇▆▅▅▆▅▅▆▆▅▅▄▅▄▅▄▅▄▄▃▃▃▂▂▂▂▂▂ ▅ 331 ns Histogram: frequency by time 382 ns < Memory estimate: 0 bytes, allocs estimate: 0.
This implementation of the function f! doesn't handle aliasing between res and x. Most cases of aliasing between two variables can be checked for with === (so res === x in this case). Using === does however not catch all possible cases of aliasing, for example it would not catch the aliasing between Arblib.realref(z) and Arblib.realref(z, prec = 2precision(z)).
Set the radius of the real part of an Acb.
julia> using Arblibjulia> z = Acb(1, 2)1.0 + 2.0imjulia> Arblib.set!(Arblib.radref(Arblib.realref(z)), 1e-10)1.0e-10julia> z[1.000000000 +/- 1.01e-10] + 2.0im
Compare a naive implementation of polynomial evaluation using mutable arithmetic with one not using using it.
julia> using Arblib, BenchmarkToolsjulia> p = ArbPoly(1:10)1.0 + 2.0⋅x + 3.0⋅x^2 + 4.0⋅x^3 + 5.0⋅x^4 + 6.0⋅x^5 + 7.0⋅x^6 + 8.0⋅x^7 + 9.0⋅x^8 + 10.0⋅x^9julia> x = Arb(1 // 3)[0.33333333333333333333333333333333333333333333333333333333333333333333333333333 +/- 4.78e-78]julia> res = zero(x)0julia> function eval(p, x) res = zero(x) xi = one(x) for i in 0:Arblib.degree(p) res += p[i] * xi xi *= x end return res endeval (generic function with 2 methods)julia> function eval!(res, p, x) Arblib.zero!(res) xi = one(x) for i in 0:Arblib.degree(p) Arblib.addmul!(res, Arblib.ref(p, i), xi) Arblib.mul!(xi, xi, x) end return res endeval! (generic function with 1 method)julia> @benchmark eval($p, $x) samples = 10000 evals = 30BenchmarkTools.Trial: 10000 samples with 30 evaluations per sample. Range (min … max): 2.480 μs … 11.156 ms ┊ GC (min … max): 0.00% … 75.73% Time (median): 2.746 μs ┊ GC (median): 0.00% Time (mean ± σ): 8.191 μs ± 200.983 μs ┊ GC (mean ± σ): 40.65% ± 1.71% ▄██▆▅▃▁ ▂▃▂▁ ▃▅▅▅▆▆▅▃▃▁ ▁▁▂▁▁ ▂ ████████▇▆▇█████▆▇▇▅▄▄▅▃▃▄▅▄████████████████████▇▇▅▅▃▃▃▂▄▅▄ █ 2.48 μs Histogram: log(frequency) by time 4.67 μs < Memory estimate: 3.72 KiB, allocs estimate: 70.julia> @benchmark eval!($res, $p, $x) samples = 10000 evals = 30BenchmarkTools.Trial: 10000 samples with 30 evaluations per sample. Range (min … max): 1.006 μs … 2.277 μs ┊ GC (min … max): 0.00% … 0.00% Time (median): 1.024 μs ┊ GC (median): 0.00% Time (mean ± σ): 1.037 μs ± 65.069 ns ┊ GC (mean ± σ): 0.00% ± 0.00% ▃▆█▇▆▅▄▃▂ ▁ ▂ ███████████▇▇▇▇▄▅▇▆▇▅▅▃▄▃▄▃▁▄▁▁▁▁▁▃▃▁▁▁▃▁▁▁▁▁▁▁▁▁▁▃▆████▇█ █ 1.01 μs Histogram: log(frequency) by time 1.34 μs < Memory estimate: 144 bytes, allocs estimate: 3.julia> @benchmark $p($x) samples = 10000 evals = 30 # Arb implementation for referenceBenchmarkTools.Trial: 10000 samples with 30 evaluations per sample. Range (min … max): 752.400 ns … 2.471 μs ┊ GC (min … max): 0.00% … 0.00% Time (median): 771.100 ns ┊ GC (median): 0.00% Time (mean ± σ): 783.308 ns ± 66.824 ns ┊ GC (mean ± σ): 0.00% ± 0.00% ▁▅██▆▅▄▃▂▁ ▁ ▂ ███████████▇▇▇▇▆▆▆▅▆▆▅▅▄▅▃▁▄▄▄▁▄▁▄▁▁▁▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▃▃▅███▇█ █ 752 ns Histogram: log(frequency) by time 1.09 μs < Memory estimate: 144 bytes, allocs estimate: 3.
Mutable arithmetic handling both Arb/Acb and ArbSeries/AcbSeries
Since Arblib.jl version 1.3.0 adjustments have been made to the low level wrapper to make it easier to write code using mutable arithmetic that can handle both the scalar types Arb and Acb as well as the series types ArbSeries and AcbSeries. For example the following implementation of sin(atan(x) / x) can support all of these types.
julia> using Arblibjulia> function f!(res, x) Arblib.atan!(res, x) Arblib.div!(res, res, x) Arblib.sin!(res, res) return res endf! (generic function with 1 method)julia> f!(Arb(prec = 53), Arb(2))[0.525731112119133 +/- 6.90e-16]julia> f!(Acb(prec = 53), Acb(2, 3))[0.277808890086608 +/- 8.76e-16] + [-0.283570391073164 +/- 7.96e-16]imjulia> f!(ArbSeries(degree = 2, prec = 53), ArbSeries((2, 1)))[0.525731112119133 +/- 6.90e-16] + [-0.150384157104163 +/- 2.78e-16]⋅x + [0.032950523196531 +/- 3.27e-16]⋅x^2 + 𝒪(x^3)julia> f!(AcbSeries(degree = 2, prec = 53), AcbSeries((2 + 3im, 1)))([0.277808890086608 +/- 8.76e-16] + [-0.283570391073164 +/- 7.96e-16]im) + ([-0.003614644237407 +/- 6.26e-16] + [0.101930812170967 +/- 5.42e-16]im)⋅x + ([-0.016354294747851 +/- 5.15e-16] + [-0.021459140853368 +/- 3.66e-16]im)⋅x^2 + 𝒪(x^3)
To make this possible there is special handling of the series functions in Arb, see Series methods for more details.