Math.fma - fused multiply–add (FMA) or fused multiply–accumulate (FMAC) in JavaScript

Math.fma.js

// a * b + c
// Math.fma - fused multiply–add (FMA) or fused multiply–accumulate (FMAC) in JavaScript
// For implementation details, see "The Handbook of Applied Cryptography"
// http://www.cacr.math.uwaterloo.ca/hac/about/chap14.pdf
// and
// https://github.com/JuliaLang/openlibm/blob/master/src/s_fma.c
//!? BUGGY

(function () {
  "use strict";

  var SPLIT = Math.pow(2, 27) + 1;
  var MIN_VALUE = Math.pow(2, -1022);
  var EPSILON = Math.pow(2, -52);
  // (1022 + 52) / 3 < C <= (1022 - 53 - 53 + 4) / 2 - ?
  var C = 416;
  var A = Math.pow(2, +C);
  var B = Math.pow(2, -C);

  var multiply = function (a, b) {
    var at = SPLIT * a;
    var ahi = at - (at - a);
    var alo = a - ahi;
    var bt = SPLIT * b;
    var bhi = bt - (bt - b);
    var blo = b - bhi;
    var p = a * b;
    var e = ((ahi * bhi - p) + ahi * blo + alo * bhi) + alo * blo;
    return {
      p: p,
      e: e
    };
  };

  var add = function (a, b) {
    var s = a + b;
    var v = s - a;
    var e = (a - (s - v)) + (b - v);
    return {
      s: s,
      e: e
    };
  };

  var adjust = function (x, y) {
    return x !== 0 && y !== 0 && SPLIT * x - (SPLIT * x - x) === x ? x * (1 + (x < 0 ? -1 : +1) * (y < 0 ? -1 : +1) * EPSILON) : x;
  };

  Math.fma = function (x, y, z) {
    x = Number(x);
    y = Number(y);
    z = Number(z);

    if (x === 0 || x !== x || x === +1 / 0 || x === -1 / 0 ||
        y === 0 || y !== y || y === +1 / 0 || y === -1 / 0) {
      return x * y + z;
    }
    if (z === 0) {
      return x * y;
    }
    if (z !== z || z === +1 / 0 || z === -1 / 0) {
      return z;
    }

    var scale = 1;
    while (Math.abs(x) > A) {
      scale *= A;
      x *= B;
    }
    while (Math.abs(y) > A) {
      scale *= A;
      y *= B;
    }
    if (scale === 1 / 0) {
      return x * y * scale;
    }
    while (Math.abs(x) < B) {
      scale *= B;
      x *= A;
    }
    while (Math.abs(y) < B) {
      scale *= B;
      y *= A;
    }
    if (scale === 0) {
      return z;
    }

    var xs = x;
    var ys = y;
    var zs = z / scale;

    if (Math.abs(zs) > Math.abs(xs * ys) * 4 / EPSILON) {
      return z;
    }
    if (Math.abs(zs) < Math.abs(xs * ys) * EPSILON / 4 * EPSILON / 4) {
      zs = (z < 0 ? -1 : +1) * MIN_VALUE;
    }

    var xy = multiply(xs, ys);
    var s = add(xy.p, zs);
    var u = add(xy.e, s.e);
    var i = add(s.s, u.s);

    var f = i.s + adjust(i.e, u.e);
    if (f === 0) {
      return f;
    }

    var fs = f * scale;
    if (Math.abs(fs) > MIN_VALUE) {
      return fs;
    }

    // It is possible that there was extra rounding for a denormalized value.
    return fs + adjust(f - fs / scale, i.e) * scale;
  };

}());

Here another implementation, but not passing overflow tests.

export function fma(a, b, c) {
  let LO;

  function twoSum(a, b) {
    let s = a + b;
    let a1  = s - b;
    LO = (a - a1) + (b - (s - a1));
    return s;
  }

  function twoProd(a, b) {
    let t = 134217729 * a;
    let ah = t + (a - t), al = a - ah;
    t = 134217729 * b;
    let bh = t + (b - t), bl = b - bh;
    t = a * b;
    LO = (((ah * bh - t) + ah * bl) + al * bh) + al * bl;
    return t;
  }

  if (!isFinite(a) || !isFinite(b) || !isFinite(c)) {
    return a * b + c;
  }

  if (a === 0 || b === 0 || c === 0) {
    return a * b + c;
  }

  let mhi = twoProd(a, b);
  let mlo = LO;

  let shi = twoSum(mhi, c);
  let slo = LO;
  
  slo += mlo;
  return shi + slo;
}

//check this tests before use, as I remember only 'overflow' fails. So Yaffle version maybe good here.

