RangeLab User’s Guide
RangeLab Version 1.0

Matthieu Martel
 
DALI - Digits Architectures et Logiciels Informatiques
Université de Perpignan Via Domitia
52 avenue Paul Alduy
66860 Perpignan Cedex 9, France
 
and
 
LIRMM - Laboratoire d’Informatique, Robotique et Microélectronique de Montpellier
CNRS: UMR 5506 - Université Montpellier 2
161 rue Ada, 34095 Montpellier Cedex 5, France

Abstract
RangeLab is a tool to automatically validate the accuracy of floating-point or fixed-point computations. Given intervals for the input variables, RangeLab computes the range of the outputs of simple functions with conditional or loops as well as ranges for the roundoff errors arising during these computations. Hence the users not only obtains the range of the result of the computation in the computer arithmetic but also a bound on the difference between the computer result and the result in infinite precision. RangeLab is based on static analysis by abstract interpretation. This manual introduces the syntax and semantics of the arithmetic expressions, operations on matrices, simple control structures as conditionals, while and for loops and functions. It is completed by a few case studies showing what results can be obtained and how to tune the system parameters.

Contents

1  Introduction

RangeLab automatically computes the ranges of the numerical outputs of simple programs and binds the roundoff errors arising during floating-point or fixed-point computations. Typically, RangeLab takes as entry a simple code, with arithmetic expressions, arrays, assignments and simple control structures. Given some ranges for the inputs, the system computes safe ranges for the outputs, accordingly to the computer arithmetic (IEEE754 floating-point arithmetic or fixed-point arithmetic). In addition, the system computes a safe range of the roundoff error on each result, i.e. a range for the difference between the result in the computer arithmetic and the result in infinite precision. Doing so, RangeLab makes it possible to assert the accuracy of some computations described by a simple piece of code.

RangeLab is an interactive systems which displays a prompt symbol ’-]’ in order to indicate that the tool is waiting for user’s commands.

> rangelab

      [ R A N G E . L A B ]

          Version 1.0.1

            June 2011

-] 0.1+0.2

ans = 
       float64: [2.999999999999999E-1,
                 3.000000000000001E-1]
         error: [-4.857225732735059E-17,
                 4.857225732735060E-17]

-] a=2*ans

a = 
     float64: [5.999999999999999E-1,6.000000000000001E-1]
       error: [-1.526556658859590E-16,1.526556658859591E-16]

-] trapeze(0.25,50,100)

ans = 
       float64: [6.557548333631875E1,1.311509666726376E2]
          error: [-1.456853305716057E-12,1.457731336808933E-12]

RangeLab always assumes that the floating-point computations are rounded to the nearest. In the above session, the system computes the numerical integral ∫0.2550 g(x)dx of some function g : R R by means of a procedure trapeze implementing the trapeze rule. The parameters require to truncate the interval [0.25,50] in 100 steps in order to compute the result. The code of the trapeze method is stored in a separate file trapeze.m whose content is:

function r = trapeze(a,b,n)
    r = 0.0; 
    xa = a; 
    h = (b-a)/n; 
    while xa<b,
        xb = xa+h;
        if xb>b, xb = b end;
        r = r + ((g(xb)+g(xa))/2)*h;
        xa = xa+h;
    end

The code of the function g is stored in a file g.m:

function y = g(x)
    y = [1.0,2.0]/(x*x*x*x)

The function g has an interval coefficient indicating that g represents the family of functions G={ xa/x4, 1≤ a ≤ 2}.

In the above session, RangeLab computes that the result of trapeze(0.25,50.0,100) is

F=[6.557548333631875· 101,1.311509666726376· 102]

with an error

E=[−1.4568… · 10−12, 1.4577…  · 10−12]

This means that, for any function in G, trapeze(0.25,50.0,100) returns a value in F and that the numerical absolute error on this result belongs to E, assuming that the errors on the coefficient a is null. We may claim several things from this execution:

Let us also remark that RangeLab does not detect the method error which affects the result of this computation: trapeze(0.25,50.0,100) differs significantly from ∫0.2550 g(x)dx even in infinite precision. Details on the parameters used to perform this analysis of the trapeze function are given in Section 6.4.

RangeLab uses abstract interpretation techniques [2] to compute safely the outputs of a program with loops and conditionals. This manual is an introduction to the system concerning its syntax and semantics. Section 2 concerns the data formats and their interpretation by the system. Section 3 introduces the operations on matrices. Section 4 presents the control structures and other instructions. Section 5 explains how to tune the system. Finally, cases studies are presented in Section 6.

2  Values in RangeLab

Values is RangeLab are somewhat non-standard. We are going to see that error terms are attached to the floating-point and fixed-point numbers and that RangeLab uses intervals instead of scalar values.

2.1  Datatypes

RangeLab assigns a datatype to any value returned by the interactive system. This type is displayed in the answer, beside the result itself:

-] 8

ans = 
       int32: 8

-] 5+4

ans = 
       int32: 9

In the above example, int32 is the type of 32-bits signed integer numbers. RangeLab offers 6 other scalar datatypes: int64 for 64-bits signed integers, float32, float64 and float128 for 32-bits, 64-bits and 128-bits floating-point numbers, mp(X) for multiple-precision floating-point numbers with X bits and fixed(X,Y) for fixed-point numbers with X bits for the integer part and Y bits for the fractional part. In addition Rangelab accepts matrices of scalar elements.

-] 4.5

ans = 
       float64: 4.500000000000000E0
         error: 0.000000000000000E-1

-] 8.25#40

ans = 
       mp(40): 8.250000000000E0
        error: 0.000000000000000E-1

-] 10_75

ans = 
       fixed(4,2): 10_75
            error: 0.000000000000000E-1

-] [1 2;3 4]

ans = 
       int32[]: 1 2 
                3 4

2.2  Error terms

In RangeLab, error terms are attached to the floating-point and fixed-point numbers. They indicate a range for the roundoff error due to the rounding to the nearest of the exact value. The exact error term being possibly not representable in finite precision, RangeLab computes an over-approximation and returns an interval with bounds made of multiple precision floating-point numbers. When the error term can be computed exactly, the error interval is a singleton and RangeLab displays it as a scalar value. In particular, when the error term is null the error displayed by the system is zero.

-] 1.0/3.0

ans = 
       float64: 3.333333333333333E-1
         error: [1.850371707708594E-17,1.850371707708595E-17]

-] 1.0/5.0

ans = 
       float64: 2.000000000000000E-1
         error: [-1.110223024625156E-17,-1.110223024625157E-17]

-] 1.0/2.0

ans = 
       float64: 5.000000000000000E-1
         error: 0.000000000000000E-1

Note that the errors can be either positive or negative. This depends on the direction of the rounding mode to the nearest which can create either an upper or a lower approximation. In the examples above, the errors on 1/3 and on 1/5 come from the rounding of the result of the division of exact values.

By default, RangeLab displays absolute errors. The command set relative makes the system display relative errors (the symbol % then appears right after the error term). The command set absolute makes the system come back to display absolute errors.

-] set relative
-] 1.0/10e8

ans = 
       float64: 1.000000000000000E-9
         error: [6.228159145777985E-17,6.228159145777986E-17]%

-] set absolute
-] 1.0/10e8

ans = 
       float64: 1.000000000000000E-9
         error: [-6.228159145777985E-26,-6.228159145777986E-26]

By default, error terms are intervals of multiple-precision numbers of 256 digits. It is possible to modify this precision and to make RangeLab display more digits for the error terms with the commands set error and set error_display (see Section 5.2).

Error terms are propagated among computations. The error on the result of some operation xy is the propagation of the errors on x and y through the operator • plus the error due to the rounding of the result of the operation itself.

Let x and y be to numbers reprensented in RangeLab by the pairs (fx,ex) and (fy,ey) where fx and fy are the floating-point or fixed-point numbers approximating x and y and ex and ey the error terms of both operands. Let ∘(v) be the rounding of the value v in the current rounding mode and let ε(v) be the rounding error, i.e. the error arising when rounding v into ∘(v). We have ε(v)=v−∘(v). RangeLab uses the following equations to compute the result of an addition, subtraction or product.

x+y = 
∘(fx+fy),ex+ey+ε(fx+fy)
    (1)
xy = 
∘(fxfy),exey+ε(fxfy)
    (2)
x× y = 
∘(fx× fy),fy× ex+fx× ey+ex× ey+ε(fx× fy)
    (3)

For an addition, the errors on the operands are added to the error due to the roundoff of the result. For a subtraction, the errors on the operands are subtracted. Finally, the semantics of the multiplication comes from the development of (fx+ex)× (fy+ey). The equations describing the propagation of errors through a division or square root are more complicated. They are defined in [8].

