This manual is for the GNU Octave interval package, version 3.2.2.

Copyright © 2015 Oliver Heimlich
Copyright © 2017 Joel Dahne
Copyright © 2008 Simone Pernice

Permission is granted to copy, distribute and/or modify this document under the terms of the GNU General Public License, Version 3 or any later version published by the Free Software Foundation. A copy of the license is included in GNU General Public License.

## Table of Contents

## Preface ¶

Welcome to the user manual of the interval package for GNU Octave. This chapter presents background information and may safely be skipped. First-time users who want to cut right to the chase should read Getting Started, which teaches basic concepts and first steps with the package. Users who are not familiar with interval arithmetic should read Introduction to Interval Arithmetic first. Still feeling undecided? Look at the Examples and see how easy you can put this software to great use!

Development of the GNU Octave Interval Package started in September 2014. The IEEE standard for interval arithmetic, IEEE Std 1788-2015, had been drafted by its working group until July 2014 and was about to enter the balloting process. In January 2015 a first package release could be made, which contained the full set of functions required by the standard’s draft. On June 11 the standard finally became approved and this interval package can be seen as the first ever completed standard conforming interval arithmetic library.

The creation of the interval package has been straightforward, although the author had no previous experience with Octave. Octave is a great environment for getting things done and its active community helps a lot. In this spirit, the interval package wants to be an easy to use tool for experimenting with and quick prototyping of interval arithmetic algorithms and applications.

Originally it was intended to only implement the operations required by the standard document, but support for fundamental concepts of Octave as well as interval vectors and interval matrices have soon been added. Today the package contains many useful interval analysis algorithms and solvers, which, together with basic arithmetic functions, form a powerful and versatile library.

In 2017 the Google Summer of Code (GSoC) program has lead to support for N-dimensional arrays of intervals.

Like Octave, the interval package has nothing to do with music. Above picture is a pun and shows an “octave interval” between the notes d’ and d”. The frequencies of these notes can be enclosed by an interval which is a subset of [293, 588] Hz.

### Acknowledgments ¶

The GNU Octave interval package is build upon great third-party software.

- Most correctly rounded arithmetic operations are based on the GNU MPFR library by Guillaume Hanrot, Vincent Lefèvre, Patrick Pélissier, Philippe Théveny and Paul Zimmermann.
- Several correctly rounded arithmetic operations are based on the correctly rounded math library by Jean-Michel Muller, Florent de Dinechin, Christoph Lauter, David Defour, Catherine Daramy-Loirat, Matthieu Gallet, and Nicolas Gast.
- Jiří Rohn has published his comprehensive verification toolbox VERSOFT as free software on July 26, 2016. On November 24 he has fully disclosed the source code and several high-level functions could be included in the interval package with minor adjustments, e. g., @infsup/chol, vereigvec, and verlinprog. Also, many thanks to Vladik Kreinovich who has helped to clear the licensing issue with VERSOFT.
- A linear system solver @infsup/mldivide and a polynomial evaluation algorithm @infsup/polyval for bare intervals are derived from routines developed for C-XSC at University of Wuppertal, Germany.
- Introduction to Interval Arithmetic is partly based on the documentation for the former SIMP package for Octave by Simone Pernice.
- A French translation of the package description has been made by Rodéric Moitié.
- In the Examples for finding root enclosures a function and code by Helmut Podhaisky has been used.
- The online course on interval analysis by Luc Jaulin and Jordan Ninin at ENSTA-Bretagne has inspired me to implement the set inversion algorithms in @infsup/fsolve.
- I have gained access to scientific literature thanks to the Weihenstephan-Triesdorf University of Applied Sciences.
- Most unit tests are written in portable ITL format and converted into GNU Octave test cases with the Interval Testing Framework for IEEE 1788. The framework has been developed by Maximilian Kiesner and Marco Nehmeier, Chair of Software Engineering, Department of Computer Science, University of Würzburg, Germany.
- Several unit tests are derived from libieeep1788, a C++ implementation of IEEE Std 1788-2015, IEEE standard for interval arithmetic. The library contains several unit tests for its IEEE 754 flavor, which is compatible with the arithmetic of the GNU Octave interval package. I have converted nearly 6000 of these test cases into portable ITL format for verification of this package.

  >     libieeep1788
  >     ============

  >     Copyright 2013 - 2015

  >     Marco Nehmeier (nehmeier@informatik.uni-wuerzburg.de)
  >     Department of Computer Science,
  >     University of Wuerzburg, Germany

  >     This product includes software developed at
  >     Chair of Software Engineering,
  >     Department of Computer Science,
  >     University of Wuerzburg, Germany
  >     http://se.informatik.uni-wuerzburg.de/

- Several unit tests are derived from MPFI, a C++ interval arithmetic library based on GNU MPFR. The library contains several unit tests for binary64 numbers, which are compatible with the arithmetic of the GNU Octave interval package. I have converted nearly 1500 of these test cases into portable ITL format for verification of this package.
- Several unit tests are derived from FI_LIB, an ANSI-C interval arithmetic library based on binary64 numbers. The library contains several unit tests, which are compatible with the arithmetic of the GNU Octave interval package. I have converted 800 of these test cases into portable ITL format for verification of this package.
- Some unit tests are derived from C-XSC, a C++ class library for interval arithmetic. The library contains some unit tests, which are compatible with the arithmetic of the GNU Octave interval package. I have converted 160 of these test cases into portable ITL format for verification of this package.
- Fast matrix multiplication (see @infsup/mtimes) as well as the linear system solver (see @infsup/mldivide) use BLAS routines with directed rounding. An OCT-file interface for setting the rounding mode has been developed by Kai Torben Ohlus, Institute for Reliable Computing, Hamburg University of Technology, Germany.

Last, but not least, many thanks to everybody who has contributed to the success of free software!

### Philosophy ¶

Features

- Free software licensed under the terms of the GNU General Public License (Version 3 or later)
- Many interval arithmetic functions with high, system-independent accuracy
- Conforming to IEEE Std 1788-2015, IEEE standard for interval arithmetic
- Support for interval vectors, interval matrices and N-dimensional interval arrays
  - very accurate vector sum, vector dot and matrix multiplication (correctly rounded)
  - fast matrix multiplication and fast solver for dense linear systems (BLAS routines)
  - vectorized function evaluation
- Easy usage
  - GNU Octave function names
  - convenient interval constructors
  - broadcasting

Limitations

- No complex numbers
- No sparse matrices (maybe in the future, if requested by users)

The interval arithmetic provided by the interval package focuses on easy usage, accuracy and correctness. It is rather slow compared to other arithmetic libraries.

If accurate type checking during compile time—a substantial feature for verified computing—is needed, the user is advised to try third-party interval libraries for strongly typed programming languages like C/C++. The interval package for GNU Octave can nonetheless be used for prototyping of interval algorithms.

Why is the interval package slow? All arithmetic interval operations are simulated in high-level Octave language using C99 or multi-precision floating-point routines, which is a lot slower than a hardware implementation. Building interval arithmetic operations from floating-point routines is easy for simple monotonic functions, e. g., addition and subtraction, but is complex for others, e. g., interval power function, atan2, or reverse functions.

For some interval operations it is not even possible to rely on floating-point routines, since not all required routines are available in C99 or BLAS. For example, accurate multiplication of matrices with many elements becomes unfeasible as it takes a lot of time.

| `plus` | `log` | `pow` | `mtimes` | `mtimes` | `inv`
--- | --- | --- | --- | --- | --- | ---
Interval matrix size | tightest accuracy | tightest accuracy | tightest accuracy | valid accuracy | tightest accuracy | valid accuracy
10 × 10 | < 0.001 | 0.001 | 0.008 | 0.001 | 0.002 | 0.025
100 × 100 | 0.003 | 0.055 | 0.61 | 0.012 | 0.53 | 0.30
500 × 500 | 0.060 | 1.3 | 15 | 0.30 | 63 | 4.2

Table 1: Approximate runtime for certain functions (wall clock time in seconds) — Results have been produced with GNU Octave 3.8.2 and Interval package 0.1.4 on an Intel Core i5-4340M CPU (2.9–3.6 GHz)

Why is the interval package accurate? The GNU Octave built-in floating-point routines are not useful for interval arithmetic: Their results depend on hardware, system libraries and compilation options. The interval package handles all arithmetic functions with the help of the GNU MPFR library. With MPFR it is possible to compute system-independent, valid and tight enclosures of the correct results for most functions. However, it should be noted that some reverse operations and matrix operations do not exists in GNU MPFR and therefore cannot be computed with the same accuracy.

It is possible to use faster (BLAS based) routines during computation of the matrix multiplication @infsup/mtimes, because correctly rounded matrix multiplication could be considered too slow for certain applications. However, this is not the default behavior and must be explicitly activated by the user.

### Distribution and Development ¶

The interval package is free software: Everyone is encouraged to use it, copy it and redistribute it, as well as to make changes under the terms of the GNU General Public License.

The interval package is part of Octave Forge, a sibling of the GNU Octave project. Official releases are published at https://octave.sourceforge.io/.

The source code repository is located at Octave Forge and contains the latest development version. Information for developers can be found on the package’s page at Octave wiki.

Bug reports and feature requests for either the software or this manual may be posted under the Octave Project at Savannah.

Contributions to the software and this manual are highly appreciated.

### Getting Help ¶

The interval package contains online help for every function, which can be accessed with the `help` command from Octave. The interval arithmetic is implemented with specialized data types, which override standard functions. Whilst the command `help function name` shows the documentation for core Octave functions, the interval variants of these functions can be requested with the command `help @infsup/function name`.

Further help can be seeked at the Octave Help mailing list.

## 1 Getting Started ¶

This chapter takes you by the hand and gives a quick overview on the interval packages basic capabilities. More detailed information can usually be found in the functions’ documentation.

### 1.1 Installation ¶

It is recommended to install the package from specialized distributors for the particular platform, e. g., Debian GNU/Linux, MacPorts (for Mac OS X), FreshPorts (for FreeBSD), and so on. Since Octave version 4.0.1 the package is included in the official installer for Microsoft Windows and is installed automatically on that platform.

In any case, the interval package can alternatively be installed with the `pkg` command from within Octave. Latest release versions are published at Octave Forge and can be automatically downloaded with the `-forge` option.

    pkg install -forge interval
      -| For information about changes from previous versions
      -| of the interval package, run 'news interval'.

During this kind of installation parts of the interval package are compiled for the target system, which requires development libraries for GNU Octave (version ≥ 4.2.0) and GNU MPFR (version ≥ 3.1.0) to be installed. It might be necessary to install packages “liboctave-dev” and “libmpfr-dev”, which are provided by most GNU distributions (names may vary).

In order to use the interval package during an Octave session, it must have been loaded, i. e., added to the path. In the following parts of the manual it is assumed that the package has been loaded, which can be accomplished with the `pkg load interval` command. It is recommended to add this command at the beginning of script files, especially if script files are published or shared. Automatic loading of the interval package can be activated by adding the line `pkg load interval` to your `.octaverc` file located in your user folder, for more information see Startup Files in GNU Octave manual.

That’s it. The package is ready to be used within Octave.

### 1.2 Set-based Interval Arithmetic ¶

The most important and fundamental concepts in the context of the interval package are:

- Intervals are closed, connected subsets of the real numbers. Intervals may be unbound (in either or both directions) or empty. In special cases `+inf` and `-inf` are used to denote boundaries of unbound intervals, but any member of the interval is a finite real number.
- Classical functions are extended to interval functions as follows: The result of function f evaluated on interval x is an interval enclosure of all possible values of f over x where the function is defined. Most interval arithmetic functions in this package manage to produce a very accurate such enclosure.
- The result of an interval arithmetic function is an interval in general. It might happen, that the mathematical range of a function consist of several intervals, but their union will be returned, e. g., 1 / [-1, 1] = [Entire].

More details can be found in Introduction to Interval Arithmetic.

### 1.3 Input and Output ¶

Before exercising interval arithmetic, interval objects must be created from non-interval data. There are interval constants empty and entire and the interval constructors @infsupdec/infsupdec (create an interval from boundaries), midrad (create an interval from midpoint and radius) and hull (create an interval enclosure for a list of mixed arguments: numbers, intervals or interval literals). The constructors are very sophisticated and can be used with several kinds of parameters: Interval boundaries can be given by numeric values or string values with decimal numbers.

