Mutable arithmetic

The high level interface can be combined with the low level wrapper to allow for efficient computations using mutable arithmetic.

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.radrefFunction
radref(x::ArbLike, prec = precision(x))

Return a MagRef referencing the radius of x.

Arblib.midrefFunction
midref(x::ArbLike, prec = precision(x))

Return an ArfRef referencing the midpoint of x.

Arblib.realrefFunction
realref(z::AcbLike, prec = precision(z))

Return an ArbRef referencing the real part of x.

Arblib.imagrefFunction
imagref(z::AcbLike, prec = precision(z))

Return an ArbRef referencing the imaginary part of x.

Arblib.refFunction
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.

!!! Note: For reading coefficients this is always safe, but if the coefficient is mutated then care has to be taken. See the comment further down for how to handle this.

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
    Arblib.set_length!(p, len)
    Arblib.normalise!(p)
    where len is 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.000000000000000000 + 2.000000000000000000⋅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.000000000000000000 +

julia> Arblib.normalise!(p) # Normalising the polynomial updates the degree
1.000000000000000000

julia> Arblib.degree(p) # This is now correct
0

julia> p # And so is printing
1.000000000000000000

Here 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.000000000000000000 + 2.000000000000000000⋅x + 𝒪(x^3)

julia> Arblib.one!(Arblib.ref(q, 2)) # Set the leading coefficient to 1
1.000000000000000000

julia> q # The leading coefficient is not picked up
1.000000000000000000 + 2.000000000000000000⋅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.000000000000000000 + 2.000000000000000000⋅x + 1.000000000000000000⋅x^2 + 𝒪(x^3)

Examples

Compare computing $\sqrt{x^2 + y^2}$ using mutable arithmetic with the default.

julia> using Arblib, BenchmarkTools
julia> x = Arb(1 // 3)[0.33333333333333333333333333333333333333333333333333333333333333333333333333333 +/- 4.78e-78]
julia> y = Arb(1 // 5)[0.20000000000000000000000000000000000000000000000000000000000000000000000000000 +/- 3.89e-78]
julia> res = zero(x)0
julia> 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: 7719 samples with 500 evaluations. Range (minmax): 600.882 ns384.005 μs GC (min … max): 0.00% … 83.55% Time (median): 705.834 ns GC (median): 0.00% Time (mean ± σ): 1.292 μs ± 9.427 μs GC (mean ± σ): 20.81% ± 2.92% ▅▃▃▂▁▂▂█▆▅▄▃▃▃▄▃▂▁▂▁▁▁▁▁▁▁▁ ▁ ███████████████████████████████▇▇▇██████▇▆▇▆▄▄▄▆▅▆▇▇▆▅▆▅▆▆▃ █ 601 ns Histogram: log(frequency) by time 1.34 μs < Memory estimate: 448 bytes, allocs estimate: 8.
julia> @benchmark f!($res, $x, $y) samples=10000 evals=500BenchmarkTools.Trial: 10000 samples with 500 evaluations. Range (minmax): 376.640 ns 1.811 μs GC (min … max): 0.00% … 0.00% Time (median): 378.926 ns GC (median): 0.00% Time (mean ± σ): 421.860 ns ± 89.071 ns GC (mean ± σ): 0.00% ± 0.00% ▂▃▃▃▃▂▂▂▁▁▂▁▂▁▂ ▁ ▁ ████████████████████▇▇▇▇▇▇▇▇▇▆▇▇▆▇▇▆▆▆▅▆▅▆▅▆▅▆▆▆▆▆▆▅▆▆▇▇▆▆ █ 377 ns Histogram: log(frequency) by time 770 ns < Memory estimate: 0 bytes, allocs estimate: 0.

Set the radius of the real part of an Acb.

julia> using Arblib
julia> z = Acb(1, 2)1.0000000000000000000000000000000000000000000000000000000000000000000000000000 + 2.0000000000000000000000000000000000000000000000000000000000000000000000000000im
julia> Arblib.set!(Arblib.radref(Arblib.realref(z)), 1e-10)(922337204 * 2^-63)
julia> z[1.000000000 +/- 1.01e-10] + 2.0000000000000000000000000000000000000000000000000000000000000000000000000000im

Compare a naive implementation of polynomial evaluation using mutable arithmetic with one not using using it.

julia> using Arblib, BenchmarkTools
julia> p = ArbPoly(1:10)1.0000000000000000000000000000000000000000000000000000000000000000000000000000 + 2.0000000000000000000000000000000000000000000000000000000000000000000000000000⋅x + 3.0000000000000000000000000000000000000000000000000000000000000000000000000000⋅x^2 + 4.0000000000000000000000000000000000000000000000000000000000000000000000000000⋅x^3 + 5.0000000000000000000000000000000000000000000000000000000000000000000000000000⋅x^4 + 6.0000000000000000000000000000000000000000000000000000000000000000000000000000⋅x^5 + 7.0000000000000000000000000000000000000000000000000000000000000000000000000000⋅x^6 + 8.0000000000000000000000000000000000000000000000000000000000000000000000000000⋅x^7 + 9.0000000000000000000000000000000000000000000000000000000000000000000000000000⋅x^8 + 10.000000000000000000000000000000000000000000000000000000000000000000000000000⋅x^9
julia> x = Arb(1 // 3)[0.33333333333333333333333333333333333333333333333333333333333333333333333333333 +/- 4.78e-78]
julia> res = zero(x)0
julia> 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. Range (minmax): 3.000 μs 6.987 ms GC (min … max): 0.00% … 82.54% Time (median): 3.161 μs GC (median): 0.00% Time (mean ± σ): 9.222 μs ± 137.947 μs GC (mean ± σ): 32.80% ± 2.26% █▇▄▃▃▃▃▂▂▂▂▄▃▃▃▃▂▁▂ ▁ ▁ ▁ ▁ ▂ ████████████████████████████████████████████▇▆▇▇▆▇▆▄▅▄▅▅▃▄ █ 3 μs Histogram: log(frequency) by time 7.49 μs < Memory estimate: 3.94 KiB, allocs estimate: 70.
julia> @benchmark eval!($res, $p, $x) samples = 10000 evals = 30BenchmarkTools.Trial: 10000 samples with 30 evaluations. Range (minmax): 1.317 μs 1.075 ms GC (min … max): 0.00% … 41.63% Time (median): 1.352 μs GC (median): 0.00% Time (mean ± σ): 1.504 μs ± 10.734 μs GC (mean ± σ): 2.97% ± 0.42% ▅█▅▄▂ ▁ ██████▇▅▅▅▆▅▇▆▇▇▆▆▆▆▅▆▅▅▄▅▄▃▅▄▅▅▄▄▅▄▅▅▄▄▅▄▄▃▄▄▄▂▄▃▄▃▄▅▄▅ █ 1.32 μs Histogram: log(frequency) by time 2.52 μs < Memory estimate: 160 bytes, allocs estimate: 3.
julia> @benchmark $p($x) samples = 10000 evals = 30 # Arb implementation for referenceBenchmarkTools.Trial: 10000 samples with 30 evaluations. Range (minmax): 889.200 ns 1.069 ms GC (min … max): 0.00% … 43.03% Time (median): 918.950 ns GC (median): 0.00% Time (mean ± σ): 1.066 μs ± 10.682 μs GC (mean ± σ): 4.32% ± 0.43% ▄█ ██▇▄▃▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▁ 889 ns Histogram: frequency by time 1.85 μs < Memory estimate: 160 bytes, allocs estimate: 3.