-] x=1.0/3.0

x = 
     float64: 3.333333333333333E-1
       error: [1.850371707708594E-17,1.850371707708595E-17]

-] y=1.0/7.0

y = 
     float64: 1.428571428571428E-1
       error: [7.930164461608261E-18,7.930164461608262E-18]

-] x+y

ans = 
       float64: 4.761904761904762E-1
         error: [2.643388153869420E-17,2.643388153869421E-17]

-] x-y

ans = 
       float64: 1.904761904761905E-1
         error: [1.057355261547768E-17,1.057355261547769E-17]

-] x*y

ans = 
       float64: 4.761904761904762E-2
         error: [2.643388153869420E-18,2.643388153869421E-18]

2.3  Intervals

In RangeLab, any value is an interval. When this interval is a singleton, as in the examples of Section 2.1 and Section 2.2, it is displayed as a single value rounded to the nearest in order to be printed with a reasonable number of digits (16 by default). Intervals arise when an input value v is not exactly representable in the current format. In this case, v is stored in an interval whose bounds are the round down and round up of the exact value. In addition, the user may directly define intervals of any width.

-] 0.1

ans = 
       float64: [9.999999999999999E-2,1.000000000000001E-1]
         error: [-6.938893903907228E-18,6.938893903907229E-18]

-] [-5,10]

ans = 
       int32: [-5,10]

-] [2.0,5.0]

ans = 
       float64: [2.000000000000000E0,5.000000000000000E0]
         error: 0.000000000000000E-1

-] [0.1,0.2]

ans = 
       float64: [9.999999999999999E-2,2.000000000000001E-1]
         error: 0.000000000000000E-1

When the bounds of an interval are not exactly representable, as in the last example above, RangeLab rounds down the lower bound and rounds up the upper bound. However, by default the error assigned to a user-defined interval is null. The user may explicitely attach error terms to intervals as follows.

-] [0.1,0.2,-5e-8,5e-8]

ans = 
       float64: [9.999999999999999E-2,2.000000000000001E-1]
         error: [-5.000000000000000E-8,5.000000000000001E-8]

Intervals can be used in expressions as ground values.

-] [3,8]/[2,5]

ans = 
       int32: [0,4]

-] a=[-1.0,1.0]

a = 
     float64: [-1.000000000000000E0,1.000000000000000E0]
       error: 0.000000000000000E-1

-] a+a

ans = 
       float64: [-2.000000000000000E0,2.000000000000000E0]
         error: [-1.110223024625156E-16,1.110223024625157E-16]

-] a-a

ans = 
       float64: 0.000000000000000E-1
         error: 0.000000000000000E-1

-] a*a-2.0*a+1.0      

ans = 
       float64: [-1.000000000000000E0,4.000000000000000E0]
         error: [-5.551115123125782E-16,5.551115123125783E-16]

In RangeLab, a pair (f,ε) made of a floating or fixed-point number f and of an error ε is used to represent the real number r=f+ε rounded to f in the current computer format. This is then extended to intervals. By definition, the pair of intervals ([f,f],[ε,ε]) represents the set of floating/fixed-point numbers belonging to [f,f] whose error with respect to the real number they aim at representing belongs to [ε,ε]:

([
f
,
f
],[
ε
,
ε
]) = 

(f,ε)∈  F× R : f∈[
f
,
f
], ε∈[
ε
,
ε
]

    (4)

In Equation (4), F denotes the current format used to represent the number in memory (float32, float64, float128, mp(X) or fixed(X,Y) and R denotes the set of real numbers.

RangeLab arithmetic, based on Equations (1) to (4), differs from the usual interval arithmetic in the following way. In the usual interval arithmetic, an interval [a,b] is used to bound both the floating-point number stored by the machine and the real number that we would have in infinite precision. In other terms, in interval arithmetic, x∈[a,b] means that if the computer stores some value f F∩ [a,b] then it would also store some value in [a,b] if it was working in infinite precision. This property is maintained through computations by rounding the operations on intervals towards outside. Doing so, interval arithmetic does not distinguish two different kinds of errors:

  1. Sensivity errors, due to the fact that a small variation of the inputs may yield a large variation of the outputs, even with real numbers.
  2. Numerical errors, due to the fact that a finite-precision calculation may diverge from a real calculation.

For example, numerical errors arise in the following program which uses simple precision floating-point numbers:

function x = oneminuseps(n,eps)

x = float32(1.0);

for i=1:n,
   x = x - eps; % if eps<<x, absorption: x-eps == x 
end

Let us set ε=10−8. The value 10−8 being subtracted n times to 1, the exact result in the reals is 1−n× 10−8. But 10−8 is less than the least significant digit of float32(1.0) and, consequently, an absorption arises: 1.0 − 10−8 = 1.0 in single precision. So, at the end of the iteration, the floating-point result is x = 1.0 and the error is n× 10−8. This error is a numerical error, there is no sensibility error in this example. RangeLab clearly indicates this matter of fact:

-] oneminuseps(100,float32(1e-8))

ans = 
       float32: 1.0000000E0
         error: [-5.960464560279433E-6,5.960464483616534E-6]

-] oneminuseps(1000,float32(1e-8))

ans = 
       float32: 1.0000000E0
         error: [-5.960464560279433E-5,5.960464483616534E-5]

The usual interval arithmetic (in which the exact result for the left bound in rounded towards −∞ and the result fot the right bound is rounded towards +∞) would return some interval [a,b] such that [1−n× 10−8,1]⊆ [a,b] since both the floating-point and exact result have to belong to [a,b]. For example, we could have [9.9999899· 10−1,1.0] and [9.9998999· 10−1,1.0] for n=100 and n=1000, respectively. Such results of the interval arithmetic are clearly less informative that RangeLab’s outputs.

As pointed out by the last examples, intervals of floating-point or fixed-point numbers and intervals of errors are propagated by RangeLab throughtout the computations. Formally, let ([f,f],[ε,ε]) be the result of some computation C whose inputs are ([f1,f1],[ε1,ε1]),…, ([fn,fn],[εn,εn]). Then, for any f1∈[f1,f1], ε1∈[ε1,ε1],… fn∈[fn,fn], εn∈[εn,εn] we have C((f11),…,(fnn))∈ ([f,f],[ε,ε]). In other words:

C



[
f1
,
f1
],[
ε1
,
ε1
]



[
fn
,
fn
],[
εn
,
εn
]



=([
f
,
f
],[
ε
,
ε
])
=⇒ 
∀ f1∈[
f1
,
f1
], ∀ε1∈[
ε1
,
ε1
],… ∀ fn∈[
fn
,
fn
], ∀εn∈[
εn
,
εn
], C((f11),…,(fnn))∈ ([
f
,
f
],[
ε
,
ε
])
    (5)

2.4  Details about the Data Formats

2.4.1  Integers

RangeLab offers to types of integers, int32 and int64, which correspond to 32 bits and 64 bits signed integers. The default type given to an integer value is int32. It is possible to convert 32 bits values into 64 values and conversely by means of the functions int32(val) and int64(val) or, alternatively, using the command cast(val,’type’).

-] 66

ans = 
       int32: 66

-] int64(ans)

ans = 
       int64: 66

-] cast(ans,'int32')

ans = 
       int32: 66

It is also possible to ask RangeLab to display values in binary. The user may switch between the decimal and binary display using the commands set binary and set decimal, or for short, set bin and set dec. For integers, the sign corresponds to the leftmost bit.

-] set bin
-] 5506

ans = 
       int32: 00000000000000000001010110000010

-] int64(ans)

ans = 
       int64: 0000000000000000000000000000000000000000000000000001010110000010

-] -5506

ans = 
       int32: 10000000000000000001010110000010

-] set dec
-] ans

ans = 
       int64: -5506

When computations involve both 32 bits and 64 bits integers, RangeLab automatically converts the 32 bits values into 64 bits ones. We say that Type int32 is coercible into Type int64.

-] int32(32)+int64(64)

ans = 
       int64: 96

2.4.2  Floating-point Numbers

RangeLab relies on the GMP library for floating-point numbers [13]. Four data formats are available: float32, float64, float128 and mp(X). They correspond to different sizes of mantissa. The representation of the exponent is left to GMP (RangeLab does not aim at detecting overflows). Values of types float32, float64 and float128 are encoded on mantissas of 24, 53 and 113 bits respectively. This corresponds to the IEEE754 formats [1]. The type mp(X) let the user specify the number X of bits of the mantissa.

By default, RangeLab considers that floating-point values are float64. As for integers, it is possible to convert a value from one type to another using the commands float32(val), float64(val), float128(val), mp(X)(val) and cast(val,’type’).