Create intervals for performing interval arithmetic

    ## Interval with a single number
    infsupdec (1)
      ⇒ ans = [1]_com

    ## Interval defined by lower and upper bound
    infsupdec (1, 2)
      ⇒ ans = [1, 2]_com

    ## Boundaries are converted from strings
    infsupdec ("3", "4")
      ⇒ ans = [3, 4]_com

    ## Decimal number
    infsupdec ("1.1")
      ⇒ ans ⊂ [1.0999, 1.1001]_com

    ## Decimal number with scientific notation
    infsupdec ("5.8e-17")
      ⇒ ans ⊂ [5.7999e-17, 5.8001e-17]_com

    ## Interval around 12 with uncertainty of 3
    midrad (12, 3)
      ⇒ ans = [9, 15]_com

    ## Again with decimal numbers
    midrad ("4.2", "1e-3")
      ⇒ ans ⊂ [4.1989, 4.2011]_com

    ## Interval members with arbitrary order
    hull (3, 42, "19.3", "-2.3")
      ⇒ ans ⊂ [-2.3001, +42]_com

    ## Symbolic numbers
    hull ("pi", "e")
      ⇒ ans ⊂ [2.7182, 3.1416]_com

Warning: In above examples decimal fractions are passed as a string to the constructor. Otherwise it is possible, that GNU Octave introduces conversion errors when the numeric literal is converted into a floating-point number before it is passed to the constructor. The interval construction is a critical process, but after this the interval package takes care of any further conversion errors, representational errors, round-off errors and inaccurate numeric functions.

Beware of the conversion pitfall

    ## The numeric constant 0.3 is an approximation of the
    ## decimal number 0.3.  An interval around this approximation
    ## will not contain the decimal number 0.3.
    sprintf("%[.17g]", infsupdec (0.3))
      ⇒ ans = [0.29999999999999998, 0.29999999999999999]_com

    ## However, passing the decimal number 0.3 as a string
    ## to the interval constructor will create an interval which
    ## actually encloses the decimal number.
    format short
    infsupdec ("0.3")
      ⇒ ans ⊂ [0.29999, 0.30001]_com

For maximum portability it is recommended to use interval literals, which are standardized by IEEE Std 1788-2015. Both interval boundaries are then given as a string in the form `[l, u]`. The output in the examples above gives examples of several interval literals.

    ## Interval literal
    infsupdec ("[20, 4.2e10]")
      ⇒ ans = [20, 4.2e+10]_com

The default text representation of intervals is not guaranteed to be exact, because this would massively spam console output. For example, the exact text representation of `realmin` would be over 700 decimal places long! However, the output is correct as it guarantees to contain the actual boundaries: a displayed lower (upper) boundary is always less (greater) than or equal to the actual boundary.

#### 1.3.1 Interval Vectors, Matrices and Arrays ¶

Vectors, matrices and arrays of intervals can be created by passing numerical arrays, string or cell arrays to the interval constructors. With cell arrays it is also possible to mix several types of boundaries.

Interval arrays behave like normal arrays in GNU Octave and can be used for broadcasting and vectorized function evaluation. Vectorized function evaluation usually is the key to create very fast programs.

Create interval arrays

    M = infsupdec (magic (3))
      ⇒ M = 3×3 interval matrix

        [8]_com   [1]_com   [6]_com
        [3]_com   [5]_com   [7]_com
        [4]_com   [9]_com   [2]_com

    infsupdec (magic (3), magic (3) + 1)
      ⇒ ans = 3×3 interval matrix

        [8, 9]_com    [1, 2]_com   [6, 7]_com
        [3, 4]_com    [5, 6]_com   [7, 8]_com
        [4, 5]_com   [9, 10]_com   [2, 3]_com

    infsupdec ("0.1; 0.2; 0.3; 0.4; 0.5")
      ⇒ ans ⊂ 5×1 interval vector

        [0.099999, 0.10001]_com
         [0.19999, 0.20001]_com
         [0.29999, 0.30001]_com
         [0.39999, 0.40001]_com
                      [0.5]_com

    infsupdec ("1 [2, 3]; 4, 5, 6")
      ⇒ ans = 2×3 interval matrix

        [1]_com   [2, 3]_com   [Empty]_trv
        [4]_com      [5]_com       [6]_com

    infsupdec ({1; eps; "4/7"; "pi"}, {2; 1; "e"; "0xff"})
      ⇒ ans ⊂ 4×1 interval vector

                  [1, 2]_com
         [2.2204e-16, 1]_com
       [0.57142, 2.7183]_com
           [3.1415, 255]_com

    infsupdec (ones (2, 2, 2))
      ⇒ ans = 2×2×2 interval array

      ans(:,:,1) =

         [1]_com   [1]_com
         [1]_com   [1]_com

      ans(:,:,2) =

         [1]_com   [1]_com
         [1]_com   [1]_com

Strings can easily be used to create vectors and matrices of intervals, `,` and `;` are used to denote the next element in the row or a new row. Octave does however not have way of representing arrays with three or more dimensions using strings in the same way. Therefore you can only create such arrays by passing numerical matrices or cells to the constructor. Alternatively you can build it up in steps.

    A = infsupdec ("1 [2, 3]; 4, 5; 6, 7")
      ⇒ A = 3×2 interval matrix

         [1]_com   [2, 3]_com
         [4]_com      [5]_com
         [6]_com      [7]_com

    A(:,:,2) = infsupdec ("0.1, 0.2; 0.3, 0.4; 0.5, 0.6")
      ⇒ A ⊂ 3×2×2 interval array

      ans(:,:,1) =

         [1]_com   [2, 3]_com
         [4]_com      [5]_com
         [6]_com      [7]_com

      ans(:,:,2) =

         [0.099999, 0.10001]_com   [0.19999, 0.20001]_com
          [0.29999, 0.30001]_com   [0.39999, 0.40001]_com
                       [0.5]_com   [0.59999, 0.60001]_com

### 1.4 Arithmetic Operations ¶

The interval package comprises many interval arithmetic operations. A complete list can be found in its function reference. Function names match GNU Octave standard functions where applicable and follow recommendations by IEEE Std 1788-2015 otherwise, see Function Names.

The interval arithmetic flavor used by this package is the “set-based” interval arithmetic and follows these rules: Intervals are sets. They are subsets of the set of real numbers. The interval version of an elementary function such as sin(x) is essentially the natural extension to sets of the corresponding point-wise function on real numbers. That is, the function is evaluated for each number in the interval where the function is defined and the result must be an enclosure of all possible values that may occur.

By default arithmetic functions are computed with best possible accuracy (which is more than what is guaranteed by GNU Octave core functions). The result will therefore be a tight and very accurate enclosure of the true mathematical value in most cases. Details on each function’s accuracy can be found in its documentation, which is accessible with GNU Octave’s `help` command.

Examples of using interval arithmetic functions

    sin (infsupdec (0.5))
      ⇒ ans ⊂ [0.47942, 0.47943]_com

    power (infsupdec (2), infsupdec (3, 4))
      ⇒ ans = [8, 16]_com

    atan2 (infsupdec (1), infsupdec (1))
      ⇒ ans ⊂ [0.78539, 0.7854]_com

    midrad (magic (3), 0.5) * pascal (3)
      ⇒ ans = 3×3 interval matrix

        [13.5, 16.5]_com   [25, 31]_com   [42, 52]_com
        [13.5, 16.5]_com   [31, 37]_com   [55, 65]_com
        [13.5, 16.5]_com   [25, 31]_com   [38, 48]_com

### 1.5 Numerical Operations ¶

Some interval functions do not return an interval enclosure, but a single number (in binary64 precision). Most important are @infsup/inf and @infsup/sup, which return the lower and upper interval boundaries.

More such operations are @infsup/mid (approximation of the interval’s midpoint), @infsup/wid (approximation of the interval’s width), @infsup/rad (approximation of the interval’s radius), @infsup/mag (interval’s magnitude) and @infsup/mig (interval’s mignitude).

    ## Enclosure of the decimal number 0.1 is not exact
    ## and results in an interval with a small uncertainty.
    wid (infsupdec ("0.1"))
      ⇒ ans =    1.3878e-17

### 1.6 Boolean Operations ¶

Interval comparison operations produce boolean results. While some comparisons are especially for intervals (@infsup/subset, @infsup/interior, @infsup/ismember, @infsup/isempty, @infsup/disjoint, …) others are interval extensions of simple numerical comparison. For example, the less-or-equal comparison is mathematically defined as ∀a ∃b a ≤ b ∧ ∀b ∃a a ≤ b.

    infsup (1, 3) <= infsup (2, 4)
      ⇒ ans = 1

### 1.7 Matrix and Array Operations ¶

Above mentioned operations can also be applied element-wise to interval vectors, matrices or arrays. Many operations use vectorization techniques.

In addition, there are operations on interval matrices and arrays. These operations comprise: dot product, matrix multiplication, vector sums (all with tightest accuracy), matrix inversion, matrix powers, and solving linear systems (the latter are less accurate). As a result of missing hardware / low-level library support and missing optimizations, these operations are relatively slow compared to familiar operations in floating-point arithmetic.

Examples of using interval matrix functions

    A = infsupdec ([1, 2, 3; 4, 0, 0; 0, 0, 1]);
    A (2, 3) = "[0, 6]"
      ⇒ A = 3×3 interval matrix

        [1]_com   [2]_com      [3]_com
        [4]_com   [0]_com   [0, 6]_com
        [0]_com   [0]_com      [1]_com

    B = inv (A)
      ⇒ B = 3×3 interval matrix

          [0]_trv     [0.25]_trv       [-1.5, 0]_trv
        [0.5]_trv   [-0.125]_trv   [-1.5, -0.75]_trv
          [0]_trv        [0]_trv             [1]_trv

    A * B
      ⇒ ans = 3×3 interval matrix

        [1]_trv   [0]_trv   [-1.5, +1.5]_trv
        [0]_trv   [1]_trv       [-6, +6]_trv
        [0]_trv   [0]_trv            [1]_trv

    A = infsupdec (magic (3))
      ⇒ A = 3×3 interval matrix

        [8]_com   [1]_com   [6]_com
        [3]_com   [5]_com   [7]_com
        [4]_com   [9]_com   [2]_com

    c = A \ [3; 4; 5]
      ⇒ c ⊂ 3×1 interval vector

        [0.18333, 0.18334]_trv
        [0.43333, 0.43334]_trv
        [0.18333, 0.18334]_trv

    A * c
      ⇒ ans ⊂ 3×1 interval vector

        [2.9999, 3.0001]_trv
        [3.9999, 4.0001]_trv
        [4.9999, 5.0001]_trv

#### 1.7.1 Notes on Linear Systems ¶

A system of linear equations in the form Ax = b with intervals can be seen as a range of classical linear systems, which can be solved simultaneously. Whereas classical algorithms compute an approximation for a single solution of a single linear system, interval algorithms compute an enclosure for all possible solutions of (possibly several) linear systems. Some characteristics should definitely be known when linear interval systems are solved:

- If the linear system is underdetermined and has infinitely many solutions, the interval solution will be unbound in at least one of its coordinates. Contrariwise, from an unbound result it can not be concluded whether the linear system is underdetermined or has solutions.
- If the interval result is empty in at least one of its coordinates, the linear system is guaranteed to be overdetermined and has no solutions. Contrariwise, from a non-empty result it can not be concluded whether all or some of the systems have solutions or not.
- Wide intervals within the matrix A can easily lead to a superposition of cases, where the rank of A is no longer unique. If the linear interval system contains cases of linear independent equations as well as linear dependent equations, the resulting enclosure of solutions will inevitably be very broad.

However, solving linear systems with interval arithmetic can produce useful results in many cases and automatically carries a guarantee for error boundaries. Additionally, it can give better information than the floating-point variants for some cases.

Standard floating point arithmetic versus interval arithmetic on ill-conditioned linear systems

    A = [1, 0; 2, 0];
    ## This linear system has no solutions
    A \ [3; 0]
      -| warning: ...matrix singular to machine precision...
      ⇒ ans =

           0.6000
                0

    ## This linear system has many solutions
    A \ [4; 8]
      ⇒ ans =

           4
           0

    ## The empty interval vector proves that there is no solution
    infsup (A) \ [3; 0]
      ⇒ ans = 2×1 interval vector

           [Empty]
           [Empty]

    ## The unbound interval vector indicates
    ## that there may be many solutions
    infsup (A) \ [4; 8]
      ⇒ ans = 2×1 interval vector

                [4]
           [Entire]

### 1.8 Plotting ¶

Plotting of intervals in 2D and 3D can be achieved with the functions @infsup/plot and @infsup/plot3 respectively. However, some differences in comparison with classical plotting in Octave shall be noted.

When plotting classical (non-interval) functions in Octave, one normally uses a vector and evaluates a function on that vector element-wise. The resulting X and Y (and possibly Z) coordinates are then drawn against each other, whilst coordinates can be connected using interpolated lines. The plot shows an approximation of the function’s graph and the accuracy (and smoothness of the graph) primarily depends on the number of coordinates where the function has been evaluated.

Evaluating the same function on a single interval (e. g. the part of the function’s domain that is of interest) yields a single interval result which covers the actual range of the function. Plotting just two intervals, input and output, against each other is boring, because the plot would only show a single rectangle. Contrariwise, evaluating the function for many individual points (e. g. using @infsup/linspace) would hardly fit in the philosophy of interval arithmetic. Individual points of evaluation are not interconnected by the interval plotting functions, because that would introduce errors.

