Naive Implementations

Credit: The Guardian
Credit: The Guardian

With ES6, the JavaScript math library gained a level in badass with the addition of seventeen new functions. While I was polishing my UsefulJS library for publication I decided I’d spend a couple of fun evenings implementing shims for them since, if I want to calculate the hypotenuse, I do not want to have to do this:

var x = 3, y = 4,
    r = ("hypot" in Math) ? Math.hypot(x, y) :
        Math.sqrt(x*x + y*y);

The C++ documentation for hypot has an interesting statement in the notes:

POSIX specifies that underflow may only occur when both arguments are subnormal and the correct result is also subnormal (this forbids naive implementations).

By naive implementations, they mean something like this:

(function() {
    if (!("hypot" in Math)) {
        Math.hypot = function(x, y) {
            return Math.sqrt(x*x + y*y);
        };
    }
})();

As it happens, this wouldn’t be up to standard anyway since Math.hypot takes any number of arguments but that’s a detail irrelevant to this discussion. Since this is the formula for Pythagoras’s theorem, what’s naive about it? Tell me, computer, what is the hypotenuse when the two shorter sides are 10155 and 1?

That’s why it’s naive: even though the correct result (which is 10155) can be represented, there’s an overflow in the intermediate calculations that causes the wrong value to be returned. Moler and Morrison published a non-naive algorithm for calculating the hypotenuse in 1983. Not only does it avoid underflow and overflow in intermediate products, it also avoids taking a square root (which was presumably an attractive quality back then). Let’s give it a whirl:

var mm = function(x, y) {
    // Use abs values
    if (x < 0) {
        x = -x;
    }
    if (y < 0) {
        y = -y;
    }
    var tmp, r, s;
    // Ensure that x >= y
    if (x < y) {
        tmp = x;
        x = y;
        y = tmp;
    }
    do {
        tmp = y / x;
        r = tmp * tmp;
        s = r / (4 + r);
        x = x + 2 * s * x;
        y *= s;
    }
    // While y is significant
    while (1 + y > 1);
    return x;
};

Tell me, computer, what, according to the Moler-Morrison algorithm, is the hypotenuse when the two shorter sides are 3 and 4?

Damn and blast! Everyone knows that’s 5! Can we do better? Let’s rearrange the formula a little starting from the familiar form:
r = \sqrt{x^2+y^2}
The sum inside the square root can be expressed as:
x^2(1 + (y/x)^2)
Why should we care? Because we can take x2 outside the square root to give:
r = |x|\sqrt{1 + (y/x)^2}
As long as we ensure that x >= y, (y/x)2 will not overflow. It might underflow, but then the result is simply |x| which to all intents and purposes is correct. Let’s give it a try:

var hypot = function(x, y) {
    var tmp, r;
    if (x < 0) {
        x = -x;
    }
    if (y < 0) {
        y = -y;
    }
    if (x < y) {
        tmp = x;
        x = y;
        y = tmp;
    }
    if (!x) {
        // hypot(0, 0)
        return 0;
    }
    r = y / x;
    return x * Math.sqrt(1 + r * r);
};

The only gotcha here is to avoid the division 0/0 which would result in NaN, but otherwise it’s straightforward. Now tell me, computer, what is the hypotenuse when the two shorter sides are:

  • 3 and 4
  • 0 and 0
  • 10155 and 10155

Works like a charm! You may well ask whether there's any point in being able to calculate the hypotenuse of extremely large numbers. I can't see myself ever needing to do it yet I went to the trouble of making it possible. The reason is simple: pride. I do not want someone looking at my code and thinking "n00b!".

Several of the new Math functions, such as expm1 and log1p, can give the wrong result if they're not carefully implemented. Tell me, computer, does your Math library contain naive implementations?

If you're seeing "Dunno", your browser doesn't support ES6 Math functions. "Looks like it" means that the area hyperbolic sine (asinh) of 10-16 was given as 0. To machine precision, this should be 10-16. In fairness to the engineers, the deficiency probably lies in the underlying C libraries. "Not as far as I can tell" means that this particular Math function at least has been carefully implemented. The difference between 0 and 10-16 may seem trivial but it's a 100% loss of accuracy, an example of what is known as catastrophic cancellation.

As to how naive or not my own math code is, see for yourself! (Math code starts around line 1250).

Leave a Reply

Your email address will not be published. Required fields are marked *