Floating point expression inspector

The following tool lets you inspect the computational flow of a floating point expression, seeing where rounding occurs, when exceptions are triggered, when precision may be lost, when special values propagate, when error accumulates, and other floating point headaches.

Example

[fpinspect]# ./fpinspect "sqrt(45.0*e+phi)/pi"
Exception: 0 (1 roundings) INEXACT (45.000000 * e)
  Trace (1 operations) MUL

Exception: 0 (1 roundings) INEXACT phi
  Trace (1 operations) MUL

Exception: 0 (2 roundings) INEXACT ((45.000000 * e) + phi)
Exception: 1 (2 roundings) INEXACT ((45.000000 * e) + phi)
  Trace (2 operations) MUL ADD

Exception: 0 (3 roundings) INEXACT sqrt(((45.000000 * e) + phi))
Exception: 1 (3 roundings) INEXACT sqrt(((45.000000 * e) + phi))
Exception: 2 (3 roundings) INEXACT sqrt(((45.000000 * e) + phi))
  Trace (3 operations) MUL ADD ADD

Exception: 0 (3 roundings) INEXACT pi
Exception: 1 (3 roundings) INEXACT pi
Exception: 2 (3 roundings) INEXACT pi
  Trace (3 operations) MUL ADD ADD

Exception: 0 (4 roundings) INEXACT (sqrt(((45.000000 * e) + phi)) / pi)
Exception: 1 (4 roundings) INEXACT (sqrt(((45.000000 * e) + phi)) / pi)
Exception: 2 (4 roundings) INEXACT (sqrt(((45.000000 * e) + phi)) / pi)
Exception: 3 (4 roundings) INEXACT (sqrt(((45.000000 * e) + phi)) / pi)
  Trace (4 operations) MUL ADD ADD DIV

(sqrt(((45.000000 * e) + phi)) / pi)
        ans: 3.54370117187500
        err: 0.00000126456894

As you can see, the expression sqrt(45.0*e+phi)/pi produces a lot of output, each empty-line-separated region is a subexpression which triggered an exception, in this case because 45 * e is an inexact value, the inexact exception is presented first. Here you can see that such an expression involved 1 operations, total and in this case the operation is just a MUL. We can also see that the resulting expression, because it's inexact, incurred one rounding.

Following down the exception list, we can see that the exception propagated to phi in a MUL (which is also an inexact value), and continued, with each new inexact subexpression resulting in several roundings. Since kernels like sqrt might themselves use operations like add, we also see the final group of exceptions contains an additional ADD in it's trace.

The final result of the expression is given in ans: and below that you will find the accumulative error err: of evaluating that expression, in this case this function is exact to five mantissa digits of precision, out of a total of seven, which means this expression has ~0.71 ULP of error.

Documentation

Run the program with no expression or -h to see the options.

Here's some constants and functions available for use in expressions.

Constants

• e
• pi
• phi

Functions

• floor
• ceil
• trunc
• sqrt
• abs
• min
• max
• copysign

How it works

This program implements IEEE-754 floating point completely in software, emulating all rounding modes, exceptions, and tininess detection methods which can be configured when evaluating an expression. With exception to transcendental functions, all floating point computation is also accurate to <= 1 ULP of error.

Currently there is support for 32-bit single-precision floating-point float32.{h,c} and 64-bit double-precision floating-point float64.{h,c}, as double-precision is necessary for 32-bit single-precision kernels kernel32.{h,c} to produce correctly rounded and truncated results to <= 1 ULP of error.

64-bit double-precision floating-point makes use of 128-bit modular arithmetic implemented in uint128.{h,c}

Accumulative error accounting is handled by real32.{h,c} and real64.{h,c} for single-precision and double-precision floating-point, respectively.

NOTE:

There are currently no 64-bit kernels, as that would require either 80-bit extended-precision floating-point, or 128-bit quadruple-precision floating-point to be implemented in software to have the precision necessary to produce correctly rounded and truncated results to <= 1 ULP of error.