-] x=0.1       

x = 
     float64: [9.999999999999999E-2,1.000000000000001E-1]
       error: [-6.938893903907228E-18,6.938893903907229E-18]

-] float32(x)

ans = 
       float32: [9.9999994E-2,1.0000001E-1]
         error: [-5.960464476151283E-9,1.490116120772545E-9]

-] cast(x,'float128')

ans = 
       float128: [9.999999999999999167332731531132594E-2,1.000000000000000055511151231257828E-1]
          error: [-6.938893903907228E-18,6.938893903907229E-18]

-] mp(12)(x)

ans = 
       mp(12): [9.997E-2,1.001E-1]
        error: [-2.441406249999861E-5,6.103515625001388E-6]

It is important to note that, by default, the floating-point numbers entered by the user are first stored in memory in the float64 format and then possibly converted into another format, if asked. For example the command float128(0.1) makes RangeLab convert the string 0.1 into a float64 value and then convert this value into float128. In this case, the error attached to the resulting float128 value corresponds to the roundoff of 0.1 in the float64 format.

-] 0.1   

ans = 
       float64: [9.999999999999999E-2,1.000000000000001E-1]
         error: [-6.938893903907228E-18,6.938893903907229E-18]

-] float128(0.1)

ans = 
       float128: [9.999999999999999167332731531132594E-2,1.000000000000000055511151231257828E-1]
          error: [-6.938893903907228E-18,6.938893903907229E-18]

In order to make RangeLab avoid the float64 intermediary representation to encode a value on a certain number of bits N, the user may append #N to the value: val#N makes RangeLab convert directly the string val into a floating-point number with N bits of mantissa.

-] 0.1#113

ans = 
       float128: [9.999999999999999999999999999999999E-2,1.000000000000000000000000000000001E-1]
          error: [-6.018531076210112E-36,6.018531076210113E-36]

-] 0.1#250

ans = 
       mp(250): [1.000000000000000000000000000000000000188079096131566001274997845955559308450E-1,
                 1.000000000000000000000000000000000000188079096131566001274997845955559308451E-1]
         error: [-3.454467422037777E-77,3.454467422037778E-77]

-] 0.1#12

ans = 
       mp(12): [9.997E-2,1.001E-1]
        error: [-1.525878906250000E-5,1.525878906250000E-5]

As for integers, floating-point values may be displayed in binary. In this case, if the value is a float32, float64 or float128 then RangeLab prints its IEEE754 representation.

-] set bin
-] -118.625

ans = 
       float64: 1 10000000101 1101101010000000000000000000000000000000000000000000
         error: 0.000000000000000E-1

-] float32(0.1)

ans = 
       float32: [0 01111011 10011001100110011001100,0 01111011 10011001100110011001101]
         error: [-5.960464476151283E-9,1.490116120772545E-9]

The sign, exponent and mantissa are separated by spaces. For multiple-precision numbers, the size of the exponent is not normalized. RangeLab uses the minimal number of bit needed to represent the value as a signed binary number, with no bias.

-] set bin
-] 0.1#12

ans = 
       mp(12): [0 -100 10011001100,0 -100 10011001101]
        error: [-1.525878906250000E-5,1.525878906250000E-5]

Recall, from Section 2.3, that it is possible to assign an initial error term to a floating-point value with the syntactic construct [f,f,ε,ε]:

-] [0.1,0.2,-5e-5,5e-4]

ans = 
       float64: [9.999999999999999E-2,2.000000000000001E-1]
         error: [-5.000000000000000E-5,5.000000000000001E-4]

-] [0.1,0.1,0.0,0.0]

ans = 
       float64: [9.999999999999999E-2,1.000000000000001E-1]
         error: 0.000000000000000E-1

RangeLab automatically performs some coercions between floating-point values of different types. We have float32 coercible to float64 coercible to float128, mp(X), coercible to mp(Y) if XY. Note that float32=mp(24), float64=mp(53) and that float128=mp(113)

-] float32(3.2)+float64(6.4)

ans = 
       float64: [9.599999809265135E0,9.600000047683717E0]
         error: [-1.907348645691087E-7,4.768371719698906E-8]

-] float32(3.2)+float128(12.8)

ans = 
       float128: [1.599999980926513565293589635984972E1,1.600000004768371653085523576010019E1]
          error: [-1.907348641250194E-7,4.768371675289985E-8]

-] float64(6.4)+mp(16)(1.6)

ans = 
       float64: [7.999975585937499E0,8.000006103515626E0]
         error: [-2.441406250131006E-5,6.103515626354473E-6]

Integer values are also coerced to float64 if they occur in the same expression.

-] 1+1.0

ans = 
       float64: 2.000000000000000E0
         error: 0.000000000000000E-1

2.4.3  Fixed-point Numbers

Syntactically, in RangeLab, fixed-point numbers are written with an underscore ’_’ in place of the dot. This symbol separates the integer part of the number from the fractional part.

-] 4_75

ans = 
       fixed(3,2): 4_75
            error: 0.000000000000000E-1

The type fixed(X,Y) corresponds to the signed fixed-point numbers whose integer part is encoded on X bits and whose fractional part is encoded on Y bits, as we can check by activating the binary display mode (the bit for the sign is not included in X and Y).

-] set bin
-] 4_75

ans = 
       fixed(3,2): 100_11
            error: 0.000000000000000E-1

As already mentioned, RangeLab is not intended to manage the overflows and, by default, it allocates enough bits to represent exactly the integer part (but no more). For the fractional part, RangeLab allocates by default 32 bits or less, if the value can be represented exactly with less than 32 bits.

-] 0_1

ans = 
       fixed(1,32): 0_0999999998603
             error: [0.000000000000000E-1,4.656612873077393E-10]

Fixed-point numbers are not rounded like floating-point numbers. When a number is not representable in the current format, it is truncated (this corresponds to a rounding towards zero). The error is then an upper-approximation of the truncated sequence.

As for multiple-precision floating-point numbers, the user may force to N the size of the fractional part of a fixed-point number by appending #N to the number.

-] 4_75#1

ans = 
       fixed(3,1): 4_5
            error: 2.500000000000000E-1

In the above example, the value 4.75 is rounded to 4.5 to fit into a fraction of size 1 and the error attached to this number is 0.25.

It is also possible to modify the default size allocated to the fractional part of the fixed-point numbers with the command set frac:

-] set frac 40
-] -0_1

ans = 
       fixed(1,40): -0_0999999999995
             error: [-1.818989403545856E-12,0.000000000000000E-1]

When fixed-point numbers of different formats occur inside an arithmetic expression, the type of the result is obtained as follows: a new size is computed for the integer part, in order to fit the result and, for the fractional part, the result takes the longuest size of the fractional parts of the operands.

-] a=2_125

a = 
     fixed(2,3): 2_125
          error: 0.000000000000000E-1

-] b=8_5

b = 
     fixed(4,1): 8_5
          error: 0.000000000000000E-1

-] c=a+b

ans = 
       fixed(4,3): 10_625
            error: 0.000000000000000E-1

-] c+c

ans = 
       fixed(5,3): 21_25
            error: 0.000000000000000E-1

Integers and floating-point numbers are also coerced to fixed-point numbers when they occur in the same expression. For integers, the result is a fixed-point number whose integer part is chosen to fit to the value. For floating-point numbers, however, it is not reasonable in general to attempt do represent the integer part without overflow. When the operands of some addition, subtraction or product involve a fixed-point and a floating-point number, the floating-point number is coerced into a fixed-point number in the same format than the other operand. This possibly yields an overflow.

-] 8+1_0

ans = 
       fixed(4,1): 9_
            error: 0.000000000000000E-1

-] 8.0+1_0

ans = 
       fixed(1,1): +oo
            error: 0.000000000000000E-1

As already mentioned, in RangeLab any value is an interval and this also holds for fixed-point numbers. An interval [f,f] with type fixed(X,Y) represents the set of fixed-point numbers of type fixed(X,Y) ranging from f to f.

-] [-0_1,0_1]

ans = 
       fixed(1,32): [-0_0999999998603,0_0999999998603]
             error: [0.000000000000000E-1,4.656612873077393E-10]

If the bounds of the interval have different precisions, the less precise bound is coerced to most precise one.

RangeLab also allows the user to attach initial errors to fixed-point numbers in the same way than for floating-point numbers. Remark that the error terms attached to the fixed-point values are (multiple-precision) floating-point numbers.

-] [-0_2,0_2,0.005,0.01]

ans = 
       fixed(1,32): [-0_199999999953,0_199999999953]
             error: [4.999999999999999E-3,1.000000000000001E-2]

3  Vectors and Matrices

