RangeLab User’s Guide |
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={ x↦ a/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:
-] trapeze(0.25,75.0,150) ans = float64: [6.567656427626613E1,1.313531285525323E2] error: [-2.172235958601054E-12,2.165505284708336E-12]The error is greater while the magnitude of the floating-point result is almost identical.
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.
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.
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
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 x• y 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) |
x−y = | ⎛ ⎝ | ∘(fx−fy),ex−ey+ε(fx−fy) | ⎞ ⎠ | (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]
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× R : 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:
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((f1,ε1),…,(fn,εn))∈ ([f,f],[ε,ε]). In other words:
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
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 X≤ Y. 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
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]
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≤ i≤ m, 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 ]
RangeLab offers a few functions which return special matrices:
-] eye(3) ans = float64[]: 1.000000000000000E0 0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1 1.000000000000000E0 0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1 1.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
-] zeros(2,3) ans = float64[]: 0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1 Errors: 0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1
-] ones(2) ans = float64[]: 1.000000000000000E0 1.000000000000000E0 1.000000000000000E0 1.000000000000000E0 Errors: 0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1
-] rand(3,2) ans = float64[]: 2.952038208008805E-1 4.290434377475177E-1 8.365036245865376E-4 6.252741394976103E-1 1.287955345420474E-1 3.091271919116317E-1 Errors: 0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1
-] rand_errors(2,3) ans = float64[]: 0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1 Errors: 9.866993724538916E-1 3.180417737722820E-1 8.605357203798341E-1 3.340488095401736E-1 4.274791997825566E-1 8.466969162065545E-1
-] linspace(0,100,5) ans = int32[]: 0 25 50 75 100
-] logspace(1,4,3) ans = float64[]: 1.000000000000000E1 5.005000000000000E3 1.000000000000000E4 Errors: 0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1
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
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≤ i≤ m, 1≤ j≤ p, 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:
-] [1 2 3]' ans = int32[]: 1 2 3
-] A=rand(3) A = float64[]: 5.355956450658870E-1 3.929013568034833E-1 3.834966181141171E-1 5.176696778695991E-1 4.457960304114821E-1 3.248969022202455E-1 5.156534544293233E-1 7.798839675067862E-1 6.856899927903758E-1 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 -] diag(A) ans = float64[]: 5.355956450658870E-1 4.457960304114821E-1 6.856899927903758E-1 Errors: 0.000000000000000E-1 0.000000000000000E-1 0.000000000000000E-1 -] diag(A,2) ans = float64[]: 3.929013568034833E-1 3.248969022202455E-1 Errors: 0.000000000000000E-1 0.000000000000000E-1
-] size(ones(3,5)) ans = int32[]: 3 5 -] size([1 1 1]) ans = int32[]: 1 3
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.
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:
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:
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
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:
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:
iteration | env1(a) | env2(a) | env3(a) |
1 | 0 | 1 | [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.
iteration | env1(a) | env2(a) | env3(a) |
1 | 0 | 1 | [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:
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 i−n. 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).
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]
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:
fun | ction 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∈ [ |
| , |
| ],…,∀ xn∈ [ |
| , |
| ], f(x1,…,xn)∈ [ |
| , |
| ] (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 ]
RangeLab offers a few other commands:
-] abs(-2.0) ans = float64: 2.000000000000000E0 error: 0.000000000000000E-1
-] sqrt(2.0) ans = float64: 1.414213562373095E0 error: [-1.110223024625156E-16,1.110223024625157E-16]
-] a=[5 4 3 2 1 8 7 6]; -] min(a) ans = int32: 1 -] max(a) ans = int32: 8
A few system commands can be used to control the system. They are described below.
-] who x a i j b y -] whos x float64 a float64[3,3] i int32 j int32 b float64[5,10] y float64
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.
-] 2_25 ans = fixed(2,2): 2_25 error: 0.000000000000000E-1 -] 2_1 ans = fixed(2,32): 2_0999999998603 error: [0.000000000000000E-1,4.656612873077393E-10]
-] set bin -] 2.5 ans = float64: 0 10000000000 0100000000000000000000000000000000000000000000000000 error: 0.000000000000000E-1 -] set dec -] ans ans = float64: 2.500000000000000E0 error: 0.000000000000000E-1
-] set abs -] 0.1 ans = float64: [9.999999999999999E-2,1.000000000000001E-1] error: [-6.938893903907228E-18,6.938893903907229E-18] -] set rel -] ans ans = float64: [9.999999999999999E-2,1.000000000000001E-1] error: [0.000000000000000E-1,6.938893903907229E-17]%
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.
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) : | ⎧ ⎨ ⎩ |
| (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.
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= | ⎛ ⎜ ⎜ ⎜ ⎝ |
| ⎞ ⎟ ⎟ ⎟ ⎠ | b= | ⎛ ⎜ ⎜ ⎜ ⎝ |
| ⎞ ⎟ ⎟ ⎟ ⎠ |
|
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.
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.
frac accuracy 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]
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).
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.
This document was translated from LATEX by HEVEA.