Math.fma.js

// This simple snippet is a good start for FMA emulation in javascript.
// It pass a LOT of accuracy tests but can't handle cases with overflows and odd rounding
// IEEE compliant version is developed here: https://github.com/munsocket/proposal-fma/

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;
}

Click here to see passed tests:

test(`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('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('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');
});