The BoundRoundingError package

This code provides tools to help analyse floating-point rounding errors in simple algorithms.

Latest version

BoundRoundingError 0.3 (May 2024).

Once downloaded, put the folder somewhere in your home directory.

If BoundRoundingError_path is that location, set the variable libname so that Maple knows how to find the package, by libname:="BoundRoundingError_path",libname:

It is convenient to add that line to your .mapleinit file (creating it if it does not exist) so that Maple will find the correct version of the package in your future sessions.

You can then check that the installation worked by asking BoundRoundinngError:-version;

This command should return the number above. (You may have to restart Maple for your changes to be taken into account.)

The source code can be read by anyone who can read Maple code, but if you want to use the package, then it’s better to download it with the link above.

Reference for the package

This package was designed to illustrate the approach in

J.-M. Muller and B. Salvy, “Effective Quadratic Error Bounds for Floating-Point Algorithms Computing the Hypotenuse Function,” arXiv, May 2024.

A demo file can be found with the arXiv version of the article.

Example

With input the following description of a program computing the hypothenuse

Algo:=[
   Input(x=0..2^16,y=0..2^16,_u=0..1/4),
   sx=RN(x^2),
   sy=RN(y^2),
   sigma=RN(sx+sy),
   rho=RN(sqrt(sigma))];

The command

BoundRoundingError(Algo);

produces the output \[ 2 \textit{_u} +\left(\frac{72}{5}-\frac{32 \sqrt{6}}{5}\right) \textit{_u}^{2} \] which is a bound on the relative error \(\rho/\sqrt{x^2+y^2}\) when \(\rho\) is computed by this program for \((x,y)\in[0,2^{16}]^2\). It is assumed that no overflow or underflow occurs, that the floating-point system is binary with precision \(p\ge2\). The value \(\textit{_u}\) in the output stands for \(\textit{_u}=2^{-p}\), so that the output is parameterized by the precision.

Syntax

The algorithms are written in the form of a list [input, assignments] where input is a statement of the form Input(var1=range1,var2=range2,...) with the names of the (real) arguments used by the program and (optionally) a range where these arguments are confined. Examples: Input(x,y) Input(x=0..1,y) Input(x=0..infinity,y=-1..1)

The assignments are of the form name = RN(args, opt_errtype) where

  • name is the name of a local variable
  • RN stands for “Round to Nearest”
  • args is expressed in terms of constants or variables defined previously. It is one of:
    • a constant;
    • a variable;
    • the square root of a variable;
    • a sum, product, difference, quotient of two constants or variables;
    • a multiply-add: a*x+b, where (a,x,b) are constants or variables.
  • opt_errtype is optional and indicates an upper bound on the error made during this rounding. It can be one of
    • 0: rounding is exact
    • absolute(expr): a bound on the absolute error, whose value is expr, usually expressed in terms of _u (see below);
    • relative(expr) with expr one of
      • “Knuth”: relative error \(\textit{_u}/(1+\textit{_u})\)
      • “JRdiv”: Jeannerod-Rump bound on relative error for division
      • “JRsqrt”: Jeannerod-Rump bound on relative error for square root.

The variable _u stands for ulp. It is expected to belong to an interval of the form

_u=0..MAXU where MAXU is a power of 2, with default value 1/64. A different value like 1/4 is obtained by adding _u=0..1/4 to the set of input variables. Because we use the bounds by Jeannerod & Rump, _u cannot get larger than 1/4.

Use

The package provides just one command BoundRoundingError. Its arguments are

  • algo: the algorithm as described above;
  • steps=… (optional): the steps of the analysis that need be performed. The possibilities are
    • “step-by-step”: performs only the analysis of the algorithm without computing any bound;
    • “linear”: computes the linear part (wrt \(\textit{_u}\)) of a bound on the error;
    • “quadratic” (default): computes a quadratic bound on the error.
  • type=… (optional): type of error to be analyzed. It is one of:
    • “relative” (default): relative error
    • “absolute”: absolute error
  • ‘optpoint’ (optional): a name. If that argument is given, the name is assigned the point where the extremum of the bound is reached. This does not mean that the bound is reached when evaluating the code there, but only that our overall algebraic bound on the error reaches its maximal value at that point.

More and more detailed information on the analysis can be obtained by giving a value between 1 and 5 to the variable infolevel[BoundRoundingError]. For instance, assigning

infolevel[BoundRoundingError]:=3;

before the analysis shows the individual bounds computed for each of the variables during the step-by-step analysis, as well as the timings of each stage of the derivation.

Bugs

This is still an early prototype. If you see any problem when using it, send me an email at bruno.salvy@inria.fr.