3.1  Specific Aspects of Vectors and Matrices in RangeLab

We have seen in Section 2.1 that, in RangeLab, any value is an interval and that, in addition, intervals of error are attached to floating-point and fixed-point values. This also holds for vectors and matrices. Let M=(aij) be a matrix of size m× n whose values are floating-point or fixed-point numbers with errors such that aij=([fij,fij],[εij,εij]), 1≤ im, 1≤ n. In order to display the matrix M, RangeLab displays successively four matrices: (fij), (fij), (εij) and finally (εij).

-] a=[0.1,0.2]

a = 
     float64: [9.999999999999999E-2,2.000000000000001E-1]
       error: [-1.387778780781445E-17,1.387778780781446E-17]

-] [a a a; a a a]

ans = 
       float64[]: [9.999999999999999E-2 9.999999999999999E-2 9.999999999999999E-2 
                   9.999999999999999E-2 9.999999999999999E-2 9.999999999999999E-2 
                  ,
                   2.000000000000001E-1 2.000000000000001E-1 2.000000000000001E-1 
                   2.000000000000001E-1 2.000000000000001E-1 2.000000000000001E-1 
                  ]
          Errors: [-1.387778780781445E-17 -1.387778780781445E-17 -1.387778780781445E-17 
                   -1.387778780781445E-17 -1.387778780781445E-17 -1.387778780781445E-17 
                  ,
                   1.387778780781446E-17 1.387778780781446E-17 1.387778780781446E-17 
                   1.387778780781446E-17 1.387778780781446E-17 1.387778780781446E-17 
                  ]

In the above example, a 2× 3 matrix is defined. The type float64[] indicates that the corresponding value is a matrix of float64 elements (in RangeLab, all the elements of a vector or matrix must have the same type).

Matrices are defined between brackets ’[’ and ’]’. The elements of the same line are separated by spaces and the lines are separated by semi-colons ’;’. Vectors are one-dimentional matrices. It is also possible to force Rangelab to omit to display the result of a command by terminating the line by a semi-colon. This is sometimes useful to avoid to display large matrices.

-] b=[-2_0,4_0]

b = 
     fixed(3,1): [-2_,4_]
          error: 0.000000000000000E-1

-] M=[b b;b b;b b];
-]

RangeLab also displays simplified forms of the matrices when either the floating/fixed-point values (fij) or the error terms (εij) or both are singletons.

-] x=1.0/3.0

x = 
     float64: 3.333333333333333E-1
       error: [1.850371707708594E-17,1.850371707708595E-17]

-] [x x x]

ans = 
       float64[]: 3.333333333333333E-1 3.333333333333333E-1 3.333333333333333E-1 
                  
          Errors: [1.850371707708594E-17 1.850371707708594E-17 1.850371707708594E-17 
                  ,
                   1.850371707708595E-17 1.850371707708595E-17 1.850371707708595E-17 
                  ]

3.2  Special Matrices

RangeLab offers a few functions which return special matrices:

3.3  Sub-Matrices

RangeLab offers a certain number of syntactic elements to describe blocks of matrices. In this section, we consider that a 2× 3 matrix m has been defined:

-] m=rand(2,3)

m = 
     float64[]: 8.256695256490975E-1 6.983298051515964E-1 4.187619622237631E-1 
                3.639862792660326E-1 9.389103586922102E-1 6.938722786960372E-1 
                
        Errors: 0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1 
                0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1 

An element is addressed by m(i,j) and m(i,:) and m(:,j) represent the ith line and jth column of m, respectively. Line and column numbering starts at 1.

-] m(2,2)

ans = 
       float64: 9.389103586922102E-1
         error: 0.000000000000000E-1

-] m(2,:)

ans = 
       float64[]: 3.639862792660326E-1 9.389103586922102E-1 6.938722786960372E-1 
                  
          Errors: 0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1 
                  

-] m(:,2)

ans = 
       float64[]: 6.983298051515964E-1 
                  9.389103586922102E-1 
                  
          Errors: 0.000000000000000E-1 
                  0.000000000000000E-1 

It is also possible to refer to a contiguous set of lines or columns by k:l.

-] m(1,2:3)

ans = 
       float64[]: 6.983298051515964E-1 4.187619622237631E-1 
                  
          Errors: 0.000000000000000E-1 0.000000000000000E-1 
                  
-] m(1:2,2)=77

m = 
     float64[]: 8.256695256490975E-1 7.700000000000000E1 4.187619622237631E-1 
                3.639862792660326E-1 7.700000000000000E1 6.938722786960372E-1 
                
        Errors: 0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1 
                0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1 
                
-] m(1:2,1:2)=m(1:2,2:3)

m = 
     float64[]: 7.700000000000000E1 4.187619622237631E-1 4.187619622237631E-1 
                7.700000000000000E1 6.938722786960372E-1 6.938722786960372E-1 
                
        Errors: 0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1 
                0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1

Finally, the expression m(k:l) denotes the lines k to l of m.

-] m(1:2)

ans = 
       float64[]: 7.700000000000000E1 4.187619622237631E-1 4.187619622237631E-1 
                  7.700000000000000E1 6.938722786960372E-1 6.938722786960372E-1 
                  
          Errors: 0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1 
                  0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1 

3.4  Operations on Matrices

In RangeLab, the operators ’+’, ’-’, ’.*’ and ’./’ perform element by element sums, subtractions, products and divisions. It is also possible to add, subtract, multiply or divise a matrix by a scalar value with the operators ’+’, ’-’, ’*’ and ’/’. The operation is performed element by element.

-] 2*eye(3)+ones(3)

ans = 
       float64[]: 3.000000000000000E0 1.000000000000000E0 1.000000000000000E0 
                  1.000000000000000E0 3.000000000000000E0 1.000000000000000E0 
                  1.000000000000000E0 1.000000000000000E0 3.000000000000000E0 
                  
          Errors: 0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1 
                  0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1 
                  0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1 
                  

-] ans .* ans

ans = 
       float64[]: 9.000000000000000E0 1.000000000000000E0 1.000000000000000E0 
                  1.000000000000000E0 9.000000000000000E0 1.000000000000000E0 
                  1.000000000000000E0 1.000000000000000E0 9.000000000000000E0 
                  
          Errors: 0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1 
                  0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1 
                  0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1

Applied to vectors or matrices, the operator ’*’ performs a scalar product or a matrix multiplication. However, in this case, RangeLab does not bound the errors on the elements of the resulting vector or matrix because there are many way to compute them (the addition being not associative in floating-point or fixed-point arithmetic). Let a=(aij) be a m× n matrix, let b=(bij) be a n× p matrix and let c=a× b. For any cij, 1≤ im, 1≤ jp, all the parsings of the expression cij=k=1n aikbkj should be considered. Instead, RangeLab only considers one parsing : the terms aikbkj are accumulated in increasing order of the index k. This corresponds to:

cij=


(ai1b1j)+(ai2b2j)
+(ai3b3j)
+…
    (6)

RangeLab warns the user about this arbitrary choice:

-] [0.1 0.2 0.3]*[0.4;0.5;0.6]

!!! Scalar product computed using left associativity.
ans = 
       float64: [3.199999999999999E-1,3.200000000000001E-1]
         error: [-1.151856388048599E-16,1.151856388048600E-16]

-] m=ones(2)*6.6

m = 
     float64[]: [6.599999999999999E0 6.599999999999999E0 
                 6.599999999999999E0 6.599999999999999E0 
                ,
                 6.600000000000001E0 6.600000000000001E0 
                 6.600000000000001E0 6.600000000000001E0 
                ]
        Errors: [-8.881784197001252E-16 -8.881784197001252E-16 
                 -8.881784197001252E-16 -8.881784197001252E-16 
                ,
                 8.881784197001253E-16 8.881784197001253E-16 
                 8.881784197001253E-16 8.881784197001253E-16 
                ]

-] m*m

!!! Matrix product computed using left associativity.
ans = 
       float64[]: [8.711999999999999E1 8.711999999999999E1 
                   8.711999999999999E1 8.711999999999999E1 
                  ,
                   8.712000000000002E1 8.712000000000002E1 
                   8.712000000000002E1 8.712000000000002E1 
                  ]
          Errors: [-4.121147867408581E-14 -4.121147867408581E-14 
                   -4.121147867408581E-14 -4.121147867408581E-14 
                  ,
                   4.121147867408582E-14 4.121147867408582E-14 
                   4.121147867408582E-14 4.121147867408582E-14 
                  ]

Finally, RangeLab offers a few additional operations on matrices:

4  Control Structures and Commands

RangeLab’s evaluation of control structures such as conditional or loops differs from the usual one. This is mainly due to the fact that RangeLab values are intervals. The system then returns a superset of the outputs that we can obtain with a standard semantics by selecting scalar entries in the ranges of the input intervals.