The solution for plotting functions with interval arithmetic is called: “mincing”. The @infsup/mince function divides an interval into many smaller adjacent subsets, which can be used for range evaluations of the function. As a result, one gets vectors of intervals, which produce a coverage of the function’s graph using rectangles. Please note, how the rectangles cover the sine function’s true range from −1 to 1 in the following example, whilst the interpolated lines make a poor approximation.

    hold on
    blue = [38 139 210] ./ 255;
    shade = [238 232 213] ./ 255;

    ## Interval plotting
    x = mince (2*infsup (0, "pi"), 6);
    plot (x, sin (x), shade)

    ## Classical plotting
    x = linspace (0, 2*pi, 7);
    plot (x, sin (x), 'linewidth', 2, 'color', blue)

    set (gca, 'XTick', 0 : pi : 2*pi)
    set (gca, 'XTickLabel', {'0', 'pi', '2 pi'})

For 3D plotting the Octave meshgrid function, as usual, becomes handy. The following example shows how two different ranges for X and Y coordinates are used to construct a grid, where the function @infsup/atan2 is evaluated. In this particular case the interval grid has gaps, because X and Y coordinates have been constructed such that intervals do not intersect.

    red = [220 50 47] ./ 255;
    shade = [238 232 213] ./ 255;

    x = midrad (1 : 6, 0.25);
    y = midrad (-3 : 3, 0.25);
    [x, y] = meshgrid (x, y);
    z = atan2 (y, x);
    plot3 (x, y, z, shade, red)

    view ([-35, 30])
    box off
    set (gca, "xgrid", "on", "ygrid", "on", "zgrid", "on")

## 2 Introduction to Interval Arithmetic ¶

> Give a digital computer a problem in arithmetic, and it will grind away methodically, tirelessly, at gigahertz speed, until ultimately it produces the wrong answer. … An interval computation yields a pair of numbers, an upper and a lower bound, which are guaranteed to enclose the exact answer. Maybe you still don’t know the truth, but at least you know how much you don’t know.

— Brian Hayes, DOI: 10.1511/2003.6.484

Interval arithmetic adds two unique features to ordinary computer arithmetic: (1) Functions can be evaluated over (connected) subsets of their domain, and (2) any computational errors are automatically considered and are accumulated in the final outcome. In conjunction they yield a verified result enclosure over a range of input values.

These possibilities of interval arithmetic enable great new possibilities, but what is wrong with the well-known computer arithmetic in the first place?

### 2.1 Motivation ¶

Floating-point arithmetic, as specified by IEEE Std 754, is available in almost every computer system today. It is wide-spread, implemented in common hardware and integral part in programming languages. For example, the binary64 format (a.k.a. double-precision) is the default numeric data type in GNU Octave. Benefits are obvious: The results of arithmetic operations are (mostly) well-defined and comparable between different systems and computation is highly efficient.

However, there are some downsides of floating-point arithmetic in practice, which will eventually produce errors in computations. Generally speaking, most of these problems occur in any arithmetic with finite precision.

- Floating-point arithmetic is often used mindlessly by developers. http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html http://www.cs.berkeley.edu/~wkahan/Mindless.pdf http://www.ima.umn.edu/~arnold/disasters/
- The binary (base-2) data types categorically are not suitable for doing financial computations. Very often representational errors are introduced when using “real world” decimal numbers. http://en.wikipedia.org/wiki/Decimal_computer
- Even if the developer would be proficient, most developing environments / technologies limit floating-point arithmetic capabilities to a very limited subset of IEEE Std 754: Only one or two data types, no rounding modes, missing functions, … http://www.cs.berkeley.edu/~wkahan/JAVAhurt.pdf
- Results are hardly predictable. All operations produce the best possible accuracy at runtime, this is how a floating point works. Contrariwise, financial computer systems typically use a fixed-point arithmetic (COBOL, PL/I, …), where overflow and rounding can be precisely predicted at compile-time.
- Results are system dependent. All but the most basic floating-point operations are not guaranteed to be accurate and produce different results depending on low level libraries and hardware.
- If you do not know the technical details (cf. first bullet) you ignore the fact that the computer lies to you in many situations. For example, when looking at numerical output and the computer says “`ans = 0.1`,” this is not absolutely correct. In fact, the value is only close enough to the value 0.1. Additionally, many functions produce limit values (∞ × −∞ = −∞, ∞ ÷ 0 = ∞, ∞ ÷ −0 = −∞, log (0) = −∞), which is sometimes (but not always!) useful when overflow and underflow occur.

Interval arithmetic addresses above problems in its very special way. It accepts the fact that numbers cannot be stored or computed with infinite precision and introduces enclosures of exact values, which can be computed on any machine with finite precision.

This introduces new possibilities for algorithms. Any errors are covered by the range of an interval during the course of computation. All members of intervals are by definition finite real numbers, which results in an exception free and mathematically well-defined arithmetic. The possibility to actually evaluate a function on a connected range of values and compute a guaranteed enclosure of all possible values is a unique selling point.

For example, the interval newton method (see @infsup/fzero) is able to find all zeros of a particular function. More precisely, the algorithm is able to reliably eliminate ranges of values where the function cannot have a root by a simple interval evaluation of either the function itself or its derivative. Global convergence can be achieved by bisecting the intermediate ranges.

### 2.2 Error bounds in real life ¶

Intervals can be used instead of simple numbers to automatically take into account the tolerance (or uncertainty) of the values used in calculation. Every day we need to compute the result of a lot of simple mathematical equations. For example the cost of the apples bought at the farmer’s market is given by:

        apple price = apple cost per kilo · kilos of apple bought

When we need the result of those mathematical expressions, we put the values on the right hand side of the equation and we compute its result for the left hand side. We usually put wrong (erroneous) numbers into the equation and therefore where is no doubt we get wrong results. There are a lot of reasons why we use incorrect values, for example

1. Most of the values are measured, therefore they are known within a given tolerance. Wikipedia sv. accuracy and precision
2. Some values have an infinite number of digits after the (decimal) point, e.g. π.
3. Some values change with time or samples (or whatever), like the weight of a person, which can change of 5 percent during the day, or the current gain of a bipolar junction transistor (BJP), which can change of 50 percent on the samples of the same series.
4. Some values are estimation or guess—something like between a minimum and a maximum.

For example, if a pipe breaks and you want to buy a new one you need its diameter. If you do not have a caliber, you may measure its circumference and divide it by π.

    diameter = circumference / π

Here are two errors: the circumference is known within the tolerance given by your meter, moreover π has an infinite number of digits while only few of them can be used in the operation. You may think the error is negligible, the result is enough accurate to buy a new pipe in a hardware shop. However, the not infinite precision of those operations avoid the use of computers as automatic theorem demonstration tools and so on.

This kind of issue is quite common in engineer design. What engineers do is to make sure their design will work in the worst case or in most of the cases (usually more than 99.9 percent). A simple example follows.

Let us say you want to repaint the walls of your living room completely messed up by your children. You need to compute how many paint cans you need to buy. The equation is quite simple:

    paint cans = 2 · (room width + room length) · room height
        / (paint per can · paint efficiency),

where “paint efficiency” is how many square meters of surface can be painted with a liter of paint. The problem here is that usually we do not have a long enough meter to measure the room width and length. It is much simpler to count the number of steps to go through it (1 step is about a meter, let us say from 0.9 to 1.1 meters). Moreover, the paint provider usually declares a paint efficiency range.

Here is the data:

- room width = 6 steps (5.4m to 6.6m)
- room length = 4 steps (3.6m to 4.4m)
- room height = 3m (it is assumed to be correct)
- paint efficiency = from 0.7 to 1.3 square meters per liter
- paint liters per can = 40 (it is assumed to be correct)

To compute the average result just put average values in. We get: paint cans = 2 · (6 + 4) · 3 / (40 · 1) = 1.5, which means two paint cans unless you are able to buy just half of the second can.

What happens in the worst case? Just put pessimistic values in the equation. We get: paint cans = 2 · (6.6 + 4.4) · 3 / (40 · 0.7) = 2.36. That is, in the worst case we would be short 0.36 cans of paint. It makes sense to buy 3 cans.

Last, consider the best case. Is it enough to only buy a single can of paint? Just put optimistic values in the equation. We get: paint cans = 2 · (5.4 + 3.6) · 3 / (40 · 1.3) = 1.04, which means one can of paint would not be enough.

You have to buy at least two cans, but probably need one more. For this result we had to go through the equation multiple times (at least twice) and carefully consider for each variable, which would be the most optimistic / pessimistic value assignment, which is not trivial. For example consider the room size versus the paint efficiency: It depends whether the highest or the lowest value takes an optimistic or pessimistic role—and this was a simple example with basic arithmetic operations.

Using interval arithmetic it is possible to compute the result in a single run with ranges as inputs. The following example demonstrates this and further below is explained how it works.

    step = midrad (1, "0.1");
    w = 6 * step;
    l = 4 * step;
    h = 3;
    eff = infsupdec ("[0.7, 1.3]");
    cansize = 40;
    cans = 2 * (w + l) * h / (eff * cansize)
      ⇒ cans ⊂ [1.0384, 2.3572]_com

    ## Since we can only buy whole cans
    ceil (cans)
      ⇒ ans = [2, 3]_def

### 2.3 Pros and Cons ¶

Interval arithmetic, introduced in the 1960s, is a young and powerful technique. Its first application has been to control errors in computations and simplify error analysis for engineers (rounding errors, truncation errors, and conversion errors). The range evaluation of functions has soon been exploited for reliably checking for certain function values and for self-verifying algorithms. Latest usage scenarios comprise root finding, function approximation, and robust pattern recognition. More useful applications are certainly left to be detected.

The major problem in interval arithmetic is that errors can easily build up, such that the final result is too wide to be useful. This is especially true when the dependency problem applies, that is, a single variable occurs several times within a computation and is represented by an interval in each occurrence. Then, the variable virtually may take different values independently, which introduces a systematic error. For example, computing `x .^ 2` will always yield a subset of `times (x, x)`, the latter considers two intervals independent of each other.

    x = infsupdec ("[-1, 3]");
    x .^ 2
      ⇒ ans = [0, 9]_com
    times (x, x)
      ⇒ ans = [-3, +9]_com

After all, it is possible to reduce overestimation errors by subdividing the function’s domain into smaller intervals, e. g., with bisection. This technique is called “mincing”. The computational errors are proportional to the interval width and a linear convergence can be achieved.

    x1 = infsupdec ("[-1, 1]");
    x2 = infsupdec ("[1, 3]");
    hull (x1 .^ 2, x2 .^ 2)
      ⇒ ans = [0, 9]_com
    hull (times (x1, x1), times (x2, x2))
      ⇒ ans = [-1, +9]_com

However, this does not help when ranges of input values are too big. For certain applications it is better to use statistical models, where infinite domains are supported.

### 2.4 Theory ¶

There are good introductions to interval arithmetic available and should be consulted for a deeper understanding of the topic. The following recommendations can make a starting point.

- Wikipedia sv. interval arithmetic
- Introduction to Interval Analysis (2009), by Ramon E. Moore, R. Baker Kearfott, and Michael J. Cloud. Cambridge University Press. ISBN 978-0898716696.
- Introduction to the methods used in interval arithmetic (French).
- Interval analysis in MATLAB Note: The INTLAB toolbox for Matlab is not entirely compatible with this interval package for GNU Octave. However, basic operations can be compared and should be compatible for common intervals.
- Interval related collection of links

## 3 Examples ¶

This chapter presents some more or less exotic use cases for the interval package.

### 3.1 Arithmetic with System-independent Accuracy ¶

According to IEEE Std 754 only the most basic floating-point operations must be provided with high accuracy. This is also true for the arithmetic functions in Octave. It is no surprise that many arithmetic functions fail to provide perfect results and their output may be system dependent.

We compute the cosecant for 100 different values.

    x = vec (1 ./ magic (10));
    sum (subset (csc (x), csc (infsupdec (x))))
      ⇒ ans =  98

Due to the general containment rule of interval arithmetic `x ∈ X ⇒ f (x) ∈ f (X)` one would expect the `csc (x)` to always be contained in the interval version of the cosecant for the same input. However, the classic cosecant is not very accurate whereas the interval version is. In 2 out 100 cases the built-in cosecant is less accurate than 1 ULP.

### 3.2 Prove the Existence of a Fixed Point ¶

A weaker formulation of Brower’s fixed-point theorem goes: If x is a bounded interval and function f is continuous and f (x) ⊂ x, then there exists a point x₀ ∈ x such that f (x₀) = x₀.

These properties can be tested automatically. Decorated intervals can even prove that the function is continuous.

    x = infsupdec ("[-1, +1]");
    f = @cos;
    subset (f (x), x)
      ⇒ ans =  1
    iscommoninterval (x)
      ⇒ ans =  1
    continuous = strcmp (decorationpart (f (x)), "com")
      ⇒ continuous =  1

