Commit cc774c603b
Changed files (4)
lib
std
special
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);
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");
+}
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);
+}
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 });