4.1  Conditionals

The syntax of a conditional statement in Rangelab is: if cond, stmt else stmt end. The else branch may be omitted, in which case no code is executed if the condition is false: if cond, stmt end.

-] if 1<2, a=1 else b=2 end
-] a

ans = 
       int32: 1

The conditions are predicates among numerical expressions, defined using the operators:

== equality = not equal
<= less than or equal         >= greater than or equal
< less than > greater than
& logical and | logical or
^ exclusive or  logical negation

When a condition containing interval variables can evaluate to both true and false, RangeLab proceeds as follows:

  1. The then branch is evaluated in an environment where the intervals corresponding to the variables occuring in the condition are reduced to the set of values which make the condition true.
  2. Similarly, the else branch is evaluated in an environment where the intervals corresponding to the variables occuring in the condition are reduced to the set of values which make the condition false.
  3. The values of the variables after execution of the then and else branches are joined.

For example, let us consider the following example:

-] x=[1,10]

x = 
     int32: [1,10]

-] clear a
-] clear b
-] if x<5, a=x else b=x end
-] a

ans = 
       int32: [1,4]

-] b

ans = 
       int32: [5,10]

-] x

ans = 
    int32: [1,10]

In the above example, x=[1,10] and the condition x<5 is interpreted as follows:

  1. The environment for the then branch is computed by intersecting the domain of x with the domain where the condition is true: xx ∩ ]−∞,4]. The evaluation of the branch then computes a=[1,4].
  2. The environment for the else branch is computed by intersecting the domain of x with the domain where the condition is false: xx ∩ [5,+∞[. The evaluation of the branch then computes b=[5,10].
  3. The resulting environments are joined, we have, after execution of the conditional statement: x=[1,10], a=[1,4] and b=[5,10].

Note that this is a superset of the results that we can obtain by chosing a scalar value for x in the interval [1,10]. It corresponds to a generalization to control structures of the property given in Equation (5) for expressions. This technique may also lead to over-approximations:

-] clear a
-] if x<5, a=0 else a=5 end
-] a

ans = 
       int32: [0,5]

In this example the result for a is the interval [0,5] while the best anwser would be the set {0,5}.

When a condition involves floating-point or fixed-point value (f,ε), only f can be reduced, the error term ε being left unchanged. This is because the comparison only concerns the floating-point or fixed-point numbers themselves, not the errors attached to them.

-] a=[0.0,10.0,-0.5,0.5]

a = 
     float64: [0.000000000000000E-1,1.000000000000000E1]
       error: [-5.000000000000000E-1,5.000000000000000E-1]


-] clear x
-] clear y
-] if a<5.0, x=a else y=a end
-] x

ans = 
       float64: [0.000000000000000E-1,5.000000000000000E0]
         error: [-5.000000000000000E-1,5.000000000000000E-1]

-] y

ans = 
       float64: [5.000000000000000E0,1.000000000000000E1]
         error: [-5.000000000000000E-1,5.000000000000000E-1]

RangeLab allows a few syntactic variant of the if statement. The then and else branches may contain sequences of instructions instead of single statements. In this case, the statements of the block are separted by semicolons ’;’. Also, the elseif keyword makes it possible to define nested conditionals.

-] if x>3.0, a=0 ; b=0 else a=1; b=1 end
-] x=[2.0,3.0]

x = 
     float64: [2.000000000000000E0,3.000000000000000E0]
       error: 0.000000000000000E-1

-] if x<1.0, a=2.0 elseif x<2.0, a=5.0, elseif x<3.0, a=10.0, else a=15.0 end
-] a

ans = 
       float64: [1.000000000000000E1,1.500000000000000E1]
         error: 0.000000000000000E-1

4.2  While Loops

The syntax of while loops in RangeLab is: while cond, stmts end, where stmts stands either for a single instruction or for a sequence of instructions separated by semicolons. The condition cond is interpreted as the conditions of the if statements and the semantics of the while in an environment env0 is:

  1. The entries are restricted to the sub-intervals which make the condition true. This yields an environment env1 which is used to execute the first iteration of the loop.
  2. The environment env2 containing the values of the variables at the end of the first iteration is joined to env1 into a new environment env3. Step 1 is then repeated with the values of env3 as entries. The process is repeated until env1=env2 (a fixed-point is reached).
  3. The environment at the exit of the loop is defined by joining env0 and env3 where env0 and env3 denote the environments env0 and env3 in which the values of the variables are restricted to the sub-intervals which make the condition false.

This is illustrated by the following example.

-] a=0; x=0;
-] while a<10, a=a+1; x=a end
-] a

ans = 
       int32: 10

-] x

ans = 
       int32: [1,10]

The values of a at each iteration are:

iterationenv1(a)env2(a)env3(a)
101[0,1]
2[0,1][1,2][0,2]
3[0,2][1,3][0,3]
9[0,8][1,9][0,9]
10[0,9][1,10][0,10]
11[0,9][1,10][0,10]

At the end of the loop, the value of a is intersected with the values which invalidate the condition, we obtain a=[0,10]∩[10,+∞[. Since x is not constrained by the condition, its value at the end of the loop corresponds to the value of the last assignment x=a.

Another particularity is that RangeLab allows the user to compute the result of a while loop without entirely executing it. In this case, the system computes a conservative extrapolation of the result using a kind of convergence acceleration technique (called widening in abstract interpretation). The results are always correct but less accurate. For example, let us consider again the preceding example with a larger number of iterations:

-] a=0; x=0;
-] while a<10000, a=a+1; x=a end
-] a

ans = 
       int32: 10000

-] x

ans = 
       int32: [1,16384]

Here, RangeLab has not executed the 10000 iterations but only a limited number (50 by default). Then, each interval growing between two consecutive iteration has been widened. The widening in RangeLab consists of substituting the next power of 2 to a growing interval bound.

iterationenv1(a)env2(a)env3(a)
101[0,1]
2[0,1][1,2][0,2]
49[0,48][1,49][0,49]
50[0,49][1,50][0,64]
51[0,64][1,65][0,128]
52[0,128][1,129][0,256]
57[0,8192][1,8193][0,16384]

In this example, after 57 iterations, the fixed-point is reached. As in the preceding example, RangeLab computes that a=10000 thanks to the constraints given by the while condition. The other variable remains to [0,16384].

When floating-point or fixed-point computations are performed inside a while loop, both the values and error terms are widened. However, as for the if construct, the error terms are not constrained by the condition of the loop. Let us first consider two examples of unstable computations.

-] x=1.0;
-] while x<1000.0, x=x*[1.2,1.3] end
-] x

ans = 
       float64: [1.000000000000000E3,1.300000000000000E3]
         error: [-oo,+oo]

-] x=0.0; i=0;  
-] while i<1000, i=i+1; x=x+0.1 end
-] x

ans = 
       float64: [9.999999999999999E-2,+oo]
         error: [-oo,+oo]

In the first case, RangeLab states that the floating-point value of x at the end of the loop belongs to [1000,1300]. In the second case, the system computes x∈[9.999999999999999· 10−2,+∞[ after the while. In both cases, RangeLab states that the roundoff error belongs to ]−∞,+∞[. This means that the difference between the computed value and the value that we would obtain in infinite precision is arbitrarily large. RangeLab has extrapolated that, in these examples, the roundoff error goes to infinity as the number of iterations goes to infinity.

Let us consider now a stable computation:

-] i=0; n=10; x=1.0; k=[0.8,0.9];
-] while i<n, i=i+1; x=x*k end
-] x

ans = 
       float64: [0.000000000000000E-1,9.000000000000001E-1]
         error: [-1.776356839400250E-15,1.776356839400251E-15]

This computation is stable, the value of x at the end of the loop being kn with k∈ [0.8,0.9] and n≥ 0. The loop performs only n=10 iterations, but a widening is invoked because x is not constrained by the condition of the loop. RangeLab computes that for any k∈[0.8,0.9] and for any 0≤ i≤ 10, kn belongs to [0.0,0.9]. Hence, by extrapolation, the system has computed that ∪i=0n[0.8,0.9]i goes towards [0.0,0.9] as n goes towards +∞. RangeLab also bounds the errors on x at any iteration by [−1.776356839400250· 10−15,1.776356839400251· 10−15]. In order to avoid the widening and to enforce the system to exactly unroll the loop, we can set the unroll parameter to a higher value than the number of iterations.

-] set unroll 20
-] i=0; n=10; x=1.0; k=[0.8,0.9];                           
-] while i<n, i=i+1; x=x*k end
-] x