Furthermore it is sometimes possible to approximate the fixed-point by repetitive evaluation of the function, although there are better methods to do so in general.

    for i = 1 : 20
        x = f (x);
    endfor
    display (x)
      ⇒ x ⊂ [0.73893, 0.73919]_com

### 3.3 Floating-point Numbers ¶

Floating-point numbers are most commonly used in binary64 format, a.k.a. double precision. Internally they are stored in the form `± m * 2 ^ e` with some integral mantissa m and exponent e. Most decimal fractions can only be stored approximately in this format.

The intervaltotext function can be used to output the approximate value up to the last decimal digit.

    intervaltotext (infsup (0.1), "[.55g]")
      ⇒ ans = [0.1000000000000000055511151231257827021181583404541015625]

It can be seen that 0.1 is converted into the most accurate floating-point number. In this case that value is greater than 0.1. The next lower value can be seen after producing an interval enclosure around 0.1 with the nearest floating-point numbers in each direction.

    intervaltotext (infsup ("0.1"), "[.55g]")
      ⇒ ans =
        [0.09999999999999999167332731531132594682276248931884765625,
         0.1000000000000000055511151231257827021181583404541015625]

The error of this approximation can be examined with the @infsup/wid function.

    wid (infsup ("0.1"))
      ⇒ ans =    1.3878e-17

With the @infsup/nextout function an interval can be enlarged in each direction up to the next floating-point number. Around zero the distance towards the next floating point number is very small, but gets bigger for numbers of higher magnitude.

    wid (nextout (infsup ([0, 1, 1e10, 1e100])))
      ⇒ ans =

           9.8813e-324    3.3307e-16    3.8147e-06    3.8853e+84

### 3.4 Root Finding ¶

#### 3.4.1 Interval Newton Method ¶

In numerical analysis, Newton’s method can find an approximation to a root of a function. Starting at a location x₀ the algorithms executes the following step to produce a better approximation:

    x₁ = x₀ − f (x₀) / f’ (x₀)

The step can be interpreted geometrically as an intersection of the graph’s tangent with the x-axis. Eventually, this may converge to a single root. In interval analysis, we start with an interval x₀ and utilize the following interval Newton step:

    x₁ = (mid (x₀) − f (mid (x₀)) / f’ (x₀)) ∩ x₀

Here we use the pivot element `mid (x₀)` and produce an enclosure of all possible tangents with the x-axis. In special cases the division with `f' (x₀)` yields two intervals and the algorithm bisects the search range. Eventually this algorithm produces enclosures for all possible roots of the function f in the interval x₀. The interval newton method is implemented by the function @infsup/fzero.

To produce the derivative of function f, the automatic differentiation from the symbolic package bears a helping hand. However, be careful since this may introduce numeric errors with coefficients.

    f = @(x) sqrt (x) + (x + 1) .* cos (x);

    pkg load symbolic
    df = function_handle (diff (formula (f (sym ("x")))))
      ⇒ df =

       @(x) -(x + 1) .* sin (x) + cos (x) + 1 ./ (2 * sqrt (x))

    fzero (f, infsup ("[0, 6]"), df)
      ⇒ ans ⊂ 2×1 interval vector

            [2.059, 2.0591]
           [4.3107, 4.3108]

We could find two roots in the interval [0, 6].

#### 3.4.2 Bisection ¶

Consider the function `f (x, y) = -(5*y - 20*y^3 + 16*y^5)^6 + (-(5*x - 20*x^3 + 16*x^5)^3 + 5*y - 20*y^3 + 16*y^5)^2`, which has several roots in the area x, y ∈ [-1, 1].

f (x, y) which shows a lot of roots for the function" />

The function is particular difficult to compute with intervals, because its variables appear several times in the expression, which benefits overestimation from the dependency problem. Computing root enclosures with the @infsup/fsolve function is unfeasible, because many bisections would be necessary until the algorithm terminates with a useful result. It is possible to reduce the overestimation with the @infsup/polyval function to some degree, but since this function is quite costly to compute, it does not speed up the bisecting algorithm.

    f = @(x,y) ...
        -(5.*y - 20.*y.^3 + 16.*y.^5).^6 + ...
            (-(5.*x - 20.*x.^3 + 16.*x.^5).^3 + ...
               5.*y - 20.*y.^3 + 16.*y.^5).^2;
    X = Y = infsup ("[-1, 1]");
    has_roots = n = 1;

    for iter = 1 : 10
        ## Bisect
        [i,j] = ind2sub ([n,n], has_roots);
        X = infsup ([X.inf,X.inf,X.mid,X.mid],[X.mid,X.mid,X.sup,X.sup]);
        Y = infsup ([Y.inf,Y.mid,Y.inf,Y.mid],[Y.mid,Y.sup,Y.mid,Y.sup]);
        ii = [2*(i-1)+1,2*(i-1)+2,2*(i-1)+1,2*(i-1)+2] ;
        jj = [2*(j-1)+1,2*(j-1)+1,2*(j-1)+2,2*(j-1)+2] ;
        has_roots = sub2ind ([2*n,2*n], ii, jj);
        n *= 2;

        ## Check if function value covers zero
        fval = f (X, Y);
        zero_contained = find (ismember (0, fval));

        ## Discard values without roots
        has_roots = has_roots(zero_contained);
        X = X(zero_contained);
        Y = Y(zero_contained);
    endfor

    colormap gray
    B = false (n);
    B(has_roots) = true;
    imagesc (B)
    axis equal
    axis off

f (x, y)" />

Now we use the same algorithm with the same number of iterations, but also utilize the mean value theorem to produce better enclosures of the function value with first order approximation of the function. The function is evaluated at the interval’s midpoint and a range evaluation of the derivative can be used to produce an enclosure of possible function values.

    f_dx = @(x,y) ...
           -6.*(5 - 60.*x.^2 + 80.*x.^4) .* ...
               (5.*x - 20.*x.^3 + 16.*x.^5).^2 .* ...
               (-(5.*x - 20.*x.^3 + 16.*x.^5).^3 + ...
                    5.*y - 20.*y.^3 + 16.*y.^5);
    f_dy = @(x,y) ...
           -6.*(5 - 60.*y.^2 + 80.*y.^4) .* ...
               (5.*y - 20.*y.^3 + 16.*y.^5).^5 + ...
           2.*(5 - 60.*y.^2 + 80.*y.^4) .* ...
              (-(5.*x - 20.*x.^3 + 16.*x.^5).^3 + ...
                    5.*y - 20.*y.^3 + 16.*y.^5);

    for iter = 1 : 10
        ...
        ## Check if function value covers zero
        fval1 = f (X, Y);
        fval2 = f (mid (X), mid (Y)) + ...
                (X - mid (X)) .* f_dx (X, Y) + ...
                (Y - mid (Y)) .* f_dy (X, Y);
        fval = intersect (fval1, fval2);
        ...
    endfor

By using the derivative, it is possible to reduce overestimation errors and achieve a much better convergence behavior.

f (x, y)" />

### 3.5 Parameter Estimation ¶

#### 3.5.1 Small Search Space ¶

Consider the model `y (t) = p1 * exp (p2 * t)`. The parameters p1 and p2 are unknown, but it is known that the model fulfills the following constraints, which have been obtained using measurements with known error bounds.

    p1, p2  ∈   [-3, 3]
    y (0.2) ∈  [1.5, 2]
    y (1)   ∈  [0.7, 0.8]
    y (2)   ∈  [0.1, 0.3]
    y (4)   ∈ [-0.1, 0.03]

A better enclosure of the parameters p1 and p2 can be estimated with the @infsup/fsolve function.

    ## Model
    y = @(p1, p2, t) p1 .* exp (p2 .* t);
    ## Observations / Constraints
    t = [0.2; 1; 2; 4];
    y_t = infsup ("[1.5, 2]; [0.7, 0.8]; [0.1, 0.3]; [-0.1, 0.03]");
    ## Estimate parameters
    f = @(p1, p2) y (p1, p2, t);
    p = fsolve (f, infsup ("[-3, 3]; [-3, 3]"), y_t)
      ⇒ p ⊂ 2×1 interval vector

             [1.9863, 2.6075]
           [-1.3243, -1.0429]

The resulting `p` guarantees to contain all parameters `[p1; p2]` which satisfy all constraints on y. It is no surprise that `f (p)` intersects the constraints for y.

    f (p(1), p(2))
      ⇒ ans ⊂ 4×1 interval vector

                [1.5241, 2.1166]
              [0.52838, 0.91888]
              [0.14055, 0.32382]
           [0.0099459, 0.040216]

#### 3.5.2 Larger Search Space ¶

Consider the function `f (x) = p1 ^ x * (p2 + p3 * x + p4 * x^2)`. Let’s say we have some known function values (measurements) and want to find matching parameters p1 through p4. The data sets (x, y) can be simulated. The parameters shall be reconstructed from the observed values on the search range p.

Using plain @infsup/fsolve would take considerably longer, because the search range has 4 dimensions. Bisecting intervals requires an exponential number of steps and can easily become inefficient. Thus we use a contractor for function f, which in addition to the function value can produce a refinement for its parameter constraints. Contractors can easily be build using interval reverse operations like @infsup/mulrev, @infsup/sqrrev, @infsup/powrev1, etc.

    ## Simulate some data sets and add uncertainty
    x = -6 : 3 : 18;
    f = @(p1, p2, p3, p4) ...
        p1 .^ x .* (p2 + p3 .* x + p4 .* x .^ 2);
    y = f (1.5, 1, -3, 0.5) .* infsup ("[0.999, 1.001]");

    function [fval, p1, p2, p3, p4] = ...
        contractor (y, p1, p2, p3, p4)
        x = -6 : 3 : 18;
        ## Forward evaluation
        a = p1 .^ x;
        b = p3 .* x;
        c = p2 + b;
        d = p4 .* x .^ 2;
        e = c + d;
        fval = a .* e;
        ## Reverse evaluation and
        ## undo broadcasting of x
        y = intersect (y, fval);
        a = mulrev (e, y, a);
        e = mulrev (a, y, e);
        p1 = powrev1 (x, a, p1);
        p1 = intersect (p1, [], 2);
        c = intersect (c, e - d);
        d = intersect (d, e - c);
        p2 = intersect (p2, c - b);
        p2 = intersect (p2, [], 2);
        b = intersect (b, c - p2);
        p3 = mulrev (x, b, p3);
        p3 = intersect (p3, [], 2);
        p4 = mulrev (x .^ 2, d, p4);
        p4 = intersect (p4, [], 2);
    endfunction

Now, search for solutions in the range of `p` and try to restore the function parameters.

    p = infsup ("[1.1, 2] [1, 5] [-5, -1] [0.1, 5]");
    p = fsolve (@contractor, ...
                p, y, ...
                struct ("Contract", true))'
      ⇒ p ⊂ 4×1 interval vector

             [1.4991, 1.5009]
                  [1, 1.0011]
           [-3.0117, -2.9915]
           [0.49772, 0.50578]

The function parameters 1.5, 1, −3, and 0.5 from above could be restored. The contractor function could significantly improve the convergence speed of the algorithm.

#### 3.5.3 Combination of Functions ¶

Sometimes it is hard to express the search range in terms of a single function and its constraints, when the preimage of the function consists of a union or intersection of different parts. Several contractor functions can be combined using ctc_union or ctc_intersect to make a contractor function for more complicated sets. The combined contractor function allows one to solve for more complicated sets in a single step.

    ## General ring contractor
    function [fval, cx1, cx2] = ctc_ring (y, c1, c2, x1, x2)
        ## Forward evaluation
        x1_c1 = x1 - c1;
        x2_c2 = x2 - c2;
        sqr_x1_c1 = x1_c1 .^ 2;
        sqr_x2_c2 = x2_c2 .^ 2;
        fval = hypot (x1_c1, x2_c2);
        ## Reverse evaluation
        y = intersect (y, fval);
        sqr_y = y .^ 2;
        sqr_x1_c1 = intersect (sqr_x1_c1, sqr_y - sqr_x2_c2);
        sqr_x2_c2 = intersect (sqr_x2_c2, sqr_y - sqr_x1_c1);
        x1_c1 = sqrrev (sqr_x1_c1, x1_c1);
        x2_c2 = sqrrev (sqr_x2_c2, x2_c2);
        cx1 = intersect (x1, x1_c1 + c1);
        cx2 = intersect (x2, x2_c2 + c2);
    endfunction

    ## Ring 1 with center at (1, 3)
    ## Ring 2 with center at (2, -1)
    ctc_ring1 = @(y, x1, x2) ctc_ring (y, 1, 3, x1, x2);
    ctc_ring2 = @(y, x1, x2) ctc_ring (y, 2, -1, x1, x2);

    ## Unite ring 1 with radius 3..4 and ring 2 with radius 5..6
    ctc_union_of_rings = ctc_union (ctc_ring1, "[3, 4]", ...
                                    ctc_ring2, "[5, 6]");

    ## Compute a paving to approximate the union of rings
    ## in the area x, y = -10..10
    [~, paving] = fsolve (ctc_union_of_rings, ...
                          infsup ("[-10, 10] [-10, 10]"), ...
                          struct ("Contract", true));
    plot (paving(1, :), paving(2, :))
    axis equal

