diff --git a/lib/std/math/acos.zig b/lib/std/math/acos.zig index 6f734dbacb..dae3f70c82 100644 --- a/lib/std/math/acos.zig +++ b/lib/std/math/acos.zig @@ -3,10 +3,13 @@ // // https://git.musl-libc.org/cgit/musl/tree/src/math/acosf.c // https://git.musl-libc.org/cgit/musl/tree/src/math/acos.c +// https://git.musl-libc.org/cgit/musl/tree/src/math/acosl.c const std = @import("../std.zig"); const math = std.math; -const expect = std.testing.expect; +const testing = std.testing; +const builtin = @import("builtin"); +const native_endian = builtin.cpu.arch.endian(); /// Returns the arc-cosine of x. /// @@ -15,71 +18,120 @@ const expect = std.testing.expect; pub fn acos(x: anytype) @TypeOf(x) { const T = @TypeOf(x); return switch (T) { - f32 => acos32(x), - f64 => acos64(x), + f16 => acosBinary16(x), + f32 => acosBinary32(x), + f64 => acosBinary64(x), + f80 => acosExtended80(x), + f128 => acosBinary128(x), else => @compileError("acos not implemented for " ++ @typeName(T)), }; } -fn r32(z: f32) f32 { - const pS0 = 1.6666586697e-01; - const pS1 = -4.2743422091e-02; - const pS2 = -8.6563630030e-03; - const qS1 = -7.0662963390e-01; +fn approxBinary16(z: f32) f32 { + const S0: f32 = 1.0000001e0; + const S1: f32 = 1.6664918e-1; + const S2: f32 = 7.55022e-2; + const S3: f32 = 3.9513987e-2; + const S4: f32 = 5.0883885e-2; + return S0 + z * (S1 + z * (S2 + z * (S3 + z * S4))); +} + +fn acosBinary16(x: f16) f16 { + const pio2: f32 = math.pi / 2.0; + + const hx: u16 = @bitCast(x); + const ix: u16 = hx & 0x7fff; + + // |x| >= 1 or nan + if (ix >= 0x3c00) { + if (ix == 0x3c00) { + if (hx >> 15 != 0) { + return @floatCast(2.0 * pio2 + 0x1p-120); + } + return 0.0; + } + return 0.0 / (x - x); + } + + const xf: f32 = @floatCast(x); + + // |x| < 0.5 + if (ix < 0x3800) { + return @floatCast(pio2 - xf * approxBinary16(xf * xf)); + } + + // x < -0.5 + if (hx >> 15 != 0) { + const z = (1.0 + xf) * 0.5; + const s = @sqrt(z); + const w = approxBinary16(z) * s; + return @floatCast(2.0 * (pio2 - w)); + } + + // x > 0.5 + const z = (1.0 - xf) * 0.5; + const s = @sqrt(z); + const w = approxBinary16(z) * s; + return @floatCast(2.0 * w); +} + +fn rationalApproxBinary32(z: f32) f32 { + const pS0: f32 = 1.6666586697e-01; + const pS1: f32 = -4.2743422091e-02; + const pS2: f32 = -8.6563630030e-03; + const qS1: f32 = -7.0662963390e-01; const p = z * (pS0 + z * (pS1 + z * pS2)); const q = 1.0 + z * qS1; return p / q; } -fn acos32(x: f32) f32 { - const pio2_hi = 1.5707962513e+00; - const pio2_lo = 7.5497894159e-08; +fn acosBinary32(x: f32) f32 { + const pio2_hi: f32 = 1.5707962513e+00; + const pio2_lo: f32 = 7.5497894159e-08; - const hx: u32 = @as(u32, @bitCast(x)); - const ix: u32 = hx & 0x7FFFFFFF; + const hx: u32 = @bitCast(x); + const ix: u32 = hx & 0x7fff_ffff; // |x| >= 1 or nan - if (ix >= 0x3F800000) { - if (ix == 0x3F800000) { + if (ix >= 0x3f800000) { + if (ix == 0x3f800000) { if (hx >> 31 != 0) { return 2.0 * pio2_hi + 0x1.0p-120; - } else { - return 0.0; } - } else { - return (x - x) / 0; + return 0.0; } + return 0.0 / (x - x); } // |x| < 0.5 - if (ix < 0x3F000000) { - if (ix <= 0x32800000) { // |x| < 2^(-26) + if (ix < 0x3f00_0000) { + // |x| < 2^(-26) + if (ix <= 0x3280_0000) { return pio2_hi + 0x1.0p-120; - } else { - return pio2_hi - (x - (pio2_lo - x * r32(x * x))); } + return pio2_hi - (x - (pio2_lo - x * rationalApproxBinary32(x * x))); } // x < -0.5 if (hx >> 31 != 0) { const z = (1 + x) * 0.5; const s = @sqrt(z); - const w = r32(z) * s - pio2_lo; - return 2 * (pio2_hi - (s + w)); + const w = rationalApproxBinary32(z) * s - pio2_lo; + return 2.0 * (pio2_hi - (s + w)); } // x > 0.5 const z = (1.0 - x) * 0.5; const s = @sqrt(z); - const jx = @as(u32, @bitCast(s)); - const df = @as(f32, @bitCast(jx & 0xFFFFF000)); + const hs: u32 = @bitCast(s); + const df: f32 = @bitCast(hs & 0xffff_f000); const c = (z - df * df) / (s + df); - const w = r32(z) * s + c; - return 2 * (df + w); + const w = rationalApproxBinary32(z) * s + c; + return 2.0 * (df + w); } -fn r64(z: f64) f64 { +fn rationalApproxBinary64(z: f64) f64 { const pS0: f64 = 1.66666666666666657415e-01; const pS1: f64 = -3.25565818622400915405e-01; const pS2: f64 = 2.01212532134862925881e-01; @@ -96,91 +148,292 @@ fn r64(z: f64) f64 { return p / q; } -fn acos64(x: f64) f64 { +fn acosBinary64(x: f64) f64 { const pio2_hi: f64 = 1.57079632679489655800e+00; const pio2_lo: f64 = 6.12323399573676603587e-17; - const ux = @as(u64, @bitCast(x)); - const hx = @as(u32, @intCast(ux >> 32)); - const ix = hx & 0x7FFFFFFF; + const hx: u32 = @intCast(@as(u64, @bitCast(x)) >> 32); + const ix: u32 = hx & 0x7fff_ffff; // |x| >= 1 or nan - if (ix >= 0x3FF00000) { - const lx = @as(u32, @intCast(ux & 0xFFFFFFFF)); - - // acos(1) = 0, acos(-1) = pi - if ((ix - 0x3FF00000) | lx == 0) { + if (ix >= 0x3ff0_0000) { + const lx: u32 = @truncate(@as(u64, @bitCast(x))); + if ((ix - 0x3ff0_0000 | lx) == 0) { if (hx >> 31 != 0) { - return 2 * pio2_hi + 0x1.0p-120; - } else { - return 0; + return 2.0 * pio2_hi + 0x1.0p-120; } + return 0.0; } - - return (x - x) / 0; + return 0.0 / (x - x); } // |x| < 0.5 - if (ix < 0x3FE00000) { + if (ix < 0x3fe0_0000) { // |x| < 2^(-57) - if (ix <= 0x3C600000) { + if (ix <= 0x3c60_0000) { return pio2_hi + 0x1.0p-120; - } else { - return pio2_hi - (x - (pio2_lo - x * r64(x * x))); } + return pio2_hi - (x - (pio2_lo - x * rationalApproxBinary64(x * x))); } // x < -0.5 if (hx >> 31 != 0) { const z = (1.0 + x) * 0.5; const s = @sqrt(z); - const w = r64(z) * s - pio2_lo; + const w = rationalApproxBinary64(z) * s - pio2_lo; return 2 * (pio2_hi - (s + w)); } // x > 0.5 const z = (1.0 - x) * 0.5; const s = @sqrt(z); - const jx = @as(u64, @bitCast(s)); - const df = @as(f64, @bitCast(jx & 0xFFFFFFFF00000000)); + const df: f64 = @bitCast(@as(u64, @bitCast(s)) & 0xffff_ffff_0000_0000); const c = (z - df * df) / (s + df); - const w = r64(z) * s + c; - return 2 * (df + w); + const w = rationalApproxBinary64(z) * s + c; + return 2.0 * (df + w); } -test acos { - try expect(acos(@as(f32, 0.0)) == acos32(0.0)); - try expect(acos(@as(f64, 0.0)) == acos64(0.0)); +fn rationalApproxExtended80(z: f80) f80 { + const pS0: f80 = 1.66666666666666666631e-01; + const pS1: f80 = -4.16313987993683104320e-01; + const pS2: f80 = 3.69068046323246813704e-01; + const pS3: f80 = -1.36213932016738603108e-01; + const pS4: f80 = 1.78324189708471965733e-02; + const pS5: f80 = -2.19216428382605211588e-04; + const pS6: f80 = -7.10526623669075243183e-06; + const qS1: f80 = -2.94788392796209867269e+00; + const qS2: f80 = 3.27309890266528636716e+00; + const qS3: f80 = -1.68285799854822427013e+00; + const qS4: f80 = 3.90699412641738801874e-01; + const qS5: f80 = -3.14365703596053263322e-02; + + const p = z * (pS0 + z * (pS1 + z * (pS2 + z * (pS3 + z * (pS4 + z * (pS5 + z * pS6)))))); + const q = 1.0 + z * (qS1 + z * (qS2 + z * (qS3 + z * (qS4 + z * qS5)))); + return p / q; } -test acos32 { - const epsilon = 0.000001; +fn acosExtended80(x: f80) f80 { + const pio2_hi: f80 = 1.57079632679489661926; + const pio2_lo: f80 = -2.50827880633416601173e-20; - try expect(math.approxEqAbs(f32, acos32(0.0), 1.570796, epsilon)); - try expect(math.approxEqAbs(f32, acos32(0.2), 1.369438, epsilon)); - try expect(math.approxEqAbs(f32, acos32(0.3434), 1.220262, epsilon)); - try expect(math.approxEqAbs(f32, acos32(0.5), 1.047198, epsilon)); - try expect(math.approxEqAbs(f32, acos32(0.8923), 0.468382, epsilon)); - try expect(math.approxEqAbs(f32, acos32(-0.2), 1.772154, epsilon)); + const hx: u80 = @bitCast(x); + const se: u16 = @truncate(hx >> 64); + const e = se & 0x7fff; + + // |x| >= 1 or nan + if (e >= 0x3fff) { + if (x == 1.0) { + return 0.0; + } + if (x == -1.0) { + return 2.0 * pio2_hi + 0x1p-120; + } + return 0.0 / (x - x); + } + // |x| < 0.5 + if (e < 0x3fff - 1) { + if (e < 0x3fff - math.floatFractionalBits(f80)) { + return pio2_hi + 0x1p-120; + } + return pio2_hi - (rationalApproxExtended80(x * x) * x - pio2_lo + x); + } + // x < -0.5 + if (se >> 15 != 0) { + const z = (1 + x) * 0.5; + const s = @sqrt(z); + return 2.0 * (pio2_hi - (rationalApproxExtended80(z) * s - pio2_lo + s)); + } + // x > 0.5 + const z = (1.0 - x) * 0.5; + const s = @sqrt(z); + const hs: u80 = @bitCast(s); + const f: f80 = @bitCast(hs & 0xffff_ffff_ffff_0000_0000); + const c = (z - f * f) / (s + f); + return 2.0 * (rationalApproxExtended80(z) * s + c + f); } -test acos64 { - const epsilon = 0.000001; +fn rationalApproxBinary128(z: f128) f128 { + const pS0: f128 = 1.66666666666666666666666666666700314e-01; + const pS1: f128 = -7.32816946414566252574527475428622708e-01; + const pS2: f128 = 1.34215708714992334609030036562143589e+00; + const pS3: f128 = -1.32483151677116409805070261790752040e+00; + const pS4: f128 = 7.61206183613632558824485341162121989e-01; + const pS5: f128 = -2.56165783329023486777386833928147375e-01; + const pS6: f128 = 4.80718586374448793411019434585413855e-02; + const pS7: f128 = -4.42523267167024279410230886239774718e-03; + const pS8: f128 = 1.44551535183911458253205638280410064e-04; + const pS9: f128 = -2.10558957916600254061591040482706179e-07; + const qS1: f128 = -4.84690167848739751544716485245697428e+00; + const qS2: f128 = 9.96619113536172610135016921140206980e+00; + const qS3: f128 = -1.13177895428973036660836798461641458e+01; + const qS4: f128 = 7.74004374389488266169304117714658761e+00; + const qS5: f128 = -3.25871986053534084709023539900339905e+00; + const qS6: f128 = 8.27830318881232209752469022352928864e-01; + const qS7: f128 = -1.18768052702942805423330715206348004e-01; + const qS8: f128 = 8.32600764660522313269101537926539470e-03; + const qS9: f128 = -1.99407384882605586705979504567947007e-04; - try expect(math.approxEqAbs(f64, acos64(0.0), 1.570796, epsilon)); - try expect(math.approxEqAbs(f64, acos64(0.2), 1.369438, epsilon)); - try expect(math.approxEqAbs(f64, acos64(0.3434), 1.220262, epsilon)); - try expect(math.approxEqAbs(f64, acos64(0.5), 1.047198, epsilon)); - try expect(math.approxEqAbs(f64, acos64(0.8923), 0.468382, epsilon)); - try expect(math.approxEqAbs(f64, acos64(-0.2), 1.772154, epsilon)); + const p = z * (pS0 + z * (pS1 + z * (pS2 + z * (pS3 + z * (pS4 + z * (pS5 + z * (pS6 + z * (pS7 + z * (pS8 + z * pS9))))))))); + const q = 1.0 + z * (qS1 + z * (qS2 + z * (qS3 + z * (qS4 + z * (qS5 + z * (qS6 + z * (qS7 + z * (qS8 + z * qS9)))))))); + return p / q; } -test "acos32.special" { - try expect(math.isNan(acos32(-2))); - try expect(math.isNan(acos32(1.5))); +fn acosBinary128(x: f128) f128 { + const pio2_hi: f128 = 1.57079632679489661923132169163975140; + const pio2_lo: f128 = 4.33590506506189051239852201302167613e-35; + + const hx: u128 = @bitCast(x); + const se: u16 = @truncate(hx >> 112); + const e = se & 0x7fff; + + // |x| >= 1 or nan + if (e >= 0x3fff) { + if (x == 1.0) { + return 0.0; + } + if (x == -1.0) { + return 2 * pio2_hi + 0x1p-120; + } + return 0.0 / (x - x); + } + // |x| < 0.5 + if (e < 0x3fff - 1) { + if (e < 0x3fff - math.floatFractionalBits(f128)) { + return pio2_hi + 0x1p-120; + } + return pio2_hi - (rationalApproxBinary128(x * x) * x - pio2_lo + x); + } + // x < -0.5 + if (se >> 15 != 0) { + const z = (1 + x) * 0.5; + const s = @sqrt(z); + return 2 * (pio2_hi - (rationalApproxBinary128(z) * s - pio2_lo + s)); + } + // x > 0.5 + const z = (1.0 - x) * 0.5; + const s = @sqrt(z); + const hs: u128 = @bitCast(s); + const f: f128 = @bitCast(hs & 0xffff_ffff_ffff_ffff_0000_0000_0000_0000); + const c = (z - f * f) / (s + f); + return 2.0 * (rationalApproxBinary128(z) * s + c + f); } -test "acos64.special" { - try expect(math.isNan(acos64(-2))); - try expect(math.isNan(acos64(1.5))); +test "acosBinary16.special" { + try testing.expectApproxEqAbs(acosBinary16(0x0p+0), 0x1.92p0, math.floatEpsAt(f16, 0x1.92p0)); + try testing.expectApproxEqAbs(acosBinary16(-0x1p+0), 0x1.92p1, math.floatEpsAt(f16, 0x1.92p1)); + try testing.expectEqual(acosBinary16(0x1p+0), 0x0p+0); + try testing.expect(math.isNan(acosBinary16(0x1.004p0))); + try testing.expect(math.isNan(acosBinary16(-0x1.004p0))); + try testing.expect(math.isNan(acosBinary16(math.inf(f16)))); + try testing.expect(math.isNan(acosBinary16(-math.inf(f16)))); + try testing.expect(math.isNan(acosBinary16(math.nan(f16)))); +} + +test "acosBinary16" { + try testing.expectApproxEqAbs(acosBinary16(0x1.db4p-5), 0x1.834p0, math.floatEpsAt(f16, 0x1.834p0)); + try testing.expectApproxEqAbs(acosBinary16(-0x1.068p-2), 0x1.d48p0, math.floatEpsAt(f16, 0x1.d48p0)); + try testing.expectApproxEqAbs(acosBinary16(-0x1.2c4p-3), 0x1.b7cp0, math.floatEpsAt(f16, 0x1.b7cp0)); + try testing.expectApproxEqAbs(acosBinary16(0x1.65p-3), 0x1.654p0, math.floatEpsAt(f16, 0x1.654p0)); + try testing.expectApproxEqAbs(acosBinary16(0x1.dfcp-1), 0x1.6d8p-2, math.floatEpsAt(f16, 0x1.6d8p-2)); + try testing.expectApproxEqAbs(acosBinary16(-0x1.764p-1), 0x1.32p1, math.floatEpsAt(f16, 0x1.32p1)); + try testing.expectApproxEqAbs(acosBinary16(0x1.b18p-3), 0x1.5b8p0, math.floatEpsAt(f16, 0x1.5b8p0)); + try testing.expectApproxEqAbs(acosBinary16(0x1.5acp-3), 0x1.668p0, math.floatEpsAt(f16, 0x1.668p0)); + try testing.expectApproxEqAbs(acosBinary16(-0x1.18cp-1), 0x1.134p1, math.floatEpsAt(f16, 0x1.134p1)); + try testing.expectApproxEqAbs(acosBinary16(-0x1.03p-1), 0x1.0dp1, math.floatEpsAt(f16, 0x1.0dp1)); +} + +test "acosBinary32.special" { + try testing.expectApproxEqAbs(acosBinary32(0x0p+0), 0x1.921fb6p+0, math.floatEpsAt(f32, 0x1.921fb6p+0)); + try testing.expectApproxEqAbs(acosBinary32(-0x1p+0), 0x1.921fb6p+1, math.floatEpsAt(f32, 0x1.921fb6p+1)); + try testing.expectEqual(acosBinary32(0x1p+0), 0x0p+0); + try testing.expect(math.isNan(acosBinary32(0x1.000002p+0))); + try testing.expect(math.isNan(acosBinary32(-0x1.000002p+0))); + try testing.expect(math.isNan(acosBinary32(math.inf(f32)))); + try testing.expect(math.isNan(acosBinary32(-math.inf(f32)))); + try testing.expect(math.isNan(acosBinary32(math.nan(f32)))); +} + +test "acosBinary32" { + try testing.expectApproxEqAbs(acosBinary32(-0x1.13284cp-2), 0x1.d7c4e6p+0, math.floatEpsAt(f32, 0x1.d7c4e6p+0)); + try testing.expectApproxEqAbs(acosBinary32(0x1.6ca8ep-1), 0x1.8e6756p-1, math.floatEpsAt(f32, 0x1.8e6756p-1)); + try testing.expectApproxEqAbs(acosBinary32(0x1.c2ca6p-1), 0x1.f9d74cp-2, math.floatEpsAt(f32, 0x1.f9d74cp-2)); + try testing.expectApproxEqAbs(acosBinary32(-0x1.55f12p-1), 0x1.26abdcp+1, math.floatEpsAt(f32, 0x1.26abdcp+1)); + try testing.expectApproxEqAbs(acosBinary32(-0x1.15679ep-2), 0x1.d85a44p+0, math.floatEpsAt(f32, 0x1.d85a44p+0)); + try testing.expectApproxEqAbs(acosBinary32(-0x1.41e132p-5), 0x1.9c2f68p+0, math.floatEpsAt(f32, 0x1.9c2f68p+0)); + try testing.expectApproxEqAbs(acosBinary32(0x1.281b0ep-1), 0x1.e881bp-1, math.floatEpsAt(f32, 0x1.e881bp-1)); + try testing.expectApproxEqAbs(acosBinary32(0x1.b5ce34p-1), 0x1.1713f6p-1, math.floatEpsAt(f32, 0x1.1713f6p-1)); + try testing.expectApproxEqAbs(acosBinary32(-0x1.583482p-3), 0x1.bd5accp+0, math.floatEpsAt(f32, 0x1.bd5accp+0)); + try testing.expectApproxEqAbs(acosBinary32(-0x1.ea8224p-1), 0x1.6ce7d8p+1, math.floatEpsAt(f32, 0x1.6ce7d8p+1)); +} + +test "acosBinary64.special" { + try testing.expectApproxEqAbs(acosBinary64(0x0p+0), 0x1.921fb54442d18p+0, math.floatEpsAt(f64, 0x1.921fb54442d18p+0)); + try testing.expectApproxEqAbs(acosBinary64(-0x1p+0), 0x1.921fb54442d18p+1, math.floatEpsAt(f64, 0x1.921fb54442d18p+1)); + try testing.expectEqual(acosBinary64(0x1p+0), 0x0p+0); + try testing.expect(math.isNan(acosBinary64(0x1.0000000000001p+0))); + try testing.expect(math.isNan(acosBinary64(-0x1.0000000000001p+0))); + try testing.expect(math.isNan(acosBinary64(math.inf(f64)))); + try testing.expect(math.isNan(acosBinary64(-math.inf(f64)))); + try testing.expect(math.isNan(acosBinary64(math.nan(f64)))); +} + +test "acosBinary64" { + try testing.expectApproxEqAbs(acosBinary64(-0x1.13284b2b5006dp-2), 0x1.d7c4e61020905p+0, math.floatEpsAt(f64, 0x1.d7c4e61020905p+0)); + try testing.expectApproxEqAbs(acosBinary64(0x1.6ca8dfb825911p-1), 0x1.8e6756e27c366p-1, math.floatEpsAt(f64, 0x1.8e6756e27c366p-1)); + try testing.expectApproxEqAbs(acosBinary64(0x1.c2ca609de7505p-1), 0x1.f9d748eaf956p-2, math.floatEpsAt(f64, 0x1.f9d748eaf956p-2)); + try testing.expectApproxEqAbs(acosBinary64(-0x1.55f11fba96889p-1), 0x1.26abdc68d07aap+1, math.floatEpsAt(f64, 0x1.26abdc68d07aap+1)); + try testing.expectApproxEqAbs(acosBinary64(-0x1.15679e27084ddp-2), 0x1.d85a44ea44fe4p+0, math.floatEpsAt(f64, 0x1.d85a44ea44fe4p+0)); + try testing.expectApproxEqAbs(acosBinary64(-0x1.41e131b093c41p-5), 0x1.9c2f688eee8abp+0, math.floatEpsAt(f64, 0x1.9c2f688eee8abp+0)); + try testing.expectApproxEqAbs(acosBinary64(0x1.281b0d18455f5p-1), 0x1.e881b1d4eb2a1p-1, math.floatEpsAt(f64, 0x1.e881b1d4eb2a1p-1)); + try testing.expectApproxEqAbs(acosBinary64(0x1.b5ce34a51b239p-1), 0x1.1713f567a87efp-1, math.floatEpsAt(f64, 0x1.1713f567a87efp-1)); + try testing.expectApproxEqAbs(acosBinary64(-0x1.583481079de4dp-3), 0x1.bd5acbe8fcc59p+0, math.floatEpsAt(f64, 0x1.bd5acbe8fcc59p+0)); + try testing.expectApproxEqAbs(acosBinary64(-0x1.ea8223103b871p-1), 0x1.6ce7d66f628e5p+1, math.floatEpsAt(f64, 0x1.6ce7d66f628e5p+1)); +} + +test "acosExtended80.special" { + try testing.expectApproxEqAbs(acosExtended80(0x0p+0), 0x1.921fb54442d1846ap+0, math.floatEpsAt(f80, 0x1.921fb54442d1846ap+0)); + try testing.expectApproxEqAbs(acosExtended80(-0x1p+0), 0x1.921fb54442d1846ap+1, math.floatEpsAt(f80, 0x1.921fb54442d1846ap+1)); + try testing.expectEqual(acosExtended80(0x1p+0), 0x0p+0); + try testing.expect(math.isNan(acosExtended80(0x1.0000000000000002p+0))); + try testing.expect(math.isNan(acosExtended80(-0x1.0000000000000002p+0))); + try testing.expect(math.isNan(acosExtended80(math.inf(f80)))); + try testing.expect(math.isNan(acosExtended80(-math.inf(f80)))); + try testing.expect(math.isNan(acosExtended80(math.nan(f80)))); +} + +test "acosExtended80" { + try testing.expectApproxEqAbs(acosExtended80(0x1.72068a321edc8804p-1), 0x1.86b349040d28f794p-1, math.floatEpsAt(f80, 0x1.86b349040d28f794p-1)); + try testing.expectApproxEqAbs(acosExtended80(-0x1.06d0a467d22977ecp-2), 0x1.d4923ade73ec379cp0, math.floatEpsAt(f80, 0x1.d4923ade73ec379cp0)); + try testing.expectApproxEqAbs(acosExtended80(0x1.77d21385faa9798ap-3), 0x1.62e0e8898c6d04f2p0, math.floatEpsAt(f80, 0x1.62e0e8898c6d04f2p0)); + try testing.expectApproxEqAbs(acosExtended80(-0x1.73ee3e8bc2a44dbep-1), 0x1.3123cbcd5dc4bd58p1, math.floatEpsAt(f80, 0x1.3123cbcd5dc4bd58p1)); + try testing.expectApproxEqAbs(acosExtended80(0x1.0a2dd1f6ffcf668ap-1), 0x1.062a6d562df2d316p0, math.floatEpsAt(f80, 0x1.062a6d562df2d316p0)); + try testing.expectApproxEqAbs(acosExtended80(0x1.8e835c490a3aff9ep-3), 0x1.5ffd68b520aa55fap0, math.floatEpsAt(f80, 0x1.5ffd68b520aa55fap0)); + try testing.expectApproxEqAbs(acosExtended80(0x1.add20cdc1565064cp-3), 0x1.5bfe6cabda700684p0, math.floatEpsAt(f80, 0x1.5bfe6cabda700684p0)); + try testing.expectApproxEqAbs(acosExtended80(0x1.21986d43727fca72p-8), 0x1.90fe1c993b571924p0, math.floatEpsAt(f80, 0x1.90fe1c993b571924p0)); + try testing.expectApproxEqAbs(acosExtended80(0x1.d61e0b3fae6a0564p-2), 0x1.18044ccc626e7f9ep0, math.floatEpsAt(f80, 0x1.18044ccc626e7f9ep0)); + try testing.expectApproxEqAbs(acosExtended80(-0x1.171e7c4a41883ccap-4), 0x1.a39513b6c16532b4p0, math.floatEpsAt(f80, 0x1.a39513b6c16532b4p0)); +} + +test "acosBinary128.special" { + try testing.expectApproxEqAbs(acosBinary128(0x0p+0), 0x1.921fb54442d18469898cc51701b8p0, math.floatEpsAt(f128, 0x1.921fb54442d18469898cc51701b8p0)); + try testing.expectApproxEqAbs(acosBinary128(-0x1p+0), 0x1.921fb54442d18469898cc51701b8p1, math.floatEpsAt(f128, 0x1.921fb54442d18469898cc51701b8p1)); + try testing.expectEqual(acosBinary128(0x1p+0), 0x0p+0); + try testing.expect(math.isNan(acosBinary128(0x1.0000000000000000000000000001p0))); + try testing.expect(math.isNan(acosBinary128(-0x1.0000000000000000000000000001p0))); + try testing.expect(math.isNan(acosBinary128(math.inf(f128)))); + try testing.expect(math.isNan(acosBinary128(-math.inf(f128)))); + try testing.expect(math.isNan(acosBinary128(math.nan(f128)))); +} + +test "acosBinary128" { + try testing.expectApproxEqAbs(acosBinary128(-0x1.511bdb99a3c4373bedf834ef4f68p-1), 0x1.250e9a58f049eeafa99db4360c88p1, math.floatEpsAt(f128, 0x1.250e9a58f049eeafa99db4360c88p1)); + try testing.expectApproxEqAbs(acosBinary128(-0x1.5879cc3ad6dfd2a52e9891c69808p-1), 0x1.2786664b1c676c99437b68590004p1, math.floatEpsAt(f128, 0x1.2786664b1c676c99437b68590004p1)); + try testing.expectApproxEqAbs(acosBinary128(0x1.3f988ba64a7eb97a751c5f0b3077p-1), 0x1.cb190cd361c7c03a09c470b4caebp-1, math.floatEpsAt(f128, 0x1.cb190cd361c7c03a09c470b4caebp-1)); + try testing.expectApproxEqAbs(acosBinary128(-0x1.3f2d96c7768e4c4fa02315727959p-1), 0x1.1f373be697880111758f582b1a96p1, math.floatEpsAt(f128, 0x1.1f373be697880111758f582b1a96p1)); + try testing.expectApproxEqAbs(acosBinary128(0x1.fad303c2e28c1f4d8f9fd0e5686fp-2), 0x1.0d92fd2a0a6ca3e4853c1de9ea6ap0, math.floatEpsAt(f128, 0x1.0d92fd2a0a6ca3e4853c1de9ea6ap0)); + try testing.expectApproxEqAbs(acosBinary128(0x1.ddde322bd1a2ee50c5ba30c9c617p-2), 0x1.15d4b306e16fbf9ea4f29e82b154p0, math.floatEpsAt(f128, 0x1.15d4b306e16fbf9ea4f29e82b154p0)); + try testing.expectApproxEqAbs(acosBinary128(-0x1.b02f6adefcbeb1d48666b827ff17p-1), 0x1.49b0a0355a5539052388e8a6dc11p1, math.floatEpsAt(f128, 0x1.49b0a0355a5539052388e8a6dc11p1)); + try testing.expectApproxEqAbs(acosBinary128(0x1.c8581cce7cd3f6efab0fc60d9b7dp-2), 0x1.1be0b757f4cef022f5d2422b9c78p0, math.floatEpsAt(f128, 0x1.1be0b757f4cef022f5d2422b9c78p0)); + try testing.expectApproxEqAbs(acosBinary128(-0x1.bf887b8c4e33cbef59993056f3dep-1), 0x1.513270e671db2d840f20b0186c2cp1, math.floatEpsAt(f128, 0x1.513270e671db2d840f20b0186c2cp1)); + try testing.expectApproxEqAbs(acosBinary128(0x1.0c0f600ab6f9c84c6102942044cep-3), 0x1.70851a509f0e8bfbe780aa8f29f9p0, math.floatEpsAt(f128, 0x1.70851a509f0e8bfbe780aa8f29f9p0)); } diff --git a/lib/std/math/asin.zig b/lib/std/math/asin.zig index cb4571c24e..113fb584bb 100644 --- a/lib/std/math/asin.zig +++ b/lib/std/math/asin.zig @@ -3,10 +3,14 @@ // // https://git.musl-libc.org/cgit/musl/tree/src/math/asinf.c // https://git.musl-libc.org/cgit/musl/tree/src/math/asin.c +// https://git.musl-libc.org/cgit/musl/tree/src/math/asinl.c const std = @import("../std.zig"); const math = std.math; -const expect = std.testing.expect; +const mem = std.mem; +const testing = std.testing; +const builtin = @import("builtin"); +const native_endian = builtin.cpu.arch.endian(); /// Returns the arc-sin of x. /// @@ -16,62 +20,101 @@ const expect = std.testing.expect; pub fn asin(x: anytype) @TypeOf(x) { const T = @TypeOf(x); return switch (T) { - f32 => asin32(x), - f64 => asin64(x), + f16 => asinBinary16(x), + f32 => asinBinary32(x), + f64 => asinBinary64(x), + f80 => asinExtended80(x), + f128 => asinBinary128(x), else => @compileError("asin not implemented for " ++ @typeName(T)), }; } -fn r32(z: f32) f32 { - const pS0 = 1.6666586697e-01; - const pS1 = -4.2743422091e-02; - const pS2 = -8.6563630030e-03; - const qS1 = -7.0662963390e-01; +fn approxBinary16(z: f32) f32 { + const S0: f32 = 1.0000001e0; + const S1: f32 = 1.6664918e-1; + const S2: f32 = 7.55022e-2; + const S3: f32 = 3.9513987e-2; + const S4: f32 = 5.0883885e-2; + return S0 + z * (S1 + z * (S2 + z * (S3 + z * S4))); +} + +fn asinBinary16(x: f16) f16 { + const pio2: f32 = math.pi / 2.0; + + const hx: u16 = @bitCast(x); + const ix = hx & 0x7fff; + + // |x| >= 1 + if (ix >= 0x3c00) { + // |x| == 1 + if (ix == 0x3c00) { + // asin(+-1) = +-pi/2 with inexact + return @floatCast(x * pio2 + 0x1.0p-120); + } + // asin(|x| > 1) is nan + return 0.0 / (x - x); + } + + // |x| < 0.5 + if (ix < 0x3800) { + return @floatCast(x * approxBinary16(x * x)); + } + + // 1 > |x| >= 0.5 + const z = (1.0 - @abs(x)) * 0.5; + const s = @sqrt(z); + const x_local = pio2 - 2.0 * s * approxBinary16(z); + if (hx >> 15 != 0) { + return @floatCast(-x_local); + } + return @floatCast(x_local); +} + +fn rationalApproxBinary32(z: f32) f32 { + const pS0: f32 = 1.6666586697e-01; + const pS1: f32 = -4.2743422091e-02; + const pS2: f32 = -8.6563630030e-03; + const qS1: f32 = -7.0662963390e-01; const p = z * (pS0 + z * (pS1 + z * pS2)); const q = 1.0 + z * qS1; return p / q; } -fn asin32(x: f32) f32 { - const pio2 = 1.570796326794896558e+00; +fn asinBinary32(x: f32) f32 { + const pio2: f64 = 1.570796326794896558e+00; - const hx: u32 = @as(u32, @bitCast(x)); - const ix: u32 = hx & 0x7FFFFFFF; + const hx: u32 = @bitCast(x); + const ix = hx & 0x7fff_ffff; // |x| >= 1 - if (ix >= 0x3F800000) { - // |x| >= 1 - if (ix == 0x3F800000) { - return x * pio2 + 0x1.0p-120; // asin(+-1) = +-pi/2 with inexact - } else { - return math.nan(f32); // asin(|x| > 1) is nan + if (ix >= 0x3f80_0000) { + // |x| == 1 + if (ix == 0x3f80_0000) { + // asin(+-1) = +-pi/2 with inexact + return @floatCast(@as(f64, @floatCast(x)) * pio2 + 0x1.0p-120); } + // asin(|x| > 1) is nan + return 0.0 / (x - x); } // |x| < 0.5 - if (ix < 0x3F000000) { + if (ix < 0x3f00_0000) { // 0x1p-126 <= |x| < 0x1p-12 - if (ix < 0x39800000 and ix >= 0x00800000) { + if (ix < 0x3980_0000 and ix >= 0x0080_0000) { return x; - } else { - return x + x * r32(x * x); } + return x + x * rationalApproxBinary32(x * x); } // 1 > |x| >= 0.5 - const z = (1 - @abs(x)) * 0.5; - const s = @sqrt(z); - const fx = pio2 - 2 * (s + s * r32(z)); - - if (hx >> 31 != 0) { - return -fx; - } else { - return fx; - } + const z = (1.0 - @abs(x)) * 0.5; + const s: f64 = @floatCast(@sqrt(z)); + const x_local: f32 = @floatCast(pio2 - 2.0 * (s + s * @as(f64, @floatCast(rationalApproxBinary32(z))))); + return if (hx >> 31 != 0) -x_local else x_local; } -fn r64(z: f64) f64 { +fn rationalApproxBinary64(z: f64) f64 { const pS0: f64 = 1.66666666666666657415e-01; const pS1: f64 = -3.25565818622400915405e-01; const pS2: f64 = 2.01212532134862925881e-01; @@ -88,96 +131,307 @@ fn r64(z: f64) f64 { return p / q; } -fn asin64(x: f64) f64 { +fn asinBinary64(x: f64) f64 { const pio2_hi: f64 = 1.57079632679489655800e+00; const pio2_lo: f64 = 6.12323399573676603587e-17; - const ux = @as(u64, @bitCast(x)); - const hx = @as(u32, @intCast(ux >> 32)); - const ix = hx & 0x7FFFFFFF; + const hx: u32 = @intCast(@as(u64, @bitCast(x)) >> 32); + const ix = hx & 0x7fffffff; // |x| >= 1 or nan - if (ix >= 0x3FF00000) { - const lx = @as(u32, @intCast(ux & 0xFFFFFFFF)); - + if (ix >= 0x3ff0_0000) { + const lx: u32 = @truncate(@as(u64, @bitCast(x))); // asin(1) = +-pi/2 with inexact - if ((ix - 0x3FF00000) | lx == 0) { + if ((ix - 0x3ff0_0000 | lx) == 0) { return x * pio2_hi + 0x1.0p-120; - } else { - return math.nan(f64); } + return 0.0 / (x - x); } // |x| < 0.5 - if (ix < 0x3FE00000) { + if (ix < 0x3fe0_0000) { // if 0x1p-1022 <= |x| < 0x1p-26 avoid raising overflow - if (ix < 0x3E500000 and ix >= 0x00100000) { + if (ix < 0x3e50_0000 and ix >= 0x0010_0000) { return x; - } else { - return x + x * r64(x * x); } + return x + x * rationalApproxBinary64(x * x); } // 1 > |x| >= 0.5 - const z = (1 - @abs(x)) * 0.5; + const z = (1.0 - @abs(x)) * 0.5; const s = @sqrt(z); - const r = r64(z); - var fx: f64 = undefined; - + const r = rationalApproxBinary64(z); // |x| > 0.975 - if (ix >= 0x3FEF3333) { - fx = pio2_hi - 2 * (s + s * r); - } else { - const jx = @as(u64, @bitCast(s)); - const df = @as(f64, @bitCast(jx & 0xFFFFFFFF00000000)); - const c = (z - df * df) / (s + df); - fx = 0.5 * pio2_hi - (2 * s * r - (pio2_lo - 2 * c) - (0.5 * pio2_hi - 2 * df)); + if (ix >= 0x3fef_3333) { + const x_local = pio2_hi - (2 * (s + s * r) - pio2_lo); + return if (hx >> 31 != 0) -x_local else x_local; + } + // f+c = sqrt(z) + const hs: u64 = @bitCast(s); + const f: f64 = @bitCast(hs & 0xffff_ffff_0000_0000); + const c: f64 = (z - f * f) / (s + f); + const x_local = 0.5 * pio2_hi - (2.0 * s * r - (pio2_lo - 2.0 * c) - (0.5 * pio2_hi - 2.0 * f)); + return if (hx >> 31 != 0) -x_local else x_local; +} + +fn rationalApproxExtended80(z: f80) f80 { + const pS0: f80 = 1.66666666666666666631e-01; + const pS1: f80 = -4.16313987993683104320e-01; + const pS2: f80 = 3.69068046323246813704e-01; + const pS3: f80 = -1.36213932016738603108e-01; + const pS4: f80 = 1.78324189708471965733e-02; + const pS5: f80 = -2.19216428382605211588e-04; + const pS6: f80 = -7.10526623669075243183e-06; + const qS1: f80 = -2.94788392796209867269e+00; + const qS2: f80 = 3.27309890266528636716e+00; + const qS3: f80 = -1.68285799854822427013e+00; + const qS4: f80 = 3.90699412641738801874e-01; + const qS5: f80 = -3.14365703596053263322e-02; + + const p = z * (pS0 + z * (pS1 + z * (pS2 + z * (pS3 + z * (pS4 + z * (pS5 + z * pS6)))))); + const q = 1.0 + z * (qS1 + z * (qS2 + z * (qS3 + z * (qS4 + z * qS5)))); + return p / q; +} + +fn asinExtended80(x: f80) f80 { + const pio2_hi: f80 = 1.57079632679489661926; + const pio2_lo: f80 = -2.50827880633416601173e-20; + + const hx: u80 = @bitCast(x); + const se: u16 = @truncate(hx >> 64); + const e = se & 0x7fff; + const sign = se >> 15 != 0; + + // |x| >= 1 or nan + if (e >= 0x3fff) { + // asin(+-1)=+-pi/2 with inexact + if (x == 1.0 or x == -1.0) { + return x * pio2_hi + 0x1p-120; + } + return 0.0 / (x - x); } - if (hx >> 31 != 0) { - return -fx; - } else { - return fx; + // |x| < 0.5 + if (e < 0x3fff - 1) { + if (e < 0x3fff - (math.floatMantissaBits(f80) + 1) / 2) { + // return x with inexact if x!=0 + mem.doNotOptimizeAway(x + 0x1p120); + return x; + } + return x + x * rationalApproxExtended80(x * x); } + + // 1 > |x| >= 0.5 + const z = (1.0 - @abs(x)) * 0.5; + const s = @sqrt(z); + const r = rationalApproxExtended80(z); + + const m: u64 = @truncate(hx & 0x0000_ffff_ffff_ffff_ffff); + if ((m >> 56) >= 0xf7) { + const x_local = pio2_hi - (2.0 * (s + s * r) - pio2_lo); + return if (sign) -x_local else x_local; + } + + const hs: u80 = @bitCast(s); + const f: f80 = @bitCast(hs & 0xffff_ffff_ffff_0000_0000); + const c = (z - f * f) / (s + f); + const x_local = 0.5 * pio2_hi - (2.0 * s * r - (pio2_lo - 2.0 * c) - (0.5 * pio2_hi - 2.0 * f)); + return if (sign) -x_local else x_local; } -test asin { - try expect(asin(@as(f32, 0.0)) == asin32(0.0)); - try expect(asin(@as(f64, 0.0)) == asin64(0.0)); +fn rationalApproxBinary128(z: f128) f128 { + const pS0: f128 = 1.66666666666666666666666666666700314e-01; + const pS1: f128 = -7.32816946414566252574527475428622708e-01; + const pS2: f128 = 1.34215708714992334609030036562143589e+00; + const pS3: f128 = -1.32483151677116409805070261790752040e+00; + const pS4: f128 = 7.61206183613632558824485341162121989e-01; + const pS5: f128 = -2.56165783329023486777386833928147375e-01; + const pS6: f128 = 4.80718586374448793411019434585413855e-02; + const pS7: f128 = -4.42523267167024279410230886239774718e-03; + const pS8: f128 = 1.44551535183911458253205638280410064e-04; + const pS9: f128 = -2.10558957916600254061591040482706179e-07; + const qS1: f128 = -4.84690167848739751544716485245697428e+00; + const qS2: f128 = 9.96619113536172610135016921140206980e+00; + const qS3: f128 = -1.13177895428973036660836798461641458e+01; + const qS4: f128 = 7.74004374389488266169304117714658761e+00; + const qS5: f128 = -3.25871986053534084709023539900339905e+00; + const qS6: f128 = 8.27830318881232209752469022352928864e-01; + const qS7: f128 = -1.18768052702942805423330715206348004e-01; + const qS8: f128 = 8.32600764660522313269101537926539470e-03; + const qS9: f128 = -1.99407384882605586705979504567947007e-04; + + const p = z * (pS0 + z * (pS1 + z * (pS2 + z * (pS3 + z * (pS4 + z * (pS5 + z * (pS6 + z * (pS7 + z * (pS8 + z * pS9))))))))); + const q = 1.0 + z * (qS1 + z * (qS2 + z * (qS3 + z * (qS4 + z * (qS5 + z * (qS6 + z * (qS7 + z * (qS8 + z * qS9)))))))); + return p / q; } -test asin32 { - const epsilon = 0.000001; +fn asinBinary128(x: f128) f128 { + const pio2_hi: f128 = 1.57079632679489661923132169163975140; + const pio2_lo: f128 = 4.33590506506189051239852201302167613e-35; - try expect(math.approxEqAbs(f32, asin32(0.0), 0.0, epsilon)); - try expect(math.approxEqAbs(f32, asin32(0.2), 0.201358, epsilon)); - try expect(math.approxEqAbs(f32, asin32(-0.2), -0.201358, epsilon)); - try expect(math.approxEqAbs(f32, asin32(0.3434), 0.350535, epsilon)); - try expect(math.approxEqAbs(f32, asin32(0.5), 0.523599, epsilon)); - try expect(math.approxEqAbs(f32, asin32(0.8923), 1.102415, epsilon)); + const hx: u128 = @bitCast(x); + const se: u16 = @truncate(hx >> 112); + const e = se & 0x7fff; + const sign = se >> 15 != 0; + + // |x| >= 1 or nan + if (e >= 0x3fff) { + // asin(+-1)=+-pi/2 with inexact + if (x == 1.0 or x == -1.0) { + return x * pio2_hi + 0x1p-120; + } + return 0.0 / (x - x); + } + + // |x| < 0.5 + if (e < 0x3fff - 1) { + if (e < 0x3fff - (math.floatMantissaBits(f128) + 2) / 2) { + // return x with inexact if x!=0 + mem.doNotOptimizeAway(x + 0x1p120); + return x; + } + return x + x * rationalApproxBinary128(x * x); + } + + // 1 > |x| >= 0.5 + const z = (1.0 - @abs(x)) * 0.5; + const s = @sqrt(z); + const r = rationalApproxBinary128(z); + + const top: u16 = @truncate((hx >> 96) & 0x0000_ffff); + if (top >= 0xee00) { + const x_local = pio2_hi - (2.0 * (s + s * r) - pio2_lo); + return if (sign) -x_local else x_local; + } + + const hs: u128 = @bitCast(s); + const f: f128 = @bitCast(hs & 0xffff_ffff_ffff_ffff_0000_0000_0000_0000); + const c = (z - f * f) / (s + f); + const x_local = 0.5 * pio2_hi - (2.0 * s * r - (pio2_lo - 2.0 * c) - (0.5 * pio2_hi - 2.0 * f)); + return if (sign) -x_local else x_local; } -test asin64 { - const epsilon = 0.000001; - - try expect(math.approxEqAbs(f64, asin64(0.0), 0.0, epsilon)); - try expect(math.approxEqAbs(f64, asin64(0.2), 0.201358, epsilon)); - try expect(math.approxEqAbs(f64, asin64(-0.2), -0.201358, epsilon)); - try expect(math.approxEqAbs(f64, asin64(0.3434), 0.350535, epsilon)); - try expect(math.approxEqAbs(f64, asin64(0.5), 0.523599, epsilon)); - try expect(math.approxEqAbs(f64, asin64(0.8923), 1.102415, epsilon)); +test "asinBinary16.special" { + try testing.expectApproxEqAbs(asinBinary16(0x1p+0), 0x1.92p0, math.floatEpsAt(f16, 0x1.92p0)); + try testing.expectApproxEqAbs(asinBinary16(-0x1p+0), -0x1.92p0, math.floatEpsAt(f16, -0x1.92p0)); + try testing.expectEqual(asinBinary16(0x0p+0), 0x0p+0); + try testing.expectEqual(asinBinary16(-0x0p+0), 0x0p+0); + try testing.expect(math.isNan(asinBinary16(0x1.004p0))); + try testing.expect(math.isNan(asinBinary16(-0x1.004p0))); + try testing.expect(math.isNan(asinBinary16(math.inf(f16)))); + try testing.expect(math.isNan(asinBinary16(-math.inf(f16)))); + try testing.expect(math.isNan(asinBinary16(math.nan(f16)))); } -test "asin32.special" { - try expect(math.isPositiveZero(asin32(0.0))); - try expect(math.isNegativeZero(asin32(-0.0))); - try expect(math.isNan(asin32(-2))); - try expect(math.isNan(asin32(1.5))); +test "asinBinary16" { + try testing.expectApproxEqAbs(asinBinary16(-0x1.e4cp-6), -0x1.e4cp-6, math.floatEpsAt(f16, -0x1.e4cp-6)); + try testing.expectApproxEqAbs(asinBinary16(0x1.d68p-1), 0x1.2a8p0, math.floatEpsAt(f16, 0x1.2a8p0)); + try testing.expectApproxEqAbs(asinBinary16(-0x1.a4cp-1), -0x1.eep-1, math.floatEpsAt(f16, -0x1.eep-1)); + try testing.expectApproxEqAbs(asinBinary16(-0x1.0a4p-2), -0x1.0d4p-2, math.floatEpsAt(f16, -0x1.0d4p-2)); + try testing.expectApproxEqAbs(asinBinary16(0x1.28cp-1), 0x1.3c8p-1, math.floatEpsAt(f16, 0x1.3c8p-1)); + try testing.expectApproxEqAbs(asinBinary16(0x1.284p-3), 0x1.298p-3, math.floatEpsAt(f16, 0x1.298p-3)); + try testing.expectApproxEqAbs(asinBinary16(-0x1.574p-1), -0x1.784p-1, math.floatEpsAt(f16, -0x1.784p-1)); + try testing.expectApproxEqAbs(asinBinary16(-0x1.4ccp-1), -0x1.6a4p-1, math.floatEpsAt(f16, -0x1.6a4p-1)); + try testing.expectApproxEqAbs(asinBinary16(0x1.a18p-1), 0x1.e84p-1, math.floatEpsAt(f16, 0x1.e84p-1)); + try testing.expectApproxEqAbs(asinBinary16(0x1.7a8p-2), 0x1.83cp-2, math.floatEpsAt(f16, 0x1.83cp-2)); } -test "asin64.special" { - try expect(math.isPositiveZero(asin64(0.0))); - try expect(math.isNegativeZero(asin64(-0.0))); - try expect(math.isNan(asin64(-2))); - try expect(math.isNan(asin64(1.5))); +test "asinBinary32.special" { + try testing.expectApproxEqAbs(asinBinary32(0x1p+0), 0x1.921fb6p+0, math.floatEpsAt(f32, 0x1.921fb6p+0)); + try testing.expectApproxEqAbs(asinBinary32(-0x1p+0), -0x1.921fb6p+0, math.floatEpsAt(f32, -0x1.921fb6p+0)); + try testing.expectEqual(asinBinary32(0x0p+0), 0x0p+0); + try testing.expectEqual(asinBinary32(-0x0p+0), 0x0p+0); + try testing.expect(math.isNan(asinBinary32(0x1.000002p+0))); + try testing.expect(math.isNan(asinBinary32(-0x1.000002p+0))); + try testing.expect(math.isNan(asinBinary32(math.inf(f32)))); + try testing.expect(math.isNan(asinBinary32(-math.inf(f32)))); + try testing.expect(math.isNan(asinBinary32(math.nan(f32)))); +} + +test "asinBinary32" { + try testing.expectApproxEqAbs(asinBinary32(-0x1.4c2906p-4), -0x1.4c868p-4, math.floatEpsAt(f32, -0x1.4c868p-4)); + try testing.expectApproxEqAbs(asinBinary32(0x1.05fcfap-1), 0x1.130648p-1, math.floatEpsAt(f32, 0x1.130648p-1)); + try testing.expectApproxEqAbs(asinBinary32(0x1.fab976p-2), 0x1.090abcp-1, math.floatEpsAt(f32, 0x1.090abcp-1)); + try testing.expectApproxEqAbs(asinBinary32(0x1.8b4b8cp-1), 0x1.c39fa2p-1, math.floatEpsAt(f32, 0x1.c39fa2p-1)); + try testing.expectApproxEqAbs(asinBinary32(0x1.7117c2p-1), 0x1.9c332p-1, math.floatEpsAt(f32, 0x1.9c332p-1)); + try testing.expectApproxEqAbs(asinBinary32(0x1.e5e112p-5), 0x1.e62a1cp-5, math.floatEpsAt(f32, 0x1.e62a1cp-5)); + try testing.expectApproxEqAbs(asinBinary32(-0x1.07673p-2), -0x1.0a65dep-2, math.floatEpsAt(f32, -0x1.0a65dep-2)); + try testing.expectApproxEqAbs(asinBinary32(-0x1.2108dep-2), -0x1.25046p-2, math.floatEpsAt(f32, -0x1.25046p-2)); + try testing.expectApproxEqAbs(asinBinary32(-0x1.4e6e6cp-1), -0x1.6c6f0cp-1, math.floatEpsAt(f32, -0x1.6c6f0cp-1)); + try testing.expectApproxEqAbs(asinBinary32(0x1.22a16ap-1), 0x1.350f7ap-1, math.floatEpsAt(f32, 0x1.350f7ap-1)); +} + +test "asinBinary64.special" { + try testing.expectApproxEqAbs(asinBinary64(0x1p+0), 0x1.921fb54442d18p+0, math.floatEpsAt(f64, 0x1.921fb54442d18p+0)); + try testing.expectApproxEqAbs(asinBinary64(-0x1p+0), -0x1.921fb54442d18p+0, math.floatEpsAt(f64, -0x1.921fb54442d18p+0)); + try testing.expectEqual(asinBinary64(0x0p+0), 0x0p+0); + try testing.expectEqual(asinBinary64(-0x0p+0), 0x0p+0); + try testing.expect(math.isNan(asinBinary64(0x1.000002p+0))); + try testing.expect(math.isNan(asinBinary64(-0x1.000002p+0))); + try testing.expect(math.isNan(asinBinary64(math.inf(f64)))); + try testing.expect(math.isNan(asinBinary64(-math.inf(f64)))); + try testing.expect(math.isNan(asinBinary64(math.nan(f64)))); +} + +test "asinBinary64" { + try testing.expectApproxEqAbs(asinBinary64(0x1.e674fba3e40d5p-2), 0x1.fae86c5941692p-2, math.floatEpsAt(f64, 0x1.fae86c5941692p-2)); + try testing.expectApproxEqAbs(asinBinary64(-0x1.30fd0566fd979p-1), -0x1.46b6ad730c93ap-1, math.floatEpsAt(f64, -0x1.46b6ad730c93ap-1)); + try testing.expectApproxEqAbs(asinBinary64(0x1.6444a25abfeaap-2), 0x1.6be0be8074eep-2, math.floatEpsAt(f64, 0x1.6be0be8074eep-2)); + try testing.expectApproxEqAbs(asinBinary64(0x1.40a53228d1a13p-1), 0x1.5a7e98f53f717p-1, math.floatEpsAt(f64, 0x1.5a7e98f53f717p-1)); + try testing.expectApproxEqAbs(asinBinary64(-0x1.ccc6d64845cfdp-1), -0x1.1ea2602d14e8p0, math.floatEpsAt(f64, -0x1.1ea2602d14e8p0)); + try testing.expectApproxEqAbs(asinBinary64(-0x1.94bd91b7fc74bp-1), -0x1.d2c2634193158p-1, math.floatEpsAt(f64, -0x1.d2c2634193158p-1)); + try testing.expectApproxEqAbs(asinBinary64(-0x1.8d741b5797fccp-2), -0x1.982d5f1895d2p-2, math.floatEpsAt(f64, -0x1.982d5f1895d2p-2)); + try testing.expectApproxEqAbs(asinBinary64(-0x1.3e8e7e15881c5p-3), -0x1.3fdaf7dfdc864p-3, math.floatEpsAt(f64, -0x1.3fdaf7dfdc864p-3)); + try testing.expectApproxEqAbs(asinBinary64(-0x1.88222d8ab8ca9p-2), -0x1.9269540735b7bp-2, math.floatEpsAt(f64, -0x1.9269540735b7bp-2)); + try testing.expectApproxEqAbs(asinBinary64(-0x1.41c0e9babcbd2p-2), -0x1.474c4c6625527p-2, math.floatEpsAt(f64, -0x1.474c4c6625527p-2)); +} + +test "asinExtended80.special" { + try testing.expectApproxEqAbs(asinExtended80(0x1p+0), 0x1.921fb54442d1846ap+0, math.floatEpsAt(f80, 0x1.921fb54442d1846ap+0)); + try testing.expectApproxEqAbs(asinExtended80(-0x1p+0), -0x1.921fb54442d1846ap+0, math.floatEpsAt(f80, -0x1.921fb54442d1846ap+0)); + try testing.expectEqual(asinExtended80(0x0p+0), 0x0p+0); + try testing.expectEqual(asinExtended80(-0x0p+0), 0x0p+0); + try testing.expect(math.isNan(asinExtended80(0x1.0000000000000002p+0))); + try testing.expect(math.isNan(asinExtended80(-0x1.0000000000000002p+0))); + try testing.expect(math.isNan(asinExtended80(math.inf(f80)))); + try testing.expect(math.isNan(asinExtended80(-math.inf(f80)))); + try testing.expect(math.isNan(asinExtended80(math.nan(f80)))); +} + +test "asinExtended80" { + try testing.expectApproxEqAbs(asinExtended80(0x1.63cf98bc52ce0da8p-9), 0x1.63cfb560149daa9p-9, math.floatEpsAt(f80, 0x1.63cfb560149daa9p-9)); + try testing.expectApproxEqAbs(asinExtended80(-0x1.0473756f7ae930dp-1), -0x1.113cbacd8cd1b96cp-1, math.floatEpsAt(f80, -0x1.113cbacd8cd1b96cp-1)); + try testing.expectApproxEqAbs(asinExtended80(-0x1.2310057e005cc288p-2), -0x1.2721b231d197b064p-2, math.floatEpsAt(f80, -0x1.2721b231d197b064p-2)); + try testing.expectApproxEqAbs(asinExtended80(0x1.f13b03bd685d96eap-1), 0x1.547c408c5d2b05aap0, math.floatEpsAt(f80, 0x1.547c408c5d2b05aap0)); + try testing.expectApproxEqAbs(asinExtended80(-0x1.d5c507e3ef84041cp-1), -0x1.296b76bfadbb5cecp0, math.floatEpsAt(f80, -0x1.296b76bfadbb5cecp0)); + try testing.expectApproxEqAbs(asinExtended80(0x1.8222cbc9147153d8p-1), 0x1.b572da8729a84f2ap-1, math.floatEpsAt(f80, 0x1.b572da8729a84f2ap-1)); + try testing.expectApproxEqAbs(asinExtended80(-0x1.42c9e6b4a088a246p-11), -0x1.42c9e80ac0524dap-11, math.floatEpsAt(f80, -0x1.42c9e80ac0524dap-11)); + try testing.expectApproxEqAbs(asinExtended80(-0x1.8f78d49deadb521cp-3), -0x1.920ca86aef6c3028p-3, math.floatEpsAt(f80, -0x1.920ca86aef6c3028p-3)); + try testing.expectApproxEqAbs(asinExtended80(-0x1.ab98792783515774p-2), -0x1.b91cb4f7204d92fp-2, math.floatEpsAt(f80, -0x1.b91cb4f7204d92fp-2)); + try testing.expectApproxEqAbs(asinExtended80(-0x1.104fe30cef6800aap-1), -0x1.1f20815fdc4c5304p-1, math.floatEpsAt(f80, -0x1.1f20815fdc4c5304p-1)); +} + +test "asinBinary128.special" { + try testing.expectApproxEqAbs(asinBinary128(0x1p+0), 0x1.921fb54442d18469898cc51701b8p0, math.floatEpsAt(f128, 0x1.921fb54442d18469898cc51701b8p0)); + try testing.expectApproxEqAbs(asinBinary128(-0x1p+0), -0x1.921fb54442d18469898cc51701b8p0, math.floatEpsAt(f128, -0x1.921fb54442d18469898cc51701b8p0)); + try testing.expectEqual(asinBinary128(0x0p+0), 0x0p+0); + try testing.expectEqual(asinBinary128(-0x0p+0), 0x0p+0); + try testing.expect(math.isNan(asinBinary128(0x1.0000000000000000000000000001p0))); + try testing.expect(math.isNan(asinBinary128(-0x1.0000000000000000000000000001p0))); + try testing.expect(math.isNan(asinBinary128(math.inf(f128)))); + try testing.expect(math.isNan(asinBinary128(-math.inf(f128)))); + try testing.expect(math.isNan(asinBinary128(math.nan(f128)))); +} + +test "asinBinary128" { + try testing.expectApproxEqAbs(asinBinary128(0x1.85868ce287ca0196b01c25fec5ffp-3), 0x1.87e9c740d7837f8e8fa667988fbep-3, math.floatEpsAt(f128, 0x1.87e9c740d7837f8e8fa667988fbep-3)); + try testing.expectApproxEqAbs(asinBinary128(0x1.8718d6d30b4daed08d04ef59f478p-1), 0x1.bd11a474e864213b48e0f005f1f4p-1, math.floatEpsAt(f128, 0x1.bd11a474e864213b48e0f005f1f4p-1)); + try testing.expectApproxEqAbs(asinBinary128(0x1.11a67640cd7f0ba5d5e362f3abfap-1), 0x1.20b56f8b42649fe72d1f8d68a378p-1, math.floatEpsAt(f128, 0x1.20b56f8b42649fe72d1f8d68a378p-1)); + try testing.expectApproxEqAbs(asinBinary128(-0x1.bd13bf14a9dce22188e52650daa7p-1), -0x1.0dc3a7ddb9736e5ad699bf338566p0, math.floatEpsAt(f128, -0x1.0dc3a7ddb9736e5ad699bf338566p0)); + try testing.expectApproxEqAbs(asinBinary128(-0x1.dee0bc217fc462af57c484eefa71p-2), -0x1.f250716038f70fa50a5826c03802p-2, math.floatEpsAt(f128, -0x1.f250716038f70fa50a5826c03802p-2)); + try testing.expectApproxEqAbs(asinBinary128(-0x1.ea7df9139371c10b9d6fd2bbccd3p-1), -0x1.47a8b4cdd327f90056722feddbabp0, math.floatEpsAt(f128, -0x1.47a8b4cdd327f90056722feddbabp0)); + try testing.expectApproxEqAbs(asinBinary128(0x1.04aaea6de3b5a616460702f26dfcp-2), 0x1.079178d52be662dec67e2cd7f6e9p-2, math.floatEpsAt(f128, 0x1.079178d52be662dec67e2cd7f6e9p-2)); + try testing.expectApproxEqAbs(asinBinary128(-0x1.c7ea85e6b61be666435a7d99444cp-1), -0x1.192df5a8d71702cf1e27014887b2p0, math.floatEpsAt(f128, -0x1.192df5a8d71702cf1e27014887b2p0)); + try testing.expectApproxEqAbs(asinBinary128(-0x1.6e210214e40edf6c8479998189d1p-1), -0x1.97f1092fd94ac0fdfddae2e1222bp-1, math.floatEpsAt(f128, -0x1.97f1092fd94ac0fdfddae2e1222bp-1)); + try testing.expectApproxEqAbs(asinBinary128(-0x1.95061bf93ed6986a45d20f0e1064p-3), -0x1.97b62bc5ae6512093828828325e1p-3, math.floatEpsAt(f128, -0x1.97b62bc5ae6512093828828325e1p-3)); }