ans = 
       float64: [1.073741823999999E-1,3.486784401000002E-1]
         error: [-4.596600877704306E-16,4.596600877704307E-16]

In this session, the body of the while loop has been executed 10 times in a standard way, the values of x at different iterations have not been joined and no widening has been performed. We can remark that the result is more accurate.

The unroll parameter makes RangeLab unroll the body of the loop and this improves the accuracy of the results in several ways [4]. If unroll=n, the loop:

while c, stmts end

is understood by the system as it was:

while c,   

        if ¬ c, exit() else stmts  end
        if ¬ c, exit() else stmts  end

end

If the unroll parameter is greater than the number of iterations then this makes it possible to execute the loop in a standard way. Otherwise, the unrolling makes the system work with n environments and the environment at iteration i is joined with the environment at iteration in. In many cases, this techniques improves significantly the accuracy of the results.

By default unroll=5 and widen=10. This makes the first widening arise after 50 iterations (10 iterations of the unrolled loop).

4.3  For Loops

In RangeLab, for loops obey to the syntax: for id = a:b, stmts end or for id = a:step:b, stmts end. By default, no widening is performed during the execution of a for.

-] a=ones(4,2);
-] for i=1:4, a(i,1)=0.1*i end
-] for i=1:2:4, a(i,2)=i / 7.0 end
-] a

ans = 
       float64[]: [9.999999999999999E-2 1.428571428571428E-1 
                   1.999999999999999E-1 1.000000000000000E0 
                   2.999999999999999E-1 4.285714285714285E-1 
                   3.999999999999999E-1 1.000000000000000E0 
                  ,
                   1.000000000000001E-1 1.428571428571429E-1 
                   2.000000000000001E-1 1.000000000000000E0 
                   3.000000000000001E-1 4.285714285714286E-1 
                   4.000000000000001E-1 1.000000000000000E0 
                  ]
          Errors: [-6.938893903907228E-18 7.930164461608261E-18 
                   -2.775557561562891E-17 0.000000000000000E-1 
                   -4.857225732735059E-17 2.379049338482478E-17 
                   -5.551115123125782E-17 0.000000000000000E-1 
                  ,
                   6.938893903907229E-18 7.930164461608262E-18 
                   2.775557561562892E-17 0.000000000000000E-1 
                   4.857225732735060E-17 2.379049338482479E-17 
                   5.551115123125783E-17 0.000000000000000E-1 
                  ]

The widening of for loops can be enabled or disabled by means of the command set widen_for. The examples below can be compared to the examples of Section 4.2.

-] x=1.0;
-] for i=1:100, x=x*[0.8,0.9] end
-] x

ans = 
       float64: [2.037035976334470E-10,2.656139888758755E-5]
         error: [-3.743573134895952E-19,3.743573134895953E-19]
-] set widen_for
-] x=1.0;
-] for i=1:100, x=x*[0.8,0.9] end
-] x

ans = 
       float64: [0.000000000000000E-1,9.000000000000001E-1]
         error: [-1.776356839400250E-15,1.776356839400251E-15]

It is also possible to define decreasing for loops. If a>b in the statement for id = a:b, stmts end then RangeLab assumes that the step is −1. It is possible to define explicitely a negative step using the syntax for id = a:step:b, stmts end. For example, compare the difference of accuracy of the two loops hereafter.

-] s1=0.0; s2=0.0;
-] for i=0:2:50, s1=s1+0.1*i end
-] s1

ans = 
       float64: [6.499999999999998E1,6.500000000000000E1]
         error: [-5.612177389480166E-14,5.612177389480167E-14]

-] for i=50:-2:0, s2=s2+0.1*i end
-] s2

ans = 
       float64: [6.500000000000000E1,6.500000000000002E1]
         error: [-9.325873406851314E-14,9.325873406851315E-14]

4.4  Functions

Functions in RangeLab are stored in files with a .m extension. The function and its file must have the same name. Inside a file, a function is defined following the syntax:

function res = func_name(arg1,,argn)
stmts

This defines a function of n arguments which returns the value of the identifier res after evaluation of stmts. For example let us consider the file mean.m:

function m = mean(a)

    n = max(size(a));
    s = 0.0;
    for i=1:n, 
     s=s+a(i)
    end

    m = s / n;

The function mean takes one argument, a. The returned value is the value of m at the end of the execution. A function is called by invoking its name from the interactive loop:

-] u=[1.1 2.2 3.3 4.4 5.5];
-] mean(u)

ans = 
       float64: 3.300000000000000E0
         error: [-3.330669073875469E-16,6.883382752675971E-16]

The .m files are read from the current directory which can be changed with the command cd(’dir’).

-] r=rand(1,10);
-] mean(r)
??? Undefined identifier 'mean'.

-] cd('Examples')
-] mean(r)

ans = 
       float64: 5.047699905524450E-1
         error: [1.221245327087672E-16,1.221245327087673E-16]

Functions may return many values. The syntax is: function [res1,,resn] = func_name(arg1,,argn). For example, let us consider the file minmax.m

function [rmin,rmax] = minmax(a)

    n = max(size(a));
    rmin = a(1);
    rmax = a(1);
 
    for i=2:n, 
        if (a(i) < rmin), rmin=a(i) end;
        if (a(i) > rmax), rmax=a(i) end
    end

An example of execution is:

-] d=[0.2,0.7]*[5 4 3 2 1 9 7 6];
-] minmax(d)

ans = 
       float64[]: [1.999999999999999E-1 1.199999999999999E0 
                  ,
                   7.000000000000001E-1 6.300000000000001E0 
                  ]
          Errors: [-4.996003610813204E-16 -9.436895709313830E-16 
                  ,
                   4.996003610813205E-16 9.436895709313831E-16 
                  ]

Note that, in RangeLab, it is not allowed to define many functions inside a .m file.

In RangeLab functions can be called with intervals as arguments. In this case, the returned values are an over-approximation of the return values that we could obtain by calling the function with simple values taken inside the interval arguments. In other words,

f([
x1
,
x1
],,[
xn
,
xn
])=[
x
,
x
] =⇒ ∀ x1∈ [
x1
,
x1
],…,∀ xn∈ [
xn
,
xn
], f(x1,…,xn)∈ [
x
,
x
]     (7)

For example, let us take the functions f.m and euler.m below:

% file f.m
function y = f(x)

    y=0.8*x+0.05

% file euler.m 
function y = euler(a,b,h,y0)

    x = [a:h:b];
    y = 0*x;

 y(1) = y0;

 for i=2:size(y),
     y(i) = y(i-1) + h*(f(y(i-1)));
 end

The function euler can be called with intervals as arguments. In this case, we obtain a range for the integral of f, computed by euler, on all the interval passed as argument. For example, euler(0.0,2.0,0.25, [1.0,2.0]) computes an over-approximation of the value of the integral of f between 0.0 and 2.0 and for any initial value y0∈[1.0,2.0].

-] euler(0.0,2.0,0.25,[1.0,2.0])

ans = 
       float64[]: [1.000000000000000E0 1.212499999999999E0 1.467499999999999E0 
                   1.773499999999999E0 2.140699999999999E0 2.581339999999999E0 
                   3.110107999999999E0 3.744629599999999E0 4.506055519999999E0 
                  ,
                   2.000000000000000E0 2.412500000000001E0 2.907500000000001E0 
                   3.501500000000001E0 4.214300000000000E0 5.069660000000000E0 
                   6.096092000000000E0 7.327810400000000E0 8.805872480000000E0 
                  ]
          Errors: [0.000000000000000E-1 -3.339342691255353E-16 -7.403799795469012E-16 
                   -1.318251063864295E-15 -2.241984375928041E-15 -3.360356437553946E-15 
                   -4.880806869778098E-15 -6.719591993942231E-15 -9.387316879383434E-15 
                  ,
                   0.000000000000000E-1 3.339342691255354E-16 7.403799795469013E-16 
                   1.318251063864296E-15 2.241984375928042E-15 3.360356437553947E-15 
                   4.880806869778099E-15 6.719591993942232E-15 9.387316879383435E-15 
                  ]

4.5  Other Commands

RangeLab offers a few other commands:

5  System Commands and Parameters

5.1  System Commands

A few system commands can be used to control the system. They are described below.

5.2  Parameters

RangeLab’s parameters concern the data formats and the accuracy of the analysis. The set command is used to display or modify the parameters. In Section 6.4, we illustrate the effects of the parameters on the analysis of the trapeze procedure introduced in Section 1.

-] set