Intersections of contractor functions are especially useful to apply several constraints at once. For example, when it is known that a particular location has a distance of a ∈ [3, 4] from object A, located at coordinates (1, 3), and a distance of b ∈ [5, 6] from object B, located at coordinates (2, -1), the intersection of both rings yields all possible locations in the search range. The combined contractor function enables fast convergence of the search algorithm.

    ## Intersect ring 1 with radius 3..4 and ring 2 with radius 5..6
    ctc_intersection_of_rings = ctc_intersect (ctc_ring1, "[3, 4]", ...
                                               ctc_ring2, "[5, 6]");

    ## Compute a paving to approximate the intersection of rings
    ## in the area x, y = -10..10
    [~, paving] = fsolve (ctc_intersection_of_rings, ...
                          infsup ("[-10, 10] [-10, 10]"), ...
                          struct ("Contract", true));
    plot (paving(1, :), paving(2, :))
    axis equal

### 3.6 Path Planning ¶

Figure 3.1: Cameleon Problem: The polygon has to be moved from the left to the right without touching any obstacles along the path.

The problem presented here is a simplified version from the paper L. Jaulin (2001). Path planning using intervals and graphs. Reliable Computing, issue 1, volume 7, 1–15.

There is an object, a simple polygon in this case, which shall be moved from a starting position to a specified target position. Along the way there are obstacles which may not be touched by the polygon. The polygon can be moved in one direction (left to right or right to left) and may be rotated around its lower left corner.

This makes a two dimensional parameter space and any feasible positions can be determined using interval arithmetic like in the examples above.

Then we use a simple path planning algorithm: We move along the centers of adjacent and feasible boxes in the parameter space until we have a closed path from the start position to the end position. The path is guaranteed to be feasible, that is, there will be no collisions if we follow the path.

    # We can build the simple polygon from interval boxes
    global polygon_x = ...
      infsup ("[18,20] [0,20] [0, 2] [ 0,14] [12,14] [10,14]")';
    global polygon_y = ...
      infsup ("[ 0,18] [0, 2] [0,14] [12,14] [ 6,14] [ 6, 8]")';
    global obstacle_x = infsup ("[ 8,11] [25,28]");
    global obstacle_y = infsup ("[10,10] [10,10]");
    color_feasible = [238 232 213] ./ 255;
    color_uncertain = [220 50 47] ./ 255;
    color_path = [133 153 0] ./ 255;

    function feasible = check_collision (obstacle_x, obstacle_y)
      global polygon_x;
      global polygon_y;

      feasible = infsup (zeros (size (obstacle_x)), ...
                         ones (size (obstacle_x)));

      # Check if the obstacle is inside the polygon
      inside = any (...
        subset (obstacle_x, polygon_x) & subset (obstacle_y, polygon_y));
      feasible(inside) = 0;

      # Check if the obstacle is outside the polygon
      outside = all (...
        disjoint (obstacle_x, polygon_x) | ...
        disjoint (obstacle_y, polygon_y));
      feasible(outside) = 1;
    endfunction

    function feasible = check_parameters (x_offset, angle)
      global obstacle_x;
      global obstacle_y;

      # Instead of rotating the polygon, we rotate the obstacles (reverse)
      s = sin (-angle);
      c = cos (-angle);
      f_x = @(x, y) (x - x_offset) .* c - y .* s;
      f_y = @(x, y) (x - x_offset) .* s + y .* c;

      # All obstacles must be considered
      feasible = 1;
      for i = 1 : numel (obstacle_x)
        feasible = min (feasible, ...
          check_collision (f_x (obstacle_x(i), obstacle_y(i)), ...
                           f_y (obstacle_x(i), obstacle_y(i))));
      endfor
    endfunction

    # Compute a paving of feasible polygon states
    [x, paving, inner] = fsolve (...
      @check_parameters, ...
      infsup ("[-28, 57] [-1.4, 2.7]"), ...
      1, ...
      struct ('MaxIter', 21, 'TolX', 0.03));
    hold on
    plot (paving(1, inner), paving(2, inner), color_feasible);
    plot (paving(1, not (inner)), paving(2, not (inner)), color_uncertain);

    # Consider only states that are guaranteed to be feasible
    paving = paving (:, inner);

    # Path search
    start_idx = find (all (ismember ([0; 0], paving)), 1);
    end_idx = find (all (ismember ([17; 0], paving)), 1);
    adjacency = not (disjoint (paving(1, :), transpose (paving(1, :))) | ...
                     disjoint (paving(2, :), transpose (paving(2, :))));

    # Do a Dijkstra search until we reach end_idx
    distance = inf (columns (paving), 1); # nan = visited
    previous = zeros (columns (paving), 1);
    distance(start_idx) = 0;
    while (not (isnan (distance(end_idx))))
      [pivot_distance, pivot_idx] = min (distance);
      visited = isnan (distance);
      neighbors_idx = adjacency(:, pivot_idx) & not (visited);
      if (not (any (neighbors_idx)))
        error ("Cannot reach target location")
      endif
      neighbors_distance = subsasgn (...
        distance, ...
        substruct ("()", {neighbors_idx}), ...
        pivot_distance + hypot (...
          # Compute distance between centers of boxes
          mid (paving(1, pivot_idx)) - mid (paving(1, neighbors_idx)), ...
          mid (paving(2, pivot_idx)) - mid (paving(2, neighbors_idx))));
      shorter_path = neighbors_idx & (neighbors_distance < distance);
      previous(shorter_path) = pivot_idx;
      distance(shorter_path) = neighbors_distance(shorter_path);
      distance(pivot_idx) = nan;
    endwhile

    # Plot the path to the target location
    last_idx = end_idx;
    while (last_idx != start_idx)
      next_idx = previous(last_idx);
      x1 = mid (paving(1, last_idx));
      y1 = mid (paving(2, last_idx));
      x2 = mid (paving(1, next_idx));
      y2 = mid (paving(2, next_idx));
      plot ([x1 x2], [y1 y2], 'linewidth', 2, 'color', color_path);
      last_idx = next_idx;
    endwhile

The script visualizes the solution in the parameter space. Unfeasible parameters are white, and uncertain combinations of parameters are red. The algorithm’s accuracy is just good enough to find a closed path, which is drawn in green color. The uncertain red area is quite big because we have used a very simple check for verification whether the polygon overlaps the obstacles. This could be improved.

The solution is not optimal, please refer to Luc Jaulin’s paper for more sophisticated approaches. However, we could find a valid solution that moves the polygon as desired without touching any obstacles.

Figure 3.2: Cameleon Problem: A possible solution which moves the polygon from the left to the right without touching obstacles.

## 4 Advanced Topics ¶

### 4.1 Error Handling ¶

Due to the nature of set-based interval arithmetic, one should not observe errors (in the sense of raised GNU Octave error messages) during computation unless operations are evaluated for incompatible data types. Arithmetic operations which are not defined for (parts of) their input, simply ignore anything that is outside of their domain.

However, the interval constructors can produce warnings depending on the input. The @infsup/infsup constructor will warn if the interval boundaries are invalid and returns empty intervals in these cases. Contrariwise, the (preferred) @infsupdec/infsupdec, midrad and hull constructors will only issue a warning and return [NaI] objects, which will propagate and survive through computations. NaI stands for “not an interval”.

Effects of set-based interval arithmetic on partial functions and the NaI object

    ## Evaluation of a function outside of its domain
    ## returns an empty interval
    infsupdec (2) / 0
      ⇒ ans = [Empty]_trv
    infsupdec (0) ^ infsupdec (0)
      ⇒ ans = [Empty]_trv

    ## Illegal interval construction creates a NaI
    infsupdec (3, 2)
      -| warning: illegal interval boundaries:
      -| infimum greater than supremum
      ⇒ ans = [NaI]

    ## NaI even survives through computations
    ans + 1
      ⇒ ans = [NaI]

There are some situations where the interval package cannot decide whether an error occurred or not and issues a warning. The user may choose to ignore these warnings or handle them as errors, see `help warning` for instructions.

`interval:PossiblyUndefined`
  Reason

    Interval construction with boundaries in decimal format, and the constructor can’t decide whether the lower boundary is smaller than the upper boundary. Both boundaries are very close and lie between subsequent binary64 numbers.

  Possible consequences

    The constructed interval is a valid and tight enclosure of both numbers. If the lower boundary was actually greater than the upper boundary, this illegal interval is not considered an error.

`interval:ImplicitPromote`
  Reason

    An interval operation has been evaluated on both, a bare and a decorated interval. The bare interval has been converted into a decorated interval in order to produce a decorated result. Note: This warning does not occur if a bare interval literal string gets promoted into a decorated interval, e. g., `infsupdec (1, 2) + "[3, 4]"` does not produce this warning whereas `infsupdec (1, 2) + infsup (3, 4)` does. A bare interval can be explicitly promoted with the @infsup/newdec function.

  Possible consequences

    The implicit conversion applies the best possible decoration for the bare interval. If the bare interval has been produced from an interval arithmetic computation, this branch of computation is not covered by the decoration information and the final decoration could be considered wrong. For example, `infsupdec (1, 2) + infsup (0, 1) ^ 0` would ignore that 0^0 is undefined.

`interval:UndefinedOperation`
  Reason

    An error has occurred during interval construction and the NaI object has been produced (an empty interval in case of the bare interval constructor). The warning text contains further details. A NaI can be explicitly created with the nai function.

  Possible consequences

    Nothing bad is going to happen, because the semantics of NaI and empty intervals are well defined by IEEE Std 1788-2015. However, the user might choose to cancel the algorithm immediately when the NaI is encountered for the first time.

Table 4.1: Warning IDs

### 4.2 Decorations ¶

The interval package provides a powerful decoration system for intervals, as specified by IEEE Std 1788-2015, IEEE standard for interval arithmetic. By default any interval carries a decoration, which collects additional information about the course of function evaluation on the interval data.

Only the (unfavored) @infsup/infsup constructor creates bare, undecorated intervals and the @infsupdec/intervalpart operation may be used to demote decorated intervals into bare, undecorated ones. It is highly recommended to always use the decorated interval arithmetic, which gives additional information about an interval result in exchange for a tiny overhead.

The following decorations are available:

Decoration | Bounded | Continuous | Defined | Definition
--- | --- | --- | --- | ---
com (common) | ✓ | ✓ | ✓ | x is a bounded, nonempty subset of Dom(f); f is continuous at each point of x; and the computed interval f(x) is bounded
dac (defined and continuous) |  | ✓ | ✓ | x is a nonempty subset of Dom(f); and the restriction of f to x is continuous
def (defined) |  |  | ✓ | x is a nonempty subset of Dom(f)
trv (trivial) |  |  |  | always true (so gives no information)
ill (ill-formed) |  |  |  | Not an interval, at least one interval constructor failed during the course of computation

The decoration information is especially useful after a very long and complicated function evaluation. For example, when the “def” decoration survives until the final result, it is proven that the overall function is actually defined for all values covered by the input intervals.

Examples of using the decoration system

    x = infsupdec (3, 4)
      ⇒ x = [3, 4]_com
    y = x - 3.5
      ⇒ y = [-0.5, +0.5]_com

    ## The square root function ignores any negative part of the input,
    ## but the decoration indicates whether this has or has not happened.
    sqrt (x)
      ⇒ ans ⊂ [1.732, 2]_com
    sqrt (y)
      ⇒ ans ⊂ [0, 0.70711]_trv

Please note that decoration information will not survive through reverse operations (see below) and set operations.

### 4.3 Specialized interval constructors ¶

Above mentioned interval construction with decimal numbers or numeric data is straightforward. Beyond that, there are more ways to define intervals or interval boundaries.

- Hexadecimal-floating-constant form: Each interval boundary may be defined by a hexadecimal number (optionally containing a point) and an exponent field with an integral power of two as defined by the C99 standard ISO/IEC9899, N1256, §6.4.4.2. This can be used as a convenient way to define interval boundaries in binary64 precision, because the hexadecimal form is much shorter than the decimal representation of many numbers.
- Rational literals: Each interval boundary may be defined as a fraction of two decimal numbers. This is especially useful if interval boundaries shall be tightest enclosures of fractions, that would be hard to write down as a decimal number.
- Uncertain form: The interval as a whole can be defined by a midpoint or upper/lower boundary and an integral number of “units in last place” (ULPs) as an uncertainty. The format is `m?ruE`, where
  `m`

    is a mantissa in decimal,

  `r`

    is either empty (which means ½ ULP) or is a non-negative decimal integral ULP count or is the `?` character (for unbounded intervals),

  `u`

    is either empty (symmetrical uncertainty of r ULPs in both directions) or is either `u` (up) or `d` (down),

  `E`

    is either empty or an exponent field comprising the character `e` followed by a decimal integer exponent (base 10).

