Fused multiply-add: Compute x * y + z with a single rounding error.
We use scaling to avoid overflow/underflow, along with the canonical precision-doubling technique adapted from:
Dekker, T. A Floating-Point Technique for Extending the Available Precision. Numer. Math. 18, 224-242 (1971).
pub fn fmaq(x: f128, y: f128, z: f128) callconv(.c) f128
pub fn fmaq(x: f128, y: f128, z: f128) callconv(.c) f128 {
if (!math.isFinite(x) or !math.isFinite(y)) {
return x * y + z;
}
if (!math.isFinite(z)) {
return z;
}
if (x == 0.0 or y == 0.0) {
return x * y + z;
}
if (z == 0.0) {
return x * y;
}
const x1 = math.frexp(x);
const ex = x1.exponent;
const xs = x1.significand;
const x2 = math.frexp(y);
const ey = x2.exponent;
const ys = x2.significand;
const x3 = math.frexp(z);
const ez = x3.exponent;
var zs = x3.significand;
var spread = ex + ey - ez;
if (spread <= 113 * 2) {
zs = math.scalbn(zs, -spread);
} else {
zs = math.copysign(math.floatMin(f128), zs);
}
const xy = dd_mul128(xs, ys);
const r = dd_add128(xy.hi, zs);
spread = ex + ey;
if (r.hi == 0.0) {
return xy.hi + zs + math.scalbn(xy.lo, spread);
}
const adj = add_adjusted128(r.lo, xy.lo);
if (spread + math.ilogb(r.hi) > -16383) {
return math.scalbn(r.hi + adj, spread);
} else {
return add_and_denorm128(r.hi, adj, spread);
}
}