Mantissa size of error terms (error) ............. 256
Mantissa size of error display (error_display) ... 53
Size of the default fractional part (frac) ....... 32
Base for data display (bin or dec)................ Decimal
Error terms display (abs or rel).................. Absolute
Loop unrolling level (unroll) .................... 5
Widening level (widen) ........................... 10
For loops widening (widen_for) ................... Disabled
Slice function calls (slice_fun) ................. Disabled
Interval slice number (slice) .................... 8
Maximal formal expression height (formal) ........ 5
Current Folder ................................... '/Users/mmartel/Documents/Rangelab'

For each line displayed by the set command, the word in parentheses is the keyword to use in the set command in order to modify the parameter. The command set default resets all the parameters to their default value. In the rest of this section, we describe the effect of RangeLab’s parameters on the results of computations.

6  Case Studies

In this section, we introduce four case studies. The first example shows how we can detect data sensitivity problems during the numerical integration of an ODE. The second example concerns the detection of numerical errors in an expensive function. It uses a gaussian elimination procedure. Next, our third example shows how RangeLab handles fixed-point arithmetic. We focus on the evaluation of a polynomial expression using Horner’s rule. Finally, the last example illustrate the effect of parameters on the analysis of the trapeze function introduced in Section 1.

6.1  Runge-Kutta and Data Sensitivity

We have seen in Section 1 how RangeLab can detect numerical errors in our implementation of the trapeze method. In this section, we aim at showing how a sensitivity error is detected in the second order Runge-Kutta method given below.

function [x_res,y_res] = rk2(x0,y0,tmax,h)
  t = [0:h:tmax]; % initialisations
  x = t * 0.0; 
  y = t * 0.0;
  x(1) = x0; 
  y(1) = y0; 
  nstep = tmax / h;
  for n = 1:nstep     % iterations 
    kx1 = h*f1(x(n),y(n));
    kx2 = h*f1(x(n)+kx1,y(n));
    x(n+1) = x(n)+(kx1+kx2)*0.5;
    ky1 = h*f2(x(n),y(n)); 
    ky2 = h*f2(x(n),y(n)+ky1); 
    y(n+1) = y(n)+((ky1+ky2)*0.5);
  end
  x_res = x(nstep); 
  y_res = y(nstep)

We use this rk2 method to solve the system (S) of first order ODEs defined by

(S) : 

=1.4x+1.6y
=1.35x+1.45y
    (8)

We take as initial condition x(0)=y(0)=2.75 and we consider two cases: in the first case we assume that the initial condition is exact while, in the second case, we attach an initial error term e corresponding to a possible rounding of the value: e=[−0.01,0.01]. To speedup the computation time, we set formal=2. We obtain:

-] set formal 2
-] c0 = 2.75;
-] rk2(c0,c0,3.0,0.01)

ans = 
       float64[]: [1.535801972726615E4
                   1.435291699379384E4 
                  ,
                   1.535801972726621E4 
                   1.435291699379390E4 
                  ]
          Errors: [-4.416249024531930E-10 
                   -4.135988688467320E-10 
                  ,
                   4.416249024531931E-10 
                   4.135988688467321E-10 
                  ]
-] c0 = [2.75,2.75,-0.01,0.01]; 
-] rk2(c0,c0,3.0,0.01)

ans = 
       float64[]: [1.535801972726615E4 
                   1.435291699379384E4 
                  ,
                   1.535801972726621E4 
                   1.435291699379390E4 
                  ]
          Errors: [-5.584734446322786E1 
                   -5.219242543239143E1 
                  ,
                   5.584734446322787E1 
                   5.219242543239144E1 
                  ]

The solution is growing but, concerning the floating-point results, we obtain the same value in both executions. However, the error terms are significantly different: RangeLab clearly detects the sensitivity to the initial condition c0. In comparison, interval arithmetic would compute that with c0=[2.74,2.76], the result is a wide interval and many reasons could be invoked like the wrapping effect or that the function is expansive.

6.2  Gaussian Elimination: Expensive Functions and Roundoff Errors

Our second example is a gaussian elimination procedure. In interval arithmetic, when a function is expansive, as it is the case here, the result is usually a large interval and the user cannot know whether numerical errors arose during the evaluation or not. It is only possible to state that the error may be as large as the width of the resulting interval, which is not quite informative. We use the following function which also gives a flavor of the matrix operations supported by RangeLab.

function [x] = gaussel(A,b)
  N = max(size(A));
  for j = 2:N,   % Gaussian elimination
    for i = j:N,
      m = A(i,j-1) / A(j-1,j-1);
      A(i,:) = A(i,:) - A(j-1,:)*m;
      b(i) = b(i) - m*b(j-1);
    end
  end

  x = zeros(N,1);  % Back substitution
  x(N) = b(N) / A(N,N);
  for j = N-1:-1:1,
    x(j) = (b(j)-A(j,j+1:N)*x(j+1:N))/A(j,j);
  end

We take the matrix A and vector b defined by:

A=



u0u
0v0v 
00uv 
0vu0




   b=







u=
([4.54
54
· 10−8,1.54
54
· 10−7], [−3.45… · 10−23,3.45… · 10−23]
v=
([−1.72
72
· 10−7,−5.45
45
· 10−8], [−3.57… · 10−23,3.57… · 10−23])

RangeLab computes the following result.

-] x=[5.0,17.0]/11.0

x = 
     float64: [4.545454545454545E-1,1.545454545454546E0]
       error: [-1.110223024625156E-16,1.110223024625157E-16]

-]  y=[-19.0,-6.0]/11.0

y = 
     float64: [-1.727272727272727E0,-5.454545454545455E-1]
       error: [-1.110223024625156E-16,1.110223024625157E-16]

-] A=[x 0.0 x 0.0; 0.0 y 0.0 y; 0.0 0.0 x y; 0.0 y x 0.0]*0.0000001

A = 
     float64[]: [4.545454545454545E-8 0.000000000000000E-1 4.545454545454545E-8 0.000000000000000E-1 
                 0.000000000000000E-1 -1.727272727272727E-7 0.000000000000000E-1 -1.727272727272727E-7 
                 0.000000000000000E-1 0.000000000000000E-1 4.545454545454545E-8 -1.727272727272727E-7 
                 0.000000000000000E-1 -1.727272727272727E-7 4.545454545454545E-8 0.000000000000000E-1 
                ,
                 1.545454545454546E-7 0.000000000000000E-1 1.545454545454546E-7 0.000000000000000E-1 
                 0.000000000000000E-1 -5.454545454545454E-8 0.000000000000000E-1 -5.454545454545454E-8 
                 0.000000000000000E-1 0.000000000000000E-1 1.545454545454546E-7 -5.454545454545454E-8 
                 0.000000000000000E-1 -5.454545454545454E-8 1.545454545454546E-7 0.000000000000000E-1 
                ]
        Errors: [-3.456408034775562E-23 0.000000000000000E-1 -3.456408034775562E-23 0.000000000000000E-1 
                 0.000000000000000E-1 -3.576725214783275E-23 0.000000000000000E-1 -3.576725214783275E-23 
                 0.000000000000000E-1 0.000000000000000E-1 -3.456408034775562E-23 -3.576725214783275E-23 
                 0.000000000000000E-1 -3.576725214783275E-23 -3.456408034775562E-23 0.000000000000000E-1 
                ,
                 3.456408034775563E-23 0.000000000000000E-1 3.456408034775563E-23 0.000000000000000E-1 
                 0.000000000000000E-1 3.576725214783276E-23 0.000000000000000E-1 3.576725214783276E-23 
                 0.000000000000000E-1 0.000000000000000E-1 3.456408034775563E-23 3.576725214783276E-23 
                 0.000000000000000E-1 3.576725214783276E-23 3.456408034775563E-23 0.000000000000000E-1 
                ]

-] B=[1.0 1.0 1.0 1.0]

B = 
     float64[]: 1.000000000000000E0 1.000000000000000E0 1.000000000000000E0 1.000000000000000E0 
                
        Errors: 0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1 
                

-] gaussel(A,B)

!!! Scalar product computed using left associativity.
ans = 
     float64[]: [-2.042984771573605E8 
                 -5.546531302876484E7 
                 -6.138527354765938E8 
                 -1.673296672306825E8 
                ,
                 2.109099300620420E9 
                 5.115439462304950E8 
                 6.655837563451781E7 
                 1.172588832487311E7 
                ]
        Errors: [-4.038941344594429E-6 
                 -6.651550726194469E-7 
                 -8.316690681233233E-7 
                 -1.223143375128125E-7 
                ,
                 4.038941344594430E-6 
                 6.651550726194470E-7 
                 8.316690681233234E-7 
                 1.223143375128126E-7 
                ]

Firstly, the system warns about the scalar product done in the back substitution. This is because the roundoff errors arising during this computation depend on the order of the operations performed to compute the scalar product. This order is not unique and the error returned by the tool is correct only if the terms of the product are accumulated in a certain order (increasing order of the indexes.)

