GitHub

Week 6 Lab: IntervalArithmetic.jl and Arblib.jl

In this lab we will look at the two Julia packages IntervalArithmetic.jl and Arblib.jl. IntervalArithmetic.jl implements "regular" interval arithmetic, where the intervals are represented by their lower and upper bounds. Arblib.jl implements ball arithmetic, where the intervals are represented by a midpoint and a radius. We will make use of both of these packages throughout the course, though likely Arblib.jl to a larger extent (mostly because I'm more used to that package).

For the first part of the lab (intro to IntervalArithmetic.jl and Arblib.jl) we will use the Julia REPL and the instructions below. For the second part (computing $\sin$) we will use the lab-6.jl Pluto notebook that you can find in the notebooks directory. The instructions for the second part are also included below, but we will use the ones in the notebook.

Intro to IntervalArithmetic.jl

Let us start by taking a closer look at IntervalArithmetic.jl.

Construction

Intervals are constructed using the interval function. With interval(a, b) we can construct the interval $[a, b]$ and with interval(a) we get the thin interval $[a, a]$.

julia> using IntervalArithmetic
julia> interval(1, 2) # We will get back to what _com means[1.0, 2.0]_com
julia> interval(0.1, 0.2)[0.1, 0.2]_com
julia> interval(2)[2.0, 2.0]_com

By default it will create intervals where the endpoints are of type Float64.

julia> typeof(interval(1))IntervalArithmetic.Interval{Float64}

We can create intervals with endpoints of different types by giving the type as first argument.

julia> interval(BigFloat, 1, 2) # Note the ₂₅₆ in the output, this is the BigFloat precision[1.0, 2.0]₂₅₆_com
julia> typeof(interval(BigFloat, 1, 2))IntervalArithmetic.Interval{BigFloat}
julia> interval(Rational{Int}, 1, 2)[1//1, 2//1]_com
julia> typeof(interval(Rational{Int}, 1, 2))IntervalArithmetic.Interval{Rational{Int64}}

In practice we will mostly use Float64 and sometimes BigFloat. Having endpoints which are rational numbers can sometimes be useful, but we won't see it too much.

Finally, you can create an interval from a string with its decimal representation using parse. This avoid issues with rounding of floating points and guarantees that the given number is contained in the interval.

julia> parse(Interval{Float64}, "0.1") # The printing of this is weird, see next section[0.1, 0.1]_com
julia> parse(Interval{Float64}, "[0.1, 0.2]")[0.1, 0.2]_com

Printing

There are a couple of different options for how intervals are printed. The default can be set with setdisplay, see its documentation for more details. The default is setdisplay(:infsup).

By default it only rounds to 6 significant digits, with no special handling of the different endpoints. This means that if the endpoints share the first 6 significant digits then they print the same.

julia> setdisplay(:infsup)Display options:
  - format: infsup
  - decorations: true
  - NG flag: true
  - significant digits: 6
julia> interval(π)[3.14159, 3.14159]_com
julia> interval(1, 1 + 1e-10)[1.0, 1.0]_com
julia> interval(Float64, 1 // 3) # Just interval(1 // 3) gives us a Rational one[0.333333, 0.333333]_com

As you can see in the above examples it also prints _com after the interval. This is a decoration. They keep track of extra information regarding the functions used to compute the interval. We won't care about these for now, but might come back to them later in the course.

To print everything you can use setdisplay(:full).

julia> setdisplay(:full)Display options:
  - format: full
  - decorations: true (ignored)
  - NG flag: true (ignored)
  - significant digits: 6 (ignored)
julia> interval(1, 2)Interval{Float64}(1.0, 2.0, com, true)
julia> interval(π)Interval{Float64}(3.141592653589793, 3.1415926535897936, com, true)
julia> interval(1, 1 + 1e-10)Interval{Float64}(1.0, 1.0000000001, com, true)
julia> interval(Float64, 1 // 3)Interval{Float64}(0.3333333333333333, 0.33333333333333337, com, true)

In this case the full endpoints are printed. It should however be noted that the printing doesn't take into account rounding. For example the printing of the following interval makes it look like it would contain $0.1 = 1 / 10$.

julia> setdisplay(:full)Display options:
  - format: full
  - decorations: true (ignored)
  - NG flag: true (ignored)
  - significant digits: 6 (ignored)
julia> interval(0.1, 0.2)Interval{Float64}(0.1, 0.2, com, true)

But this is an artefact of $0.1$ not being exactly representable in Float64. We can print more digits by creating an interval of type BigFloat, and then we can see that the interval does in fact not contain $0.1$ (it does contain $0.2$ though).

julia> setdisplay(:full)Display options:
  - format: full
  - decorations: true (ignored)
  - NG flag: true (ignored)
  - significant digits: 6 (ignored)
julia> interval(BigFloat, 0.1, 0.2)Interval{BigFloat}(0.1000000000000000055511151231257827021181583404541015625, 0.200000000000000011102230246251565404236316680908203125, com, true)

You should be careful with reading information of intervals from their printed representation. It is usually better to use a predicate to check something directly

julia> in_interval(1 // 10, interval(0.1, 0.2)) # 1 // 10 is exactly 0.1false

Useful tools

You can get the lower and upper bounds of the interval using inf and sup.

julia> x = interval(1, 2)Interval{Float64}(1.0, 2.0, com, true)
julia> inf(x)1.0
julia> sup(x)2.0

You can get the midpoint and radius using mid and radius, or both using midradius. You can also get the diameter with diam.

julia> x = interval(1, 2)Interval{Float64}(1.0, 2.0, com, true)
julia> mid(x)1.5
julia> radius(x)0.5
julia> midradius(x)(1.5, 0.5)
julia> diam(x)1.0

Note that these functions return floating points and that they therefore round. For radius and diam the rounding is outwards, so they should always be at least as large as the true value. For the midpoint the rounding is to nearest. Since they round you should however be a bit careful with using these.

julia> x = interval(1e-100, 1)Interval{Float64}(1.0e-100, 1.0, com, true)
julia> mid(x)0.5
julia> radius(x)0.5
julia> midradius(x)(0.5, 0.5)
julia> diam(x)1.0

Intro to Arblib.jl

Construction

For Arblib.jl balls are constructed using the Arb constructor. With Arb(x) we get a ball enclosing the value of $x$, with Arb((a, b)) we get a ball enclosing the interval $[a, b]$.

julia> using Arblib
julia> Arb(1)1.0
julia> Arb(1 // 3)[0.33333333333333333333333333333333333333 +/- 3.83e-39]
julia> Arb(π)[3.1415926535897932384626433832795028842 +/- 1.06e-38]
julia> Arb("0.1")[0.10000000000000000000000000000000000000 +/- 1.48e-40]
julia> Arb("[0.1 +/- 1e-10]")[0.100000000 +/- 1.01e-10]
julia> Arb((1, 2)) # Note that this prints like [+/- 2.01], we'll get back to this[+/- 2.01]
julia> Arb((-1, 1))[+/- 1.01]
julia> Arb((1 // 3, π))[+/- 3.15]

An alternative, slightly lower level, constructor is setball. We have that setball(Arb, m, r) creates a ball with midpoint m and radius r.

julia> setball(Arb, 0, 1)[+/- 1.01]
julia> setball(Arb, 1, 1e-10)[1.000000000 +/- 1.01e-10]

With setball the midpoint m is first rounded to a floating point. This for example means that setball(Arb, 1 // 3, 0) will not actually contain the number 1 / 3.

julia> x = setball(Arb, 1 // 3, 0)[0.33333333333333333333333333333333333333 +/- 3.83e-39]
julia> radius(x)0

If rounding for the midpoint needs to be taken into account add_error is more useful, add_error(x::Arb, err) will return the ball x with err added to its radius.

julia> x = Arb(1 // 3)[0.33333333333333333333333333333333333333 +/- 3.83e-39]
julia> y = add_error(x, Arb("1e-50"))[0.33333333333333333333333333333333333333 +/- 3.83e-39]

Printing

Similar to for IntervalArithmetic.jl there are a couple of options when printing Arb values. For Arblib.jl there is however no global setting, instead the string function is used with different arguments. See the documentation of string for all the options.

Arb values are printed on the form [m +/- r] and by default the printed value is guaranteed to enclose the true interval. The output for the midpoint is rounded so that the value is correct up to 1 ulp (unit in the last decimal place). The following example from the documentation of string shows that this behavior in some cases can be slightly confusing

julia> x = Arb((1, 2))[+/- 2.01]
julia> string(x)"[+/- 2.01]"
julia> string(x, more = true)"[1.5000000000000000000000000000000000000 +/- 0.501]"
julia> string(x, digits = 5, more = true)"[1.5000 +/- 0.501]"

Compared to IntervalArithmetic.jl you can hence trust that the true value is contained in the printed output. As the above example shows there are however often large overestimations in the printed value, in particular when no or very few significant digits of the output can be determined.

Useful tools

You can get the midpoint and radius of a ball using midpoint and radius or getball.

julia> x = Arb(1 // 3)[0.33333333333333333333333333333333333333 +/- 3.83e-39]
julia> midpoint(x)0.33333333333333333333333333333333333333
julia> radius(x)1.46936794e-39
julia> getball(x)(0.33333333333333333333333333333333333333, 1.46936794e-39)

By default the midpoint and radius are returned as floating points (Arf for the midpoint, Mag for the radius). You can get them as Arb values (with radius zero)

julia> x = Arb(1 // 3)[0.33333333333333333333333333333333333333 +/- 3.83e-39]
julia> midpoint(Arb, x)[0.33333333333333333333333333333333333333 +/- 2.36e-39]
julia> radius(Arb, x)[1.4693679385278593849609206715278070973e-39 +/- 2.67e-77]
julia> getball(Arb, x)([0.33333333333333333333333333333333333333 +/- 2.36e-39], [1.4693679385278593849609206715278070973e-39 +/- 2.67e-77])

This is often more useful if you want to later do calculations with them.

You can get lower and upper bounds with lbound, ubound and getinterval.

julia> x = Arb(1 // 3)[0.33333333333333333333333333333333333333 +/- 3.83e-39]
julia> lbound(x)0.33333333333333333333333333333333333333
julia> ubound(x)0.33333333333333333333333333333333333333
julia> getinterval(x)(0.33333333333333333333333333333333333333, 0.33333333333333333333333333333333333333)

By default these are returned as floating points, you can get them as Arb values (with radius zero)

julia> x = Arb(1 // 3)[0.33333333333333333333333333333333333333 +/- 3.83e-39]
julia> lbound(Arb, x)[0.33333333333333333333333333333333333333 +/- 8.85e-40]
julia> ubound(Arb, x)[0.33333333333333333333333333333333333333 +/- 3.83e-39]
julia> getinterval(Arb, x)([0.33333333333333333333333333333333333333 +/- 8.85e-40], [0.33333333333333333333333333333333333333 +/- 3.83e-39])

Note that these values are not the exact lower and upper bounds of the interval, they are floating point values given by rounding the true values outwards.

The Arblib.rel_accuracy_bits function can be used to get the relative accuracy of a ball.

julia> x = Arb(1 // 3)[0.33333333333333333333333333333333333333 +/- 3.83e-39]
julia> Arblib.rel_accuracy_bits(x)126
julia> y = add_error(x, Arb(1e-16))[0.333333333333333 +/- 4.34e-16]
julia> Arblib.rel_accuracy_bits(y)51

Computing $\sin$

We are now ready to compute $\sin$ using both IntervalArithmetic.jl and Arblib.jl!

As before we will use a Taylor expansion around zero. Compared to before we will however not use a fixed degree for the expansion, but allow it to be easily changed. Recall that

\[\sin(x) = \sum_{n = 0}^\infty (-1)^n\frac{x^{2n + 1}}{(2n + 1)!}.\]

Since the series is alternating the error we get if we only sum up to term $N$ is bounded by term $N + 1$, i.e.

\[\left|\sin(x) - \sum_{n = 0}^N (-1)^n\frac{x^{2n + 1}}{(2n + 1)!}\right| \leq \frac{|x|^{2N + 3}}{(2N + 3)!}.\]

This bound can be improved, but it is good enough for our purposes. Note that we don't put any restrictions on the value of $x$ here. For large values of $x$ the error bound will of course be very large, but it is still valid.

Here is an implementation of this approach for Float64. Note that it doesn't actually do anything with the computed error bound.

function my_sin(x::Float64; N::Integer = 6)
    y = 0.0
    for n in 0:N
        y += (-1)^n * x^(2n + 1) / factorial(2n + 1)
    end
    # For Float64 we can't actually do anything with this value...
    err = abs(x)^(2N + 3) / factorial(2N + 3)
    return y
end
my_sin (generic function with 1 method)

Let us check that it seems to work

julia> my_sin(1.0)0.8414709848086585
julia> sin(1.0)0.8414709848078965
julia> my_sin(1.0) ≈ sin(1.0)true
julia> my_sin(10.0) # For large values of x we would need more terms548.9651500762614
julia> sin(10.0)-0.5440211108893698

Let us now implement Interval{Float64} and Arb versions of this function! Here we do want to correctly handle the error bounds as well.

function my_sin(x::Interval{Float64}; N::Integer = 6)
    # TASK: Implement this
end

function my_sin(x::Arb; N::Integer = 6)
    # TASK: Implement this
end
Solution
function my_sin(x::Interval{Float64}; N::Integer = 6)
    y = zero(x)
    for n in 0:N
        y += (-1)^n * x^(2n + 1) / factorial(2n + 1)
    end
    err = abs(x)^(2N + 3) / factorial(2N + 3)
    return y + interval(-sup(err), sup(err))
end

function my_sin(x::Arb; N::Integer = 6)
    y = zero(x)
    for n in 0:N
        y += (-1)^n * x^(2n + 1) / factorial(2n + 1)
    end
    err = abs(x)^(2N + 3) / factorial(2N + 3)
    return add_error(y, err)
end
my_sin (generic function with 3 methods)
julia> setdisplay(:full) # To see the full outputDisplay options:
  - format: full
  - decorations: true (ignored)
  - NG flag: true (ignored)
  - significant digits: 6 (ignored)
julia> my_sin(interval(1))Interval{Float64}(0.8414709848078934, 0.8414709848094235, com, false)
julia> sin(interval(1))Interval{Float64}(0.8414709848078965, 0.8414709848078966, com, true)
julia> issubset_interval(sin(interval(1)), my_sin(interval(1)))true
julia> my_sin(Arb(1))[0.84147098481 +/- 2.11e-12]
julia> sin(Arb(1))[0.84147098480789650665250232163029899962 +/- 3.95e-39]
julia> Arblib.contains(my_sin(Arb(1)), sin(Arb(1)))true
julia> my_sin(interval(10)) # For large values of x we get enormous error boundsInterval{Float64}(-215.75122310572146, 1313.6815232582435, com, false)
julia> my_sin(Arb(10))[+/- 1.32e+3]

With these two implementations, some of the things we can look at are:

  • How does the radius of the Interval and Arb versions compare? For example for $x = 1$.
  • What happens if you increase $N$? Hint: At some point you will get an error and need to adjust the code.
  • What happens for wide input? $x = [0.999, 1.001],\ [0.9, 1.1],\ [0, 1],\ [0, 2]$? Or even larger!
  • Can we do better for wide input? There are two reasonable approaches, one better for Interval and one for Arb.
  • For large values of $x$ you first want to reduce the argument to a smaller value using the periodicity. How could this be done?