Examples of different formats during interval construction

    infsupdec ("0x1.999999999999Ap-4") # hex-form
      ⇒ ans ⊂ [0.1, 0.10001]_com

    infsupdec ("1/3", "7/9") # rational form
      ⇒ ans ⊂ [0.33333, 0.77778]_com

    infsupdec ("121.2?") # uncertain form
      ⇒ ans ⊂ [121.14, 121.25]_com

    infsupdec ("5?32e2") # uncertain form with ulp count
      ⇒ ans = [-2700, +3700]_com

    infsupdec ("-42??u") # unbound uncertain form
      ⇒ ans = [-42, +inf]_dac

The hex-form can be set for output with the `format hex` command.

### 4.4 Reverse Arithmetic Operations ¶

Some arithmetic functions also provide reverse mode operations. That is inverse functions with interval constraints. For example the @infsup/sqrrev function can compute the inverse of the `x .^ 2` function on intervals. The syntax is `sqrrev (C, X)` and will compute the enclosure of all numbers x ∈ X that fulfill the constraint x² ∈ C.

In the following example, we compute the constraints for base and exponent of the power function pow as shown in the figure.

Figure 4.1: Reverse power operations. A relevant subset of the function’s domain is outlined and hatched. In this example we use x^y ∈ [2, 3].

    x = powrev1 (infsupdec ("[1.1, 1.45]"), infsupdec (2, 3))
      ⇒ x ⊂ [1.6128, 2.7149]_trv
    y = powrev2 (infsupdec ("[2.14, 2.5]"), infsupdec (2, 3))
      ⇒ y ⊂ [0.75647, 1.4441]_trv

### 4.5 Tips and Tricks ¶

For convenience it is possible to implicitly call the interval constructor during all interval operations if at least one input already is an interval object.

    infsupdec ("17.7") + 1
      ⇒ ans ⊂ [18.699, 18.701]_com
    ans + "[0, 2]"
      ⇒ ans ⊂ [18.699, 20.701]_com

Interval functions with only one argument can be called by using property syntax, e. g. `x.inf`, `x.sup` or even `x.sqrt`.

Whilst most functions (@infsup/size, @infsup/isvector, @infsup/ismatrix, …) work as expected on interval data types, the function @infsup/isempty is evaluated element-wise and checks if an interval equals the empty set.

    builtin ("isempty", empty ())
      ⇒ ans = 0
    isempty (empty ())
      ⇒ ans =  1

### 4.6 Validation ¶

The interval package contains an extensive test suite, which can be run with the command `__run_test_suite__ ({pkg("list", "interval"){}.dir}, {})` to verify correct functionality for a particular system.

In addition, examples from the package documentation can be verified using the doctest package:

    pkg load doctest
    doctest (pkg ("list", "interval"){}.dir)

## Appendix A IEEE Std 1788-2015 ¶

The IEEE standard for interval arithmetic is an important asset for the general use of interval arithmetic. Several interval arithmetic libraries have been created (most popular for the language C++), which vary greatly in their philosophy, completeness and—most important—mathematical definition of certain functions and arithmetic evaluation. The standard grants support for several interval arithmetic flavors, but fights incompatibilities on many layers: Interval arithmetic applications shall be portable, predictable, and reproducible. This is especially important since interval arithmetic shall lead to reliable results. Also a common standard is necessary to catalyze the availability of (fast) interval operations in hardware.

For all conforming implementations certain accuracy constraints must be satisfied and a good amount of interval functions must be implemented. It is defined how to handle functions that are not globally defined or have limiting values. Also such basic things like interval representation, many useful constructors, and interchange encoding are addressed.

The interval package for GNU Octave is the first complete implementation that claims to be standard conforming.

### A.1 Function Names ¶

In terms of a better integration into the GNU Octave language, several operations use a function name which is different from the name proposed in the standard document. The following table translates and lists the implemented function names of the IEEE standard for interval arithmetic.

The implementation provides several additional functions, but this section lists only functions that are mentioned in IEEE Std 1788-2015.

#### A.1.1 Interval constants ¶

See empty and entire.

#### A.1.2 Constructors ¶

The operations textToInterval (S), numsToInterval (l, u), and setDec (x) are implemented by the class constructors @infsup/infsup for bare intervals and @infsupdec/infsupdec for decorated intervals.

#### A.1.3 Required functions ¶

Operation | Implementation | Tightness
--- | --- | ---
| Basic operations |
neg (x) | @infsup/uminus | tightest
add (x, y) | @infsup/plus | tightest
sub (x, y) | @infsup/minus | tightest
mul (x, y) | @infsup/times | tightest
div (x, y) | @infsup/rdivide | tightest
recip (x) | `1 ./ x` | tightest
sqr (x) | `x .^ 2` | tightest
sqrt (x) | @infsup/realsqrt | tightest
fma (x, y, z) | @infsup/fma | tightest
| Power functions |
pown (x, p) | @infsup/pown | tightest
pow (x, y) | @infsup/pow | tightest
exp (x) | @infsup/exp | tightest
exp2 (x) | @infsup/pow2 | tightest
exp10 (x) | @infsup/pow10 | tightest
log (x) | @infsup/log | tightest
log2 (x) | @infsup/log2 | tightest
log10 (x) | @infsup/log10 | tightest
| Trigonometric / hyperbolic |
sin (x) | @infsup/sin | tightest
cos (x) | @infsup/cos | tightest
tan (x) | @infsup/tan | tightest
asin (x) | @infsup/asin | tightest
acos (x) | @infsup/acos | tightest
atan (x) | @infsup/atan | tightest
atan2 (y, x) | @infsup/atan2 | tightest
sinh (x) | @infsup/sinh | tightest
cosh (x) | @infsup/cosh | tightest
tanh (x) | @infsup/tanh | tightest
asinh (x) | @infsup/asinh | tightest
acosh (x) | @infsup/acosh | tightest
atanh (x) | @infsup/atanh | tightest
| Integer functions |
sign (x) | @infsup/sign | tightest
ceil (x) | @infsup/ceil | tightest
floor (x) | @infsup/floor | tightest
trunc (x) | @infsup/fix | tightest
roundTiesToEven (x) | @infsup/roundb | tightest
roundTiesToAway (x) | @infsup/round | tightest
| Absmax functions |
abs (x) | @infsup/abs | tightest
min (x, y) | @infsup/min | tightest
max (x, y) | @infsup/max | tightest

Table A.1: Required forward elementary functions

Operation | Implementation | Tightness
--- | --- | ---
| From unary functions |
sqrRev (c, x) | @infsup/sqrrev | tightest
absRev (c, x) | @infsup/absrev | tightest
pownRev (c, x, p) | @infsup/pownrev | valid (tightest for p ≥ −2)
sinRev (c, x) | @infsup/sinrev | valid
cosRev (c, x) | @infsup/cosrev | valid
tanRev (c, x) | @infsup/tanrev | valid
coshRev (c, x) | @infsup/coshrev | tightest
| From binary functions |
mulRev (b, c, x) | @infsup/mulrev | tightest
powRev1 (b, c, x) | @infsup/powrev1 | valid
powRev2 (a, c, x) | @infsup/powrev2 | valid
atan2Rev1 (b, c, x) | @infsup/atan2rev1 | valid
atan2Rev2 (a, c, x) | @infsup/atan2rev2 | valid
| Two-output division |
mulRevToPair (b, c) | @infsup/mulrev | tightest
| Cancellative addition and subtraction |
cancelMinus (x, y) | @infsup/cancelminus | tightest
cancelPlus (x, y) | @infsup/cancelplus | tightest

Table A.2: Required reverse functions

Operation | Implementation | Tightness
--- | --- | ---
intersection (x, y) | @infsup/intersect | tightest
convexHull (x, y) | @infsup/union | tightest

Table A.3: Required set operations

Operation | Implementation | Rounding mode
--- | --- | ---
inf (x) | @infsup/inf |
sup (x) | @infsup/sup |
mid (x) | @infsup/mid | to nearest, ties to even
wid (x) | @infsup/wid | toward +∞
rad (x) | @infsup/rad | toward +∞
mag (x) | @infsup/mag |
mig (x) | @infsup/mig |

Table A.4: Required numeric functions of intervals

Operation | Implementation | Description
--- | --- | ---
isEmpty (x) | @infsup/isempty | x is the empty set
isEntire (x) | @infsup/isentire | x is the whole line
equal (a, b) | @infsup/eq | a equals b
subset (a, b) | @infsup/subset | a is a subset of b
less (a, b) | @infsup/le | a is weakly less than b
precedes (a, b) | @infsup/precedes | a is left of but may touch b
interior (a, b) | @infsup/interior | a is interior to b
strictLess (a, b) | @infsup/lt | a is strictly less than b
strictPrecedes (a, b) | @infsup/strictprecedes | a is strictly left of b
disjoint (a, b) | @infsup/disjoint | a and b are disjoint

Table A.5: Required boolean functions of intervals

#### A.1.4 Recommended functions ¶

Operation | Implementation | Tightness / Comments
--- | --- | ---
| Elementary functions |
rootn (x, q) | @infsup/nthroot | valid (tightest for q ≥ −2)
expm1 (x) | @infsup/expm1 | tightest
logp1 (x) | @infsup/log1p | tightest
hypot (x, y) | @infsup/hypot | tightest
rSqrt (x) | @infsup/rsqrt | tightest
| Boolean functions |
isCommonInterval (x) | @infsup/iscommoninterval | (=bound and non-empty)
isSingleton (x) | @infsup/issingleton | (=single real)
isMember (m, x) | @infsup/ismember |
| Extended comparison |
overlap (a, b) | @infsup/overlap |

Table A.6: Recommended functions

#### A.1.5 Operations on/with decorations ¶

See @infsup/newdec, @infsupdec/intervalpart, and @infsupdec/decorationpart. The operation setDec is implemented by @infsupdec/infsupdec.

For comparison of decorations with respect to the propagation order `com > dac > def > trv > ill` use the numeric value returned by `decorationpart (x, "uint8")`.

#### A.1.6 Reduction operations ¶

See mpfr_vector_sum_d and mpfr_vector_dot_d. The operations `sumAbs` and `sumSquare` can be computed with `mpfr_vector_sum_d (rounding mode, abs (x))` and `mpfr_vector_dot_d (rounding mode, x, x)` respectively.

#### A.1.7 Input ¶

The operation `textToInterval` is implemented by the class constructors @infsup/infsup for bare intervals and @infsupdec/infsupdec for decorated intervals. Both are able to operate on interval literals provided as strings.

#### A.1.8 Output ¶

See intervaltotext.

#### A.1.9 Exact text representation ¶

See exacttointerval and @infsup/intervaltoexact.

#### A.1.10 Interchange representation and encoding ¶

See interval_bitpack and @infsup/bitunpack.

### A.2 Conformance Claim ¶

The inverval package version 3.2.2 for GNU Octave is conforming to IEEE Std 1788-2015, IEEE Standard for Interval Arithmetic. It is conforming to the set-based flavor with IEEE 754 conformance for the infsup binary64 interval type and without compressed arithmetic. Additionally it provides no further flavors.

### A.3 Conformance Questionnaire ¶

1. Implementation-defined behavior
   1. What status flags or other means to signal the occurrence of certain decoration values in computations does the implementation provide if any?

      The implementation does not signal the occurrence of decoration values.

2. Documentation of behavior
   1. If the implementation supports implicit interval types, how is the interval hull operation realized?

      The implementation supports explicit interval types only.

   2. What accuracy is achieved (i.e., tightest, accurate, or valid) for each of the implementation’s interval operations?

      The accuracy requirements of IEEE Std 1788-2015 are fulfilled. Most operations achieve tightest accuracy, some operations (especially reverse operations) do not. The tightness of each operation is documented in the function’s documentation string and can be displayed with the `help` command.

   3. Under what conditions is a constructor unable to determine whether a Level 1 value exists that corresponds to the supplied inputs?

      When two different string boundaries for an interval both lie between the same two subsequent binary64 numbers, a PossiblyUndefined warning is created.

   4. How are ties broken in rounding numbers if multiple numbers qualify as the rounded result?

      Tie-breaking uses the IEEE Std 754 default: round ties to even.

   5. How are interval datums converted to their exact text representations?

      The binary64 boundaries are converted into hexadecimal-significand form as required by the standard. It is also possible to convert interval datums to exact interval literals in decimal form, see optional arguments of function intervaltotext for that purpose.