Secondly, RangeLab indicates that while the result of the computation belongs to a large interval (for instance the width of the interval is about 2.3· 109 in the first dimension), the numerical error is always very small (always less than e=4.04· 10−6 in absolute value). This means that for any matrix with floating-point coefficients taken in the intervals of A, the numerical error on the result of the gaussian elimination is always less than or equal to e in absolute value, which corresponds to the error precision in the IEEE754 double format. In other words, for any matrix with floating-point coefficients taken in the intervals of A, the system ensures that the roundoff error on the result of the gaussian elimination are negligible. This result could not be obtained by the usual interval arithmetic.

6.3  Fixed-point Arithmetic

A fixed number [14] is made of an integer part i and of a fractional part f. Recall that it is written i_f in RangeLab’s syntax (see Section 2.4.3). For an input value, RangeLab always adjusts the size of i to the minimal size needed to represent the value. The size of the fractional part can be specified by the user (a default size is used instead). Concerning the result of an elementary operation, the integer part is extended if needed, in order to avoid overflows, and the size of the fractional part corresponds to the maximal size of the fractional parts of the operands. Operations are rounded towards zero and introduce roundoff errors, just like floating-point numbers do. For example, the following command adds the values 1.1 and 1.2 represented with a fractional part of 16 bits.

-] 1_1#16+1_2#16

ans = 
       fixed(2,16): 2_299987792969
             error: 9.155273438254952E-6

The tool states that the result is encoded on a 2-bits integer part and on a 16-bits fractional part and that the roundoff error on the result is 9.15… · 10−6. RangeLab also handles intervals of fixed-point numbers in the same way than intervals of floating-points numbers.

Let us consider now the implementation of Horner’s method for the evaluation of a polynomial. The function below takes as arguments an array with the coefficients of the polynomial p and the point x and computes p(x).

function y = horner(a,x)
  y = a(1); 
  n = max(size(a));
  for i = 2:n
    y = a(i) + x * y
  end

If we set the default size of the fractions to 16, we obtain:

-] set frac 16
-] a = [[0_8,1_2] [-6_2,-5_8] [10_8,11_2] [-6_2,5_8]];
-] horner(a,[0_0,4_0])

ans = 
 fixed(6,16): [-42_199996948242,54_599975585938]
       error: [-3.051757812500000E-5,6.103515625000000E-5]

The system computes that for the given coefficients the result is always representable without overflow on 6 bits for the integer part and that, with 16 bits of fraction the roundoff error on the result of the evaluation is bounded by 6.10… · 10−5. With various set frac commands we obtain the accuracies given in Table 1.


fracaccuracy
8[−7.8125· 10−3,1.5625· 10−2]
16[−3.0517578125· 10−5,6.103515625· 10−5]
24[−1.192092895507812· 10−7,2.384185791015625· 10−7]
32[−4.656612873077392· 10−10,9.313225746154786· 10−10]
Table 1: Accuracy of Horner’s method for various fraction sizes.

6.4  Tuning the Parameters for the Trapeze Function

In this section, we illustrate the effect of RangeLab’s parameters on the analysis of the trapeze function introduced in Section 1. First we analyze the function with all the default parameters. After a few minuts, RangeLab keeps computing and still do not answer. We interrupt it (control-c), we decrease the formal parameter in order to speed up the runtime and we restart the analysis of the trapeze function.

-] set

Mantissa size of error terms (error) ............. 256
Mantissa size of error display (error_display) ... 53
Size of the default fractional part (frac) ....... 32
Base for data display (bin or dec)................ Decimal
Error terms display (abs or rel).................. Absolute
Loop unrolling level (unroll) .................... 5
Widening level (widen) ........................... 10
For loops widening (widen_for) ................... Disabled
Slice function calls (slice_fun) ................. Disabled
Interval slice number (slice) .................... 8
Maximal formal expression height (formal) ........ 4
Current Folder ................................... '/Users/mmartel/Documents/Rangelab/Examples'

-] trapeze(0.25,50.0,100)
^C
!!! Interrupted.

-] set formal 2
-] set

Mantissa size of error terms (error) ............. 256
Mantissa size of error display (error_display) ... 53
Size of the default fractional part (frac) ....... 32
Base for data display (bin or dec)................ Decimal
Error terms display (abs or rel).................. Absolute
Loop unrolling level (unroll) .................... 5
Widening level (widen) ........................... 10
For loops widening (widen_for) ................... Disabled
Slice function calls (slice_fun) ................. Disabled
Interval slice number (slice) .................... 8
Maximal formal expression height (formal) ........ 2
Current Folder ................................... '/Users/mmartel/Documents/Rangelab/Examples'

-] trapeze(0.25,50,100)

ans = 
       float64: [6.447674308233614E1,+oo]
         error: [-oo,+oo]

We obtain a result where the floating-point value and the error term are infinite. This is correct while very inaccurate. Such over-approximations often come from the analysis of loops that we can tune thanks to the unroll and widen parameters (see Section 4.2). For example, let us ask RangeLab to fully unroll the loop.

-] set unroll 110
-] trapeze(0.25,50.0,100)

 float64: [6.557548333631875E1,1.311509666726376E2]
   error: [-1.456853305716057E-12,1.457731336808933E-12]

The execution lasts a few minuts but we obtain a satisfying result (see Section 1).

7  Conclusion

This manual gives an overview of the features of RangeLab, an interactive tool to detect numerical errors in floating-point and fixed-point computations. More details on the techniques implemented in RangeLab can be found in [8].

RangeLab 1.0 handles the elementary arithmetic operations and it would be interesting that RangeLab 2.0 handles more advanced features. However this raises some theoretical issues. Firstly, to add the usual elementary mathematical functions (sine, logarithms, etc.) the way roundoff errors are propagated inside these functions must be defined carefully. Secondly, we would like solvers, for example for ODEs. However, in order to estimate the error on the results, guaranteed integration techniques should be employed [11]. The same remark holds for other possible extensions, for example concerning matrix computations. But despite the huge work that would represent the implementation of these features, let us insist on the fact that RangeLab is open source and that we encourage other programmers to contribute to the development of the tool.

References

[1]
ANSI/IEEE. IEEE Standard for Binary Floating-point Arithmetic, std 754-2008 edition, 2008.
[2]
P. Cousot and R. Cousot. Abstract interpretation: A unified lattice model for static analysis of programs by construction of approximations of fixed points. In Principles of Programming Languages 4, pages 238–252. ACM Press, 1977.
[3]
D. Goldberg. What every computer scientist should know about floating-point arithmetic. ACM Computing Surveys, 23(1):5–48, 1991.
[4]
E. Goubault, M. Martel, and S. Putot. Static analysis-based validation of floating-point computations. In follow-up of the seminary on Numerical Software with Result Verification, at Dagstuhl, Germany, number 2991 in Lecture Notes in Computer Science, 2004.
[5]
L. Jaulin, M. Kieffer, O. Didrit, and E. Walter. Applied Interval Arithmetics with Examples in Parameter and State Estimation, Robust Control and Robotics. Springer-Verlag, 2001.
[6]
M. Martel. An overview of semantics for the validation of numerical programs. In Verification, Model Checking and Abstract Interpretation, VMCAI’05, number 3385 in Lecture Notes in Computer Science, pages 59–77. Springer-Verlag, 2005.
[7]
M. Martel. Semantics of roundoff error propagation in finite precision calculations. Journal of Higher Order and Symbolic Computation, 19:7–30, 2006.
[8]
M. Martel. RangeLab: a static-analyzer to bound the accuracy of finite-precision computations. 2011.
[9]
R. E. Moore, R. B. Kearfott, and M. J. Cloud. Introduction To Interval Analysis. Cambridge Uni Press (CUP), 2009.
[10]
Jean-Michel Muller. On the definition of ulp (x). Technical Report 5504, INRIA-Rhones-Alpes, 2005.
[11]
N. S. Nedialkov and J. D. Pryce. Solving differential-algebraic equations by taylor series (i): Computing taylor coefficients. BIT, 45(3):561–591, 2005.
[12]
S.M. Rump. INTLAB - INTerval LABoratory. In Tibor Csendes, editor, Developments in Reliable Computing, pages 77–104. Kluwer Academic Publishers, Dordrecht, 1999.
[13]
Torbjorn Granlund and the GMP development team. The GNU Multiple Precision Arithmetic Library, 5.0.2 edition, 2011. http://gmplib.org.
[14]
Randy Yates. Fixed-Point Arithmetic: An Introduction. Digital Signal Labs, July 2009.

This document was translated from LATEX by HEVEA.