From 2fe7b06f3df06b7a442ffc3d9b951c5d52a11a24 Mon Sep 17 00:00:00 2001 From: Andrew Kelley Date: Tue, 5 Oct 2021 12:32:26 -0700 Subject: [PATCH] add support for f128 `@mulAdd` std: add f128 implementations of fma, frexp, and ilogb. Expose `fmal` in zig's freestanding libc. This makes `@mulAdd` work correctly for f128. Fixes a CI regression from yesterday, where I added a usage of f128 `@mulAdd` into the self-hosted compiler. --- lib/std/math.zig | 3 +- lib/std/math/fma.zig | 173 +++++++++++++++++++++++++++++++++-- lib/std/math/frexp.zig | 109 ++++++++++++++++++---- lib/std/math/ilogb.zig | 63 +++++++++++-- lib/std/special/c_stage1.zig | 4 + test/behavior/muladd.zig | 13 ++- 6 files changed, 326 insertions(+), 39 deletions(-) diff --git a/lib/std/math.zig b/lib/std/math.zig index d913265bde..0f3a22d2b1 100644 --- a/lib/std/math.zig +++ b/lib/std/math.zig @@ -229,8 +229,7 @@ pub const floor = @import("math/floor.zig").floor; pub const trunc = @import("math/trunc.zig").trunc; pub const round = @import("math/round.zig").round; pub const frexp = @import("math/frexp.zig").frexp; -pub const frexp32_result = @import("math/frexp.zig").frexp32_result; -pub const frexp64_result = @import("math/frexp.zig").frexp64_result; +pub const Frexp = @import("math/frexp.zig").Frexp; pub const modf = @import("math/modf.zig").modf; pub const modf32_result = @import("math/modf.zig").modf32_result; pub const modf64_result = @import("math/modf.zig").modf64_result; diff --git a/lib/std/math/fma.zig b/lib/std/math/fma.zig index ea6e15538f..54d4d071a2 100644 --- a/lib/std/math/fma.zig +++ b/lib/std/math/fma.zig @@ -1,6 +1,7 @@ -// Ported from musl, which is licensed under the MIT license: +// Ported from musl, which is MIT licensed: // https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT // +// https://git.musl-libc.org/cgit/musl/tree/src/math/fmal.c // https://git.musl-libc.org/cgit/musl/tree/src/math/fmaf.c // https://git.musl-libc.org/cgit/musl/tree/src/math/fma.c @@ -13,6 +14,7 @@ pub fn fma(comptime T: type, x: T, y: T, z: T) T { return switch (T) { f32 => fma32(x, y, z), f64 => fma64(x, y, z), + f128 => fma128(x, y, z), else => @compileError("fma not implemented for " ++ @typeName(T)), }; } @@ -142,12 +144,159 @@ fn add_and_denorm(a: f64, b: f64, scale: i32) f64 { return math.scalbn(sum.hi, scale); } -test "math.fma" { - try expect(fma(f32, 0.0, 1.0, 1.0) == fma32(0.0, 1.0, 1.0)); - try expect(fma(f64, 0.0, 1.0, 1.0) == fma64(0.0, 1.0, 1.0)); +/// A struct that represents a floating-point number with twice the precision +/// of f128. We maintain the invariant that "hi" stores the high-order +/// bits of the result. +const dd128 = struct { + hi: f128, + lo: f128, +}; + +/// Compute a+b exactly, returning the exact result in a struct dd. We assume +/// that both a and b are finite, but make no assumptions about their relative +/// magnitudes. +fn dd_add128(a: f128, b: f128) dd128 { + var ret: dd128 = undefined; + ret.hi = a + b; + const s = ret.hi - a; + ret.lo = (a - (ret.hi - s)) + (b - s); + return ret; } -test "math.fma32" { +/// Compute a+b, with a small tweak: The least significant bit of the +/// result is adjusted into a sticky bit summarizing all the bits that +/// were lost to rounding. This adjustment negates the effects of double +/// rounding when the result is added to another number with a higher +/// exponent. For an explanation of round and sticky bits, see any reference +/// on FPU design, e.g., +/// +/// J. Coonen. An Implementation Guide to a Proposed Standard for +/// Floating-Point Arithmetic. Computer, vol. 13, no. 1, Jan 1980. +fn add_adjusted128(a: f128, b: f128) f128 { + var sum = dd_add128(a, b); + if (sum.lo != 0) { + var uhii = @bitCast(u128, sum.hi); + if (uhii & 1 == 0) { + // hibits += copysign(1.0, sum.hi, sum.lo) + const uloi = @bitCast(u128, sum.lo); + uhii += 1 - ((uhii ^ uloi) >> 126); + sum.hi = @bitCast(f128, uhii); + } + } + return sum.hi; +} + +/// Compute ldexp(a+b, scale) with a single rounding error. It is assumed +/// that the result will be subnormal, and care is taken to ensure that +/// double rounding does not occur. +fn add_and_denorm128(a: f128, b: f128, scale: i32) f128 { + var sum = dd_add128(a, b); + // If we are losing at least two bits of accuracy to denormalization, + // then the first lost bit becomes a round bit, and we adjust the + // lowest bit of sum.hi to make it a sticky bit summarizing all the + // bits in sum.lo. With the sticky bit adjusted, the hardware will + // break any ties in the correct direction. + // + // If we are losing only one bit to denormalization, however, we must + // break the ties manually. + if (sum.lo != 0) { + var uhii = @bitCast(u128, sum.hi); + const bits_lost = -@intCast(i32, (uhii >> 112) & 0x7FFF) - scale + 1; + if ((bits_lost != 1) == (uhii & 1 != 0)) { + const uloi = @bitCast(u128, sum.lo); + uhii += 1 - (((uhii ^ uloi) >> 126) & 2); + sum.hi = @bitCast(f128, uhii); + } + } + return math.scalbn(sum.hi, scale); +} + +/// Compute a*b exactly, returning the exact result in a struct dd. We assume +/// that both a and b are normalized, so no underflow or overflow will occur. +/// The current rounding mode must be round-to-nearest. +fn dd_mul128(a: f128, b: f128) dd128 { + var ret: dd128 = undefined; + const split: f128 = 0x1.0p57 + 1.0; + + var p = a * split; + var ha = a - p; + ha += p; + var la = a - ha; + + p = b * split; + var hb = b - p; + hb += p; + var lb = b - hb; + + p = ha * hb; + var q = ha * lb + la * hb; + + ret.hi = p + q; + ret.lo = p - ret.hi + q + la * lb; + return ret; +} + +/// 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). +fn fma128(x: f128, y: f128, z: f128) 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); + var ex = x1.exponent; + var xs = x1.significand; + const x2 = math.frexp(y); + var ey = x2.exponent; + var ys = x2.significand; + const x3 = math.frexp(z); + var 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(f128, math.f128_min, 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); + } +} + +test "type dispatch" { + try expect(fma(f32, 0.0, 1.0, 1.0) == fma32(0.0, 1.0, 1.0)); + try expect(fma(f64, 0.0, 1.0, 1.0) == fma64(0.0, 1.0, 1.0)); + try expect(fma(f128, 0.0, 1.0, 1.0) == fma128(0.0, 1.0, 1.0)); +} + +test "32" { const epsilon = 0.000001; try expect(math.approxEqAbs(f32, fma32(0.0, 5.0, 9.124), 9.124, epsilon)); @@ -159,7 +308,7 @@ test "math.fma32" { try expect(math.approxEqAbs(f32, fma32(123123.234375, 5.0, 9.124), 615625.295875, epsilon)); } -test "math.fma64" { +test "64" { const epsilon = 0.000001; try expect(math.approxEqAbs(f64, fma64(0.0, 5.0, 9.124), 9.124, epsilon)); @@ -170,3 +319,15 @@ test "math.fma64" { try expect(math.approxEqAbs(f64, fma64(89.123, 5.0, 9.124), 454.739, epsilon)); try expect(math.approxEqAbs(f64, fma64(123123.234375, 5.0, 9.124), 615625.295875, epsilon)); } + +test "128" { + const epsilon = 0.000001; + + try expect(math.approxEqAbs(f128, fma128(0.0, 5.0, 9.124), 9.124, epsilon)); + try expect(math.approxEqAbs(f128, fma128(0.2, 5.0, 9.124), 10.124, epsilon)); + try expect(math.approxEqAbs(f128, fma128(0.8923, 5.0, 9.124), 13.5855, epsilon)); + try expect(math.approxEqAbs(f128, fma128(1.5, 5.0, 9.124), 16.624, epsilon)); + try expect(math.approxEqAbs(f128, fma128(37.45, 5.0, 9.124), 196.374, epsilon)); + try expect(math.approxEqAbs(f128, fma128(89.123, 5.0, 9.124), 454.739, epsilon)); + try expect(math.approxEqAbs(f128, fma128(123123.234375, 5.0, 9.124), 615625.295875, epsilon)); +} diff --git a/lib/std/math/frexp.zig b/lib/std/math/frexp.zig index f0d0c57af5..31168d28d4 100644 --- a/lib/std/math/frexp.zig +++ b/lib/std/math/frexp.zig @@ -1,6 +1,7 @@ -// Ported from musl, which is licensed under the MIT license: +// Ported from musl, which is MIT licensed: // https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT // +// https://git.musl-libc.org/cgit/musl/tree/src/math/frexpl.c // https://git.musl-libc.org/cgit/musl/tree/src/math/frexpf.c // https://git.musl-libc.org/cgit/musl/tree/src/math/frexp.c @@ -8,14 +9,12 @@ const std = @import("../std.zig"); const math = std.math; const expect = std.testing.expect; -fn frexp_result(comptime T: type) type { +pub fn Frexp(comptime T: type) type { return struct { significand: T, exponent: i32, }; } -pub const frexp32_result = frexp_result(f32); -pub const frexp64_result = frexp_result(f64); /// Breaks x into a normalized fraction and an integral power of two. /// f == frac * 2^exp, with |frac| in the interval [0.5, 1). @@ -24,17 +23,20 @@ pub const frexp64_result = frexp_result(f64); /// - frexp(+-0) = +-0, 0 /// - frexp(+-inf) = +-inf, 0 /// - frexp(nan) = nan, undefined -pub fn frexp(x: anytype) frexp_result(@TypeOf(x)) { +pub fn frexp(x: anytype) Frexp(@TypeOf(x)) { const T = @TypeOf(x); return switch (T) { f32 => frexp32(x), f64 => frexp64(x), + f128 => frexp128(x), else => @compileError("frexp not implemented for " ++ @typeName(T)), }; } -fn frexp32(x: f32) frexp32_result { - var result: frexp32_result = undefined; +// TODO: unify all these implementations using generics + +fn frexp32(x: f32) Frexp(f32) { + var result: Frexp(f32) = undefined; var y = @bitCast(u32, x); const e = @intCast(i32, y >> 23) & 0xFF; @@ -70,8 +72,8 @@ fn frexp32(x: f32) frexp32_result { return result; } -fn frexp64(x: f64) frexp64_result { - var result: frexp64_result = undefined; +fn frexp64(x: f64) Frexp(f64) { + var result: Frexp(f64) = undefined; var y = @bitCast(u64, x); const e = @intCast(i32, y >> 52) & 0x7FF; @@ -107,7 +109,44 @@ fn frexp64(x: f64) frexp64_result { return result; } -test "math.frexp" { +fn frexp128(x: f128) Frexp(f128) { + var result: Frexp(f128) = undefined; + + var y = @bitCast(u128, x); + const e = @intCast(i32, y >> 112) & 0x7FFF; + + if (e == 0) { + if (x != 0) { + // subnormal + result = frexp128(x * 0x1.0p120); + result.exponent -= 120; + } else { + // frexp(+-0) = (+-0, 0) + result.significand = x; + result.exponent = 0; + } + return result; + } else if (e == 0x7FFF) { + // frexp(nan) = (nan, undefined) + result.significand = x; + result.exponent = undefined; + + // frexp(+-inf) = (+-inf, 0) + if (math.isInf(x)) { + result.exponent = 0; + } + + return result; + } + + result.exponent = e - 0x3FFE; + y &= 0x8000FFFFFFFFFFFFFFFFFFFFFFFFFFFF; + y |= 0x3FFE0000000000000000000000000000; + result.significand = @bitCast(f128, y); + return result; +} + +test "type dispatch" { const a = frexp(@as(f32, 1.3)); const b = frexp32(1.3); try expect(a.significand == b.significand and a.exponent == b.exponent); @@ -115,11 +154,15 @@ test "math.frexp" { const c = frexp(@as(f64, 1.3)); const d = frexp64(1.3); try expect(c.significand == d.significand and c.exponent == d.exponent); + + const e = frexp(@as(f128, 1.3)); + const f = frexp128(1.3); + try expect(e.significand == f.significand and e.exponent == f.exponent); } -test "math.frexp32" { +test "32" { const epsilon = 0.000001; - var r: frexp32_result = undefined; + var r: Frexp(f32) = undefined; r = frexp32(1.3); try expect(math.approxEqAbs(f32, r.significand, 0.65, epsilon) and r.exponent == 1); @@ -128,9 +171,9 @@ test "math.frexp32" { try expect(math.approxEqAbs(f32, r.significand, 0.609558, epsilon) and r.exponent == 7); } -test "math.frexp64" { +test "64" { const epsilon = 0.000001; - var r: frexp64_result = undefined; + var r: Frexp(f64) = undefined; r = frexp64(1.3); try expect(math.approxEqAbs(f64, r.significand, 0.65, epsilon) and r.exponent == 1); @@ -139,8 +182,19 @@ test "math.frexp64" { try expect(math.approxEqAbs(f64, r.significand, 0.609558, epsilon) and r.exponent == 7); } -test "math.frexp32.special" { - var r: frexp32_result = undefined; +test "128" { + const epsilon = 0.000001; + var r: Frexp(f128) = undefined; + + r = frexp128(1.3); + try expect(math.approxEqAbs(f128, r.significand, 0.65, epsilon) and r.exponent == 1); + + r = frexp128(78.0234); + try expect(math.approxEqAbs(f128, r.significand, 0.609558, epsilon) and r.exponent == 7); +} + +test "32 special" { + var r: Frexp(f32) = undefined; r = frexp32(0.0); try expect(r.significand == 0.0 and r.exponent == 0); @@ -158,8 +212,8 @@ test "math.frexp32.special" { try expect(math.isNan(r.significand)); } -test "math.frexp64.special" { - var r: frexp64_result = undefined; +test "64 special" { + var r: Frexp(f64) = undefined; r = frexp64(0.0); try expect(r.significand == 0.0 and r.exponent == 0); @@ -176,3 +230,22 @@ test "math.frexp64.special" { r = frexp64(math.nan(f64)); try expect(math.isNan(r.significand)); } + +test "128 special" { + var r: Frexp(f128) = undefined; + + r = frexp128(0.0); + try expect(r.significand == 0.0 and r.exponent == 0); + + r = frexp128(-0.0); + try expect(r.significand == -0.0 and r.exponent == 0); + + r = frexp128(math.inf(f128)); + try expect(math.isPositiveInf(r.significand) and r.exponent == 0); + + r = frexp128(-math.inf(f128)); + try expect(math.isNegativeInf(r.significand) and r.exponent == 0); + + r = frexp128(math.nan(f128)); + try expect(math.isNan(r.significand)); +} diff --git a/lib/std/math/ilogb.zig b/lib/std/math/ilogb.zig index e8fc56d1f8..5056ba30c2 100644 --- a/lib/std/math/ilogb.zig +++ b/lib/std/math/ilogb.zig @@ -1,6 +1,7 @@ -// Ported from musl, which is licensed under the MIT license: +// Ported from musl, which is MIT licensed. // https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT // +// https://git.musl-libc.org/cgit/musl/tree/src/math/ilogbl.c // https://git.musl-libc.org/cgit/musl/tree/src/math/ilogbf.c // https://git.musl-libc.org/cgit/musl/tree/src/math/ilogb.c @@ -21,10 +22,13 @@ pub fn ilogb(x: anytype) i32 { return switch (T) { f32 => ilogb32(x), f64 => ilogb64(x), + f128 => ilogb128(x), else => @compileError("ilogb not implemented for " ++ @typeName(T)), }; } +// TODO: unify these implementations with generics + // NOTE: Should these be exposed publicly? const fp_ilogbnan = -1 - @as(i32, maxInt(u32) >> 1); const fp_ilogb0 = fp_ilogbnan; @@ -100,12 +104,43 @@ fn ilogb64(x: f64) i32 { return e - 0x3FF; } -test "math.ilogb" { +fn ilogb128(x: f128) i32 { + var u = @bitCast(u128, x); + var e = @intCast(i32, (u >> 112) & 0x7FFF); + + if (math.isNan(x)) { + return maxInt(i32); + } + + if (e == 0) { + u <<= 16; + if (u == 0) { + math.raiseInvalid(); + return fp_ilogb0; + } + + // subnormal x + return ilogb128(x * 0x1p120) - 120; + } + + if (e == 0x7FFF) { + math.raiseInvalid(); + if (u << 16 != 0) { + return fp_ilogbnan; + } else { + return maxInt(i32); + } + } + + return e - 0x3FFF; +} + +test "type dispatch" { try expect(ilogb(@as(f32, 0.2)) == ilogb32(0.2)); try expect(ilogb(@as(f64, 0.2)) == ilogb64(0.2)); } -test "math.ilogb32" { +test "32" { try expect(ilogb32(0.0) == fp_ilogb0); try expect(ilogb32(0.5) == -1); try expect(ilogb32(0.8923) == -1); @@ -114,7 +149,7 @@ test "math.ilogb32" { try expect(ilogb32(2398.23) == 11); } -test "math.ilogb64" { +test "64" { try expect(ilogb64(0.0) == fp_ilogb0); try expect(ilogb64(0.5) == -1); try expect(ilogb64(0.8923) == -1); @@ -123,16 +158,32 @@ test "math.ilogb64" { try expect(ilogb64(2398.23) == 11); } -test "math.ilogb32.special" { +test "128" { + try expect(ilogb128(0.0) == fp_ilogb0); + try expect(ilogb128(0.5) == -1); + try expect(ilogb128(0.8923) == -1); + try expect(ilogb128(10.0) == 3); + try expect(ilogb128(-123984) == 16); + try expect(ilogb128(2398.23) == 11); +} + +test "32 special" { try expect(ilogb32(math.inf(f32)) == maxInt(i32)); try expect(ilogb32(-math.inf(f32)) == maxInt(i32)); try expect(ilogb32(0.0) == minInt(i32)); try expect(ilogb32(math.nan(f32)) == maxInt(i32)); } -test "math.ilogb64.special" { +test "64 special" { try expect(ilogb64(math.inf(f64)) == maxInt(i32)); try expect(ilogb64(-math.inf(f64)) == maxInt(i32)); try expect(ilogb64(0.0) == minInt(i32)); try expect(ilogb64(math.nan(f64)) == maxInt(i32)); } + +test "128 special" { + try expect(ilogb128(math.inf(f128)) == maxInt(i32)); + try expect(ilogb128(-math.inf(f128)) == maxInt(i32)); + try expect(ilogb128(0.0) == minInt(i32)); + try expect(ilogb128(math.nan(f128)) == maxInt(i32)); +} diff --git a/lib/std/special/c_stage1.zig b/lib/std/special/c_stage1.zig index 7ea2b95704..57eaab48b9 100644 --- a/lib/std/special/c_stage1.zig +++ b/lib/std/special/c_stage1.zig @@ -656,6 +656,10 @@ export fn ceil(x: f64) f64 { return math.ceil(x); } +export fn fmal(a: f128, b: f128, c: f128) f128 { + return math.fma(f128, a, b, c); +} + export fn fma(a: f64, b: f64, c: f64) f64 { return math.fma(f64, a, b, c); } diff --git a/test/behavior/muladd.zig b/test/behavior/muladd.zig index 87fbf57fd3..eaa30324df 100644 --- a/test/behavior/muladd.zig +++ b/test/behavior/muladd.zig @@ -24,11 +24,10 @@ fn testMulAdd() !void { var c: f64 = 6.25; try expect(@mulAdd(f64, a, b, c) == 20); } - // Awaits implementation in libm.zig - //{ - // var a: f16 = 5.5; - // var b: f128 = 2.5; - // var c: f128 = 6.25; - //try expect(@mulAdd(f128, a, b, c) == 20); - //} + { + var a: f16 = 5.5; + var b: f128 = 2.5; + var c: f128 = 6.25; + try expect(@mulAdd(f128, a, b, c) == 20); + } }