3. Implementation-defined behavior
   1. Does the implementation include the interval overlapping function? If so, how is it made available to the user?

      Yes, the interval overlapping function is implemented under the name @infsup/overlap.

   2. Does the implementation store additional information in a NaI? What functions are provided for the user to set and read this information?

      No additional information is stored in a NaI.

   3. What means if any does the implementation provide for an exception to be signaled when a NaI is produced?

      The creation of a NaI is signaled with GNU Octave’s warning mechanism.

   4. What interval types are supported besides the required ones?

      None.

   5. What mechanisms of exception handling are used in exception handlers provided by the implementation? What additional exception handling is provided by the implementation?

      The exceptions described by IEEE Std 1788-2015 raise a warning, which can be handled with GNU Octave’s warning mechanism. The warning may be customized to produce an error instead and interrupt computation. The implementation provides no additional exception handling.

   6. [Question does not apply to IEEE 754 conforming types.]
   7. Does the implementation include different versions of the same operation for a given type and how are these provided to the user?

      The interval matrix multiplication (see @infsup/mtimes) offers two implementations, with either `tightest` or `valid` accuracy. The user may chose the desired version with an optional argument during the function call.

   8. What combinations of formats are supported in interval constructors?

      Any reasonable combination of the formats described in the standard document is supported.

   9. [Question does not apply to IEEE 754 conforming types.]
   10. What methods are used to read or write strings from or to character streams? Does the implementation employ variations in locales (such as specific character case matching)? This includes the syntax used in the strings for reading and writing.

       Input and output is implemented with GNU Octave string variables in UTF-8 encoding. There is no discrimination between different locales. Character case is ignored during input. The syntax for interval literals is used as described by the standard document.

   11. What is the tightness for the interval to string conversion for all interval types?

       The general-purpose interval to string conversion produces the tightest decimal infsup form which has no more digits than are necessary to separate two binary64 numbers.

   12. What is the result of Level 3 operations for invalid inputs?

       Interval constructors prevent the creation of invalid Level 3 interval datums.

       Any non-interval input to Level 3 operations is implicitly converted into an interval and the operation silently continues on interval inputs.

       If at least one input is a decorated interval, bare interval inputs are implicitly decorated as described by `newDec` in the standard document. Implicit promotion from a bare interval to a decorated interval is signaled with GNU Octave’s warning mechanism.

       If implicit conversion fails (e.g., illegal interval literals), bare interval operations produce empty intervals, whereas the decorated interval operations continue on NaI inputs.

   13. [Question does not apply to IEEE 754 conforming types.]
   14. What decorations does the implementation provide and what is their mathematical definition? How are these decorations mapped when converting an interval to the interchange format?

       The implementation provides the decorations com, dac, def, trv, and ill as described by the standard document.

   15. [Question does not apply to IEEE 754 conforming types.]
4. [Question applies to compressed arithmetic only, which is not supported.]
5. [Questions apply to non-standard flavors only, which are not supported.]

## Appendix B GNU General Public License ¶

Version 3, 29 June 2007

    Copyright © 2007 Free Software Foundation, Inc. http://fsf.org/

    Everyone is permitted to copy and distribute verbatim copies of this
    license document, but changing it is not allowed.

### Preamble ¶

The GNU General Public License is a free, copyleft license for software and other kinds of works.

The licenses for most software and other practical works are designed to take away your freedom to share and change the works. By contrast, the GNU General Public License is intended to guarantee your freedom to share and change all versions of a program—to make sure it remains free software for all its users. We, the Free Software Foundation, use the GNU General Public License for most of our software; it applies also to any other work released this way by its authors. You can apply it to your programs, too.

When we speak of free software, we are referring to freedom, not price. Our General Public Licenses are designed to make sure that you have the freedom to distribute copies of free software (and charge for them if you wish), that you receive source code or can get it if you want it, that you can change the software or use pieces of it in new free programs, and that you know you can do these things.

To protect your rights, we need to prevent others from denying you these rights or asking you to surrender the rights. Therefore, you have certain responsibilities if you distribute copies of the software, or if you modify it: responsibilities to respect the freedom of others.

For example, if you distribute copies of such a program, whether gratis or for a fee, you must pass on to the recipients the same freedoms that you received. You must make sure that they, too, receive or can get the source code. And you must show them these terms so they know their rights.

Developers that use the GNU GPL protect your rights with two steps: (1) assert copyright on the software, and (2) offer you this License giving you legal permission to copy, distribute and/or modify it.

For the developers’ and authors’ protection, the GPL clearly explains that there is no warranty for this free software. For both users’ and authors’ sake, the GPL requires that modified versions be marked as changed, so that their problems will not be attributed erroneously to authors of previous versions.

Some devices are designed to deny users access to install or run modified versions of the software inside them, although the manufacturer can do so. This is fundamentally incompatible with the aim of protecting users’ freedom to change the software. The systematic pattern of such abuse occurs in the area of products for individuals to use, which is precisely where it is most unacceptable. Therefore, we have designed this version of the GPL to prohibit the practice for those products. If such problems arise substantially in other domains, we stand ready to extend this provision to those domains in future versions of the GPL, as needed to protect the freedom of users.

Finally, every program is threatened constantly by software patents. States should not allow patents to restrict development and use of software on general-purpose computers, but in those that do, we wish to avoid the special danger that patents applied to a free program could make it effectively proprietary. To prevent this, the GPL assures that patents cannot be used to render the program non-free.

The precise terms and conditions for copying, distribution and modification follow.

### TERMS AND CONDITIONS ¶

0. Definitions.

   “This License” refers to version 3 of the GNU General Public License.

   “Copyright” also means copyright-like laws that apply to other kinds of works, such as semiconductor masks.

   “The Program” refers to any copyrightable work licensed under this License. Each licensee is addressed as “you”. “Licensees” and “recipients” may be individuals or organizations.

   To “modify” a work means to copy from or adapt all or part of the work in a fashion requiring copyright permission, other than the making of an exact copy. The resulting work is called a “modified version” of the earlier work or a work “based on” the earlier work.

   A “covered work” means either the unmodified Program or a work based on the Program.

   To “propagate” a work means to do anything with it that, without permission, would make you directly or secondarily liable for infringement under applicable copyright law, except executing it on a computer or modifying a private copy. Propagation includes copying, distribution (with or without modification), making available to the public, and in some countries other activities as well.

   To “convey” a work means any kind of propagation that enables other parties to make or receive copies. Mere interaction with a user through a computer network, with no transfer of a copy, is not conveying.

   An interactive user interface displays “Appropriate Legal Notices” to the extent that it includes a convenient and prominently visible feature that (1) displays an appropriate copyright notice, and (2) tells the user that there is no warranty for the work (except to the extent that warranties are provided), that licensees may convey the work under this License, and how to view a copy of this License. If the interface presents a list of user commands or options, such as a menu, a prominent item in the list meets this criterion.

1. Source Code.

   The “source code” for a work means the preferred form of the work for making modifications to it. “Object code” means any non-source form of a work.

   A “Standard Interface” means an interface that either is an official standard defined by a recognized standards body, or, in the case of interfaces specified for a particular programming language, one that is widely used among developers working in that language.

   The “System Libraries” of an executable work include anything, other than the work as a whole, that (a) is included in the normal form of packaging a Major Component, but which is not part of that Major Component, and (b) serves only to enable use of the work with that Major Component, or to implement a Standard Interface for which an implementation is available to the public in source code form. A “Major Component”, in this context, means a major essential component (kernel, window system, and so on) of the specific operating system (if any) on which the executable work runs, or a compiler used to produce the work, or an object code interpreter used to run it.

   The “Corresponding Source” for a work in object code form means all the source code needed to generate, install, and (for an executable work) run the object code and to modify the work, including scripts to control those activities. However, it does not include the work’s System Libraries, or general-purpose tools or generally available free programs which are used unmodified in performing those activities but which are not part of the work. For example, Corresponding Source includes interface definition files associated with source files for the work, and the source code for shared libraries and dynamically linked subprograms that the work is specifically designed to require, such as by intimate data communication or control flow between those subprograms and other parts of the work.

   The Corresponding Source need not include anything that users can regenerate automatically from other parts of the Corresponding Source.

   The Corresponding Source for a work in source code form is that same work.

2. Basic Permissions.

   All rights granted under this License are granted for the term of copyright on the Program, and are irrevocable provided the stated conditions are met. This License explicitly affirms your unlimited permission to run the unmodified Program. The output from running a covered work is covered by this License only if the output, given its content, constitutes a covered work. This License acknowledges your rights of fair use or other equivalent, as provided by copyright law.

   You may make, run and propagate covered works that you do not convey, without conditions so long as your license otherwise remains in force. You may convey covered works to others for the sole purpose of having them make modifications exclusively for you, or provide you with facilities for running those works, provided that you comply with the terms of this License in conveying all material for which you do not control copyright. Those thus making or running the covered works for you must do so exclusively on your behalf, under your direction and control, on terms that prohibit them from making any copies of your copyrighted material outside their relationship with you.

   Conveying under any other circumstances is permitted solely under the conditions stated below. Sublicensing is not allowed; section 10 makes it unnecessary.

3. Protecting Users’ Legal Rights From Anti-Circumvention Law.

   No covered work shall be deemed part of an effective technological measure under any applicable law fulfilling obligations under article 11 of the WIPO copyright treaty adopted on 20 December 1996, or similar laws prohibiting or restricting circumvention of such measures.

   When you convey a covered work, you waive any legal power to forbid circumvention of technological measures to the extent such circumvention is effected by exercising rights under this License with respect to the covered work, and you disclaim any intention to limit operation or modification of the work as a means of enforcing, against the work’s users, your or third parties’ legal rights to forbid circumvention of technological measures.

4. Conveying Verbatim Copies.

   You may convey verbatim copies of the Program’s source code as you receive it, in any medium, provided that you conspicuously and appropriately publish on each copy an appropriate copyright notice; keep intact all notices stating that this License and any non-permissive terms added in accord with section 7 apply to the code; keep intact all notices of the absence of any warranty; and give all recipients a copy of this License along with the Program.

   You may charge any price or no price for each copy that you convey, and you may offer support or warranty protection for a fee.

5. Conveying Modified Source Versions.

   You may convey a work based on the Program, or the modifications to produce it from the Program, in the form of source code under the terms of section 4, provided that you also meet all of these conditions:

   1. The work must carry prominent notices stating that you modified it, and giving a relevant date.
   2. The work must carry prominent notices stating that it is released under this License and any conditions added under section 7. This requirement modifies the requirement in section 4 to “keep intact all notices”.
   3. You must license the entire work, as a whole, under this License to anyone who comes into possession of a copy. This License will therefore apply, along with any applicable section 7 additional terms, to the whole of the work, and all its parts, regardless of how they are packaged. This License gives no permission to license the work in any other way, but it does not invalidate such permission if you have separately received it.
   4. If the work has interactive user interfaces, each must display Appropriate Legal Notices; however, if the Program has interactive interfaces that do not display Appropriate Legal Notices, your work need not make them do so.

   A compilation of a covered work with other separate and independent works, which are not by their nature extensions of the covered work, and which are not combined with it such as to form a larger program, in or on a volume of a storage or distribution medium, is called an “aggregate” if the compilation and its resulting copyright are not used to limit the access or legal rights of the compilation’s users beyond what the individual works permit. Inclusion of a covered work in an aggregate does not cause this License to apply to the other parts of the aggregate.

