Herbie Tutorial, Part 2

See the main page for more info on Herbie.

Part 1 of this tutorial described how Herbie can be used to automatically rewrite floating point expressions, to make them more accurate. Part 1 focused on running Herbie and reading its results; this Part 2 will instead work through applying Herbie to a realistic program.

Finding numerical expressions

As an example realistic program, we'll use math.js, an extensive math library for JavaScript. In particular, we'll walk through bug 208, which found inaccuracy in the implementation of complex square root; for a full write-up of the bug itself, check out this blog post by one of the Herbie authors.

To use Herbie, you first need to find some floating-point expressions to feed to Herbie. In the case of math.js, the floating-point expressions of interest are the various functions that compute mathematical functions; in your code, there's a good chance that a small core of your application does the mathematical computations, and the rest sets up parameters, handles control flow, visualizes or print results, and so on. The mathematical core is what Herbie will be interested in.

For example, in the case of math.js, the mathematical core is in lib/function/. Each file in each subdirectory contains a collection of mathematical functions, each of which is potentially inaccurate. You can start by sending all of them into Herbie, or only the most important ones. Here, let's look at just the file arithmetic/sqrt.js, which contains real and complex square roots. In full, the code of interest is:

math.sqrt = function sqrt (x) {
  if (arguments.length != 1) {
    throw new math.error.ArgumentsError('sqrt', arguments.length, 1);
  }

  if (isNumber(x)) {
    if (x >= 0) {
      return Math.sqrt(x);
    }
    else {
      return sqrt(new Complex(x, 0));
    }
  }

  if (isComplex(x)) {
    var r = Math.sqrt(x.re * x.re + x.im * x.im);
    if (x.im >= 0) {
      return new Complex(
          0.5 * Math.sqrt(2.0 * (r + x.re)),
          0.5 * Math.sqrt(2.0 * (r - x.re))
      );
    }
    else {
      return new Complex(
          0.5 * Math.sqrt(2.0 * (r + x.re)),
          -0.5 * Math.sqrt(2.0 * (r - x.re))
      );
    }
  }

  if (x instanceof BigNumber) {
    if (x.isNegative()) {
      // negative value -> downgrade to number to do complex value computation
      return sqrt(x.toNumber());
    }
    else {
      return x.sqrt();
    }
  }

  if (isCollection(x)) {
    return collection.deepMap(x, sqrt);
  }

  if (isBoolean(x) || x === null) {
    return sqrt(+x);
  }

  throw new math.error.UnsupportedTypeError('sqrt', math['typeof'](x));
};

Extracting expressions

The code above is complex, with argument checks, dispatching over five possible types, and error handling. Herbie does not handle complex data structures (only floating-point values), so we'll want to break up the code above into multiple inputs, one for each type of data structure. Let's look at the isComplex(x) case:

var r = Math.sqrt(x.re * x.re + x.im * x.im);
if (x.im >= 0) {
  return new Complex(
      0.5 * Math.sqrt(2.0 * (r + x.re)),
      0.5 * Math.sqrt(2.0 * (r - x.re))
  );
}
else {
  return new Complex(
      0.5 * Math.sqrt(2.0 * (r + x.re)),
      -0.5 * Math.sqrt(2.0 * (r - x.re))
  );
}

This code contains a branch: one option for non-negative x.im, and one for positive x.im. While Herbie supports an if construct, it's usually better to encode branches as separate inputs to Herbie.

Finally, each branch access fields of a data structure (x is of type Complex) and constructs new data structures. Since Herbie does not understand complex data structures, we must write each floating-point value used in constructing the final output as its own test case.

So, this isComplex(x) case would become four inputs to Herbie: a real and an imaginary part, for each negative or non-negative x.im. Note that the r variable is computed outside the branch; it will have to be duplicated in each input. Each input is a single floating-point expression that computes a single floating-point output without branches, loops, or data structures:

var r = Math.sqrt(xre * xre + xim * xim);
0.5 * Math.sqrt(2.0 * (r + xre)),
// xim ≥ 0
var r = Math.sqrt(xre * xre + xim * xim);
0.5 * Math.sqrt(2.0 * (r - xre)),
// xim ≥ 0
var r = Math.sqrt(xre * xre + xim * xim);
0.5 * Math.sqrt(2.0 * (r + xre)),
// xim < 0
var r = Math.sqrt(xre * xre + xim * xim);
-0.5 * Math.sqrt(2.0 * (r - xre)),
// xim < 0

Note that x.im and x.re have changed to xim and xre. This is to emphasize that there is no longer an x structure, just its floating-point fields. The comment below each case reminds us what the bounds on the input variables are.

Translating to Herbie's input language

Now that we have simple floating-point expressions, we can translate them to Herbie's input language. The input language is a variant of Scheme. For each input, we will write a herbie-test declaration; each declaration will have a list of input variables, a description of the input, and the floating-point expression itself. Here's how we'd start for the first input above:

(herbie-test (xim xre)
    "arithmetic/sqrt.js, isComplex(x), x>=0"
    ?)

The question mark is what we will fill in with the expression. But before we do that, take a look at the other fields. The variables are specified as (xim xre), which tells Herbie that there are two input variables named xim and xre; the parentheses are mandatory. r isn't on that list, because even though it is a variable in the code above, it's not an input variable: it's not an argument to our code, but just a value computed internally.

Now, we must translate the expression itself. We can define r with a let* expression:

(herbie-test (xim xre)
    "arithmetic/sqrt.js, isComplex(x), x>=0"
    (let* ([r (sqrt (+ (sqr xre) (sqr xim)))])
      ?))

Note the peculiar syntax of let*; the first argument is a list of square-bracketed binders, each of which has a variable name (like r) and an expression to bind that variable to (here, (sqrt (+ (sqr xre) (sqr xim)))). There's only one binder here, but you could have more if you wanted.

Inside the body of the let*, which is its second argument (the question mark), you can write another expression which can refer to any of the bound variable names. We'll translate the second line there:

(herbie-test (xim xre)
    "arithmetic/sqrt.js, isComplex(x), x>=0"
    (let* ([r (sqrt (+ (sqr xre) (sqr xim)))])
      (* 0.5 (sqrt (* 2.0 (+ r xre))))))

Translating expressions is not too hard—Herbie understands many common mathematical functions, and even has shortcuts, such as sqr for squaring numbers.

The final step is to add our input bound, xim ≥ 0. You do this by changing the variable declaration in the first argument to herbie-test. Instead of just writing xim, write a binder, which has the variable name xim and a distribution to sample xim from:

(herbie-test ([xim (>= default 0)] xre)
    "arithmetic/sqrt.js, isComplex(x), x>=0"
    (let* ([r (sqrt (+ (sqr xre) (sqr xim)))])
      (* 0.5 (sqrt (* 2.0 (+ r xre))))))

We use the distribution (>= default 0), which means to sample values from default and only keep them if they are greater than 0. You usually want to use default as the input distribution, since it specifically tries very large and very small inputs in an effort to find inaccurate inputs; but there are other distributions as well, including integer for 32-bit integers and (uniform a b) for uniformly-distributed real values.

This finishes our first input to Herbie. We can translate the other four cases at this point, or go ahead with the first case. For the sake of the tutorial, let's move ahead with one input for now.

Running Herbie

Running Herbie is exactly like before: