From 4e95662a4a75817d7c1a35d7eba11b708ea43a30 Mon Sep 17 00:00:00 2001 From: LemonBoy Date: Tue, 24 Mar 2020 11:09:25 +0100 Subject: [PATCH 1/2] std: Add qNaN constants --- lib/std/math.zig | 25 ++++++++++++++++++------- 1 file changed, 18 insertions(+), 7 deletions(-) diff --git a/lib/std/math.zig b/lib/std/math.zig index ed13e3cd27..df1a522046 100644 --- a/lib/std/math.zig +++ b/lib/std/math.zig @@ -61,24 +61,36 @@ pub const f16_toint = 1.0 / f16_epsilon; pub const nan_u16 = @as(u16, 0x7C01); pub const nan_f16 = @bitCast(f16, nan_u16); +pub const qnan_u16 = @as(u16, 0x7E00); +pub const qnan_f16 = @bitCast(f16, qnan_u16); + pub const inf_u16 = @as(u16, 0x7C00); pub const inf_f16 = @bitCast(f16, inf_u16); pub const nan_u32 = @as(u32, 0x7F800001); pub const nan_f32 = @bitCast(f32, nan_u32); +pub const qnan_u32 = @as(u32, 0x7FC00000); +pub const qnan_f32 = @bitCast(f32, qnan_u32); + pub const inf_u32 = @as(u32, 0x7F800000); pub const inf_f32 = @bitCast(f32, inf_u32); pub const nan_u64 = @as(u64, 0x7FF << 52) | 1; pub const nan_f64 = @bitCast(f64, nan_u64); +pub const qnan_u64 = @as(u64, 0x7ff8000000000000); +pub const qnan_f64 = @bitCast(f64, qnan_u64); + pub const inf_u64 = @as(u64, 0x7FF << 52); pub const inf_f64 = @bitCast(f64, inf_u64); pub const nan_u128 = @as(u128, 0x7fff0000000000000000000000000001); pub const nan_f128 = @bitCast(f128, nan_u128); +pub const qnan_u128 = @as(u128, 0x7fff8000000000000000000000000000); +pub const qnan_f128 = @bitCast(f128, qnan_u128); + pub const inf_u128 = @as(u128, 0x7fff0000000000000000000000000000); pub const inf_f128 = @bitCast(f128, inf_u128); @@ -670,13 +682,12 @@ fn testRem() void { /// Returns the absolute value of the integer parameter. /// Result is an unsigned integer. -pub fn absCast(x: var) switch(@typeInfo(@TypeOf(x))) { - .ComptimeInt => comptime_int, - .Int => |intInfo| std.meta.IntType(false, intInfo.bits), - else => @compileError("absCast only accepts integers"), - } -{ - switch(@typeInfo(@TypeOf(x))) { +pub fn absCast(x: var) switch (@typeInfo(@TypeOf(x))) { + .ComptimeInt => comptime_int, + .Int => |intInfo| std.meta.IntType(false, intInfo.bits), + else => @compileError("absCast only accepts integers"), +} { + switch (@typeInfo(@TypeOf(x))) { .ComptimeInt => { if (x < 0) { return -x; From cc774c603bbba1cee9499a556a71a732a4a887f6 Mon Sep 17 00:00:00 2001 From: LemonBoy Date: Tue, 24 Mar 2020 11:09:47 +0100 Subject: [PATCH 2/2] compiler-rt: Add __divtf3 --- lib/std/special/compiler_rt.zig | 1 + lib/std/special/compiler_rt/divdf3.zig | 4 +- lib/std/special/compiler_rt/divtf3.zig | 228 ++++++++++++++++++++ lib/std/special/compiler_rt/divtf3_test.zig | 46 ++++ 4 files changed, 277 insertions(+), 2 deletions(-) create mode 100644 lib/std/special/compiler_rt/divtf3.zig create mode 100644 lib/std/special/compiler_rt/divtf3_test.zig diff --git a/lib/std/special/compiler_rt.zig b/lib/std/special/compiler_rt.zig index 4a9468fe5f..f0cf85a7a9 100644 --- a/lib/std/special/compiler_rt.zig +++ b/lib/std/special/compiler_rt.zig @@ -67,6 +67,7 @@ comptime { @export(@import("compiler_rt/divsf3.zig").__divsf3, .{ .name = "__divsf3", .linkage = linkage }); @export(@import("compiler_rt/divdf3.zig").__divdf3, .{ .name = "__divdf3", .linkage = linkage }); + @export(@import("compiler_rt/divtf3.zig").__divtf3, .{ .name = "__divtf3", .linkage = linkage }); @export(@import("compiler_rt/ashlti3.zig").__ashlti3, .{ .name = "__ashlti3", .linkage = linkage }); @export(@import("compiler_rt/lshrti3.zig").__lshrti3, .{ .name = "__lshrti3", .linkage = linkage }); diff --git a/lib/std/special/compiler_rt/divdf3.zig b/lib/std/special/compiler_rt/divdf3.zig index 7e870feec9..ed571e6024 100644 --- a/lib/std/special/compiler_rt/divdf3.zig +++ b/lib/std/special/compiler_rt/divdf3.zig @@ -203,7 +203,7 @@ pub fn __divdf3(a: f64, b: f64) callconv(.C) f64 { } } -fn wideMultiply(comptime Z: type, a: Z, b: Z, hi: *Z, lo: *Z) void { +pub fn wideMultiply(comptime Z: type, a: Z, b: Z, hi: *Z, lo: *Z) void { @setRuntimeSafety(builtin.is_test); switch (Z) { u32 => { @@ -312,7 +312,7 @@ fn wideMultiply(comptime Z: type, a: Z, b: Z, hi: *Z, lo: *Z) void { } } -fn normalize(comptime T: type, significand: *std.meta.IntType(false, T.bit_count)) i32 { +pub fn normalize(comptime T: type, significand: *std.meta.IntType(false, T.bit_count)) i32 { @setRuntimeSafety(builtin.is_test); const Z = std.meta.IntType(false, T.bit_count); const significandBits = std.math.floatMantissaBits(T); diff --git a/lib/std/special/compiler_rt/divtf3.zig b/lib/std/special/compiler_rt/divtf3.zig new file mode 100644 index 0000000000..b95f0dba66 --- /dev/null +++ b/lib/std/special/compiler_rt/divtf3.zig @@ -0,0 +1,228 @@ +const std = @import("std"); +const builtin = @import("builtin"); + +const normalize = @import("divdf3.zig").normalize; +const wideMultiply = @import("divdf3.zig").wideMultiply; + +pub fn __divtf3(a: f128, b: f128) callconv(.C) f128 { + @setRuntimeSafety(builtin.is_test); + const Z = std.meta.IntType(false, f128.bit_count); + const SignedZ = std.meta.IntType(true, f128.bit_count); + + const typeWidth = f128.bit_count; + const significandBits = std.math.floatMantissaBits(f128); + const exponentBits = std.math.floatExponentBits(f128); + + const signBit = (@as(Z, 1) << (significandBits + exponentBits)); + const maxExponent = ((1 << exponentBits) - 1); + const exponentBias = (maxExponent >> 1); + + const implicitBit = (@as(Z, 1) << significandBits); + const quietBit = implicitBit >> 1; + const significandMask = implicitBit - 1; + + const absMask = signBit - 1; + const exponentMask = absMask ^ significandMask; + const qnanRep = exponentMask | quietBit; + const infRep = @bitCast(Z, std.math.inf(f128)); + + const aExponent = @truncate(u32, (@bitCast(Z, a) >> significandBits) & maxExponent); + const bExponent = @truncate(u32, (@bitCast(Z, b) >> significandBits) & maxExponent); + const quotientSign: Z = (@bitCast(Z, a) ^ @bitCast(Z, b)) & signBit; + + var aSignificand: Z = @bitCast(Z, a) & significandMask; + var bSignificand: Z = @bitCast(Z, b) & significandMask; + var scale: i32 = 0; + + // Detect if a or b is zero, denormal, infinity, or NaN. + if (aExponent -% 1 >= maxExponent -% 1 or bExponent -% 1 >= maxExponent -% 1) { + const aAbs: Z = @bitCast(Z, a) & absMask; + const bAbs: Z = @bitCast(Z, b) & absMask; + + // NaN / anything = qNaN + if (aAbs > infRep) return @bitCast(f128, @bitCast(Z, a) | quietBit); + // anything / NaN = qNaN + if (bAbs > infRep) return @bitCast(f128, @bitCast(Z, b) | quietBit); + + if (aAbs == infRep) { + // infinity / infinity = NaN + if (bAbs == infRep) { + return @bitCast(f128, qnanRep); + } + // infinity / anything else = +/- infinity + else { + return @bitCast(f128, aAbs | quotientSign); + } + } + + // anything else / infinity = +/- 0 + if (bAbs == infRep) return @bitCast(f128, quotientSign); + + if (aAbs == 0) { + // zero / zero = NaN + if (bAbs == 0) { + return @bitCast(f128, qnanRep); + } + // zero / anything else = +/- zero + else { + return @bitCast(f128, quotientSign); + } + } + // anything else / zero = +/- infinity + if (bAbs == 0) return @bitCast(f128, infRep | quotientSign); + + // one or both of a or b is denormal, the other (if applicable) is a + // normal number. Renormalize one or both of a and b, and set scale to + // include the necessary exponent adjustment. + if (aAbs < implicitBit) scale +%= normalize(f128, &aSignificand); + if (bAbs < implicitBit) scale -%= normalize(f128, &bSignificand); + } + + // Set the implicit significand bit. If we fell through from the + // denormal path it was already set by normalize( ), but setting it twice + // won't hurt anything. + aSignificand |= implicitBit; + bSignificand |= implicitBit; + var quotientExponent: i32 = @bitCast(i32, aExponent -% bExponent) +% scale; + + // Align the significand of b as a Q63 fixed-point number in the range + // [1, 2.0) and get a Q64 approximate reciprocal using a small minimax + // polynomial approximation: reciprocal = 3/4 + 1/sqrt(2) - b/2. This + // is accurate to about 3.5 binary digits. + const q63b = @truncate(u64, bSignificand >> 49); + var recip64 = @as(u64, 0x7504f333F9DE6484) -% q63b; + // 0x7504f333F9DE6484 / 2^64 + 1 = 3/4 + 1/sqrt(2) + + // Now refine the reciprocal estimate using a Newton-Raphson iteration: + // + // x1 = x0 * (2 - x0 * b) + // + // This doubles the number of correct binary digits in the approximation + // with each iteration. + var correction64: u64 = undefined; + correction64 = @truncate(u64, ~(@as(u128, recip64) *% q63b >> 64) +% 1); + recip64 = @truncate(u64, @as(u128, recip64) *% correction64 >> 63); + correction64 = @truncate(u64, ~(@as(u128, recip64) *% q63b >> 64) +% 1); + recip64 = @truncate(u64, @as(u128, recip64) *% correction64 >> 63); + correction64 = @truncate(u64, ~(@as(u128, recip64) *% q63b >> 64) +% 1); + recip64 = @truncate(u64, @as(u128, recip64) *% correction64 >> 63); + correction64 = @truncate(u64, ~(@as(u128, recip64) *% q63b >> 64) +% 1); + recip64 = @truncate(u64, @as(u128, recip64) *% correction64 >> 63); + correction64 = @truncate(u64, ~(@as(u128, recip64) *% q63b >> 64) +% 1); + recip64 = @truncate(u64, @as(u128, recip64) *% correction64 >> 63); + + // The reciprocal may have overflowed to zero if the upper half of b is + // exactly 1.0. This would sabatoge the full-width final stage of the + // computation that follows, so we adjust the reciprocal down by one bit. + recip64 -%= 1; + + // We need to perform one more iteration to get us to 112 binary digits; + // The last iteration needs to happen with extra precision. + const q127blo: u64 = @truncate(u64, bSignificand << 15); + var correction: u128 = undefined; + var reciprocal: u128 = undefined; + + // NOTE: This operation is equivalent to __multi3, which is not implemented + // in some architechure + var r64q63: u128 = undefined; + var r64q127: u128 = undefined; + var r64cH: u128 = undefined; + var r64cL: u128 = undefined; + var dummy: u128 = undefined; + wideMultiply(u128, recip64, q63b, &dummy, &r64q63); + wideMultiply(u128, recip64, q127blo, &dummy, &r64q127); + + correction = -%(r64q63 + (r64q127 >> 64)); + + const cHi = @truncate(u64, correction >> 64); + const cLo = @truncate(u64, correction); + + wideMultiply(u128, recip64, cHi, &dummy, &r64cH); + wideMultiply(u128, recip64, cLo, &dummy, &r64cL); + + reciprocal = r64cH + (r64cL >> 64); + + // Adjust the final 128-bit reciprocal estimate downward to ensure that it + // is strictly smaller than the infinitely precise exact reciprocal. Because + // the computation of the Newton-Raphson step is truncating at every step, + // this adjustment is small; most of the work is already done. + reciprocal -%= 2; + + // The numerical reciprocal is accurate to within 2^-112, lies in the + // interval [0.5, 1.0), and is strictly smaller than the true reciprocal + // of b. Multiplying a by this reciprocal thus gives a numerical q = a/b + // in Q127 with the following properties: + // + // 1. q < a/b + // 2. q is in the interval [0.5, 2.0) + // 3. The error in q is bounded away from 2^-113 (actually, we have a + // couple of bits to spare, but this is all we need). + + // We need a 128 x 128 multiply high to compute q. + var quotient: u128 = undefined; + var quotientLo: u128 = undefined; + wideMultiply(u128, aSignificand << 2, reciprocal, "ient, "ientLo); + + // Two cases: quotient is in [0.5, 1.0) or quotient is in [1.0, 2.0). + // In either case, we are going to compute a residual of the form + // + // r = a - q*b + // + // We know from the construction of q that r satisfies: + // + // 0 <= r < ulp(q)*b + // + // If r is greater than 1/2 ulp(q)*b, then q rounds up. Otherwise, we + // already have the correct result. The exact halfway case cannot occur. + // We also take this time to right shift quotient if it falls in the [1,2) + // range and adjust the exponent accordingly. + var residual: u128 = undefined; + var qb: u128 = undefined; + + if (quotient < (implicitBit << 1)) { + wideMultiply(u128, quotient, bSignificand, &dummy, &qb); + residual = (aSignificand << 113) -% qb; + quotientExponent -%= 1; + } else { + quotient >>= 1; + wideMultiply(u128, quotient, bSignificand, &dummy, &qb); + residual = (aSignificand << 112) -% qb; + } + + const writtenExponent = quotientExponent +% exponentBias; + + if (writtenExponent >= maxExponent) { + // If we have overflowed the exponent, return infinity. + return @bitCast(f128, infRep | quotientSign); + } else if (writtenExponent < 1) { + if (writtenExponent == 0) { + // Check whether the rounded result is normal. + const round = @boolToInt((residual << 1) > bSignificand); + // Clear the implicit bit. + var absResult = quotient & significandMask; + // Round. + absResult += round; + if ((absResult & ~significandMask) > 0) { + // The rounded result is normal; return it. + return @bitCast(f128, absResult | quotientSign); + } + } + // Flush denormals to zero. In the future, it would be nice to add + // code to round them correctly. + return @bitCast(f128, quotientSign); + } else { + const round = @boolToInt((residual << 1) >= bSignificand); + // Clear the implicit bit + var absResult = quotient & significandMask; + // Insert the exponent + absResult |= @intCast(Z, writtenExponent) << significandBits; + // Round + absResult +%= round; + // Insert the sign and return + return @bitCast(f128, absResult | quotientSign); + } +} + +test "import divtf3" { + _ = @import("divtf3_test.zig"); +} diff --git a/lib/std/special/compiler_rt/divtf3_test.zig b/lib/std/special/compiler_rt/divtf3_test.zig new file mode 100644 index 0000000000..4940b2c403 --- /dev/null +++ b/lib/std/special/compiler_rt/divtf3_test.zig @@ -0,0 +1,46 @@ +const std = @import("std"); +const math = std.math; +const testing = std.testing; + +const __divtf3 = @import("divtf3.zig").__divtf3; + +fn compareResultLD(result: f128, expectedHi: u64, expectedLo: u64) bool { + const rep = @bitCast(u128, result); + const hi = @truncate(u64, rep >> 64); + const lo = @truncate(u64, rep); + + if (hi == expectedHi and lo == expectedLo) { + return true; + } + // test other possible NaN representation(signal NaN) + else if (expectedHi == 0x7fff800000000000 and expectedLo == 0) { + if ((hi & 0x7fff000000000000) == 0x7fff000000000000 and + ((hi & 0xffffffffffff) > 0 or lo > 0)) + { + return true; + } + } + return false; +} + +fn test__divtf3(a: f128, b: f128, expectedHi: u64, expectedLo: u64) void { + const x = __divtf3(a, b); + const ret = compareResultLD(x, expectedHi, expectedLo); + testing.expect(ret == true); +} + +test "divtf3" { + // qNaN / any = qNaN + test__divtf3(math.qnan_f128, 0x1.23456789abcdefp+5, 0x7fff800000000000, 0); + // NaN / any = NaN + test__divtf3(math.nan_f128, 0x1.23456789abcdefp+5, 0x7fff800000000000, 0); + // inf / any = inf + test__divtf3(math.inf_f128, 0x1.23456789abcdefp+5, 0x7fff000000000000, 0); + + test__divtf3(0x1.a23b45362464523375893ab4cdefp+5, 0x1.eedcbaba3a94546558237654321fp-1, 0x4004b0b72924d407, 0x0717e84356c6eba2); + test__divtf3(0x1.a2b34c56d745382f9abf2c3dfeffp-50, 0x1.ed2c3ba15935332532287654321fp-9, 0x3fd5b2af3f828c9b, 0x40e51f64cde8b1f2); + test__divtf3(0x1.2345f6aaaa786555f42432abcdefp+456, 0x1.edacbba9874f765463544dd3621fp+6400, 0x28c62e15dc464466, 0xb5a07586348557ac); + test__divtf3(0x1.2d3456f789ba6322bc665544edefp-234, 0x1.eddcdba39f3c8b7a36564354321fp-4455, 0x507b38442b539266, 0x22ce0f1d024e1252); + test__divtf3(0x1.2345f6b77b7a8953365433abcdefp+234, 0x1.edcba987d6bb3aa467754354321fp-4055, 0x50bf2e02f0798d36, 0x5e6fcb6b60044078); + test__divtf3(6.72420628622418701252535563464350521E-4932, 2.0, 0x0001000000000000, 0); +}