6. Conveying Non-Source Forms.

   You may convey a covered work in object code form under the terms of sections 4 and 5, provided that you also convey the machine-readable Corresponding Source under the terms of this License, in one of these ways:

   1. Convey the object code in, or embodied in, a physical product (including a physical distribution medium), accompanied by the Corresponding Source fixed on a durable physical medium customarily used for software interchange.
   2. Convey the object code in, or embodied in, a physical product (including a physical distribution medium), accompanied by a written offer, valid for at least three years and valid for as long as you offer spare parts or customer support for that product model, to give anyone who possesses the object code either (1) a copy of the Corresponding Source for all the software in the product that is covered by this License, on a durable physical medium customarily used for software interchange, for a price no more than your reasonable cost of physically performing this conveying of source, or (2) access to copy the Corresponding Source from a network server at no charge.
   3. Convey individual copies of the object code with a copy of the written offer to provide the Corresponding Source. This alternative is allowed only occasionally and noncommercially, and only if you received the object code with such an offer, in accord with subsection 6b.
   4. Convey the object code by offering access from a designated place (gratis or for a charge), and offer equivalent access to the Corresponding Source in the same way through the same place at no further charge. You need not require recipients to copy the Corresponding Source along with the object code. If the place to copy the object code is a network server, the Corresponding Source may be on a different server (operated by you or a third party) that supports equivalent copying facilities, provided you maintain clear directions next to the object code saying where to find the Corresponding Source. Regardless of what server hosts the Corresponding Source, you remain obligated to ensure that it is available for as long as needed to satisfy these requirements.
   5. Convey the object code using peer-to-peer transmission, provided you inform other peers where the object code and Corresponding Source of the work are being offered to the general public at no charge under subsection 6d.

   A separable portion of the object code, whose source code is excluded from the Corresponding Source as a System Library, need not be included in conveying the object code work.

   A “User Product” is either (1) a “consumer product”, which means any tangible personal property which is normally used for personal, family, or household purposes, or (2) anything designed or sold for incorporation into a dwelling. In determining whether a product is a consumer product, doubtful cases shall be resolved in favor of coverage. For a particular product received by a particular user, “normally used” refers to a typical or common use of that class of product, regardless of the status of the particular user or of the way in which the particular user actually uses, or expects or is expected to use, the product. A product is a consumer product regardless of whether the product has substantial commercial, industrial or non-consumer uses, unless such uses represent the only significant mode of use of the product.

   “Installation Information” for a User Product means any methods, procedures, authorization keys, or other information required to install and execute modified versions of a covered work in that User Product from a modified version of its Corresponding Source. The information must suffice to ensure that the continued functioning of the modified object code is in no case prevented or interfered with solely because modification has been made.

   If you convey an object code work under this section in, or with, or specifically for use in, a User Product, and the conveying occurs as part of a transaction in which the right of possession and use of the User Product is transferred to the recipient in perpetuity or for a fixed term (regardless of how the transaction is characterized), the Corresponding Source conveyed under this section must be accompanied by the Installation Information. But this requirement does not apply if neither you nor any third party retains the ability to install modified object code on the User Product (for example, the work has been installed in ROM).

   The requirement to provide Installation Information does not include a requirement to continue to provide support service, warranty, or updates for a work that has been modified or installed by the recipient, or for the User Product in which it has been modified or installed. Access to a network may be denied when the modification itself materially and adversely affects the operation of the network or violates the rules and protocols for communication across the network.

   Corresponding Source conveyed, and Installation Information provided, in accord with this section must be in a format that is publicly documented (and with an implementation available to the public in source code form), and must require no special password or key for unpacking, reading or copying.

7. Additional Terms.

   “Additional permissions” are terms that supplement the terms of this License by making exceptions from one or more of its conditions. Additional permissions that are applicable to the entire Program shall be treated as though they were included in this License, to the extent that they are valid under applicable law. If additional permissions apply only to part of the Program, that part may be used separately under those permissions, but the entire Program remains governed by this License without regard to the additional permissions.

   When you convey a copy of a covered work, you may at your option remove any additional permissions from that copy, or from any part of it. (Additional permissions may be written to require their own removal in certain cases when you modify the work.) You may place additional permissions on material, added by you to a covered work, for which you have or can give appropriate copyright permission.

   Notwithstanding any other provision of this License, for material you add to a covered work, you may (if authorized by the copyright holders of that material) supplement the terms of this License with terms:

   1. Disclaiming warranty or limiting liability differently from the terms of sections 15 and 16 of this License; or
   2. Requiring preservation of specified reasonable legal notices or author attributions in that material or in the Appropriate Legal Notices displayed by works containing it; or
   3. Prohibiting misrepresentation of the origin of that material, or requiring that modified versions of such material be marked in reasonable ways as different from the original version; or
   4. Limiting the use for publicity purposes of names of licensors or authors of the material; or
   5. Declining to grant rights under trademark law for use of some trade names, trademarks, or service marks; or
   6. Requiring indemnification of licensors and authors of that material by anyone who conveys the material (or modified versions of it) with contractual assumptions of liability to the recipient, for any liability that these contractual assumptions directly impose on those licensors and authors.

   All other non-permissive additional terms are considered “further restrictions” within the meaning of section 10. If the Program as you received it, or any part of it, contains a notice stating that it is governed by this License along with a term that is a further restriction, you may remove that term. If a license document contains a further restriction but permits relicensing or conveying under this License, you may add to a covered work material governed by the terms of that license document, provided that the further restriction does not survive such relicensing or conveying.

   If you add terms to a covered work in accord with this section, you must place, in the relevant source files, a statement of the additional terms that apply to those files, or a notice indicating where to find the applicable terms.

   Additional terms, permissive or non-permissive, may be stated in the form of a separately written license, or stated as exceptions; the above requirements apply either way.

8. Termination.

   You may not propagate or modify a covered work except as expressly provided under this License. Any attempt otherwise to propagate or modify it is void, and will automatically terminate your rights under this License (including any patent licenses granted under the third paragraph of section 11).

   However, if you cease all violation of this License, then your license from a particular copyright holder is reinstated (a) provisionally, unless and until the copyright holder explicitly and finally terminates your license, and (b) permanently, if the copyright holder fails to notify you of the violation by some reasonable means prior to 60 days after the cessation.

   Moreover, your license from a particular copyright holder is reinstated permanently if the copyright holder notifies you of the violation by some reasonable means, this is the first time you have received notice of violation of this License (for any work) from that copyright holder, and you cure the violation prior to 30 days after your receipt of the notice.

   Termination of your rights under this section does not terminate the licenses of parties who have received copies or rights from you under this License. If your rights have been terminated and not permanently reinstated, you do not qualify to receive new licenses for the same material under section 10.

9. Acceptance Not Required for Having Copies.

   You are not required to accept this License in order to receive or run a copy of the Program. Ancillary propagation of a covered work occurring solely as a consequence of using peer-to-peer transmission to receive a copy likewise does not require acceptance. However, nothing other than this License grants you permission to propagate or modify any covered work. These actions infringe copyright if you do not accept this License. Therefore, by modifying or propagating a covered work, you indicate your acceptance of this License to do so.

10. Automatic Licensing of Downstream Recipients.

    Each time you convey a covered work, the recipient automatically receives a license from the original licensors, to run, modify and propagate that work, subject to this License. You are not responsible for enforcing compliance by third parties with this License.

    An “entity transaction” is a transaction transferring control of an organization, or substantially all assets of one, or subdividing an organization, or merging organizations. If propagation of a covered work results from an entity transaction, each party to that transaction who receives a copy of the work also receives whatever licenses to the work the party’s predecessor in interest had or could give under the previous paragraph, plus a right to possession of the Corresponding Source of the work from the predecessor in interest, if the predecessor has it or can get it with reasonable efforts.

    You may not impose any further restrictions on the exercise of the rights granted or affirmed under this License. For example, you may not impose a license fee, royalty, or other charge for exercise of rights granted under this License, and you may not initiate litigation (including a cross-claim or counterclaim in a lawsuit) alleging that any patent claim is infringed by making, using, selling, offering for sale, or importing the Program or any portion of it.

11. Patents.

    A “contributor” is a copyright holder who authorizes use under this License of the Program or a work on which the Program is based. The work thus licensed is called the contributor’s “contributor version”.

    A contributor’s “essential patent claims” are all patent claims owned or controlled by the contributor, whether already acquired or hereafter acquired, that would be infringed by some manner, permitted by this License, of making, using, or selling its contributor version, but do not include claims that would be infringed only as a consequence of further modification of the contributor version. For purposes of this definition, “control” includes the right to grant patent sublicenses in a manner consistent with the requirements of this License.

    Each contributor grants you a non-exclusive, worldwide, royalty-free patent license under the contributor’s essential patent claims, to make, use, sell, offer for sale, import and otherwise run, modify and propagate the contents of its contributor version.

    In the following three paragraphs, a “patent license” is any express agreement or commitment, however denominated, not to enforce a patent (such as an express permission to practice a patent or covenant not to sue for patent infringement). To “grant” such a patent license to a party means to make such an agreement or commitment not to enforce a patent against the party.

    If you convey a covered work, knowingly relying on a patent license, and the Corresponding Source of the work is not available for anyone to copy, free of charge and under the terms of this License, through a publicly available network server or other readily accessible means, then you must either (1) cause the Corresponding Source to be so available, or (2) arrange to deprive yourself of the benefit of the patent license for this particular work, or (3) arrange, in a manner consistent with the requirements of this License, to extend the patent license to downstream recipients. “Knowingly relying” means you have actual knowledge that, but for the patent license, your conveying the covered work in a country, or your recipient’s use of the covered work in a country, would infringe one or more identifiable patents in that country that you have reason to believe are valid.

    If, pursuant to or in connection with a single transaction or arrangement, you convey, or propagate by procuring conveyance of, a covered work, and grant a patent license to some of the parties receiving the covered work authorizing them to use, propagate, modify or convey a specific copy of the covered work, then the patent license you grant is automatically extended to all recipients of the covered work and works based on it.

    A patent license is “discriminatory” if it does not include within the scope of its coverage, prohibits the exercise of, or is conditioned on the non-exercise of one or more of the rights that are specifically granted under this License. You may not convey a covered work if you are a party to an arrangement with a third party that is in the business of distributing software, under which you make payment to the third party based on the extent of your activity of conveying the work, and under which the third party grants, to any of the parties who would receive the covered work from you, a discriminatory patent license (a) in connection with copies of the covered work conveyed by you (or copies made from those copies), or (b) primarily for and in connection with specific products or compilations that contain the covered work, unless you entered into that arrangement, or that patent license was granted, prior to 28 March 2007.

    Nothing in this License shall be construed as excluding or limiting any implied license or other defenses to infringement that may otherwise be available to you under applicable patent law.

12. No Surrender of Others’ Freedom.

    If conditions are imposed on you (whether by court order, agreement or otherwise) that contradict the conditions of this License, they do not excuse you from the conditions of this License. If you cannot convey a covered work so as to satisfy simultaneously your obligations under this License and any other pertinent obligations, then as a consequence you may not convey it at all. For example, if you agree to terms that obligate you to collect a royalty for further conveying from those to whom you convey the Program, the only way you could satisfy both those terms and this License would be to refrain entirely from conveying the Program.

13. Use with the GNU Affero General Public License.

    Notwithstanding any other provision of this License, you have permission to link or combine any covered work with a work licensed under version 3 of the GNU Affero General Public License into a single combined work, and to convey the resulting work. The terms of this License will continue to apply to the part which is the covered work, but the special requirements of the GNU Affero General Public License, section 13, concerning interaction through a network will apply to the combination as such.

14. Revised Versions of this License.

    The Free Software Foundation may publish revised and/or new versions of the GNU General Public License from time to time. Such new versions will be similar in spirit to the present version, but may differ in detail to address new problems or concerns.

    Each version is given a distinguishing version number. If the Program specifies that a certain numbered version of the GNU General Public License “or any later version” applies to it, you have the option of following the terms and conditions either of that numbered version or of any later version published by the Free Software Foundation. If the Program does not specify a version number of the GNU General Public License, you may choose any version ever published by the Free Software Foundation.

    If the Program specifies that a proxy can decide which future versions of the GNU General Public License can be used, that proxy’s public statement of acceptance of a version permanently authorizes you to choose that version for the Program.

    Later license versions may give you additional or different permissions. However, no additional obligations are imposed on any author or copyright holder as a result of your choosing to follow a later version.

15. Disclaimer of Warranty.

    THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM “AS IS” WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING, REPAIR OR CORRECTION.

16. Limitation of Liability.

    IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES.

17. Interpretation of Sections 15 and 16.

    If the disclaimer of warranty and limitation of liability provided above cannot be given local legal effect according to their terms, reviewing courts shall apply local law that most closely approximates an absolute waiver of all civil liability in connection with the Program, unless a warranty or assumption of liability accompanies a copy of the Program in return for a fee.

### END OF TERMS AND CONDITIONS ¶

### How to Apply These Terms to Your New Programs ¶

If you develop a new program, and you want it to be of the greatest possible use to the public, the best way to achieve this is to make it free software which everyone can redistribute and change under these terms.

To do so, attach the following notices to the program. It is safest to attach them to the start of each source file to most effectively state the exclusion of warranty; and each file should have at least the “copyright” line and a pointer to where the full notice is found.

    one line to give the program's name and a brief idea of what it does.
    Copyright (C) year name of author

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or (at
    your option) any later version.

    This program is distributed in the hope that it will be useful, but
    WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see http://www.gnu.org/licenses/.

Also add information on how to contact you by electronic and paper mail.

If the program does terminal interaction, make it output a short notice like this when it starts in an interactive mode:

    program Copyright (C) year name of author
    This program comes with ABSOLUTELY NO WARRANTY; for details type ‘show w’.
    This is free software, and you are welcome to redistribute it
    under certain conditions; type ‘show c’ for details.

The hypothetical commands ‘`show w`’ and ‘`show c`’ should show the appropriate parts of the General Public License. Of course, your program’s commands might be different; for a GUI interface, you would use an “about box”.

You should also get your employer (if you work as a programmer) or school, if any, to sign a “copyright disclaimer” for the program, if necessary. For more information on this, and how to apply and follow the GNU GPL, see http://www.gnu.org/licenses/.

The GNU General Public License does not permit incorporating your program into proprietary programs. If your program is a subroutine library, you may consider it more useful to permit linking proprietary applications with the library. If this is what you want to do, use the GNU Lesser General Public License instead of this License. But first, please read http://www.gnu.org/philosophy/why-not-lgpl.html.