Precision

Similarly to BigFloat, most types in Arblib.jl support computations with arbitrary precision. However, the way the precision is handled differ slightly from how BigFloat handles it. The main points, which are expanded upon below, are:

1. Precision can be specified when constructing a new value using the prec keyword argument.
2. The default precision that is used if no precision is specified can be changed, either globally or within a dynamical scope, using setprecision.
3. Computations use the precision of the input(s). This is different to BigFloat, which always uses the default precision for all computations.
4. For Arb and Acb it is safe to mix values computed with different precisions, the error bounds keep track of the actual precision of the computations. This is not true for Arf and Acf (nor is it true for BigFloat), where mixing precisions can make the final result appear to be computed at a higher precision than it actually is.

Specifying precision in constructor

Precision can be specified when constructing a new value using the prec keyword argument.

julia> x = Arb(1 // 3, prec = 64)[0.3333333333333333333 +/- 4.24e-20]
julia> precision(x)64
julia> Arb(1 // 3, prec = 128)[0.33333333333333333333333333333333333333 +/- 3.83e-39]
julia> AcbSeries((π, 1, 1 // 3), prec = 53)[3.14159265358979 +/- 3.34e-15] + 1.0⋅x + [0.333333333333333 +/- 3.71e-16]⋅x^2 + 𝒪(x^3)
julia> ArbVector([1, 2, 3], prec = 80)3-element ArbVector:
 1.0
 2.0
 3.0

Using setprecision

Similarly to BigFloat there is a default precision that is used if no precision is specified. The default precision can be changed, either globally or within a dynamical scope, using setprecision.

setprecision(::Type{<:ArbTypes}, precision::Integer; base::Integer = 2)
setprecision(f::Function, ::Type{<:ArbTypes}, precision::Integer; base=2)

old = precision(Arb)
setprecision(Arb, precision)
f()
setprecision(Arb, old)

julia> x = Arb(0) # Default precision is 256 bits
0

julia> acos(x)
[1.5707963267948966192313216916397514420985846996875529104874722961539082031431 +/- 9.63e-78]

julia> setprecision(Arb, 128) do
           acos(x) # Still computed using 256 bits since that is the precision of x
       end
[1.5707963267948966192313216916397514420985846996875529104874722961539082031431 +/- 9.63e-78]

julia> y = BigFloat(0) # Default precision is 256 bits
0.0

julia> acos(y)
1.570796326794896619231321691639751442098584699687552910487472296153908203143099

julia> setprecision(BigFloat, 128) do
           acos(y)
       end
1.570796326794896619231321691639751442098

setprecision(x::ArbTypes, precision::Int; base::Integer = 2)

Handling of precision in computations

Contrary to BigFloat, most computations use the precision of the arguments. For functions taking multiple arguments the maximum precision of the arguments is used.

julia> x = Arb(1 // 3, prec = 64)[0.3333333333333333333 +/- 4.24e-20]
julia> y = sin(x) # Result uses 64 bits of precision[0.3271946967961522441 +/- 8.51e-20]
julia> precision(y)64
julia> a = Arb(π, prec = 128)[3.1415926535897932384626433832795028842 +/- 1.06e-38]
julia> b = Arb(1 // 3, prec = 64)[0.3333333333333333333 +/- 4.24e-20]
julia> c = a + b # Result uses 128 bits of precision[3.4749259869231265718 +/- 4.92e-20]
julia> precision(c)128

This makes it easy to use either a lower or higher precision for a part of the computation. It can be useful to use lower precision when evaluating terms that are significantly smaller in magnitude. Higher precision can be helpful when evaluating close to a removable singularity. Consider the artificial example of computing $e^x + \sin(x)$ for $x$ near zero, the $\sin(x)$ term can then be computed with significantly lower precision than the $e^x$ term.

julia> x = Arb(1e-70) # Input close to zero[9.9999999999999999566603349622936327208651652641680501722443601799944492674166e-71 +/- 1.13e-148]
julia> exp(x) + sin(x) # Compute everything in full precision[1.0000000000000000000000000000000000000000000000000000000000000000000002000000 +/- 5.06e-77]
julia> exp(x) + sin(setprecision(x, 64)) # Compute sin term in lower precision[1.0000000000000000000000000000000000000000000000000000000000000000000002000000 +/- 5.06e-77]

Common pitfalls

While the library does its best to preserve precision through computations, there are still ways to accidentally cause some operations to use the default precision instead.

The most common situation is methods that only see the type of a value. Since the type doesn't carry any information about the precision these methods fall back to the default precision. This for example applies to zero(T) and convert(T, x) (and also to oftype(y, x) which is equivalent to convert(typeof(y), x)). For example, if we are interested in computing the function $f(x) = x + 2\pi$, and we implement it using 2oftype(x, π) to compute $2\pi$, that part of the computation will use the currently set default precision and not the precision of x.

julia> f(x) = x + 2oftype(x, π)f (generic function with 1 method)
julia> y = f(Arb(1, prec = 53))[7.2831853071795864769252867665590057683943387987502116419498891846156328125724 +/- 3.85e-77]
julia> precision(y) # Results uses 256 bits of precision instead of 53256
julia> z = f(Arb(1, prec = 512))[7.2831853071795864769252867665590057683943387987502116419498891846156328125724 +/- 3.85e-77]
julia> precision(z) # Result uses 512 bits of precision as wanted512
julia> Arblib.rel_accuracy_bits(z) # But error bounds are only accurate to 256 bits255

In the above case one solution would be to define f(x::Arb) = x + 2Arb(π, prec = precision(x)). In general there is however no generic solution to this problem that doesn't require implementing a special case for f(x::Arb).