test(FMA accuracy tests, t => {
actual = Math.fma(10.20, 20.91, 30.12);
expected = 243.402;
t.ok(actual === expected, ${actual - expected}, some value);

actual = Math.fma(0.1, 10, -1);
expected = 5.5511151231257827e-17;
t.ok(actual === expected, ${actual - expected}, non zero);

actual = Math.fma(1 + eps, 1 + eps, -1 - 2 * eps);
expected = eps * eps;
t.ok(actual === expected, ${actual - expected}, Nelson H. F. Beebe's test);

actual = Math.fma(1 + Math.pow(2, -30), 1 + Math.pow(2, -23), -(1 + Math.pow(2, -23) + Math.pow(2, -30)));
expected = Math.pow(2, -53);
t.ok(actual === expected, ${actual - expected}, Accurate Algorithms test1);

actual = Math.fma(1 + Math.pow(2, -30), 1 + Math.pow(2, -52), -(1 + Math.pow(2, -30)));
expected = Math.pow(2, -52) + Math.pow(2, -82);
t.ok(actual === expected, ${actual - expected}, Accurate Algorithms test2);
});

test('FMA infinity tests', t => {
t.ok(Math.fma(Infinity, 1, 1) === Infinity, 'is Infinity');
t.ok(Math.fma(1, Infinity, 1) === Infinity, 'is Infinity');
t.ok(Math.fma(1, 1, Infinity) === Infinity, 'is Infinity');

t.ok(Math.fma(Infinity, 1, 0) === Infinity, 'is Infinity');
t.ok(Math.fma(1, Infinity, 0) === Infinity, 'is Infinity');
t.ok(Math.fma(0, 1, Infinity) === Infinity, 'is Infinity');

t.ok(Math.fma(1, Infinity, Infinity) === Infinity, 'is Infinity');
t.ok(Math.fma(Infinity, 1, Infinity) === Infinity, 'is Infinity');
t.ok(Math.fma(Infinity, Infinity, 1) === Infinity, 'is Infinity');

t.ok(Math.fma(-Infinity, 1, 1) === -Infinity, 'is Infinity');
t.ok(Math.fma(1, -Infinity, 1) === -Infinity, 'is Infinity');
t.ok(Math.fma(1, 1, -Infinity) === -Infinity, 'is Infinity');

t.ok(Math.fma(-Infinity, 1, 0) === -Infinity, 'is Infinity');
t.ok(Math.fma(1, -Infinity, 0) === -Infinity, 'is Infinity');
t.ok(Math.fma(0, 1, -Infinity) === -Infinity, 'is Infinity');

t.ok(Math.fma(1, -Infinity, -Infinity) === -Infinity, 'is Infinity');
t.ok(Math.fma(-Infinity, 1, -Infinity) === -Infinity, 'is Infinity');
t.ok(Math.fma(-Infinity, -Infinity, 1) === Infinity, 'is Infinity');
});

test('FMA NaN tests', t => {
t.ok(isNaN(Math.fma(NaN, 1, 1)), 'is NaN');
t.ok(isNaN(Math.fma(1, NaN, 1)), 'is NaN');
t.ok(isNaN(Math.fma(1, 1, NaN)), 'is NaN');

t.ok(isNaN(Math.fma(Infinity, 0, 1)), 'is NaN');
t.ok(isNaN(Math.fma(0, Infinity, 1)), 'is NaN');
t.ok(isNaN(Math.fma(Infinity, 0, NaN)), 'is NaN');
t.ok(isNaN(Math.fma(0, Infinity, NaN)), 'is NaN');

t.ok(isNaN(Math.fma(-Infinity, 0, 1)), 'is NaN');
t.ok(isNaN(Math.fma(0, -Infinity, 1)), 'is NaN');
t.ok(isNaN(Math.fma(-Infinity, 0, NaN)), 'is NaN');
t.ok(isNaN(Math.fma(0, -Infinity, NaN)), 'is NaN');

t.ok(isNaN(Math.fma(Infinity, 0, -Infinity)), 'is NaN');
t.ok(isNaN(Math.fma(0, -Infinity, Infinity)), 'is NaN');
t.ok(isNaN(Math.fma(Infinity, 1, -Infinity)), 'is NaN');
t.ok(isNaN(Math.fma(1, -Infinity, Infinity)), 'is NaN');

t.ok(isNaN(Math.fma(-Infinity, 0, Infinity)), 'is NaN');
t.ok(isNaN(Math.fma(0, Infinity, -Infinity)), 'is NaN');
t.ok(isNaN(Math.fma(-Infinity, 1, Infinity)), 'is NaN');
t.ok(isNaN(Math.fma(1, Infinity, -Infinity)), 'is NaN');

t.ok(isNaN(Math.fma(Infinity, Infinity, -Infinity)), 'is NaN');
t.ok(isNaN(Math.fma(-Infinity, Infinity, Infinity)), 'is NaN');
t.ok(isNaN(Math.fma(Infinity, -Infinity, Infinity)), 'is NaN');
t.ok(isNaN(Math.fma(-Infinity, -Infinity, -Infinity)), 'is NaN');
});

test('FMA overflow tests', t => {
actual = Math.fma(Number.MAX_VALUE, 2., -Number.MAX_VALUE);
expected = Number.MAX_VALUE;
t.ok(actual === expected, ${actual - expected}, overflow test);
});