Skip to content

Instantly share code, notes, and snippets.

@Yaffle
Last active Mar 19, 2017
Embed
What would you like to do?
Math.fma - fused multiply–add (FMA) or fused multiply–accumulate (FMAC) in JavaScript
// 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;
};
}());
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment