diff options
| author | Andrew Kelley <andrew@ziglang.org> | 2022-04-26 10:13:55 -0700 |
|---|---|---|
| committer | Andrew Kelley <andrew@ziglang.org> | 2022-04-27 12:20:44 -0700 |
| commit | 41dd2beaacade94c5c98400a4a655aea07b9e2f3 (patch) | |
| tree | d7cd75c3ded0e8517e801f62dbb883d93f3cd585 /lib/std | |
| parent | 6f4343b61afe36a709e713735947561a2b76bce8 (diff) | |
| download | zig-41dd2beaacade94c5c98400a4a655aea07b9e2f3.tar.gz zig-41dd2beaacade94c5c98400a4a655aea07b9e2f3.zip | |
compiler-rt: math functions reorg
* unify the logic for exporting math functions from compiler-rt,
with the appropriate suffixes and prefixes.
- add all missing f128 and f80 exports. Functions with missing
implementations call other functions and have TODO comments.
- also add f16 functions
* move math functions from freestanding libc to compiler-rt (#7265)
* enable all the f128 and f80 code in the stage2 compiler and behavior
tests (#11161).
* update std lib to use builtins rather than `std.math`.
Diffstat (limited to 'lib/std')
61 files changed, 2519 insertions, 2755 deletions
diff --git a/lib/std/fmt/errol.zig b/lib/std/fmt/errol.zig index 29dd2b7a63..1ce72de0fc 100644 --- a/lib/std/fmt/errol.zig +++ b/lib/std/fmt/errol.zig @@ -113,7 +113,7 @@ fn errolSlow(val: f64, buffer: []u8) FloatDecimal { // normalize the midpoint const e = math.frexp(val).exponent; - var exp = @floatToInt(i16, math.floor(307 + @intToFloat(f64, e) * 0.30103)); + var exp = @floatToInt(i16, @floor(307 + @intToFloat(f64, e) * 0.30103)); if (exp < 20) { exp = 20; } else if (@intCast(usize, exp) >= lookup_table.len) { @@ -170,10 +170,10 @@ fn errolSlow(val: f64, buffer: []u8) FloatDecimal { // digit generation var buf_index: usize = 0; while (true) { - var hdig = @floatToInt(u8, math.floor(high.val)); + var hdig = @floatToInt(u8, @floor(high.val)); if ((high.val == @intToFloat(f64, hdig)) and (high.off < 0)) hdig -= 1; - var ldig = @floatToInt(u8, math.floor(low.val)); + var ldig = @floatToInt(u8, @floor(low.val)); if ((low.val == @intToFloat(f64, ldig)) and (low.off < 0)) ldig -= 1; if (ldig != hdig) break; @@ -187,7 +187,7 @@ fn errolSlow(val: f64, buffer: []u8) FloatDecimal { } const tmp = (high.val + low.val) / 2.0; - var mdig = @floatToInt(u8, math.floor(tmp + 0.5)); + var mdig = @floatToInt(u8, @floor(tmp + 0.5)); if ((@intToFloat(f64, mdig) - tmp) == 0.5 and (mdig & 0x1) != 0) mdig -= 1; buffer[buf_index] = mdig + '0'; diff --git a/lib/std/math.zig b/lib/std/math.zig index b229c8973e..214ade39ce 100644 --- a/lib/std/math.zig +++ b/lib/std/math.zig @@ -138,7 +138,7 @@ pub fn approxEqAbs(comptime T: type, x: T, y: T, tolerance: T) bool { if (isNan(x) or isNan(y)) return false; - return fabs(x - y) <= tolerance; + return @fabs(x - y) <= tolerance; } /// Performs an approximate comparison of two floating point values `x` and `y`. @@ -166,7 +166,7 @@ pub fn approxEqRel(comptime T: type, x: T, y: T, tolerance: T) bool { if (isNan(x) or isNan(y)) return false; - return fabs(x - y) <= max(fabs(x), fabs(y)) * tolerance; + return @fabs(x - y) <= max(@fabs(x), @fabs(y)) * tolerance; } pub fn approxEq(comptime T: type, x: T, y: T, tolerance: T) bool { @@ -233,11 +233,6 @@ pub fn raiseDivByZero() void { pub const isNan = @import("math/isnan.zig").isNan; pub const isSignalNan = @import("math/isnan.zig").isSignalNan; -pub const fabs = @import("math/fabs.zig").fabs; -pub const ceil = @import("math/ceil.zig").ceil; -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 Frexp = @import("math/frexp.zig").Frexp; pub const modf = @import("math/modf.zig").modf; @@ -261,8 +256,6 @@ pub const asin = @import("math/asin.zig").asin; pub const atan = @import("math/atan.zig").atan; pub const atan2 = @import("math/atan2.zig").atan2; pub const hypot = @import("math/hypot.zig").hypot; -pub const exp = @import("math/exp.zig").exp; -pub const exp2 = @import("math/exp2.zig").exp2; pub const expm1 = @import("math/expm1.zig").expm1; pub const ilogb = @import("math/ilogb.zig").ilogb; pub const ln = @import("math/ln.zig").ln; @@ -270,16 +263,12 @@ pub const log = @import("math/log.zig").log; pub const log2 = @import("math/log2.zig").log2; pub const log10 = @import("math/log10.zig").log10; pub const log1p = @import("math/log1p.zig").log1p; -pub const fma = @import("math/fma.zig").fma; pub const asinh = @import("math/asinh.zig").asinh; pub const acosh = @import("math/acosh.zig").acosh; pub const atanh = @import("math/atanh.zig").atanh; pub const sinh = @import("math/sinh.zig").sinh; pub const cosh = @import("math/cosh.zig").cosh; pub const tanh = @import("math/tanh.zig").tanh; -pub const cos = @import("math/cos.zig").cos; -pub const sin = @import("math/sin.zig").sin; -pub const tan = @import("math/tan.zig").tan; pub const complex = @import("math/complex.zig"); pub const Complex = complex.Complex; @@ -716,17 +705,6 @@ fn testAbsInt() !void { try testing.expect((absInt(@as(i32, 10)) catch unreachable) == 10); } -pub const absFloat = fabs; - -test "absFloat" { - try testAbsFloat(); - comptime try testAbsFloat(); -} -fn testAbsFloat() !void { - try testing.expect(absFloat(@as(f32, -10.05)) == 10.05); - try testing.expect(absFloat(@as(f32, 10.05)) == 10.05); -} - /// Divide numerator by denominator, rounding toward zero. Returns an /// error on overflow or when denominator is zero. pub fn divTrunc(comptime T: type, numerator: T, denominator: T) !T { @@ -1400,11 +1378,6 @@ test "order.compare" { try testing.expect(order(1, 0).compare(.neq)); } -test "comptime sin and ln" { - const v = comptime (sin(@as(f32, 1)) + ln(@as(f32, 5))); - try testing.expect(v == sin(@as(f32, 1)) + ln(@as(f32, 5))); -} - /// Returns a mask of all ones if value is true, /// and a mask of all zeroes if value is false. /// Compiles to one instruction for register sized integers. diff --git a/lib/std/math/acos.zig b/lib/std/math/acos.zig index b90ba9c78e..e88bed7227 100644 --- a/lib/std/math/acos.zig +++ b/lib/std/math/acos.zig @@ -64,14 +64,14 @@ fn acos32(x: f32) f32 { // x < -0.5 if (hx >> 31 != 0) { const z = (1 + x) * 0.5; - const s = math.sqrt(z); + const s = @sqrt(z); const w = r32(z) * s - pio2_lo; return 2 * (pio2_hi - (s + w)); } // x > 0.5 const z = (1.0 - x) * 0.5; - const s = math.sqrt(z); + const s = @sqrt(z); const jx = @bitCast(u32, s); const df = @bitCast(f32, jx & 0xFFFFF000); const c = (z - df * df) / (s + df); @@ -133,14 +133,14 @@ fn acos64(x: f64) f64 { // x < -0.5 if (hx >> 31 != 0) { const z = (1.0 + x) * 0.5; - const s = math.sqrt(z); + const s = @sqrt(z); const w = r64(z) * s - pio2_lo; return 2 * (pio2_hi - (s + w)); } // x > 0.5 const z = (1.0 - x) * 0.5; - const s = math.sqrt(z); + const s = @sqrt(z); const jx = @bitCast(u64, s); const df = @bitCast(f64, jx & 0xFFFFFFFF00000000); const c = (z - df * df) / (s + df); diff --git a/lib/std/math/acosh.zig b/lib/std/math/acosh.zig index e42f4fd5d3..a78130d2ef 100644 --- a/lib/std/math/acosh.zig +++ b/lib/std/math/acosh.zig @@ -29,15 +29,15 @@ fn acosh32(x: f32) f32 { // |x| < 2, invalid if x < 1 or nan if (i < 0x3F800000 + (1 << 23)) { - return math.log1p(x - 1 + math.sqrt((x - 1) * (x - 1) + 2 * (x - 1))); + return math.log1p(x - 1 + @sqrt((x - 1) * (x - 1) + 2 * (x - 1))); } // |x| < 0x1p12 else if (i < 0x3F800000 + (12 << 23)) { - return math.ln(2 * x - 1 / (x + math.sqrt(x * x - 1))); + return @log(2 * x - 1 / (x + @sqrt(x * x - 1))); } // |x| >= 0x1p12 else { - return math.ln(x) + 0.693147180559945309417232121458176568; + return @log(x) + 0.693147180559945309417232121458176568; } } @@ -47,15 +47,15 @@ fn acosh64(x: f64) f64 { // |x| < 2, invalid if x < 1 or nan if (e < 0x3FF + 1) { - return math.log1p(x - 1 + math.sqrt((x - 1) * (x - 1) + 2 * (x - 1))); + return math.log1p(x - 1 + @sqrt((x - 1) * (x - 1) + 2 * (x - 1))); } // |x| < 0x1p26 else if (e < 0x3FF + 26) { - return math.ln(2 * x - 1 / (x + math.sqrt(x * x - 1))); + return @log(2 * x - 1 / (x + @sqrt(x * x - 1))); } // |x| >= 0x1p26 or nan else { - return math.ln(x) + 0.693147180559945309417232121458176568; + return @log(x) + 0.693147180559945309417232121458176568; } } diff --git a/lib/std/math/asin.zig b/lib/std/math/asin.zig index 0849fac72e..48ad04c579 100644 --- a/lib/std/math/asin.zig +++ b/lib/std/math/asin.zig @@ -60,8 +60,8 @@ fn asin32(x: f32) f32 { } // 1 > |x| >= 0.5 - const z = (1 - math.fabs(x)) * 0.5; - const s = math.sqrt(z); + const z = (1 - @fabs(x)) * 0.5; + const s = @sqrt(z); const fx = pio2 - 2 * (s + s * r32(z)); if (hx >> 31 != 0) { @@ -119,8 +119,8 @@ fn asin64(x: f64) f64 { } // 1 > |x| >= 0.5 - const z = (1 - math.fabs(x)) * 0.5; - const s = math.sqrt(z); + const z = (1 - @fabs(x)) * 0.5; + const s = @sqrt(z); const r = r64(z); var fx: f64 = undefined; diff --git a/lib/std/math/asinh.zig b/lib/std/math/asinh.zig index 8717ebbb66..65028ef5d9 100644 --- a/lib/std/math/asinh.zig +++ b/lib/std/math/asinh.zig @@ -39,15 +39,15 @@ fn asinh32(x: f32) f32 { // |x| >= 0x1p12 or inf or nan if (i >= 0x3F800000 + (12 << 23)) { - rx = math.ln(rx) + 0.69314718055994530941723212145817656; + rx = @log(rx) + 0.69314718055994530941723212145817656; } // |x| >= 2 else if (i >= 0x3F800000 + (1 << 23)) { - rx = math.ln(2 * x + 1 / (math.sqrt(x * x + 1) + x)); + rx = @log(2 * x + 1 / (@sqrt(x * x + 1) + x)); } // |x| >= 0x1p-12, up to 1.6ulp error else if (i >= 0x3F800000 - (12 << 23)) { - rx = math.log1p(x + x * x / (math.sqrt(x * x + 1) + 1)); + rx = math.log1p(x + x * x / (@sqrt(x * x + 1) + 1)); } // |x| < 0x1p-12, inexact if x != 0 else { @@ -70,15 +70,15 @@ fn asinh64(x: f64) f64 { // |x| >= 0x1p26 or inf or nan if (e >= 0x3FF + 26) { - rx = math.ln(rx) + 0.693147180559945309417232121458176568; + rx = @log(rx) + 0.693147180559945309417232121458176568; } // |x| >= 2 else if (e >= 0x3FF + 1) { - rx = math.ln(2 * x + 1 / (math.sqrt(x * x + 1) + x)); + rx = @log(2 * x + 1 / (@sqrt(x * x + 1) + x)); } // |x| >= 0x1p-12, up to 1.6ulp error else if (e >= 0x3FF - 26) { - rx = math.log1p(x + x * x / (math.sqrt(x * x + 1) + 1)); + rx = math.log1p(x + x * x / (@sqrt(x * x + 1) + 1)); } // |x| < 0x1p-12, inexact if x != 0 else { diff --git a/lib/std/math/atan.zig b/lib/std/math/atan.zig index c67e6fe8e0..3a13d943e8 100644 --- a/lib/std/math/atan.zig +++ b/lib/std/math/atan.zig @@ -73,7 +73,7 @@ fn atan32(x_: f32) f32 { } id = null; } else { - x = math.fabs(x); + x = @fabs(x); // |x| < 1.1875 if (ix < 0x3F980000) { // 7/16 <= |x| < 11/16 @@ -171,7 +171,7 @@ fn atan64(x_: f64) f64 { } id = null; } else { - x = math.fabs(x); + x = @fabs(x); // |x| < 1.1875 if (ix < 0x3FF30000) { // 7/16 <= |x| < 11/16 diff --git a/lib/std/math/atan2.zig b/lib/std/math/atan2.zig index d440d65e04..b9b37e7da4 100644 --- a/lib/std/math/atan2.zig +++ b/lib/std/math/atan2.zig @@ -108,7 +108,7 @@ fn atan2_32(y: f32, x: f32) f32 { if ((m & 2) != 0 and iy + (26 << 23) < ix) { break :z 0.0; } else { - break :z math.atan(math.fabs(y / x)); + break :z math.atan(@fabs(y / x)); } }; @@ -198,7 +198,7 @@ fn atan2_64(y: f64, x: f64) f64 { if ((m & 2) != 0 and iy +% (64 << 20) < ix) { break :z 0.0; } else { - break :z math.atan(math.fabs(y / x)); + break :z math.atan(@fabs(y / x)); } }; diff --git a/lib/std/math/complex.zig b/lib/std/math/complex.zig index 42342faa3e..2fd1cf15a1 100644 --- a/lib/std/math/complex.zig +++ b/lib/std/math/complex.zig @@ -115,7 +115,7 @@ pub fn Complex(comptime T: type) type { /// Returns the magnitude of a complex number. pub fn magnitude(self: Self) T { - return math.sqrt(self.re * self.re + self.im * self.im); + return @sqrt(self.re * self.re + self.im * self.im); } }; } diff --git a/lib/std/math/complex/atan.zig b/lib/std/math/complex/atan.zig index 484b41edf5..929b98aebd 100644 --- a/lib/std/math/complex/atan.zig +++ b/lib/std/math/complex/atan.zig @@ -66,7 +66,7 @@ fn atan32(z: Complex(f32)) Complex(f32) { t = y + 1.0; a = (x2 + (t * t)) / a; - return Complex(f32).init(w, 0.25 * math.ln(a)); + return Complex(f32).init(w, 0.25 * @log(a)); } fn redupif64(x: f64) f64 { @@ -115,7 +115,7 @@ fn atan64(z: Complex(f64)) Complex(f64) { t = y + 1.0; a = (x2 + (t * t)) / a; - return Complex(f64).init(w, 0.25 * math.ln(a)); + return Complex(f64).init(w, 0.25 * @log(a)); } const epsilon = 0.0001; diff --git a/lib/std/math/complex/cosh.zig b/lib/std/math/complex/cosh.zig index 46f7a714a2..719d0f28cd 100644 --- a/lib/std/math/complex/cosh.zig +++ b/lib/std/math/complex/cosh.zig @@ -44,12 +44,12 @@ fn cosh32(z: Complex(f32)) Complex(f32) { // |x|>= 9, so cosh(x) ~= exp(|x|) if (ix < 0x42b17218) { // x < 88.7: exp(|x|) won't overflow - const h = math.exp(math.fabs(x)) * 0.5; + const h = @exp(@fabs(x)) * 0.5; return Complex(f32).init(math.copysign(f32, h, x) * math.cos(y), h * math.sin(y)); } // x < 192.7: scale to avoid overflow else if (ix < 0x4340b1e7) { - const v = Complex(f32).init(math.fabs(x), y); + const v = Complex(f32).init(@fabs(x), y); const r = ldexp_cexp(v, -1); return Complex(f32).init(r.re, r.im * math.copysign(f32, 1, x)); } @@ -112,12 +112,12 @@ fn cosh64(z: Complex(f64)) Complex(f64) { // |x|>= 22, so cosh(x) ~= exp(|x|) if (ix < 0x40862e42) { // x < 710: exp(|x|) won't overflow - const h = math.exp(math.fabs(x)) * 0.5; + const h = @exp(@fabs(x)) * 0.5; return Complex(f64).init(h * math.cos(y), math.copysign(f64, h, x) * math.sin(y)); } // x < 1455: scale to avoid overflow else if (ix < 0x4096bbaa) { - const v = Complex(f64).init(math.fabs(x), y); + const v = Complex(f64).init(@fabs(x), y); const r = ldexp_cexp(v, -1); return Complex(f64).init(r.re, r.im * math.copysign(f64, 1, x)); } diff --git a/lib/std/math/complex/exp.zig b/lib/std/math/complex/exp.zig index ce25025ded..4ed731d85c 100644 --- a/lib/std/math/complex/exp.zig +++ b/lib/std/math/complex/exp.zig @@ -33,7 +33,7 @@ fn exp32(z: Complex(f32)) Complex(f32) { const hy = @bitCast(u32, y) & 0x7fffffff; // cexp(x + i0) = exp(x) + i0 if (hy == 0) { - return Complex(f32).init(math.exp(x), y); + return Complex(f32).init(@exp(x), y); } const hx = @bitCast(u32, x); @@ -63,7 +63,7 @@ fn exp32(z: Complex(f32)) Complex(f32) { // - x = +-inf // - x = nan else { - const exp_x = math.exp(x); + const exp_x = @exp(x); return Complex(f32).init(exp_x * math.cos(y), exp_x * math.sin(y)); } } @@ -81,7 +81,7 @@ fn exp64(z: Complex(f64)) Complex(f64) { // cexp(x + i0) = exp(x) + i0 if (hy | ly == 0) { - return Complex(f64).init(math.exp(x), y); + return Complex(f64).init(@exp(x), y); } const fx = @bitCast(u64, x); @@ -114,13 +114,13 @@ fn exp64(z: Complex(f64)) Complex(f64) { // - x = +-inf // - x = nan else { - const exp_x = math.exp(x); + const exp_x = @exp(x); return Complex(f64).init(exp_x * math.cos(y), exp_x * math.sin(y)); } } test "complex.cexp32" { - const tolerance_f32 = math.sqrt(math.floatEps(f32)); + const tolerance_f32 = @sqrt(math.floatEps(f32)); { const a = Complex(f32).init(5, 3); @@ -140,7 +140,7 @@ test "complex.cexp32" { } test "complex.cexp64" { - const tolerance_f64 = math.sqrt(math.floatEps(f64)); + const tolerance_f64 = @sqrt(math.floatEps(f64)); { const a = Complex(f64).init(5, 3); diff --git a/lib/std/math/complex/ldexp.zig b/lib/std/math/complex/ldexp.zig index db710a0438..1c2d06b858 100644 --- a/lib/std/math/complex/ldexp.zig +++ b/lib/std/math/complex/ldexp.zig @@ -26,7 +26,7 @@ fn frexp_exp32(x: f32, expt: *i32) f32 { const k = 235; // reduction constant const kln2 = 162.88958740; // k * ln2 - const exp_x = math.exp(x - kln2); + const exp_x = @exp(x - kln2); const hx = @bitCast(u32, exp_x); // TODO zig should allow this cast implicitly because it should know the value is in range expt.* = @intCast(i32, hx >> 23) - (0x7f + 127) + k; @@ -54,7 +54,7 @@ fn frexp_exp64(x: f64, expt: *i32) f64 { const k = 1799; // reduction constant const kln2 = 1246.97177782734161156; // k * ln2 - const exp_x = math.exp(x - kln2); + const exp_x = @exp(x - kln2); const fx = @bitCast(u64, exp_x); const hx = @intCast(u32, fx >> 32); diff --git a/lib/std/math/complex/log.zig b/lib/std/math/complex/log.zig index 90c51058cf..6d1b06d272 100644 --- a/lib/std/math/complex/log.zig +++ b/lib/std/math/complex/log.zig @@ -10,7 +10,7 @@ pub fn log(z: anytype) Complex(@TypeOf(z.re)) { const r = cmath.abs(z); const phi = cmath.arg(z); - return Complex(T).init(math.ln(r), phi); + return Complex(T).init(@log(r), phi); } const epsilon = 0.0001; diff --git a/lib/std/math/complex/sinh.zig b/lib/std/math/complex/sinh.zig index 851af3e62e..b21f6e59eb 100644 --- a/lib/std/math/complex/sinh.zig +++ b/lib/std/math/complex/sinh.zig @@ -44,12 +44,12 @@ fn sinh32(z: Complex(f32)) Complex(f32) { // |x|>= 9, so cosh(x) ~= exp(|x|) if (ix < 0x42b17218) { // x < 88.7: exp(|x|) won't overflow - const h = math.exp(math.fabs(x)) * 0.5; + const h = @exp(@fabs(x)) * 0.5; return Complex(f32).init(math.copysign(f32, h, x) * math.cos(y), h * math.sin(y)); } // x < 192.7: scale to avoid overflow else if (ix < 0x4340b1e7) { - const v = Complex(f32).init(math.fabs(x), y); + const v = Complex(f32).init(@fabs(x), y); const r = ldexp_cexp(v, -1); return Complex(f32).init(r.re * math.copysign(f32, 1, x), r.im); } @@ -111,12 +111,12 @@ fn sinh64(z: Complex(f64)) Complex(f64) { // |x|>= 22, so cosh(x) ~= exp(|x|) if (ix < 0x40862e42) { // x < 710: exp(|x|) won't overflow - const h = math.exp(math.fabs(x)) * 0.5; + const h = @exp(@fabs(x)) * 0.5; return Complex(f64).init(math.copysign(f64, h, x) * math.cos(y), h * math.sin(y)); } // x < 1455: scale to avoid overflow else if (ix < 0x4096bbaa) { - const v = Complex(f64).init(math.fabs(x), y); + const v = Complex(f64).init(@fabs(x), y); const r = ldexp_cexp(v, -1); return Complex(f64).init(r.re * math.copysign(f64, 1, x), r.im); } diff --git a/lib/std/math/complex/sqrt.zig b/lib/std/math/complex/sqrt.zig index 4f16e631b8..ab24e2d60d 100644 --- a/lib/std/math/complex/sqrt.zig +++ b/lib/std/math/complex/sqrt.zig @@ -43,7 +43,7 @@ fn sqrt32(z: Complex(f32)) Complex(f32) { // sqrt(-inf + i nan) = nan +- inf i // sqrt(-inf + iy) = 0 + inf i if (math.signbit(x)) { - return Complex(f32).init(math.fabs(x - y), math.copysign(f32, x, y)); + return Complex(f32).init(@fabs(x - y), math.copysign(f32, x, y)); } else { return Complex(f32).init(x, math.copysign(f32, y - y, y)); } @@ -56,15 +56,15 @@ fn sqrt32(z: Complex(f32)) Complex(f32) { const dy = @as(f64, y); if (dx >= 0) { - const t = math.sqrt((dx + math.hypot(f64, dx, dy)) * 0.5); + const t = @sqrt((dx + math.hypot(f64, dx, dy)) * 0.5); return Complex(f32).init( @floatCast(f32, t), @floatCast(f32, dy / (2.0 * t)), ); } else { - const t = math.sqrt((-dx + math.hypot(f64, dx, dy)) * 0.5); + const t = @sqrt((-dx + math.hypot(f64, dx, dy)) * 0.5); return Complex(f32).init( - @floatCast(f32, math.fabs(y) / (2.0 * t)), + @floatCast(f32, @fabs(y) / (2.0 * t)), @floatCast(f32, math.copysign(f64, t, y)), ); } @@ -94,7 +94,7 @@ fn sqrt64(z: Complex(f64)) Complex(f64) { // sqrt(-inf + i nan) = nan +- inf i // sqrt(-inf + iy) = 0 + inf i if (math.signbit(x)) { - return Complex(f64).init(math.fabs(x - y), math.copysign(f64, x, y)); + return Complex(f64).init(@fabs(x - y), math.copysign(f64, x, y)); } else { return Complex(f64).init(x, math.copysign(f64, y - y, y)); } @@ -104,7 +104,7 @@ fn sqrt64(z: Complex(f64)) Complex(f64) { // scale to avoid overflow var scale = false; - if (math.fabs(x) >= threshold or math.fabs(y) >= threshold) { + if (@fabs(x) >= threshold or @fabs(y) >= threshold) { x *= 0.25; y *= 0.25; scale = true; @@ -112,11 +112,11 @@ fn sqrt64(z: Complex(f64)) Complex(f64) { var result: Complex(f64) = undefined; if (x >= 0) { - const t = math.sqrt((x + math.hypot(f64, x, y)) * 0.5); + const t = @sqrt((x + math.hypot(f64, x, y)) * 0.5); result = Complex(f64).init(t, y / (2.0 * t)); } else { - const t = math.sqrt((-x + math.hypot(f64, x, y)) * 0.5); - result = Complex(f64).init(math.fabs(y) / (2.0 * t), math.copysign(f64, t, y)); + const t = @sqrt((-x + math.hypot(f64, x, y)) * 0.5); + result = Complex(f64).init(@fabs(y) / (2.0 * t), math.copysign(f64, t, y)); } if (scale) { diff --git a/lib/std/math/complex/tanh.zig b/lib/std/math/complex/tanh.zig index 0960c66679..e61ec1e95b 100644 --- a/lib/std/math/complex/tanh.zig +++ b/lib/std/math/complex/tanh.zig @@ -44,7 +44,7 @@ fn tanh32(z: Complex(f32)) Complex(f32) { // x >= 11 if (ix >= 0x41300000) { - const exp_mx = math.exp(-math.fabs(x)); + const exp_mx = @exp(-@fabs(x)); return Complex(f32).init(math.copysign(f32, 1, x), 4 * math.sin(y) * math.cos(y) * exp_mx * exp_mx); } @@ -52,7 +52,7 @@ fn tanh32(z: Complex(f32)) Complex(f32) { const t = math.tan(y); const beta = 1.0 + t * t; const s = math.sinh(x); - const rho = math.sqrt(1 + s * s); + const rho = @sqrt(1 + s * s); const den = 1 + beta * s * s; return Complex(f32).init((beta * rho * s) / den, t / den); @@ -87,7 +87,7 @@ fn tanh64(z: Complex(f64)) Complex(f64) { // x >= 22 if (ix >= 0x40360000) { - const exp_mx = math.exp(-math.fabs(x)); + const exp_mx = @exp(-@fabs(x)); return Complex(f64).init(math.copysign(f64, 1, x), 4 * math.sin(y) * math.cos(y) * exp_mx * exp_mx); } @@ -95,7 +95,7 @@ fn tanh64(z: Complex(f64)) Complex(f64) { const t = math.tan(y); const beta = 1.0 + t * t; const s = math.sinh(x); - const rho = math.sqrt(1 + s * s); + const rho = @sqrt(1 + s * s); const den = 1 + beta * s * s; return Complex(f64).init((beta * rho * s) / den, t / den); diff --git a/lib/std/math/cosh.zig b/lib/std/math/cosh.zig index c71e82ea1c..d633f2fa0c 100644 --- a/lib/std/math/cosh.zig +++ b/lib/std/math/cosh.zig @@ -45,7 +45,7 @@ fn cosh32(x: f32) f32 { // |x| < log(FLT_MAX) if (ux < 0x42B17217) { - const t = math.exp(ax); + const t = @exp(ax); return 0.5 * (t + 1 / t); } @@ -77,7 +77,7 @@ fn cosh64(x: f64) f64 { // |x| < log(DBL_MAX) if (w < 0x40862E42) { - const t = math.exp(ax); + const t = @exp(ax); // NOTE: If x > log(0x1p26) then 1/t is not required. return 0.5 * (t + 1 / t); } diff --git a/lib/std/math/expo2.zig b/lib/std/math/expo2.zig index f404570fb6..4345233173 100644 --- a/lib/std/math/expo2.zig +++ b/lib/std/math/expo2.zig @@ -22,7 +22,7 @@ fn expo2f(x: f32) f32 { const u = (0x7F + k / 2) << 23; const scale = @bitCast(f32, u); - return math.exp(x - kln2) * scale * scale; + return @exp(x - kln2) * scale * scale; } fn expo2d(x: f64) f64 { @@ -31,5 +31,5 @@ fn expo2d(x: f64) f64 { const u = (0x3FF + k / 2) << 20; const scale = @bitCast(f64, @as(u64, u) << 32); - return math.exp(x - kln2) * scale * scale; + return @exp(x - kln2) * scale * scale; } diff --git a/lib/std/math/fabs.zig b/lib/std/math/fabs.zig deleted file mode 100644 index 44918e75d9..0000000000 --- a/lib/std/math/fabs.zig +++ /dev/null @@ -1,45 +0,0 @@ -const std = @import("../std.zig"); -const math = std.math; -const expect = std.testing.expect; - -/// Returns the absolute value of x. -/// -/// Special Cases: -/// - fabs(+-inf) = +inf -/// - fabs(nan) = nan -pub fn fabs(x: anytype) @TypeOf(x) { - const T = @TypeOf(x); - const TBits = std.meta.Int(.unsigned, @bitSizeOf(T)); - if (@typeInfo(T) != .Float) { - @compileError("fabs not implemented for " ++ @typeName(T)); - } - - const float_bits = @bitCast(TBits, x); - const remove_sign = ~@as(TBits, 0) >> 1; - - return @bitCast(T, float_bits & remove_sign); -} - -test "math.fabs" { - // TODO add support for c_longdouble here - inline for ([_]type{ f16, f32, f64, f80, f128 }) |T| { - // normals - try expect(fabs(@as(T, 1.0)) == 1.0); - try expect(fabs(@as(T, -1.0)) == 1.0); - try expect(fabs(math.floatMin(T)) == math.floatMin(T)); - try expect(fabs(-math.floatMin(T)) == math.floatMin(T)); - try expect(fabs(math.floatMax(T)) == math.floatMax(T)); - try expect(fabs(-math.floatMax(T)) == math.floatMax(T)); - - // subnormals - try expect(fabs(@as(T, 0.0)) == 0.0); - try expect(fabs(@as(T, -0.0)) == 0.0); - try expect(fabs(math.floatTrueMin(T)) == math.floatTrueMin(T)); - try expect(fabs(-math.floatTrueMin(T)) == math.floatTrueMin(T)); - - // non-finite numbers - try expect(math.isPositiveInf(fabs(math.inf(T)))); - try expect(math.isPositiveInf(fabs(-math.inf(T)))); - try expect(math.isNan(fabs(math.nan(T)))); - } -} diff --git a/lib/std/math/hypot.zig b/lib/std/math/hypot.zig index e47a191892..981f6143fe 100644 --- a/lib/std/math/hypot.zig +++ b/lib/std/math/hypot.zig @@ -56,7 +56,7 @@ fn hypot32(x: f32, y: f32) f32 { yy *= 0x1.0p-90; } - return z * math.sqrt(@floatCast(f32, @as(f64, x) * x + @as(f64, y) * y)); + return z * @sqrt(@floatCast(f32, @as(f64, x) * x + @as(f64, y) * y)); } fn sq(hi: *f64, lo: *f64, x: f64) void { @@ -117,7 +117,7 @@ fn hypot64(x: f64, y: f64) f64 { sq(&hx, &lx, x); sq(&hy, &ly, y); - return z * math.sqrt(ly + lx + hy + hx); + return z * @sqrt(ly + lx + hy + hx); } test "math.hypot" { diff --git a/lib/std/math/ln.zig b/lib/std/math/ln.zig index bb352cd6e1..65db861587 100644 --- a/lib/std/math/ln.zig +++ b/lib/std/math/ln.zig @@ -1,12 +1,6 @@ -// Ported from musl, which is licensed under the MIT license: -// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT -// -// https://git.musl-libc.org/cgit/musl/tree/src/math/lnf.c -// https://git.musl-libc.org/cgit/musl/tree/src/math/ln.c - const std = @import("../std.zig"); const math = std.math; -const expect = std.testing.expect; +const testing = std.testing; /// Returns the natural logarithm of x. /// @@ -15,175 +9,26 @@ const expect = std.testing.expect; /// - ln(0) = -inf /// - ln(x) = nan if x < 0 /// - ln(nan) = nan +/// TODO remove this in favor of `@log`. pub fn ln(x: anytype) @TypeOf(x) { const T = @TypeOf(x); switch (@typeInfo(T)) { .ComptimeFloat => { - return @as(comptime_float, ln_64(x)); - }, - .Float => { - return switch (T) { - f32 => ln_32(x), - f64 => ln_64(x), - else => @compileError("ln not implemented for " ++ @typeName(T)), - }; + return @as(comptime_float, @log(x)); }, + .Float => return @log(x), .ComptimeInt => { - return @as(comptime_int, math.floor(ln_64(@as(f64, x)))); + return @as(comptime_int, @floor(@log(@as(f64, x)))); }, .Int => |IntType| switch (IntType.signedness) { .signed => @compileError("ln not implemented for signed integers"), - .unsigned => return @as(T, math.floor(ln_64(@as(f64, x)))), + .unsigned => return @as(T, @floor(@log(@as(f64, x)))), }, else => @compileError("ln not implemented for " ++ @typeName(T)), } } -pub fn ln_32(x_: f32) f32 { - const ln2_hi: f32 = 6.9313812256e-01; - const ln2_lo: f32 = 9.0580006145e-06; - const Lg1: f32 = 0xaaaaaa.0p-24; - const Lg2: f32 = 0xccce13.0p-25; - const Lg3: f32 = 0x91e9ee.0p-25; - const Lg4: f32 = 0xf89e26.0p-26; - - var x = x_; - var ix = @bitCast(u32, x); - var k: i32 = 0; - - // x < 2^(-126) - if (ix < 0x00800000 or ix >> 31 != 0) { - // log(+-0) = -inf - if (ix << 1 == 0) { - return -math.inf(f32); - } - // log(-#) = nan - if (ix >> 31 != 0) { - return math.nan(f32); - } - - // subnormal, scale x - k -= 25; - x *= 0x1.0p25; - ix = @bitCast(u32, x); - } else if (ix >= 0x7F800000) { - return x; - } else if (ix == 0x3F800000) { - return 0; - } - - // x into [sqrt(2) / 2, sqrt(2)] - ix += 0x3F800000 - 0x3F3504F3; - k += @intCast(i32, ix >> 23) - 0x7F; - ix = (ix & 0x007FFFFF) + 0x3F3504F3; - x = @bitCast(f32, ix); - - const f = x - 1.0; - const s = f / (2.0 + f); - const z = s * s; - const w = z * z; - const t1 = w * (Lg2 + w * Lg4); - const t2 = z * (Lg1 + w * Lg3); - const R = t2 + t1; - const hfsq = 0.5 * f * f; - const dk = @intToFloat(f32, k); - - return s * (hfsq + R) + dk * ln2_lo - hfsq + f + dk * ln2_hi; -} - -pub fn ln_64(x_: f64) f64 { - const ln2_hi: f64 = 6.93147180369123816490e-01; - const ln2_lo: f64 = 1.90821492927058770002e-10; - const Lg1: f64 = 6.666666666666735130e-01; - const Lg2: f64 = 3.999999999940941908e-01; - const Lg3: f64 = 2.857142874366239149e-01; - const Lg4: f64 = 2.222219843214978396e-01; - const Lg5: f64 = 1.818357216161805012e-01; - const Lg6: f64 = 1.531383769920937332e-01; - const Lg7: f64 = 1.479819860511658591e-01; - - var x = x_; - var ix = @bitCast(u64, x); - var hx = @intCast(u32, ix >> 32); - var k: i32 = 0; - - if (hx < 0x00100000 or hx >> 31 != 0) { - // log(+-0) = -inf - if (ix << 1 == 0) { - return -math.inf(f64); - } - // log(-#) = nan - if (hx >> 31 != 0) { - return math.nan(f64); - } - - // subnormal, scale x - k -= 54; - x *= 0x1.0p54; - hx = @intCast(u32, @bitCast(u64, ix) >> 32); - } else if (hx >= 0x7FF00000) { - return x; - } else if (hx == 0x3FF00000 and ix << 32 == 0) { - return 0; - } - - // x into [sqrt(2) / 2, sqrt(2)] - hx += 0x3FF00000 - 0x3FE6A09E; - k += @intCast(i32, hx >> 20) - 0x3FF; - hx = (hx & 0x000FFFFF) + 0x3FE6A09E; - ix = (@as(u64, hx) << 32) | (ix & 0xFFFFFFFF); - x = @bitCast(f64, ix); - - const f = x - 1.0; - const hfsq = 0.5 * f * f; - const s = f / (2.0 + f); - const z = s * s; - const w = z * z; - const t1 = w * (Lg2 + w * (Lg4 + w * Lg6)); - const t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7))); - const R = t2 + t1; - const dk = @intToFloat(f64, k); - - return s * (hfsq + R) + dk * ln2_lo - hfsq + f + dk * ln2_hi; -} - test "math.ln" { - try expect(ln(@as(f32, 0.2)) == ln_32(0.2)); - try expect(ln(@as(f64, 0.2)) == ln_64(0.2)); -} - -test "math.ln32" { - const epsilon = 0.000001; - - try expect(math.approxEqAbs(f32, ln_32(0.2), -1.609438, epsilon)); - try expect(math.approxEqAbs(f32, ln_32(0.8923), -0.113953, epsilon)); - try expect(math.approxEqAbs(f32, ln_32(1.5), 0.405465, epsilon)); - try expect(math.approxEqAbs(f32, ln_32(37.45), 3.623007, epsilon)); - try expect(math.approxEqAbs(f32, ln_32(89.123), 4.490017, epsilon)); - try expect(math.approxEqAbs(f32, ln_32(123123.234375), 11.720941, epsilon)); -} - -test "math.ln64" { - const epsilon = 0.000001; - - try expect(math.approxEqAbs(f64, ln_64(0.2), -1.609438, epsilon)); - try expect(math.approxEqAbs(f64, ln_64(0.8923), -0.113953, epsilon)); - try expect(math.approxEqAbs(f64, ln_64(1.5), 0.405465, epsilon)); - try expect(math.approxEqAbs(f64, ln_64(37.45), 3.623007, epsilon)); - try expect(math.approxEqAbs(f64, ln_64(89.123), 4.490017, epsilon)); - try expect(math.approxEqAbs(f64, ln_64(123123.234375), 11.720941, epsilon)); -} - -test "math.ln32.special" { - try expect(math.isPositiveInf(ln_32(math.inf(f32)))); - try expect(math.isNegativeInf(ln_32(0.0))); - try expect(math.isNan(ln_32(-1.0))); - try expect(math.isNan(ln_32(math.nan(f32)))); -} - -test "math.ln64.special" { - try expect(math.isPositiveInf(ln_64(math.inf(f64)))); - try expect(math.isNegativeInf(ln_64(0.0))); - try expect(math.isNan(ln_64(-1.0))); - try expect(math.isNan(ln_64(math.nan(f64)))); + try testing.expect(ln(@as(f32, 0.2)) == @log(0.2)); + try testing.expect(ln(@as(f64, 0.2)) == @log(0.2)); } diff --git a/lib/std/math/log.zig b/lib/std/math/log.zig index 6336726b39..ad2763fa54 100644 --- a/lib/std/math/log.zig +++ b/lib/std/math/log.zig @@ -15,28 +15,28 @@ pub fn log(comptime T: type, base: T, x: T) T { } else if (base == 10) { return math.log10(x); } else if ((@typeInfo(T) == .Float or @typeInfo(T) == .ComptimeFloat) and base == math.e) { - return math.ln(x); + return @log(x); } const float_base = math.lossyCast(f64, base); switch (@typeInfo(T)) { .ComptimeFloat => { - return @as(comptime_float, math.ln(@as(f64, x)) / math.ln(float_base)); + return @as(comptime_float, @log(@as(f64, x)) / @log(float_base)); }, .ComptimeInt => { - return @as(comptime_int, math.floor(math.ln(@as(f64, x)) / math.ln(float_base))); + return @as(comptime_int, @floor(@log(@as(f64, x)) / @log(float_base))); }, // TODO implement integer log without using float math .Int => |IntType| switch (IntType.signedness) { .signed => @compileError("log not implemented for signed integers"), - .unsigned => return @floatToInt(T, math.floor(math.ln(@intToFloat(f64, x)) / math.ln(float_base))), + .unsigned => return @floatToInt(T, @floor(@log(@intToFloat(f64, x)) / @log(float_base))), }, .Float => { switch (T) { - f32 => return @floatCast(f32, math.ln(@as(f64, x)) / math.ln(float_base)), - f64 => return math.ln(x) / math.ln(float_base), + f32 => return @floatCast(f32, @log(@as(f64, x)) / @log(float_base)), + f64 => return @log(x) / @log(float_base), else => @compileError("log not implemented for " ++ @typeName(T)), } }, diff --git a/lib/std/math/log10.zig b/lib/std/math/log10.zig index 84eced85f0..4f13426079 100644 --- a/lib/std/math/log10.zig +++ b/lib/std/math/log10.zig @@ -1,9 +1,3 @@ -// Ported from musl, which is licensed under the MIT license: -// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT -// -// https://git.musl-libc.org/cgit/musl/tree/src/math/log10f.c -// https://git.musl-libc.org/cgit/musl/tree/src/math/log10.c - const std = @import("../std.zig"); const math = std.math; const testing = std.testing; @@ -20,198 +14,16 @@ pub fn log10(x: anytype) @TypeOf(x) { const T = @TypeOf(x); switch (@typeInfo(T)) { .ComptimeFloat => { - return @as(comptime_float, log10_64(x)); - }, - .Float => { - return switch (T) { - f32 => log10_32(x), - f64 => log10_64(x), - else => @compileError("log10 not implemented for " ++ @typeName(T)), - }; + return @as(comptime_float, @log10(x)); }, + .Float => return @log10(x), .ComptimeInt => { - return @as(comptime_int, math.floor(log10_64(@as(f64, x)))); + return @as(comptime_int, @floor(@log10(@as(f64, x)))); }, .Int => |IntType| switch (IntType.signedness) { .signed => @compileError("log10 not implemented for signed integers"), - .unsigned => return @floatToInt(T, math.floor(log10_64(@intToFloat(f64, x)))), + .unsigned => return @floatToInt(T, @floor(@log10(@intToFloat(f64, x)))), }, else => @compileError("log10 not implemented for " ++ @typeName(T)), } } - -pub fn log10_32(x_: f32) f32 { - const ivln10hi: f32 = 4.3432617188e-01; - const ivln10lo: f32 = -3.1689971365e-05; - const log10_2hi: f32 = 3.0102920532e-01; - const log10_2lo: f32 = 7.9034151668e-07; - const Lg1: f32 = 0xaaaaaa.0p-24; - const Lg2: f32 = 0xccce13.0p-25; - const Lg3: f32 = 0x91e9ee.0p-25; - const Lg4: f32 = 0xf89e26.0p-26; - - var x = x_; - var u = @bitCast(u32, x); - var ix = u; - var k: i32 = 0; - - // x < 2^(-126) - if (ix < 0x00800000 or ix >> 31 != 0) { - // log(+-0) = -inf - if (ix << 1 == 0) { - return -math.inf(f32); - } - // log(-#) = nan - if (ix >> 31 != 0) { - return math.nan(f32); - } - - k -= 25; - x *= 0x1.0p25; - ix = @bitCast(u32, x); - } else if (ix >= 0x7F800000) { - return x; - } else if (ix == 0x3F800000) { - return 0; - } - - // x into [sqrt(2) / 2, sqrt(2)] - ix += 0x3F800000 - 0x3F3504F3; - k += @intCast(i32, ix >> 23) - 0x7F; - ix = (ix & 0x007FFFFF) + 0x3F3504F3; - x = @bitCast(f32, ix); - - const f = x - 1.0; - const s = f / (2.0 + f); - const z = s * s; - const w = z * z; - const t1 = w * (Lg2 + w * Lg4); - const t2 = z * (Lg1 + w * Lg3); - const R = t2 + t1; - const hfsq = 0.5 * f * f; - - var hi = f - hfsq; - u = @bitCast(u32, hi); - u &= 0xFFFFF000; - hi = @bitCast(f32, u); - const lo = f - hi - hfsq + s * (hfsq + R); - const dk = @intToFloat(f32, k); - - return dk * log10_2lo + (lo + hi) * ivln10lo + lo * ivln10hi + hi * ivln10hi + dk * log10_2hi; -} - -pub fn log10_64(x_: f64) f64 { - const ivln10hi: f64 = 4.34294481878168880939e-01; - const ivln10lo: f64 = 2.50829467116452752298e-11; - const log10_2hi: f64 = 3.01029995663611771306e-01; - const log10_2lo: f64 = 3.69423907715893078616e-13; - const Lg1: f64 = 6.666666666666735130e-01; - const Lg2: f64 = 3.999999999940941908e-01; - const Lg3: f64 = 2.857142874366239149e-01; - const Lg4: f64 = 2.222219843214978396e-01; - const Lg5: f64 = 1.818357216161805012e-01; - const Lg6: f64 = 1.531383769920937332e-01; - const Lg7: f64 = 1.479819860511658591e-01; - - var x = x_; - var ix = @bitCast(u64, x); - var hx = @intCast(u32, ix >> 32); - var k: i32 = 0; - - if (hx < 0x00100000 or hx >> 31 != 0) { - // log(+-0) = -inf - if (ix << 1 == 0) { - return -math.inf(f32); - } - // log(-#) = nan - if (hx >> 31 != 0) { - return math.nan(f32); - } - - // subnormal, scale x - k -= 54; - x *= 0x1.0p54; - hx = @intCast(u32, @bitCast(u64, x) >> 32); - } else if (hx >= 0x7FF00000) { - return x; - } else if (hx == 0x3FF00000 and ix << 32 == 0) { - return 0; - } - - // x into [sqrt(2) / 2, sqrt(2)] - hx += 0x3FF00000 - 0x3FE6A09E; - k += @intCast(i32, hx >> 20) - 0x3FF; - hx = (hx & 0x000FFFFF) + 0x3FE6A09E; - ix = (@as(u64, hx) << 32) | (ix & 0xFFFFFFFF); - x = @bitCast(f64, ix); - - const f = x - 1.0; - const hfsq = 0.5 * f * f; - const s = f / (2.0 + f); - const z = s * s; - const w = z * z; - const t1 = w * (Lg2 + w * (Lg4 + w * Lg6)); - const t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7))); - const R = t2 + t1; - - // hi + lo = f - hfsq + s * (hfsq + R) ~ log(1 + f) - var hi = f - hfsq; - var hii = @bitCast(u64, hi); - hii &= @as(u64, maxInt(u64)) << 32; - hi = @bitCast(f64, hii); - const lo = f - hi - hfsq + s * (hfsq + R); - - // val_hi + val_lo ~ log10(1 + f) + k * log10(2) - var val_hi = hi * ivln10hi; - const dk = @intToFloat(f64, k); - const y = dk * log10_2hi; - var val_lo = dk * log10_2lo + (lo + hi) * ivln10lo + lo * ivln10hi; - - // Extra precision multiplication - const ww = y + val_hi; - val_lo += (y - ww) + val_hi; - val_hi = ww; - - return val_lo + val_hi; -} - -test "math.log10" { - try testing.expect(log10(@as(f32, 0.2)) == log10_32(0.2)); - try testing.expect(log10(@as(f64, 0.2)) == log10_64(0.2)); -} - -test "math.log10_32" { - const epsilon = 0.000001; - - try testing.expect(math.approxEqAbs(f32, log10_32(0.2), -0.698970, epsilon)); - try testing.expect(math.approxEqAbs(f32, log10_32(0.8923), -0.049489, epsilon)); - try testing.expect(math.approxEqAbs(f32, log10_32(1.5), 0.176091, epsilon)); - try testing.expect(math.approxEqAbs(f32, log10_32(37.45), 1.573452, epsilon)); - try testing.expect(math.approxEqAbs(f32, log10_32(89.123), 1.94999, epsilon)); - try testing.expect(math.approxEqAbs(f32, log10_32(123123.234375), 5.09034, epsilon)); -} - -test "math.log10_64" { - const epsilon = 0.000001; - - try testing.expect(math.approxEqAbs(f64, log10_64(0.2), -0.698970, epsilon)); - try testing.expect(math.approxEqAbs(f64, log10_64(0.8923), -0.049489, epsilon)); - try testing.expect(math.approxEqAbs(f64, log10_64(1.5), 0.176091, epsilon)); - try testing.expect(math.approxEqAbs(f64, log10_64(37.45), 1.573452, epsilon)); - try testing.expect(math.approxEqAbs(f64, log10_64(89.123), 1.94999, epsilon)); - try testing.expect(math.approxEqAbs(f64, log10_64(123123.234375), 5.09034, epsilon)); -} - -test "math.log10_32.special" { - try testing.expect(math.isPositiveInf(log10_32(math.inf(f32)))); - try testing.expect(math.isNegativeInf(log10_32(0.0))); - try testing.expect(math.isNan(log10_32(-1.0))); - try testing.expect(math.isNan(log10_32(math.nan(f32)))); -} - -test "math.log10_64.special" { - try testing.expect(math.isPositiveInf(log10_64(math.inf(f64)))); - try testing.expect(math.isNegativeInf(log10_64(0.0))); - try testing.expect(math.isNan(log10_64(-1.0))); - try testing.expect(math.isNan(log10_64(math.nan(f64)))); -} diff --git a/lib/std/math/log2.zig b/lib/std/math/log2.zig index 556c16f5cf..c83b170208 100644 --- a/lib/std/math/log2.zig +++ b/lib/std/math/log2.zig @@ -1,13 +1,6 @@ -// Ported from musl, which is licensed under the MIT license: -// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT -// -// https://git.musl-libc.org/cgit/musl/tree/src/math/log2f.c -// https://git.musl-libc.org/cgit/musl/tree/src/math/log2.c - const std = @import("../std.zig"); const math = std.math; const expect = std.testing.expect; -const maxInt = std.math.maxInt; /// Returns the base-2 logarithm of x. /// @@ -20,15 +13,9 @@ pub fn log2(x: anytype) @TypeOf(x) { const T = @TypeOf(x); switch (@typeInfo(T)) { .ComptimeFloat => { - return @as(comptime_float, log2_64(x)); - }, - .Float => { - return switch (T) { - f32 => log2_32(x), - f64 => log2_64(x), - else => @compileError("log2 not implemented for " ++ @typeName(T)), - }; + return @as(comptime_float, @log2(x)); }, + .Float => return @log2(x), .ComptimeInt => comptime { var result = 0; var x_shifted = x; @@ -46,168 +33,7 @@ pub fn log2(x: anytype) @TypeOf(x) { } } -pub fn log2_32(x_: f32) f32 { - const ivln2hi: f32 = 1.4428710938e+00; - const ivln2lo: f32 = -1.7605285393e-04; - const Lg1: f32 = 0xaaaaaa.0p-24; - const Lg2: f32 = 0xccce13.0p-25; - const Lg3: f32 = 0x91e9ee.0p-25; - const Lg4: f32 = 0xf89e26.0p-26; - - var x = x_; - var u = @bitCast(u32, x); - var ix = u; - var k: i32 = 0; - - // x < 2^(-126) - if (ix < 0x00800000 or ix >> 31 != 0) { - // log(+-0) = -inf - if (ix << 1 == 0) { - return -math.inf(f32); - } - // log(-#) = nan - if (ix >> 31 != 0) { - return math.nan(f32); - } - - k -= 25; - x *= 0x1.0p25; - ix = @bitCast(u32, x); - } else if (ix >= 0x7F800000) { - return x; - } else if (ix == 0x3F800000) { - return 0; - } - - // x into [sqrt(2) / 2, sqrt(2)] - ix += 0x3F800000 - 0x3F3504F3; - k += @intCast(i32, ix >> 23) - 0x7F; - ix = (ix & 0x007FFFFF) + 0x3F3504F3; - x = @bitCast(f32, ix); - - const f = x - 1.0; - const s = f / (2.0 + f); - const z = s * s; - const w = z * z; - const t1 = w * (Lg2 + w * Lg4); - const t2 = z * (Lg1 + w * Lg3); - const R = t2 + t1; - const hfsq = 0.5 * f * f; - - var hi = f - hfsq; - u = @bitCast(u32, hi); - u &= 0xFFFFF000; - hi = @bitCast(f32, u); - const lo = f - hi - hfsq + s * (hfsq + R); - return (lo + hi) * ivln2lo + lo * ivln2hi + hi * ivln2hi + @intToFloat(f32, k); -} - -pub fn log2_64(x_: f64) f64 { - const ivln2hi: f64 = 1.44269504072144627571e+00; - const ivln2lo: f64 = 1.67517131648865118353e-10; - const Lg1: f64 = 6.666666666666735130e-01; - const Lg2: f64 = 3.999999999940941908e-01; - const Lg3: f64 = 2.857142874366239149e-01; - const Lg4: f64 = 2.222219843214978396e-01; - const Lg5: f64 = 1.818357216161805012e-01; - const Lg6: f64 = 1.531383769920937332e-01; - const Lg7: f64 = 1.479819860511658591e-01; - - var x = x_; - var ix = @bitCast(u64, x); - var hx = @intCast(u32, ix >> 32); - var k: i32 = 0; - - if (hx < 0x00100000 or hx >> 31 != 0) { - // log(+-0) = -inf - if (ix << 1 == 0) { - return -math.inf(f64); - } - // log(-#) = nan - if (hx >> 31 != 0) { - return math.nan(f64); - } - - // subnormal, scale x - k -= 54; - x *= 0x1.0p54; - hx = @intCast(u32, @bitCast(u64, x) >> 32); - } else if (hx >= 0x7FF00000) { - return x; - } else if (hx == 0x3FF00000 and ix << 32 == 0) { - return 0; - } - - // x into [sqrt(2) / 2, sqrt(2)] - hx += 0x3FF00000 - 0x3FE6A09E; - k += @intCast(i32, hx >> 20) - 0x3FF; - hx = (hx & 0x000FFFFF) + 0x3FE6A09E; - ix = (@as(u64, hx) << 32) | (ix & 0xFFFFFFFF); - x = @bitCast(f64, ix); - - const f = x - 1.0; - const hfsq = 0.5 * f * f; - const s = f / (2.0 + f); - const z = s * s; - const w = z * z; - const t1 = w * (Lg2 + w * (Lg4 + w * Lg6)); - const t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7))); - const R = t2 + t1; - - // hi + lo = f - hfsq + s * (hfsq + R) ~ log(1 + f) - var hi = f - hfsq; - var hii = @bitCast(u64, hi); - hii &= @as(u64, maxInt(u64)) << 32; - hi = @bitCast(f64, hii); - const lo = f - hi - hfsq + s * (hfsq + R); - - var val_hi = hi * ivln2hi; - var val_lo = (lo + hi) * ivln2lo + lo * ivln2hi; - - // spadd(val_hi, val_lo, y) - const y = @intToFloat(f64, k); - const ww = y + val_hi; - val_lo += (y - ww) + val_hi; - val_hi = ww; - - return val_lo + val_hi; -} - -test "math.log2" { - try expect(log2(@as(f32, 0.2)) == log2_32(0.2)); - try expect(log2(@as(f64, 0.2)) == log2_64(0.2)); -} - -test "math.log2_32" { - const epsilon = 0.000001; - - try expect(math.approxEqAbs(f32, log2_32(0.2), -2.321928, epsilon)); - try expect(math.approxEqAbs(f32, log2_32(0.8923), -0.164399, epsilon)); - try expect(math.approxEqAbs(f32, log2_32(1.5), 0.584962, epsilon)); - try expect(math.approxEqAbs(f32, log2_32(37.45), 5.226894, epsilon)); - try expect(math.approxEqAbs(f32, log2_32(123123.234375), 16.909744, epsilon)); -} - -test "math.log2_64" { - const epsilon = 0.000001; - - try expect(math.approxEqAbs(f64, log2_64(0.2), -2.321928, epsilon)); - try expect(math.approxEqAbs(f64, log2_64(0.8923), -0.164399, epsilon)); - try expect(math.approxEqAbs(f64, log2_64(1.5), 0.584962, epsilon)); - try expect(math.approxEqAbs(f64, log2_64(37.45), 5.226894, epsilon)); - try expect(math.approxEqAbs(f64, log2_64(123123.234375), 16.909744, epsilon)); -} - -test "math.log2_32.special" { - try expect(math.isPositiveInf(log2_32(math.inf(f32)))); - try expect(math.isNegativeInf(log2_32(0.0))); - try expect(math.isNan(log2_32(-1.0))); - try expect(math.isNan(log2_32(math.nan(f32)))); -} - -test "math.log2_64.special" { - try expect(math.isPositiveInf(log2_64(math.inf(f64)))); - try expect(math.isNegativeInf(log2_64(0.0))); - try expect(math.isNan(log2_64(-1.0))); - try expect(math.isNan(log2_64(math.nan(f64)))); +test "log2" { + try expect(log2(@as(f32, 0.2)) == @log2(0.2)); + try expect(log2(@as(f64, 0.2)) == @log2(0.2)); } diff --git a/lib/std/math/nan.zig b/lib/std/math/nan.zig index 634af1f0d6..329f67b74e 100644 --- a/lib/std/math/nan.zig +++ b/lib/std/math/nan.zig @@ -2,13 +2,13 @@ const math = @import("../math.zig"); /// Returns the nan representation for type T. pub fn nan(comptime T: type) T { - return switch (T) { - f16 => math.nan_f16, - f32 => math.nan_f32, - f64 => math.nan_f64, - f80 => math.nan_f80, - f128 => math.nan_f128, - else => @compileError("nan not implemented for " ++ @typeName(T)), + return switch (@typeInfo(T).Float.bits) { + 16 => math.nan_f16, + 32 => math.nan_f32, + 64 => math.nan_f64, + 80 => math.nan_f80, + 128 => math.nan_f128, + else => @compileError("unreachable"), }; } @@ -16,12 +16,12 @@ pub fn nan(comptime T: type) T { pub fn snan(comptime T: type) T { // Note: A signalling nan is identical to a standard right now by may have a different bit // representation in the future when required. - return switch (T) { - f16 => @bitCast(f16, math.nan_u16), - f32 => @bitCast(f32, math.nan_u32), - f64 => @bitCast(f64, math.nan_u64), - f80 => @bitCast(f80, math.nan_u80), - f128 => @bitCast(f128, math.nan_u128), - else => @compileError("snan not implemented for " ++ @typeName(T)), + return switch (@typeInfo(T).Float.bits) { + 16 => math.nan_u16, + 32 => math.nan_u32, + 64 => math.nan_u64, + 80 => math.nan_u80, + 128 => math.nan_u128, + else => @compileError("unreachable"), }; } diff --git a/lib/std/math/pow.zig b/lib/std/math/pow.zig index 040abf9a44..48c6636926 100644 --- a/lib/std/math/pow.zig +++ b/lib/std/math/pow.zig @@ -82,7 +82,7 @@ pub fn pow(comptime T: type, x: T, y: T) T { } // pow(x, +inf) = +0 for |x| < 1 // pow(x, -inf) = +0 for |x| > 1 - else if ((math.fabs(x) < 1) == math.isPositiveInf(y)) { + else if ((@fabs(x) < 1) == math.isPositiveInf(y)) { return 0; } // pow(x, -inf) = +inf for |x| < 1 @@ -108,14 +108,14 @@ pub fn pow(comptime T: type, x: T, y: T) T { // special case sqrt if (y == 0.5) { - return math.sqrt(x); + return @sqrt(x); } if (y == -0.5) { - return 1 / math.sqrt(x); + return 1 / @sqrt(x); } - const r1 = math.modf(math.fabs(y)); + const r1 = math.modf(@fabs(y)); var yi = r1.ipart; var yf = r1.fpart; @@ -123,7 +123,7 @@ pub fn pow(comptime T: type, x: T, y: T) T { return math.nan(T); } if (yi >= 1 << (@typeInfo(T).Float.bits - 1)) { - return math.exp(y * math.ln(x)); + return @exp(y * @log(x)); } // a = a1 * 2^ae @@ -136,7 +136,7 @@ pub fn pow(comptime T: type, x: T, y: T) T { yf -= 1; yi += 1; } - a1 = math.exp(yf * math.ln(x)); + a1 = @exp(yf * @log(x)); } // a *= x^yi diff --git a/lib/std/math/round.zig b/lib/std/math/round.zig deleted file mode 100644 index be33a9cfbd..0000000000 --- a/lib/std/math/round.zig +++ /dev/null @@ -1,185 +0,0 @@ -// Ported from musl, which is licensed under the MIT license: -// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT -// -// https://git.musl-libc.org/cgit/musl/tree/src/math/roundf.c -// https://git.musl-libc.org/cgit/musl/tree/src/math/round.c - -const expect = std.testing.expect; -const std = @import("../std.zig"); -const math = std.math; - -/// Returns x rounded to the nearest integer, rounding half away from zero. -/// -/// Special Cases: -/// - round(+-0) = +-0 -/// - round(+-inf) = +-inf -/// - round(nan) = nan -pub fn round(x: anytype) @TypeOf(x) { - const T = @TypeOf(x); - return switch (T) { - f32 => round32(x), - f64 => round64(x), - f128 => round128(x), - - // TODO this is not correct for some targets - c_longdouble => @floatCast(c_longdouble, round128(x)), - - else => @compileError("round not implemented for " ++ @typeName(T)), - }; -} - -fn round32(x_: f32) f32 { - const f32_toint = 1.0 / math.floatEps(f32); - - var x = x_; - const u = @bitCast(u32, x); - const e = (u >> 23) & 0xFF; - var y: f32 = undefined; - - if (e >= 0x7F + 23) { - return x; - } - if (u >> 31 != 0) { - x = -x; - } - if (e < 0x7F - 1) { - math.doNotOptimizeAway(x + f32_toint); - return 0 * @bitCast(f32, u); - } - - y = x + f32_toint - f32_toint - x; - if (y > 0.5) { - y = y + x - 1; - } else if (y <= -0.5) { - y = y + x + 1; - } else { - y = y + x; - } - - if (u >> 31 != 0) { - return -y; - } else { - return y; - } -} - -fn round64(x_: f64) f64 { - const f64_toint = 1.0 / math.floatEps(f64); - - var x = x_; - const u = @bitCast(u64, x); - const e = (u >> 52) & 0x7FF; - var y: f64 = undefined; - - if (e >= 0x3FF + 52) { - return x; - } - if (u >> 63 != 0) { - x = -x; - } - if (e < 0x3ff - 1) { - math.doNotOptimizeAway(x + f64_toint); - return 0 * @bitCast(f64, u); - } - - y = x + f64_toint - f64_toint - x; - if (y > 0.5) { - y = y + x - 1; - } else if (y <= -0.5) { - y = y + x + 1; - } else { - y = y + x; - } - - if (u >> 63 != 0) { - return -y; - } else { - return y; - } -} - -fn round128(x_: f128) f128 { - const f128_toint = 1.0 / math.floatEps(f128); - - var x = x_; - const u = @bitCast(u128, x); - const e = (u >> 112) & 0x7FFF; - var y: f128 = undefined; - - if (e >= 0x3FFF + 112) { - return x; - } - if (u >> 127 != 0) { - x = -x; - } - if (e < 0x3FFF - 1) { - math.doNotOptimizeAway(x + f128_toint); - return 0 * @bitCast(f128, u); - } - - y = x + f128_toint - f128_toint - x; - if (y > 0.5) { - y = y + x - 1; - } else if (y <= -0.5) { - y = y + x + 1; - } else { - y = y + x; - } - - if (u >> 127 != 0) { - return -y; - } else { - return y; - } -} - -test "math.round" { - try expect(round(@as(f32, 1.3)) == round32(1.3)); - try expect(round(@as(f64, 1.3)) == round64(1.3)); - try expect(round(@as(f128, 1.3)) == round128(1.3)); -} - -test "math.round32" { - try expect(round32(1.3) == 1.0); - try expect(round32(-1.3) == -1.0); - try expect(round32(0.2) == 0.0); - try expect(round32(1.8) == 2.0); -} - -test "math.round64" { - try expect(round64(1.3) == 1.0); - try expect(round64(-1.3) == -1.0); - try expect(round64(0.2) == 0.0); - try expect(round64(1.8) == 2.0); -} - -test "math.round128" { - try expect(round128(1.3) == 1.0); - try expect(round128(-1.3) == -1.0); - try expect(round128(0.2) == 0.0); - try expect(round128(1.8) == 2.0); -} - -test "math.round32.special" { - try expect(round32(0.0) == 0.0); - try expect(round32(-0.0) == -0.0); - try expect(math.isPositiveInf(round32(math.inf(f32)))); - try expect(math.isNegativeInf(round32(-math.inf(f32)))); - try expect(math.isNan(round32(math.nan(f32)))); -} - -test "math.round64.special" { - try expect(round64(0.0) == 0.0); - try expect(round64(-0.0) == -0.0); - try expect(math.isPositiveInf(round64(math.inf(f64)))); - try expect(math.isNegativeInf(round64(-math.inf(f64)))); - try expect(math.isNan(round64(math.nan(f64)))); -} - -test "math.round128.special" { - try expect(round128(0.0) == 0.0); - try expect(round128(-0.0) == -0.0); - try expect(math.isPositiveInf(round128(math.inf(f128)))); - try expect(math.isNegativeInf(round128(-math.inf(f128)))); - try expect(math.isNan(round128(math.nan(f128)))); -} diff --git a/lib/std/math/trunc.zig b/lib/std/math/trunc.zig deleted file mode 100644 index 32bd7fb0aa..0000000000 --- a/lib/std/math/trunc.zig +++ /dev/null @@ -1,141 +0,0 @@ -// Ported from musl, which is licensed under the MIT license: -// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT -// -// https://git.musl-libc.org/cgit/musl/tree/src/math/truncf.c -// https://git.musl-libc.org/cgit/musl/tree/src/math/trunc.c - -const std = @import("../std.zig"); -const math = std.math; -const expect = std.testing.expect; -const maxInt = std.math.maxInt; - -/// Returns the integer value of x. -/// -/// Special Cases: -/// - trunc(+-0) = +-0 -/// - trunc(+-inf) = +-inf -/// - trunc(nan) = nan -pub fn trunc(x: anytype) @TypeOf(x) { - const T = @TypeOf(x); - return switch (T) { - f32 => trunc32(x), - f64 => trunc64(x), - f128 => trunc128(x), - - // TODO this is not correct for some targets - c_longdouble => @floatCast(c_longdouble, trunc128(x)), - - else => @compileError("trunc not implemented for " ++ @typeName(T)), - }; -} - -fn trunc32(x: f32) f32 { - const u = @bitCast(u32, x); - var e = @intCast(i32, ((u >> 23) & 0xFF)) - 0x7F + 9; - var m: u32 = undefined; - - if (e >= 23 + 9) { - return x; - } - if (e < 9) { - e = 1; - } - - m = @as(u32, maxInt(u32)) >> @intCast(u5, e); - if (u & m == 0) { - return x; - } else { - math.doNotOptimizeAway(x + 0x1p120); - return @bitCast(f32, u & ~m); - } -} - -fn trunc64(x: f64) f64 { - const u = @bitCast(u64, x); - var e = @intCast(i32, ((u >> 52) & 0x7FF)) - 0x3FF + 12; - var m: u64 = undefined; - - if (e >= 52 + 12) { - return x; - } - if (e < 12) { - e = 1; - } - - m = @as(u64, maxInt(u64)) >> @intCast(u6, e); - if (u & m == 0) { - return x; - } else { - math.doNotOptimizeAway(x + 0x1p120); - return @bitCast(f64, u & ~m); - } -} - -fn trunc128(x: f128) f128 { - const u = @bitCast(u128, x); - var e = @intCast(i32, ((u >> 112) & 0x7FFF)) - 0x3FFF + 16; - var m: u128 = undefined; - - if (e >= 112 + 16) { - return x; - } - if (e < 16) { - e = 1; - } - - m = @as(u128, maxInt(u128)) >> @intCast(u7, e); - if (u & m == 0) { - return x; - } else { - math.doNotOptimizeAway(x + 0x1p120); - return @bitCast(f128, u & ~m); - } -} - -test "math.trunc" { - try expect(trunc(@as(f32, 1.3)) == trunc32(1.3)); - try expect(trunc(@as(f64, 1.3)) == trunc64(1.3)); - try expect(trunc(@as(f128, 1.3)) == trunc128(1.3)); -} - -test "math.trunc32" { - try expect(trunc32(1.3) == 1.0); - try expect(trunc32(-1.3) == -1.0); - try expect(trunc32(0.2) == 0.0); -} - -test "math.trunc64" { - try expect(trunc64(1.3) == 1.0); - try expect(trunc64(-1.3) == -1.0); - try expect(trunc64(0.2) == 0.0); -} - -test "math.trunc128" { - try expect(trunc128(1.3) == 1.0); - try expect(trunc128(-1.3) == -1.0); - try expect(trunc128(0.2) == 0.0); -} - -test "math.trunc32.special" { - try expect(trunc32(0.0) == 0.0); // 0x3F800000 - try expect(trunc32(-0.0) == -0.0); - try expect(math.isPositiveInf(trunc32(math.inf(f32)))); - try expect(math.isNegativeInf(trunc32(-math.inf(f32)))); - try expect(math.isNan(trunc32(math.nan(f32)))); -} - -test "math.trunc64.special" { - try expect(trunc64(0.0) == 0.0); - try expect(trunc64(-0.0) == -0.0); - try expect(math.isPositiveInf(trunc64(math.inf(f64)))); - try expect(math.isNegativeInf(trunc64(-math.inf(f64)))); - try expect(math.isNan(trunc64(math.nan(f64)))); -} - -test "math.trunc128.special" { - try expect(trunc128(0.0) == 0.0); - try expect(trunc128(-0.0) == -0.0); - try expect(math.isPositiveInf(trunc128(math.inf(f128)))); - try expect(math.isNegativeInf(trunc128(-math.inf(f128)))); - try expect(math.isNan(trunc128(math.nan(f128)))); -} diff --git a/lib/std/rand/ziggurat.zig b/lib/std/rand/ziggurat.zig index 5c18d0023b..b05ce7fd73 100644 --- a/lib/std/rand/ziggurat.zig +++ b/lib/std/rand/ziggurat.zig @@ -33,7 +33,7 @@ pub fn next_f64(random: Random, comptime tables: ZigTable) f64 { }; const x = u * tables.x[i]; - const test_x = if (tables.is_symmetric) math.fabs(x) else x; + const test_x = if (tables.is_symmetric) @fabs(x) else x; // equivalent to |u| < tables.x[i+1] / tables.x[i] (or u < tables.x[i+1] / tables.x[i]) if (test_x < tables.x[i + 1]) { @@ -106,18 +106,18 @@ const norm_r = 3.6541528853610088; const norm_v = 0.00492867323399; fn norm_f(x: f64) f64 { - return math.exp(-x * x / 2.0); + return @exp(-x * x / 2.0); } fn norm_f_inv(y: f64) f64 { - return math.sqrt(-2.0 * math.ln(y)); + return @sqrt(-2.0 * @log(y)); } fn norm_zero_case(random: Random, u: f64) f64 { var x: f64 = 1; var y: f64 = 0; while (-2.0 * y < x * x) { - x = math.ln(random.float(f64)) / norm_r; - y = math.ln(random.float(f64)); + x = @log(random.float(f64)) / norm_r; + y = @log(random.float(f64)); } if (u < 0) { @@ -151,13 +151,13 @@ const exp_r = 7.69711747013104972; const exp_v = 0.0039496598225815571993; fn exp_f(x: f64) f64 { - return math.exp(-x); + return @exp(-x); } fn exp_f_inv(y: f64) f64 { - return -math.ln(y); + return -@log(y); } fn exp_zero_case(random: Random, _: f64) f64 { - return exp_r - math.ln(random.float(f64)); + return exp_r - @log(random.float(f64)); } test "exp dist sanity" { diff --git a/lib/std/special/c.zig b/lib/std/special/c.zig index dfc2020334..525bdd267d 100644 --- a/lib/std/special/c.zig +++ b/lib/std/special/c.zig @@ -12,7 +12,6 @@ const maxInt = std.math.maxInt; const native_os = builtin.os.tag; const native_arch = builtin.cpu.arch; const native_abi = builtin.abi; -const long_double_is_f128 = builtin.target.longDoubleIs(f128); const is_wasm = switch (native_arch) { .wasm32, .wasm64 => true, @@ -55,53 +54,6 @@ comptime { } else if (is_msvc) { @export(_fltused, .{ .name = "_fltused", .linkage = .Strong }); } - - @export(trunc, .{ .name = "trunc", .linkage = .Strong }); - @export(truncf, .{ .name = "truncf", .linkage = .Strong }); - @export(truncl, .{ .name = "truncl", .linkage = .Strong }); - - @export(log, .{ .name = "log", .linkage = .Strong }); - @export(logf, .{ .name = "logf", .linkage = .Strong }); - - @export(sin, .{ .name = "sin", .linkage = .Strong }); - @export(sinf, .{ .name = "sinf", .linkage = .Strong }); - - @export(cos, .{ .name = "cos", .linkage = .Strong }); - @export(cosf, .{ .name = "cosf", .linkage = .Strong }); - - @export(exp, .{ .name = "exp", .linkage = .Strong }); - @export(expf, .{ .name = "expf", .linkage = .Strong }); - - @export(exp2, .{ .name = "exp2", .linkage = .Strong }); - @export(exp2f, .{ .name = "exp2f", .linkage = .Strong }); - - @export(log2, .{ .name = "log2", .linkage = .Strong }); - @export(log2f, .{ .name = "log2f", .linkage = .Strong }); - - @export(log10, .{ .name = "log10", .linkage = .Strong }); - @export(log10f, .{ .name = "log10f", .linkage = .Strong }); - - @export(fmod, .{ .name = "fmod", .linkage = .Strong }); - @export(fmodf, .{ .name = "fmodf", .linkage = .Strong }); - - @export(sincos, .{ .name = "sincos", .linkage = .Strong }); - @export(sincosf, .{ .name = "sincosf", .linkage = .Strong }); - - @export(fabs, .{ .name = "fabs", .linkage = .Strong }); - @export(fabsf, .{ .name = "fabsf", .linkage = .Strong }); - - @export(round, .{ .name = "round", .linkage = .Strong }); - @export(roundf, .{ .name = "roundf", .linkage = .Strong }); - @export(roundl, .{ .name = "roundl", .linkage = .Strong }); - - @export(fmin, .{ .name = "fmin", .linkage = .Strong }); - @export(fminf, .{ .name = "fminf", .linkage = .Strong }); - - @export(fmax, .{ .name = "fmax", .linkage = .Strong }); - @export(fmaxf, .{ .name = "fmaxf", .linkage = .Strong }); - - @export(sqrt, .{ .name = "sqrt", .linkage = .Strong }); - @export(sqrtf, .{ .name = "sqrtf", .linkage = .Strong }); } // Avoid dragging in the runtime safety mechanisms into this .o file, @@ -352,549 +304,6 @@ test "strncmp" { try std.testing.expect(strncmp("\xff", "\x02", 1) == 253); } -fn trunc(a: f64) callconv(.C) f64 { - return math.trunc(a); -} - -fn truncf(a: f32) callconv(.C) f32 { - return math.trunc(a); -} - -fn truncl(a: c_longdouble) callconv(.C) c_longdouble { - if (!long_double_is_f128) { - @panic("TODO implement this"); - } - return math.trunc(a); -} - -fn log(a: f64) callconv(.C) f64 { - return math.ln(a); -} - -fn logf(a: f32) callconv(.C) f32 { - return math.ln(a); -} - -fn sin(a: f64) callconv(.C) f64 { - return math.sin(a); -} - -fn sinf(a: f32) callconv(.C) f32 { - return math.sin(a); -} - -fn cos(a: f64) callconv(.C) f64 { - return math.cos(a); -} - -fn cosf(a: f32) callconv(.C) f32 { - return math.cos(a); -} - -fn exp(a: f64) callconv(.C) f64 { - return math.exp(a); -} - -fn expf(a: f32) callconv(.C) f32 { - return math.exp(a); -} - -fn exp2(a: f64) callconv(.C) f64 { - return math.exp2(a); -} - -fn exp2f(a: f32) callconv(.C) f32 { - return math.exp2(a); -} - -fn log2(a: f64) callconv(.C) f64 { - return math.log2(a); -} - -fn log2f(a: f32) callconv(.C) f32 { - return math.log2(a); -} - -fn log10(a: f64) callconv(.C) f64 { - return math.log10(a); -} - -fn log10f(a: f32) callconv(.C) f32 { - return math.log10(a); -} - -fn fmodf(x: f32, y: f32) callconv(.C) f32 { - return generic_fmod(f32, x, y); -} -fn fmod(x: f64, y: f64) callconv(.C) f64 { - return generic_fmod(f64, x, y); -} - -fn generic_fmod(comptime T: type, x: T, y: T) T { - @setRuntimeSafety(false); - - const bits = @typeInfo(T).Float.bits; - const uint = std.meta.Int(.unsigned, bits); - const log2uint = math.Log2Int(uint); - const digits = if (T == f32) 23 else 52; - const exp_bits = if (T == f32) 9 else 12; - const bits_minus_1 = bits - 1; - const mask = if (T == f32) 0xff else 0x7ff; - var ux = @bitCast(uint, x); - var uy = @bitCast(uint, y); - var ex = @intCast(i32, (ux >> digits) & mask); - var ey = @intCast(i32, (uy >> digits) & mask); - const sx = if (T == f32) @intCast(u32, ux & 0x80000000) else @intCast(i32, ux >> bits_minus_1); - var i: uint = undefined; - - if (uy << 1 == 0 or isNan(@bitCast(T, uy)) or ex == mask) - return (x * y) / (x * y); - - if (ux << 1 <= uy << 1) { - if (ux << 1 == uy << 1) - return 0 * x; - return x; - } - - // normalize x and y - if (ex == 0) { - i = ux << exp_bits; - while (i >> bits_minus_1 == 0) : ({ - ex -= 1; - i <<= 1; - }) {} - ux <<= @intCast(log2uint, @bitCast(u32, -ex + 1)); - } else { - ux &= maxInt(uint) >> exp_bits; - ux |= 1 << digits; - } - if (ey == 0) { - i = uy << exp_bits; - while (i >> bits_minus_1 == 0) : ({ - ey -= 1; - i <<= 1; - }) {} - uy <<= @intCast(log2uint, @bitCast(u32, -ey + 1)); - } else { - uy &= maxInt(uint) >> exp_bits; - uy |= 1 << digits; - } - - // x mod y - while (ex > ey) : (ex -= 1) { - i = ux -% uy; - if (i >> bits_minus_1 == 0) { - if (i == 0) - return 0 * x; - ux = i; - } - ux <<= 1; - } - i = ux -% uy; - if (i >> bits_minus_1 == 0) { - if (i == 0) - return 0 * x; - ux = i; - } - while (ux >> digits == 0) : ({ - ux <<= 1; - ex -= 1; - }) {} - - // scale result up - if (ex > 0) { - ux -%= 1 << digits; - ux |= @as(uint, @bitCast(u32, ex)) << digits; - } else { - ux >>= @intCast(log2uint, @bitCast(u32, -ex + 1)); - } - if (T == f32) { - ux |= sx; - } else { - ux |= @intCast(uint, sx) << bits_minus_1; - } - return @bitCast(T, ux); -} - -test "fmod, fmodf" { - inline for ([_]type{ f32, f64 }) |T| { - const nan_val = math.nan(T); - const inf_val = math.inf(T); - - try std.testing.expect(isNan(generic_fmod(T, nan_val, 1.0))); - try std.testing.expect(isNan(generic_fmod(T, 1.0, nan_val))); - try std.testing.expect(isNan(generic_fmod(T, inf_val, 1.0))); - try std.testing.expect(isNan(generic_fmod(T, 0.0, 0.0))); - try std.testing.expect(isNan(generic_fmod(T, 1.0, 0.0))); - - try std.testing.expectEqual(@as(T, 0.0), generic_fmod(T, 0.0, 2.0)); - try std.testing.expectEqual(@as(T, -0.0), generic_fmod(T, -0.0, 2.0)); - - try std.testing.expectEqual(@as(T, -2.0), generic_fmod(T, -32.0, 10.0)); - try std.testing.expectEqual(@as(T, -2.0), generic_fmod(T, -32.0, -10.0)); - try std.testing.expectEqual(@as(T, 2.0), generic_fmod(T, 32.0, 10.0)); - try std.testing.expectEqual(@as(T, 2.0), generic_fmod(T, 32.0, -10.0)); - } -} - -fn sincos(a: f64, r_sin: *f64, r_cos: *f64) callconv(.C) void { - r_sin.* = math.sin(a); - r_cos.* = math.cos(a); -} - -fn sincosf(a: f32, r_sin: *f32, r_cos: *f32) callconv(.C) void { - r_sin.* = math.sin(a); - r_cos.* = math.cos(a); -} - -fn fabs(a: f64) callconv(.C) f64 { - return math.fabs(a); -} - -fn fabsf(a: f32) callconv(.C) f32 { - return math.fabs(a); -} - -fn roundf(a: f32) callconv(.C) f32 { - return math.round(a); -} - -fn round(a: f64) callconv(.C) f64 { - return math.round(a); -} - -fn roundl(a: c_longdouble) callconv(.C) c_longdouble { - if (!long_double_is_f128) { - @panic("TODO implement this"); - } - return math.round(a); -} - -fn fminf(x: f32, y: f32) callconv(.C) f32 { - return generic_fmin(f32, x, y); -} - -fn fmin(x: f64, y: f64) callconv(.C) f64 { - return generic_fmin(f64, x, y); -} - -fn generic_fmin(comptime T: type, x: T, y: T) T { - if (isNan(x)) - return y; - if (isNan(y)) - return x; - return if (x < y) x else y; -} - -test "fmin, fminf" { - inline for ([_]type{ f32, f64 }) |T| { - const nan_val = math.nan(T); - - try std.testing.expect(isNan(generic_fmin(T, nan_val, nan_val))); - try std.testing.expectEqual(@as(T, 1.0), generic_fmin(T, nan_val, 1.0)); - try std.testing.expectEqual(@as(T, 1.0), generic_fmin(T, 1.0, nan_val)); - - try std.testing.expectEqual(@as(T, 1.0), generic_fmin(T, 1.0, 10.0)); - try std.testing.expectEqual(@as(T, -1.0), generic_fmin(T, 1.0, -1.0)); - } -} - -fn fmaxf(x: f32, y: f32) callconv(.C) f32 { - return generic_fmax(f32, x, y); -} - -fn fmax(x: f64, y: f64) callconv(.C) f64 { - return generic_fmax(f64, x, y); -} - -fn generic_fmax(comptime T: type, x: T, y: T) T { - if (isNan(x)) - return y; - if (isNan(y)) - return x; - return if (x < y) y else x; -} - -test "fmax, fmaxf" { - inline for ([_]type{ f32, f64 }) |T| { - const nan_val = math.nan(T); - - try std.testing.expect(isNan(generic_fmax(T, nan_val, nan_val))); - try std.testing.expectEqual(@as(T, 1.0), generic_fmax(T, nan_val, 1.0)); - try std.testing.expectEqual(@as(T, 1.0), generic_fmax(T, 1.0, nan_val)); - - try std.testing.expectEqual(@as(T, 10.0), generic_fmax(T, 1.0, 10.0)); - try std.testing.expectEqual(@as(T, 1.0), generic_fmax(T, 1.0, -1.0)); - } -} - -// NOTE: The original code is full of implicit signed -> unsigned assumptions and u32 wraparound -// behaviour. Most intermediate i32 values are changed to u32 where appropriate but there are -// potentially some edge cases remaining that are not handled in the same way. -fn sqrt(x: f64) callconv(.C) f64 { - const tiny: f64 = 1.0e-300; - const sign: u32 = 0x80000000; - const u = @bitCast(u64, x); - - var ix0 = @intCast(u32, u >> 32); - var ix1 = @intCast(u32, u & 0xFFFFFFFF); - - // sqrt(nan) = nan, sqrt(+inf) = +inf, sqrt(-inf) = nan - if (ix0 & 0x7FF00000 == 0x7FF00000) { - return x * x + x; - } - - // sqrt(+-0) = +-0 - if (x == 0.0) { - return x; - } - // sqrt(-ve) = snan - if (ix0 & sign != 0) { - return math.snan(f64); - } - - // normalize x - var m = @intCast(i32, ix0 >> 20); - if (m == 0) { - // subnormal - while (ix0 == 0) { - m -= 21; - ix0 |= ix1 >> 11; - ix1 <<= 21; - } - - // subnormal - var i: u32 = 0; - while (ix0 & 0x00100000 == 0) : (i += 1) { - ix0 <<= 1; - } - m -= @intCast(i32, i) - 1; - ix0 |= ix1 >> @intCast(u5, 32 - i); - ix1 <<= @intCast(u5, i); - } - - // unbias exponent - m -= 1023; - ix0 = (ix0 & 0x000FFFFF) | 0x00100000; - if (m & 1 != 0) { - ix0 += ix0 + (ix1 >> 31); - ix1 = ix1 +% ix1; - } - m >>= 1; - - // sqrt(x) bit by bit - ix0 += ix0 + (ix1 >> 31); - ix1 = ix1 +% ix1; - - var q: u32 = 0; - var q1: u32 = 0; - var s0: u32 = 0; - var s1: u32 = 0; - var r: u32 = 0x00200000; - var t: u32 = undefined; - var t1: u32 = undefined; - - while (r != 0) { - t = s0 +% r; - if (t <= ix0) { - s0 = t + r; - ix0 -= t; - q += r; - } - ix0 = ix0 +% ix0 +% (ix1 >> 31); - ix1 = ix1 +% ix1; - r >>= 1; - } - - r = sign; - while (r != 0) { - t1 = s1 +% r; - t = s0; - if (t < ix0 or (t == ix0 and t1 <= ix1)) { - s1 = t1 +% r; - if (t1 & sign == sign and s1 & sign == 0) { - s0 += 1; - } - ix0 -= t; - if (ix1 < t1) { - ix0 -= 1; - } - ix1 = ix1 -% t1; - q1 += r; - } - ix0 = ix0 +% ix0 +% (ix1 >> 31); - ix1 = ix1 +% ix1; - r >>= 1; - } - - // rounding direction - if (ix0 | ix1 != 0) { - var z = 1.0 - tiny; // raise inexact - if (z >= 1.0) { - z = 1.0 + tiny; - if (q1 == 0xFFFFFFFF) { - q1 = 0; - q += 1; - } else if (z > 1.0) { - if (q1 == 0xFFFFFFFE) { - q += 1; - } - q1 += 2; - } else { - q1 += q1 & 1; - } - } - } - - ix0 = (q >> 1) + 0x3FE00000; - ix1 = q1 >> 1; - if (q & 1 != 0) { - ix1 |= 0x80000000; - } - - // NOTE: musl here appears to rely on signed twos-complement wraparound. +% has the same - // behaviour at least. - var iix0 = @intCast(i32, ix0); - iix0 = iix0 +% (m << 20); - - const uz = (@intCast(u64, iix0) << 32) | ix1; - return @bitCast(f64, uz); -} - -test "sqrt" { - const V = [_]f64{ - 0.0, - 4.089288054930154, - 7.538757127071935, - 8.97780793672623, - 5.304443821913729, - 5.682408965311888, - 0.5846878579110049, - 3.650338664297043, - 0.3178091951800732, - 7.1505232436382835, - 3.6589165881946464, - }; - - // Note that @sqrt will either generate the sqrt opcode (if supported by the - // target ISA) or a call to `sqrtf` otherwise. - for (V) |val| - try std.testing.expectEqual(@sqrt(val), sqrt(val)); -} - -test "sqrt special" { - try std.testing.expect(std.math.isPositiveInf(sqrt(std.math.inf(f64)))); - try std.testing.expect(sqrt(0.0) == 0.0); - try std.testing.expect(sqrt(-0.0) == -0.0); - try std.testing.expect(isNan(sqrt(-1.0))); - try std.testing.expect(isNan(sqrt(std.math.nan(f64)))); -} - -fn sqrtf(x: f32) callconv(.C) f32 { - const tiny: f32 = 1.0e-30; - const sign: i32 = @bitCast(i32, @as(u32, 0x80000000)); - var ix: i32 = @bitCast(i32, x); - - if ((ix & 0x7F800000) == 0x7F800000) { - return x * x + x; // sqrt(nan) = nan, sqrt(+inf) = +inf, sqrt(-inf) = snan - } - - // zero - if (ix <= 0) { - if (ix & ~sign == 0) { - return x; // sqrt (+-0) = +-0 - } - if (ix < 0) { - return math.snan(f32); - } - } - - // normalize - var m = ix >> 23; - if (m == 0) { - // subnormal - var i: i32 = 0; - while (ix & 0x00800000 == 0) : (i += 1) { - ix <<= 1; - } - m -= i - 1; - } - - m -= 127; // unbias exponent - ix = (ix & 0x007FFFFF) | 0x00800000; - - if (m & 1 != 0) { // odd m, double x to even - ix += ix; - } - - m >>= 1; // m = [m / 2] - - // sqrt(x) bit by bit - ix += ix; - var q: i32 = 0; // q = sqrt(x) - var s: i32 = 0; - var r: i32 = 0x01000000; // r = moving bit right -> left - - while (r != 0) { - const t = s + r; - if (t <= ix) { - s = t + r; - ix -= t; - q += r; - } - ix += ix; - r >>= 1; - } - - // floating add to find rounding direction - if (ix != 0) { - var z = 1.0 - tiny; // inexact - if (z >= 1.0) { - z = 1.0 + tiny; - if (z > 1.0) { - q += 2; - } else { - if (q & 1 != 0) { - q += 1; - } - } - } - } - - ix = (q >> 1) + 0x3f000000; - ix += m << 23; - return @bitCast(f32, ix); -} - -test "sqrtf" { - const V = [_]f32{ - 0.0, - 4.089288054930154, - 7.538757127071935, - 8.97780793672623, - 5.304443821913729, - 5.682408965311888, - 0.5846878579110049, - 3.650338664297043, - 0.3178091951800732, - 7.1505232436382835, - 3.6589165881946464, - }; - - // Note that @sqrt will either generate the sqrt opcode (if supported by the - // target ISA) or a call to `sqrtf` otherwise. - for (V) |val| - try std.testing.expectEqual(@sqrt(val), sqrtf(val)); -} - -test "sqrtf special" { - try std.testing.expect(std.math.isPositiveInf(sqrtf(std.math.inf(f32)))); - try std.testing.expect(sqrtf(0.0) == 0.0); - try std.testing.expect(sqrtf(-0.0) == -0.0); - try std.testing.expect(isNan(sqrtf(-1.0))); - try std.testing.expect(isNan(sqrtf(std.math.nan(f32)))); -} - // TODO we should be able to put this directly in std/linux/x86_64.zig but // it causes a segfault in release mode. this is a workaround of calling it // across .o file boundaries. fix comptime @ptrCast of nakedcc functions. diff --git a/lib/std/special/compiler_rt.zig b/lib/std/special/compiler_rt.zig index 93e0ffbe1a..dccb9264bd 100644 --- a/lib/std/special/compiler_rt.zig +++ b/lib/std/special/compiler_rt.zig @@ -19,9 +19,6 @@ const strong_linkage = if (is_test) else std.builtin.GlobalLinkage.Strong; -const long_double_is_f80 = builtin.target.longDoubleIs(f80); -const long_double_is_f128 = builtin.target.longDoubleIs(f128); - comptime { // These files do their own comptime exporting logic. _ = @import("compiler_rt/atomics.zig"); @@ -726,42 +723,25 @@ comptime { @export(_aullrem, .{ .name = "\x01__aullrem", .linkage = strong_linkage }); } - if (!is_test) { - if (long_double_is_f80) { - @export(fmodx, .{ .name = "fmodl", .linkage = linkage }); - } else if (long_double_is_f128) { - @export(fmodq, .{ .name = "fmodl", .linkage = linkage }); - } else { - @export(fmodl, .{ .name = "fmodl", .linkage = linkage }); - } - if (long_double_is_f80 or builtin.zig_backend == .stage1) { - // TODO: https://github.com/ziglang/zig/issues/11161 - @export(fmodx, .{ .name = "fmodx", .linkage = linkage }); - } - @export(fmodq, .{ .name = "fmodq", .linkage = linkage }); - - @export(floorf, .{ .name = "floorf", .linkage = linkage }); - @export(floor, .{ .name = "floor", .linkage = linkage }); - @export(floorl, .{ .name = "floorl", .linkage = linkage }); - - @export(ceilf, .{ .name = "ceilf", .linkage = linkage }); - @export(ceil, .{ .name = "ceil", .linkage = linkage }); - @export(ceill, .{ .name = "ceill", .linkage = linkage }); - - @export(fma, .{ .name = "fma", .linkage = linkage }); - @export(fmaf, .{ .name = "fmaf", .linkage = linkage }); - @export(fmal, .{ .name = "fmal", .linkage = linkage }); - if (long_double_is_f80) { - @export(fmal, .{ .name = "__fmax", .linkage = linkage }); - } else { - @export(__fmax, .{ .name = "__fmax", .linkage = linkage }); - } - if (long_double_is_f128) { - @export(fmal, .{ .name = "fmaq", .linkage = linkage }); - } else { - @export(fmaq, .{ .name = "fmaq", .linkage = linkage }); - } - } + mathExport("ceil", @import("./compiler_rt/ceil.zig")); + mathExport("cos", @import("./compiler_rt/cos.zig")); + mathExport("exp", @import("./compiler_rt/exp.zig")); + mathExport("exp2", @import("./compiler_rt/exp2.zig")); + mathExport("fabs", @import("./compiler_rt/fabs.zig")); + mathExport("floor", @import("./compiler_rt/floor.zig")); + mathExport("fma", @import("./compiler_rt/fma.zig")); + mathExport("fmax", @import("./compiler_rt/fmax.zig")); + mathExport("fmin", @import("./compiler_rt/fmin.zig")); + mathExport("fmod", @import("./compiler_rt/fmod.zig")); + mathExport("log", @import("./compiler_rt/log.zig")); + mathExport("log10", @import("./compiler_rt/log10.zig")); + mathExport("log2", @import("./compiler_rt/log2.zig")); + mathExport("round", @import("./compiler_rt/round.zig")); + mathExport("sin", @import("./compiler_rt/sin.zig")); + mathExport("sincos", @import("./compiler_rt/sincos.zig")); + mathExport("sqrt", @import("./compiler_rt/sqrt.zig")); + mathExport("tan", @import("./compiler_rt/tan.zig")); + mathExport("trunc", @import("./compiler_rt/trunc.zig")); if (arch.isSPARC()) { // SPARC systems use a different naming scheme @@ -842,63 +822,44 @@ comptime { @export(__unordtf2, .{ .name = "__unordkf2", .linkage = linkage }); // LLVM PPC backend lowers f128 fma to `fmaf128`. - @export(fmal, .{ .name = "fmaf128", .linkage = linkage }); + const fmaq = @import("./compiler_rt/fma.zig").fmaq; + @export(fmaq, .{ .name = "fmaf128", .linkage = linkage }); } } -const math = std.math; - -fn fmaf(a: f32, b: f32, c: f32) callconv(.C) f32 { - return math.fma(f32, a, b, c); -} -fn fma(a: f64, b: f64, c: f64) callconv(.C) f64 { - return math.fma(f64, a, b, c); -} -fn __fmax(a: f80, b: f80, c: f80) callconv(.C) f80 { - return math.fma(f80, a, b, c); -} -fn fmaq(a: f128, b: f128, c: f128) callconv(.C) f128 { - return math.fma(f128, a, b, c); -} -fn fmal(a: c_longdouble, b: c_longdouble, c: c_longdouble) callconv(.C) c_longdouble { - return math.fma(c_longdouble, a, b, c); -} - -// TODO add intrinsics for these (and probably the double version too) -// and have the math stuff use the intrinsic. same as @mod and @rem -fn floorf(x: f32) callconv(.C) f32 { - return math.floor(x); -} -fn floor(x: f64) callconv(.C) f64 { - return math.floor(x); -} -fn floorl(x: c_longdouble) callconv(.C) c_longdouble { - if (!long_double_is_f128) { - @panic("TODO implement this"); - } - return math.floor(x); -} - -fn ceilf(x: f32) callconv(.C) f32 { - return math.ceil(x); -} -fn ceil(x: f64) callconv(.C) f64 { - return math.ceil(x); -} -fn ceill(x: c_longdouble) callconv(.C) c_longdouble { - if (!long_double_is_f128) { - @panic("TODO implement this"); - } - return math.ceil(x); -} - -const fmodq = @import("compiler_rt/fmodq.zig").fmodq; -const fmodx = @import("compiler_rt/fmodx.zig").fmodx; -fn fmodl(x: c_longdouble, y: c_longdouble) callconv(.C) c_longdouble { - if (!long_double_is_f128) { - @panic("TODO implement this"); +inline fn mathExport(double_name: []const u8, comptime import: type) void { + const half_name = "__" ++ double_name ++ "h"; + const half_fn = @field(import, half_name); + const float_name = double_name ++ "f"; + const float_fn = @field(import, float_name); + const double_fn = @field(import, double_name); + const long_double_name = double_name ++ "l"; + const xf80_name = "__" ++ double_name ++ "x"; + const xf80_fn = @field(import, xf80_name); + const quad_name = double_name ++ "q"; + const quad_fn = @field(import, quad_name); + + @export(half_fn, .{ .name = half_name, .linkage = linkage }); + @export(float_fn, .{ .name = float_name, .linkage = linkage }); + @export(double_fn, .{ .name = double_name, .linkage = linkage }); + @export(xf80_fn, .{ .name = xf80_name, .linkage = linkage }); + @export(quad_fn, .{ .name = quad_name, .linkage = linkage }); + + const pairs = .{ + .{ f16, half_fn }, + .{ f32, float_fn }, + .{ f64, double_fn }, + .{ f80, xf80_fn }, + .{ f128, quad_fn }, + }; + + inline for (pairs) |pair| { + const F = pair[0]; + const func = pair[1]; + if (builtin.target.longDoubleIs(F)) { + @export(func, .{ .name = long_double_name, .linkage = linkage }); + } } - return @floatCast(c_longdouble, fmodq(x, y)); } // Avoid dragging in the runtime safety mechanisms into this .o file, diff --git a/lib/std/math/ceil.zig b/lib/std/special/compiler_rt/ceil.zig index 686be8e58d..c7087a2c3a 100644 --- a/lib/std/math/ceil.zig +++ b/lib/std/special/compiler_rt/ceil.zig @@ -4,31 +4,16 @@ // https://git.musl-libc.org/cgit/musl/tree/src/math/ceilf.c // https://git.musl-libc.org/cgit/musl/tree/src/math/ceil.c -const std = @import("../std.zig"); +const std = @import("std"); const math = std.math; const expect = std.testing.expect; -/// Returns the least integer value greater than of equal to x. -/// -/// Special Cases: -/// - ceil(+-0) = +-0 -/// - ceil(+-inf) = +-inf -/// - ceil(nan) = nan -pub fn ceil(x: anytype) @TypeOf(x) { - const T = @TypeOf(x); - return switch (T) { - f32 => ceil32(x), - f64 => ceil64(x), - f128 => ceil128(x), - - // TODO this is not correct for some targets - c_longdouble => @floatCast(c_longdouble, ceil128(x)), - - else => @compileError("ceil not implemented for " ++ @typeName(T)), - }; +pub fn __ceilh(x: f16) callconv(.C) f16 { + // TODO: more efficient implementation + return @floatCast(f16, ceilf(x)); } -fn ceil32(x: f32) f32 { +pub fn ceilf(x: f32) callconv(.C) f32 { var u = @bitCast(u32, x); var e = @intCast(i32, (u >> 23) & 0xFF) - 0x7F; var m: u32 = undefined; @@ -61,7 +46,7 @@ fn ceil32(x: f32) f32 { } } -fn ceil64(x: f64) f64 { +pub fn ceil(x: f64) callconv(.C) f64 { const f64_toint = 1.0 / math.floatEps(f64); const u = @bitCast(u64, x); @@ -92,7 +77,12 @@ fn ceil64(x: f64) f64 { } } -fn ceil128(x: f128) f128 { +pub fn __ceilx(x: f80) callconv(.C) f80 { + // TODO: more efficient implementation + return @floatCast(f80, ceilq(x)); +} + +pub fn ceilq(x: f128) callconv(.C) f128 { const f128_toint = 1.0 / math.floatEps(f128); const u = @bitCast(u128, x); @@ -121,50 +111,44 @@ fn ceil128(x: f128) f128 { } } -test "math.ceil" { - try expect(ceil(@as(f32, 0.0)) == ceil32(0.0)); - try expect(ceil(@as(f64, 0.0)) == ceil64(0.0)); - try expect(ceil(@as(f128, 0.0)) == ceil128(0.0)); -} - -test "math.ceil32" { - try expect(ceil32(1.3) == 2.0); - try expect(ceil32(-1.3) == -1.0); - try expect(ceil32(0.2) == 1.0); +test "ceil32" { + try expect(ceilf(1.3) == 2.0); + try expect(ceilf(-1.3) == -1.0); + try expect(ceilf(0.2) == 1.0); } -test "math.ceil64" { - try expect(ceil64(1.3) == 2.0); - try expect(ceil64(-1.3) == -1.0); - try expect(ceil64(0.2) == 1.0); +test "ceil64" { + try expect(ceil(1.3) == 2.0); + try expect(ceil(-1.3) == -1.0); + try expect(ceil(0.2) == 1.0); } -test "math.ceil128" { - try expect(ceil128(1.3) == 2.0); - try expect(ceil128(-1.3) == -1.0); - try expect(ceil128(0.2) == 1.0); +test "ceil128" { + try expect(ceilq(1.3) == 2.0); + try expect(ceilq(-1.3) == -1.0); + try expect(ceilq(0.2) == 1.0); } -test "math.ceil32.special" { - try expect(ceil32(0.0) == 0.0); - try expect(ceil32(-0.0) == -0.0); - try expect(math.isPositiveInf(ceil32(math.inf(f32)))); - try expect(math.isNegativeInf(ceil32(-math.inf(f32)))); - try expect(math.isNan(ceil32(math.nan(f32)))); +test "ceil32.special" { + try expect(ceilf(0.0) == 0.0); + try expect(ceilf(-0.0) == -0.0); + try expect(math.isPositiveInf(ceilf(math.inf(f32)))); + try expect(math.isNegativeInf(ceilf(-math.inf(f32)))); + try expect(math.isNan(ceilf(math.nan(f32)))); } -test "math.ceil64.special" { - try expect(ceil64(0.0) == 0.0); - try expect(ceil64(-0.0) == -0.0); - try expect(math.isPositiveInf(ceil64(math.inf(f64)))); - try expect(math.isNegativeInf(ceil64(-math.inf(f64)))); - try expect(math.isNan(ceil64(math.nan(f64)))); +test "ceil64.special" { + try expect(ceil(0.0) == 0.0); + try expect(ceil(-0.0) == -0.0); + try expect(math.isPositiveInf(ceil(math.inf(f64)))); + try expect(math.isNegativeInf(ceil(-math.inf(f64)))); + try expect(math.isNan(ceil(math.nan(f64)))); } -test "math.ceil128.special" { - try expect(ceil128(0.0) == 0.0); - try expect(ceil128(-0.0) == -0.0); - try expect(math.isPositiveInf(ceil128(math.inf(f128)))); - try expect(math.isNegativeInf(ceil128(-math.inf(f128)))); - try expect(math.isNan(ceil128(math.nan(f128)))); +test "ceil128.special" { + try expect(ceilq(0.0) == 0.0); + try expect(ceilq(-0.0) == -0.0); + try expect(math.isPositiveInf(ceilq(math.inf(f128)))); + try expect(math.isNegativeInf(ceilq(-math.inf(f128)))); + try expect(math.isNan(ceilq(math.nan(f128)))); } diff --git a/lib/std/math/cos.zig b/lib/std/special/compiler_rt/cos.zig index 22bae0daee..295f6a47ea 100644 --- a/lib/std/math/cos.zig +++ b/lib/std/special/compiler_rt/cos.zig @@ -1,32 +1,17 @@ -// Ported from musl, which is licensed under the MIT license: -// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT -// -// https://git.musl-libc.org/cgit/musl/tree/src/math/cosf.c -// https://git.musl-libc.org/cgit/musl/tree/src/math/cos.c - -const std = @import("../std.zig"); +const std = @import("std"); const math = std.math; const expect = std.testing.expect; -const kernel = @import("__trig.zig"); -const __rem_pio2 = @import("__rem_pio2.zig").__rem_pio2; -const __rem_pio2f = @import("__rem_pio2f.zig").__rem_pio2f; - -/// Returns the cosine of the radian value x. -/// -/// Special Cases: -/// - cos(+-inf) = nan -/// - cos(nan) = nan -pub fn cos(x: anytype) @TypeOf(x) { - const T = @TypeOf(x); - return switch (T) { - f32 => cos32(x), - f64 => cos64(x), - else => @compileError("cos not implemented for " ++ @typeName(T)), - }; +const kernel = @import("trig.zig"); +const rem_pio2 = @import("rem_pio2.zig").rem_pio2; +const rem_pio2f = @import("rem_pio2f.zig").rem_pio2f; + +pub fn __cosh(a: f16) callconv(.C) f16 { + // TODO: more efficient implementation + return @floatCast(f16, cosf(a)); } -fn cos32(x: f32) f32 { +pub fn cosf(x: f32) callconv(.C) f32 { // Small multiples of pi/2 rounded to double precision. const c1pio2: f64 = 1.0 * math.pi / 2.0; // 0x3FF921FB, 0x54442D18 const c2pio2: f64 = 2.0 * math.pi / 2.0; // 0x400921FB, 0x54442D18 @@ -74,7 +59,7 @@ fn cos32(x: f32) f32 { } var y: f64 = undefined; - const n = __rem_pio2f(x, &y); + const n = rem_pio2f(x, &y); return switch (n & 3) { 0 => kernel.__cosdf(y), 1 => kernel.__sindf(-y), @@ -83,7 +68,7 @@ fn cos32(x: f32) f32 { }; } -fn cos64(x: f64) f64 { +pub fn cos(x: f64) callconv(.C) f64 { var ix = @bitCast(u64, x) >> 32; ix &= 0x7fffffff; @@ -103,7 +88,7 @@ fn cos64(x: f64) f64 { } var y: [2]f64 = undefined; - const n = __rem_pio2(x, &y); + const n = rem_pio2(x, &y); return switch (n & 3) { 0 => kernel.__cos(y[0], y[1]), 1 => -kernel.__sin(y[0], y[1], 1), @@ -112,43 +97,48 @@ fn cos64(x: f64) f64 { }; } -test "math.cos" { - try expect(cos(@as(f32, 0.0)) == cos32(0.0)); - try expect(cos(@as(f64, 0.0)) == cos64(0.0)); +pub fn __cosx(a: f80) callconv(.C) f80 { + // TODO: more efficient implementation + return @floatCast(f80, cosq(a)); +} + +pub fn cosq(a: f128) callconv(.C) f128 { + // TODO: more correct implementation + return cos(@floatCast(f64, a)); } -test "math.cos32" { +test "cos32" { const epsilon = 0.00001; - try expect(math.approxEqAbs(f32, cos32(0.0), 1.0, epsilon)); - try expect(math.approxEqAbs(f32, cos32(0.2), 0.980067, epsilon)); - try expect(math.approxEqAbs(f32, cos32(0.8923), 0.627623, epsilon)); - try expect(math.approxEqAbs(f32, cos32(1.5), 0.070737, epsilon)); - try expect(math.approxEqAbs(f32, cos32(-1.5), 0.070737, epsilon)); - try expect(math.approxEqAbs(f32, cos32(37.45), 0.969132, epsilon)); - try expect(math.approxEqAbs(f32, cos32(89.123), 0.400798, epsilon)); + try expect(math.approxEqAbs(f32, cosf(0.0), 1.0, epsilon)); + try expect(math.approxEqAbs(f32, cosf(0.2), 0.980067, epsilon)); + try expect(math.approxEqAbs(f32, cosf(0.8923), 0.627623, epsilon)); + try expect(math.approxEqAbs(f32, cosf(1.5), 0.070737, epsilon)); + try expect(math.approxEqAbs(f32, cosf(-1.5), 0.070737, epsilon)); + try expect(math.approxEqAbs(f32, cosf(37.45), 0.969132, epsilon)); + try expect(math.approxEqAbs(f32, cosf(89.123), 0.400798, epsilon)); } -test "math.cos64" { +test "cos64" { const epsilon = 0.000001; - try expect(math.approxEqAbs(f64, cos64(0.0), 1.0, epsilon)); - try expect(math.approxEqAbs(f64, cos64(0.2), 0.980067, epsilon)); - try expect(math.approxEqAbs(f64, cos64(0.8923), 0.627623, epsilon)); - try expect(math.approxEqAbs(f64, cos64(1.5), 0.070737, epsilon)); - try expect(math.approxEqAbs(f64, cos64(-1.5), 0.070737, epsilon)); - try expect(math.approxEqAbs(f64, cos64(37.45), 0.969132, epsilon)); - try expect(math.approxEqAbs(f64, cos64(89.123), 0.40080, epsilon)); + try expect(math.approxEqAbs(f64, cos(0.0), 1.0, epsilon)); + try expect(math.approxEqAbs(f64, cos(0.2), 0.980067, epsilon)); + try expect(math.approxEqAbs(f64, cos(0.8923), 0.627623, epsilon)); + try expect(math.approxEqAbs(f64, cos(1.5), 0.070737, epsilon)); + try expect(math.approxEqAbs(f64, cos(-1.5), 0.070737, epsilon)); + try expect(math.approxEqAbs(f64, cos(37.45), 0.969132, epsilon)); + try expect(math.approxEqAbs(f64, cos(89.123), 0.40080, epsilon)); } -test "math.cos32.special" { - try expect(math.isNan(cos32(math.inf(f32)))); - try expect(math.isNan(cos32(-math.inf(f32)))); - try expect(math.isNan(cos32(math.nan(f32)))); +test "cos32.special" { + try expect(math.isNan(cosf(math.inf(f32)))); + try expect(math.isNan(cosf(-math.inf(f32)))); + try expect(math.isNan(cosf(math.nan(f32)))); } -test "math.cos64.special" { - try expect(math.isNan(cos64(math.inf(f64)))); - try expect(math.isNan(cos64(-math.inf(f64)))); - try expect(math.isNan(cos64(math.nan(f64)))); +test "cos64.special" { + try expect(math.isNan(cos(math.inf(f64)))); + try expect(math.isNan(cos(-math.inf(f64)))); + try expect(math.isNan(cos(math.nan(f64)))); } diff --git a/lib/std/special/compiler_rt/divxf3_test.zig b/lib/std/special/compiler_rt/divxf3_test.zig index b79b90c6fb..0ed2b74217 100644 --- a/lib/std/special/compiler_rt/divxf3_test.zig +++ b/lib/std/special/compiler_rt/divxf3_test.zig @@ -30,9 +30,9 @@ fn test__divxf3(a: f80, b: f80) !void { const x_minus_eps = @bitCast(f80, (@bitCast(u80, x) - 1) | integerBit); // Make sure result is more accurate than the adjacent floats - const err_x = std.math.fabs(@mulAdd(f80, x, b, -a)); - const err_x_plus_eps = std.math.fabs(@mulAdd(f80, x_plus_eps, b, -a)); - const err_x_minus_eps = std.math.fabs(@mulAdd(f80, x_minus_eps, b, -a)); + const err_x = @fabs(@mulAdd(f80, x, b, -a)); + const err_x_plus_eps = @fabs(@mulAdd(f80, x_plus_eps, b, -a)); + const err_x_minus_eps = @fabs(@mulAdd(f80, x_minus_eps, b, -a)); try testing.expect(err_x_minus_eps > err_x); try testing.expect(err_x_plus_eps > err_x); diff --git a/lib/std/math/exp.zig b/lib/std/special/compiler_rt/exp.zig index 71a492c7ad..0f129dfd4c 100644 --- a/lib/std/math/exp.zig +++ b/lib/std/special/compiler_rt/exp.zig @@ -4,25 +4,16 @@ // https://git.musl-libc.org/cgit/musl/tree/src/math/expf.c // https://git.musl-libc.org/cgit/musl/tree/src/math/exp.c -const std = @import("../std.zig"); +const std = @import("std"); const math = std.math; const expect = std.testing.expect; -/// Returns e raised to the power of x (e^x). -/// -/// Special Cases: -/// - exp(+inf) = +inf -/// - exp(nan) = nan -pub fn exp(x: anytype) @TypeOf(x) { - const T = @TypeOf(x); - return switch (T) { - f32 => exp32(x), - f64 => exp64(x), - else => @compileError("exp not implemented for " ++ @typeName(T)), - }; +pub fn __exph(a: f16) callconv(.C) f16 { + // TODO: more efficient implementation + return @floatCast(f16, expf(a)); } -fn exp32(x_: f32) f32 { +pub fn expf(x_: f32) callconv(.C) f32 { const half = [_]f32{ 0.5, -0.5 }; const ln2hi = 6.9314575195e-1; const ln2lo = 1.4286067653e-6; @@ -97,7 +88,7 @@ fn exp32(x_: f32) f32 { } } -fn exp64(x_: f64) f64 { +pub fn exp(x_: f64) callconv(.C) f64 { const half = [_]f64{ 0.5, -0.5 }; const ln2hi: f64 = 6.93147180369123816490e-01; const ln2lo: f64 = 1.90821492927058770002e-10; @@ -181,37 +172,42 @@ fn exp64(x_: f64) f64 { } } -test "math.exp" { - try expect(exp(@as(f32, 0.0)) == exp32(0.0)); - try expect(exp(@as(f64, 0.0)) == exp64(0.0)); +pub fn __expx(a: f80) callconv(.C) f80 { + // TODO: more efficient implementation + return @floatCast(f80, expq(a)); } -test "math.exp32" { +pub fn expq(a: f128) callconv(.C) f128 { + // TODO: more correct implementation + return exp(@floatCast(f64, a)); +} + +test "exp32" { const epsilon = 0.000001; - try expect(exp32(0.0) == 1.0); - try expect(math.approxEqAbs(f32, exp32(0.0), 1.0, epsilon)); - try expect(math.approxEqAbs(f32, exp32(0.2), 1.221403, epsilon)); - try expect(math.approxEqAbs(f32, exp32(0.8923), 2.440737, epsilon)); - try expect(math.approxEqAbs(f32, exp32(1.5), 4.481689, epsilon)); + try expect(expf(0.0) == 1.0); + try expect(math.approxEqAbs(f32, expf(0.0), 1.0, epsilon)); + try expect(math.approxEqAbs(f32, expf(0.2), 1.221403, epsilon)); + try expect(math.approxEqAbs(f32, expf(0.8923), 2.440737, epsilon)); + try expect(math.approxEqAbs(f32, expf(1.5), 4.481689, epsilon)); } -test "math.exp64" { +test "exp64" { const epsilon = 0.000001; - try expect(exp64(0.0) == 1.0); - try expect(math.approxEqAbs(f64, exp64(0.0), 1.0, epsilon)); - try expect(math.approxEqAbs(f64, exp64(0.2), 1.221403, epsilon)); - try expect(math.approxEqAbs(f64, exp64(0.8923), 2.440737, epsilon)); - try expect(math.approxEqAbs(f64, exp64(1.5), 4.481689, epsilon)); + try expect(exp(0.0) == 1.0); + try expect(math.approxEqAbs(f64, exp(0.0), 1.0, epsilon)); + try expect(math.approxEqAbs(f64, exp(0.2), 1.221403, epsilon)); + try expect(math.approxEqAbs(f64, exp(0.8923), 2.440737, epsilon)); + try expect(math.approxEqAbs(f64, exp(1.5), 4.481689, epsilon)); } -test "math.exp32.special" { - try expect(math.isPositiveInf(exp32(math.inf(f32)))); - try expect(math.isNan(exp32(math.nan(f32)))); +test "exp32.special" { + try expect(math.isPositiveInf(expf(math.inf(f32)))); + try expect(math.isNan(expf(math.nan(f32)))); } -test "math.exp64.special" { - try expect(math.isPositiveInf(exp64(math.inf(f64)))); - try expect(math.isNan(exp64(math.nan(f64)))); +test "exp64.special" { + try expect(math.isPositiveInf(exp(math.inf(f64)))); + try expect(math.isNan(exp(math.nan(f64)))); } diff --git a/lib/std/math/exp2.zig b/lib/std/special/compiler_rt/exp2.zig index 76530ec61f..53432a831d 100644 --- a/lib/std/math/exp2.zig +++ b/lib/std/special/compiler_rt/exp2.zig @@ -4,44 +4,16 @@ // https://git.musl-libc.org/cgit/musl/tree/src/math/exp2f.c // https://git.musl-libc.org/cgit/musl/tree/src/math/exp2.c -const std = @import("../std.zig"); +const std = @import("std"); const math = std.math; const expect = std.testing.expect; -/// Returns 2 raised to the power of x (2^x). -/// -/// Special Cases: -/// - exp2(+inf) = +inf -/// - exp2(nan) = nan -pub fn exp2(x: anytype) @TypeOf(x) { - const T = @TypeOf(x); - return switch (T) { - f32 => exp2_32(x), - f64 => exp2_64(x), - else => @compileError("exp2 not implemented for " ++ @typeName(T)), - }; +pub fn __exp2h(x: f16) callconv(.C) f16 { + // TODO: more efficient implementation + return @floatCast(f16, exp2f(x)); } -const exp2ft = [_]f64{ - 0x1.6a09e667f3bcdp-1, - 0x1.7a11473eb0187p-1, - 0x1.8ace5422aa0dbp-1, - 0x1.9c49182a3f090p-1, - 0x1.ae89f995ad3adp-1, - 0x1.c199bdd85529cp-1, - 0x1.d5818dcfba487p-1, - 0x1.ea4afa2a490dap-1, - 0x1.0000000000000p+0, - 0x1.0b5586cf9890fp+0, - 0x1.172b83c7d517bp+0, - 0x1.2387a6e756238p+0, - 0x1.306fe0a31b715p+0, - 0x1.3dea64c123422p+0, - 0x1.4bfdad5362a27p+0, - 0x1.5ab07dd485429p+0, -}; - -fn exp2_32(x: f32) f32 { +pub fn exp2f(x: f32) callconv(.C) f32 { const tblsiz = @intCast(u32, exp2ft.len); const redux: f32 = 0x1.8p23 / @intToFloat(f32, tblsiz); const P1: f32 = 0x1.62e430p-1; @@ -98,6 +70,104 @@ fn exp2_32(x: f32) f32 { return @floatCast(f32, r * uk); } +pub fn exp2(x: f64) callconv(.C) f64 { + const tblsiz: u32 = @intCast(u32, exp2dt.len / 2); + const redux: f64 = 0x1.8p52 / @intToFloat(f64, tblsiz); + const P1: f64 = 0x1.62e42fefa39efp-1; + const P2: f64 = 0x1.ebfbdff82c575p-3; + const P3: f64 = 0x1.c6b08d704a0a6p-5; + const P4: f64 = 0x1.3b2ab88f70400p-7; + const P5: f64 = 0x1.5d88003875c74p-10; + + const ux = @bitCast(u64, x); + const ix = @intCast(u32, ux >> 32) & 0x7FFFFFFF; + + // TODO: This should be handled beneath. + if (math.isNan(x)) { + return math.nan(f64); + } + + // |x| >= 1022 or nan + if (ix >= 0x408FF000) { + // x >= 1024 or nan + if (ix >= 0x40900000 and ux >> 63 == 0) { + math.raiseOverflow(); + return math.inf(f64); + } + // -inf or -nan + if (ix >= 0x7FF00000) { + return -1 / x; + } + // x <= -1022 + if (ux >> 63 != 0) { + // underflow + if (x <= -1075 or x - 0x1.0p52 + 0x1.0p52 != x) { + math.doNotOptimizeAway(@floatCast(f32, -0x1.0p-149 / x)); + } + if (x <= -1075) { + return 0; + } + } + } + // |x| < 0x1p-54 + else if (ix < 0x3C900000) { + return 1.0 + x; + } + + // NOTE: musl relies on unsafe behaviours which are replicated below + // (addition overflow, division truncation, casting). Appears that this + // produces the intended result but should confirm how GCC/Clang handle this + // to ensure. + + // reduce x + var uf: f64 = x + redux; + // NOTE: musl performs an implicit 64-bit to 32-bit u32 truncation here + var i_0: u32 = @truncate(u32, @bitCast(u64, uf)); + i_0 +%= tblsiz / 2; + + const k: u32 = i_0 / tblsiz * tblsiz; + const ik: i32 = @divTrunc(@bitCast(i32, k), tblsiz); + i_0 %= tblsiz; + uf -= redux; + + // r = exp2(y) = exp2t[i_0] * p(z - eps[i]) + var z: f64 = x - uf; + const t: f64 = exp2dt[@intCast(usize, 2 * i_0)]; + z -= exp2dt[@intCast(usize, 2 * i_0 + 1)]; + const r: f64 = t + t * z * (P1 + z * (P2 + z * (P3 + z * (P4 + z * P5)))); + + return math.scalbn(r, ik); +} + +pub fn __exp2x(x: f80) callconv(.C) f80 { + // TODO: more efficient implementation + return @floatCast(f80, exp2q(x)); +} + +pub fn exp2q(x: f128) callconv(.C) f128 { + // TODO: more correct implementation + return exp2(@floatCast(f64, x)); +} + +const exp2ft = [_]f64{ + 0x1.6a09e667f3bcdp-1, + 0x1.7a11473eb0187p-1, + 0x1.8ace5422aa0dbp-1, + 0x1.9c49182a3f090p-1, + 0x1.ae89f995ad3adp-1, + 0x1.c199bdd85529cp-1, + 0x1.d5818dcfba487p-1, + 0x1.ea4afa2a490dap-1, + 0x1.0000000000000p+0, + 0x1.0b5586cf9890fp+0, + 0x1.172b83c7d517bp+0, + 0x1.2387a6e756238p+0, + 0x1.306fe0a31b715p+0, + 0x1.3dea64c123422p+0, + 0x1.4bfdad5362a27p+0, + 0x1.5ab07dd485429p+0, +}; + const exp2dt = [_]f64{ // exp2(z + eps) eps 0x1.6a09e667f3d5dp-1, 0x1.9880p-44, @@ -358,108 +428,34 @@ const exp2dt = [_]f64{ 0x1.690f4b19e9471p+0, -0x1.9780p-45, }; -fn exp2_64(x: f64) f64 { - const tblsiz: u32 = @intCast(u32, exp2dt.len / 2); - const redux: f64 = 0x1.8p52 / @intToFloat(f64, tblsiz); - const P1: f64 = 0x1.62e42fefa39efp-1; - const P2: f64 = 0x1.ebfbdff82c575p-3; - const P3: f64 = 0x1.c6b08d704a0a6p-5; - const P4: f64 = 0x1.3b2ab88f70400p-7; - const P5: f64 = 0x1.5d88003875c74p-10; - - const ux = @bitCast(u64, x); - const ix = @intCast(u32, ux >> 32) & 0x7FFFFFFF; - - // TODO: This should be handled beneath. - if (math.isNan(x)) { - return math.nan(f64); - } - - // |x| >= 1022 or nan - if (ix >= 0x408FF000) { - // x >= 1024 or nan - if (ix >= 0x40900000 and ux >> 63 == 0) { - math.raiseOverflow(); - return math.inf(f64); - } - // -inf or -nan - if (ix >= 0x7FF00000) { - return -1 / x; - } - // x <= -1022 - if (ux >> 63 != 0) { - // underflow - if (x <= -1075 or x - 0x1.0p52 + 0x1.0p52 != x) { - math.doNotOptimizeAway(@floatCast(f32, -0x1.0p-149 / x)); - } - if (x <= -1075) { - return 0; - } - } - } - // |x| < 0x1p-54 - else if (ix < 0x3C900000) { - return 1.0 + x; - } - - // NOTE: musl relies on unsafe behaviours which are replicated below - // (addition overflow, division truncation, casting). Appears that this - // produces the intended result but should confirm how GCC/Clang handle this - // to ensure. - - // reduce x - var uf: f64 = x + redux; - // NOTE: musl performs an implicit 64-bit to 32-bit u32 truncation here - var i_0: u32 = @truncate(u32, @bitCast(u64, uf)); - i_0 +%= tblsiz / 2; - - const k: u32 = i_0 / tblsiz * tblsiz; - const ik: i32 = @divTrunc(@bitCast(i32, k), tblsiz); - i_0 %= tblsiz; - uf -= redux; - - // r = exp2(y) = exp2t[i_0] * p(z - eps[i]) - var z: f64 = x - uf; - const t: f64 = exp2dt[@intCast(usize, 2 * i_0)]; - z -= exp2dt[@intCast(usize, 2 * i_0 + 1)]; - const r: f64 = t + t * z * (P1 + z * (P2 + z * (P3 + z * (P4 + z * P5)))); - - return math.scalbn(r, ik); -} - -test "math.exp2" { - try expect(exp2(@as(f32, 0.8923)) == exp2_32(0.8923)); - try expect(exp2(@as(f64, 0.8923)) == exp2_64(0.8923)); -} - -test "math.exp2_32" { +test "exp2_32" { const epsilon = 0.000001; - try expect(exp2_32(0.0) == 1.0); - try expect(math.approxEqAbs(f32, exp2_32(0.2), 1.148698, epsilon)); - try expect(math.approxEqAbs(f32, exp2_32(0.8923), 1.856133, epsilon)); - try expect(math.approxEqAbs(f32, exp2_32(1.5), 2.828427, epsilon)); - try expect(math.approxEqAbs(f32, exp2_32(37.45), 187747237888, epsilon)); - try expect(math.approxEqAbs(f32, exp2_32(-1), 0.5, epsilon)); + try expect(exp2f(0.0) == 1.0); + try expect(math.approxEqAbs(f32, exp2f(0.2), 1.148698, epsilon)); + try expect(math.approxEqAbs(f32, exp2f(0.8923), 1.856133, epsilon)); + try expect(math.approxEqAbs(f32, exp2f(1.5), 2.828427, epsilon)); + try expect(math.approxEqAbs(f32, exp2f(37.45), 187747237888, epsilon)); + try expect(math.approxEqAbs(f32, exp2f(-1), 0.5, epsilon)); } -test "math.exp2_64" { +test "exp2_64" { const epsilon = 0.000001; - try expect(exp2_64(0.0) == 1.0); - try expect(math.approxEqAbs(f64, exp2_64(0.2), 1.148698, epsilon)); - try expect(math.approxEqAbs(f64, exp2_64(0.8923), 1.856133, epsilon)); - try expect(math.approxEqAbs(f64, exp2_64(1.5), 2.828427, epsilon)); - try expect(math.approxEqAbs(f64, exp2_64(-1), 0.5, epsilon)); - try expect(math.approxEqAbs(f64, exp2_64(-0x1.a05cc754481d1p-2), 0x1.824056efc687cp-1, epsilon)); + try expect(exp2(0.0) == 1.0); + try expect(math.approxEqAbs(f64, exp2(0.2), 1.148698, epsilon)); + try expect(math.approxEqAbs(f64, exp2(0.8923), 1.856133, epsilon)); + try expect(math.approxEqAbs(f64, exp2(1.5), 2.828427, epsilon)); + try expect(math.approxEqAbs(f64, exp2(-1), 0.5, epsilon)); + try expect(math.approxEqAbs(f64, exp2(-0x1.a05cc754481d1p-2), 0x1.824056efc687cp-1, epsilon)); } -test "math.exp2_32.special" { - try expect(math.isPositiveInf(exp2_32(math.inf(f32)))); - try expect(math.isNan(exp2_32(math.nan(f32)))); +test "exp2_32.special" { + try expect(math.isPositiveInf(exp2f(math.inf(f32)))); + try expect(math.isNan(exp2f(math.nan(f32)))); } -test "math.exp2_64.special" { - try expect(math.isPositiveInf(exp2_64(math.inf(f64)))); - try expect(math.isNan(exp2_64(math.nan(f64)))); +test "exp2_64.special" { + try expect(math.isPositiveInf(exp2(math.inf(f64)))); + try expect(math.isNan(exp2(math.nan(f64)))); } diff --git a/lib/std/special/compiler_rt/fabs.zig b/lib/std/special/compiler_rt/fabs.zig new file mode 100644 index 0000000000..fbef81fc9a --- /dev/null +++ b/lib/std/special/compiler_rt/fabs.zig @@ -0,0 +1,29 @@ +const std = @import("std"); + +pub fn __fabsh(a: f16) callconv(.C) f16 { + return generic_fabs(a); +} + +pub fn fabsf(a: f32) callconv(.C) f32 { + return generic_fabs(a); +} + +pub fn fabs(a: f64) callconv(.C) f64 { + return generic_fabs(a); +} + +pub fn __fabsx(a: f80) callconv(.C) f80 { + return generic_fabs(a); +} + +pub fn fabsq(a: f128) callconv(.C) f128 { + return generic_fabs(a); +} + +inline fn generic_fabs(x: anytype) @TypeOf(x) { + const T = @TypeOf(x); + const TBits = std.meta.Int(.unsigned, @typeInfo(T).Float.bits); + const float_bits = @bitCast(TBits, x); + const remove_sign = ~@as(TBits, 0) >> 1; + return @bitCast(T, float_bits & remove_sign); +} diff --git a/lib/std/math/floor.zig b/lib/std/special/compiler_rt/floor.zig index ab5ca3583b..f6df164b58 100644 --- a/lib/std/math/floor.zig +++ b/lib/std/special/compiler_rt/floor.zig @@ -4,32 +4,11 @@ // https://git.musl-libc.org/cgit/musl/tree/src/math/floorf.c // https://git.musl-libc.org/cgit/musl/tree/src/math/floor.c -const expect = std.testing.expect; -const std = @import("../std.zig"); +const std = @import("std"); const math = std.math; +const expect = std.testing.expect; -/// Returns the greatest integer value less than or equal to x. -/// -/// Special Cases: -/// - floor(+-0) = +-0 -/// - floor(+-inf) = +-inf -/// - floor(nan) = nan -pub fn floor(x: anytype) @TypeOf(x) { - const T = @TypeOf(x); - return switch (T) { - f16 => floor16(x), - f32 => floor32(x), - f64 => floor64(x), - f128 => floor128(x), - - // TODO this is not correct for some targets - c_longdouble => @floatCast(c_longdouble, floor128(x)), - - else => @compileError("floor not implemented for " ++ @typeName(T)), - }; -} - -fn floor16(x: f16) f16 { +pub fn __floorh(x: f16) callconv(.C) f16 { var u = @bitCast(u16, x); const e = @intCast(i16, (u >> 10) & 31) - 15; var m: u16 = undefined; @@ -63,7 +42,7 @@ fn floor16(x: f16) f16 { } } -fn floor32(x: f32) f32 { +pub fn floorf(x: f32) callconv(.C) f32 { var u = @bitCast(u32, x); const e = @intCast(i32, (u >> 23) & 0xFF) - 0x7F; var m: u32 = undefined; @@ -97,7 +76,7 @@ fn floor32(x: f32) f32 { } } -fn floor64(x: f64) f64 { +pub fn floor(x: f64) callconv(.C) f64 { const f64_toint = 1.0 / math.floatEps(f64); const u = @bitCast(u64, x); @@ -128,7 +107,12 @@ fn floor64(x: f64) f64 { } } -fn floor128(x: f128) f128 { +pub fn __floorx(x: f80) callconv(.C) f80 { + // TODO: more efficient implementation + return @floatCast(f80, floorq(x)); +} + +pub fn floorq(x: f128) callconv(.C) f128 { const f128_toint = 1.0 / math.floatEps(f128); const u = @bitCast(u128, x); @@ -157,65 +141,58 @@ fn floor128(x: f128) f128 { } } -test "math.floor" { - try expect(floor(@as(f16, 1.3)) == floor16(1.3)); - try expect(floor(@as(f32, 1.3)) == floor32(1.3)); - try expect(floor(@as(f64, 1.3)) == floor64(1.3)); - try expect(floor(@as(f128, 1.3)) == floor128(1.3)); -} - -test "math.floor16" { - try expect(floor16(1.3) == 1.0); - try expect(floor16(-1.3) == -2.0); - try expect(floor16(0.2) == 0.0); +test "floor16" { + try expect(__floorh(1.3) == 1.0); + try expect(__floorh(-1.3) == -2.0); + try expect(__floorh(0.2) == 0.0); } -test "math.floor32" { - try expect(floor32(1.3) == 1.0); - try expect(floor32(-1.3) == -2.0); - try expect(floor32(0.2) == 0.0); +test "floor32" { + try expect(floorf(1.3) == 1.0); + try expect(floorf(-1.3) == -2.0); + try expect(floorf(0.2) == 0.0); } -test "math.floor64" { - try expect(floor64(1.3) == 1.0); - try expect(floor64(-1.3) == -2.0); - try expect(floor64(0.2) == 0.0); +test "floor64" { + try expect(floor(1.3) == 1.0); + try expect(floor(-1.3) == -2.0); + try expect(floor(0.2) == 0.0); } -test "math.floor128" { - try expect(floor128(1.3) == 1.0); - try expect(floor128(-1.3) == -2.0); - try expect(floor128(0.2) == 0.0); +test "floor128" { + try expect(floorq(1.3) == 1.0); + try expect(floorq(-1.3) == -2.0); + try expect(floorq(0.2) == 0.0); } -test "math.floor16.special" { - try expect(floor16(0.0) == 0.0); - try expect(floor16(-0.0) == -0.0); - try expect(math.isPositiveInf(floor16(math.inf(f16)))); - try expect(math.isNegativeInf(floor16(-math.inf(f16)))); - try expect(math.isNan(floor16(math.nan(f16)))); +test "floor16.special" { + try expect(__floorh(0.0) == 0.0); + try expect(__floorh(-0.0) == -0.0); + try expect(math.isPositiveInf(__floorh(math.inf(f16)))); + try expect(math.isNegativeInf(__floorh(-math.inf(f16)))); + try expect(math.isNan(__floorh(math.nan(f16)))); } -test "math.floor32.special" { - try expect(floor32(0.0) == 0.0); - try expect(floor32(-0.0) == -0.0); - try expect(math.isPositiveInf(floor32(math.inf(f32)))); - try expect(math.isNegativeInf(floor32(-math.inf(f32)))); - try expect(math.isNan(floor32(math.nan(f32)))); +test "floor32.special" { + try expect(floorf(0.0) == 0.0); + try expect(floorf(-0.0) == -0.0); + try expect(math.isPositiveInf(floorf(math.inf(f32)))); + try expect(math.isNegativeInf(floorf(-math.inf(f32)))); + try expect(math.isNan(floorf(math.nan(f32)))); } -test "math.floor64.special" { - try expect(floor64(0.0) == 0.0); - try expect(floor64(-0.0) == -0.0); - try expect(math.isPositiveInf(floor64(math.inf(f64)))); - try expect(math.isNegativeInf(floor64(-math.inf(f64)))); - try expect(math.isNan(floor64(math.nan(f64)))); +test "floor64.special" { + try expect(floor(0.0) == 0.0); + try expect(floor(-0.0) == -0.0); + try expect(math.isPositiveInf(floor(math.inf(f64)))); + try expect(math.isNegativeInf(floor(-math.inf(f64)))); + try expect(math.isNan(floor(math.nan(f64)))); } -test "math.floor128.special" { - try expect(floor128(0.0) == 0.0); - try expect(floor128(-0.0) == -0.0); - try expect(math.isPositiveInf(floor128(math.inf(f128)))); - try expect(math.isNegativeInf(floor128(-math.inf(f128)))); - try expect(math.isNan(floor128(math.nan(f128)))); +test "floor128.special" { + try expect(floorq(0.0) == 0.0); + try expect(floorq(-0.0) == -0.0); + try expect(math.isPositiveInf(floorq(math.inf(f128)))); + try expect(math.isNegativeInf(floorq(-math.inf(f128)))); + try expect(math.isNan(floorq(math.nan(f128)))); } diff --git a/lib/std/math/fma.zig b/lib/std/special/compiler_rt/fma.zig index 7afc6e557e..4c603bf095 100644 --- a/lib/std/math/fma.zig +++ b/lib/std/special/compiler_rt/fma.zig @@ -5,27 +5,16 @@ // https://git.musl-libc.org/cgit/musl/tree/src/math/fmaf.c // https://git.musl-libc.org/cgit/musl/tree/src/math/fma.c -const std = @import("../std.zig"); +const std = @import("std"); const math = std.math; const expect = std.testing.expect; -/// Returns x * y + z with a single rounding error. -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), - - // TODO this is not correct for some targets - c_longdouble => @floatCast(c_longdouble, fma128(x, y, z)), - - f80 => @floatCast(f80, fma128(x, y, z)), - - else => @compileError("fma not implemented for " ++ @typeName(T)), - }; +pub fn __fmah(x: f16, y: f16, z: f16) callconv(.C) f16 { + // TODO: more efficient implementation + return @floatCast(f16, fmaf(x, y, z)); } -fn fma32(x: f32, y: f32, z: f32) f32 { +pub fn fmaf(x: f32, y: f32, z: f32) callconv(.C) f32 { const xy = @as(f64, x) * y; const xy_z = xy + z; const u = @bitCast(u64, xy_z); @@ -39,8 +28,8 @@ fn fma32(x: f32, y: f32, z: f32) f32 { } } -// NOTE: Upstream fma.c has been rewritten completely to raise fp exceptions more accurately. -fn fma64(x: f64, y: f64, z: f64) f64 { +/// NOTE: Upstream fma.c has been rewritten completely to raise fp exceptions more accurately. +pub fn fma(x: f64, y: f64, z: f64) callconv(.C) f64 { if (!math.isFinite(x) or !math.isFinite(y)) { return x * y + z; } @@ -87,6 +76,65 @@ fn fma64(x: f64, y: f64, z: f64) f64 { } } +pub fn __fmax(a: f80, b: f80, c: f80) callconv(.C) f80 { + // TODO: more efficient implementation + return @floatCast(f80, fmaq(a, b, c)); +} + +/// 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 { + 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.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); + } +} + const dd = struct { hi: f64, lo: f64, @@ -242,98 +290,38 @@ fn dd_mul128(a: f128, b: f128) dd128 { 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.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); - } -} - -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)); - try expect(math.approxEqAbs(f32, fma32(0.2, 5.0, 9.124), 10.124, epsilon)); - try expect(math.approxEqAbs(f32, fma32(0.8923, 5.0, 9.124), 13.5855, epsilon)); - try expect(math.approxEqAbs(f32, fma32(1.5, 5.0, 9.124), 16.624, epsilon)); - try expect(math.approxEqAbs(f32, fma32(37.45, 5.0, 9.124), 196.374004, epsilon)); - try expect(math.approxEqAbs(f32, fma32(89.123, 5.0, 9.124), 454.739005, epsilon)); - try expect(math.approxEqAbs(f32, fma32(123123.234375, 5.0, 9.124), 615625.295875, epsilon)); + try expect(math.approxEqAbs(f32, fmaf(0.0, 5.0, 9.124), 9.124, epsilon)); + try expect(math.approxEqAbs(f32, fmaf(0.2, 5.0, 9.124), 10.124, epsilon)); + try expect(math.approxEqAbs(f32, fmaf(0.8923, 5.0, 9.124), 13.5855, epsilon)); + try expect(math.approxEqAbs(f32, fmaf(1.5, 5.0, 9.124), 16.624, epsilon)); + try expect(math.approxEqAbs(f32, fmaf(37.45, 5.0, 9.124), 196.374004, epsilon)); + try expect(math.approxEqAbs(f32, fmaf(89.123, 5.0, 9.124), 454.739005, epsilon)); + try expect(math.approxEqAbs(f32, fmaf(123123.234375, 5.0, 9.124), 615625.295875, epsilon)); } test "64" { const epsilon = 0.000001; - try expect(math.approxEqAbs(f64, fma64(0.0, 5.0, 9.124), 9.124, epsilon)); - try expect(math.approxEqAbs(f64, fma64(0.2, 5.0, 9.124), 10.124, epsilon)); - try expect(math.approxEqAbs(f64, fma64(0.8923, 5.0, 9.124), 13.5855, epsilon)); - try expect(math.approxEqAbs(f64, fma64(1.5, 5.0, 9.124), 16.624, epsilon)); - try expect(math.approxEqAbs(f64, fma64(37.45, 5.0, 9.124), 196.374, epsilon)); - 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)); + try expect(math.approxEqAbs(f64, fma(0.0, 5.0, 9.124), 9.124, epsilon)); + try expect(math.approxEqAbs(f64, fma(0.2, 5.0, 9.124), 10.124, epsilon)); + try expect(math.approxEqAbs(f64, fma(0.8923, 5.0, 9.124), 13.5855, epsilon)); + try expect(math.approxEqAbs(f64, fma(1.5, 5.0, 9.124), 16.624, epsilon)); + try expect(math.approxEqAbs(f64, fma(37.45, 5.0, 9.124), 196.374, epsilon)); + try expect(math.approxEqAbs(f64, fma(89.123, 5.0, 9.124), 454.739, epsilon)); + try expect(math.approxEqAbs(f64, fma(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)); + try expect(math.approxEqAbs(f128, fmaq(0.0, 5.0, 9.124), 9.124, epsilon)); + try expect(math.approxEqAbs(f128, fmaq(0.2, 5.0, 9.124), 10.124, epsilon)); + try expect(math.approxEqAbs(f128, fmaq(0.8923, 5.0, 9.124), 13.5855, epsilon)); + try expect(math.approxEqAbs(f128, fmaq(1.5, 5.0, 9.124), 16.624, epsilon)); + try expect(math.approxEqAbs(f128, fmaq(37.45, 5.0, 9.124), 196.374, epsilon)); + try expect(math.approxEqAbs(f128, fmaq(89.123, 5.0, 9.124), 454.739, epsilon)); + try expect(math.approxEqAbs(f128, fmaq(123123.234375, 5.0, 9.124), 615625.295875, epsilon)); } diff --git a/lib/std/special/compiler_rt/fmax.zig b/lib/std/special/compiler_rt/fmax.zig new file mode 100644 index 0000000000..a5bd68cd74 --- /dev/null +++ b/lib/std/special/compiler_rt/fmax.zig @@ -0,0 +1,43 @@ +const std = @import("std"); +const math = std.math; + +pub fn __fmaxh(x: f16, y: f16) callconv(.C) f16 { + return generic_fmax(f16, x, y); +} + +pub fn fmaxf(x: f32, y: f32) callconv(.C) f32 { + return generic_fmax(f32, x, y); +} + +pub fn fmax(x: f64, y: f64) callconv(.C) f64 { + return generic_fmax(f64, x, y); +} + +pub fn __fmaxx(x: f80, y: f80) callconv(.C) f80 { + return generic_fmax(f80, x, y); +} + +pub fn fmaxq(x: f128, y: f128) callconv(.C) f128 { + return generic_fmax(f128, x, y); +} + +inline fn generic_fmax(comptime T: type, x: T, y: T) T { + if (math.isNan(x)) + return y; + if (math.isNan(y)) + return x; + return if (x < y) y else x; +} + +test "generic_fmax" { + inline for ([_]type{ f32, f64, c_longdouble, f80, f128 }) |T| { + const nan_val = math.nan(T); + + try std.testing.expect(math.isNan(generic_fmax(T, nan_val, nan_val))); + try std.testing.expectEqual(@as(T, 1.0), generic_fmax(T, nan_val, 1.0)); + try std.testing.expectEqual(@as(T, 1.0), generic_fmax(T, 1.0, nan_val)); + + try std.testing.expectEqual(@as(T, 10.0), generic_fmax(T, 1.0, 10.0)); + try std.testing.expectEqual(@as(T, 1.0), generic_fmax(T, 1.0, -1.0)); + } +} diff --git a/lib/std/special/compiler_rt/fmin.zig b/lib/std/special/compiler_rt/fmin.zig new file mode 100644 index 0000000000..cc4dbf082b --- /dev/null +++ b/lib/std/special/compiler_rt/fmin.zig @@ -0,0 +1,43 @@ +const std = @import("std"); +const math = std.math; + +pub fn __fminh(x: f16, y: f16) callconv(.C) f16 { + return generic_fmin(f16, x, y); +} + +pub fn fminf(x: f32, y: f32) callconv(.C) f32 { + return generic_fmin(f32, x, y); +} + +pub fn fmin(x: f64, y: f64) callconv(.C) f64 { + return generic_fmin(f64, x, y); +} + +pub fn __fminx(x: f80, y: f80) callconv(.C) f80 { + return generic_fmin(f80, x, y); +} + +pub fn fminq(x: f128, y: f128) callconv(.C) f128 { + return generic_fmin(f128, x, y); +} + +inline fn generic_fmin(comptime T: type, x: T, y: T) T { + if (math.isNan(x)) + return y; + if (math.isNan(y)) + return x; + return if (x < y) x else y; +} + +test "generic_fmin" { + inline for ([_]type{ f32, f64, c_longdouble, f80, f128 }) |T| { + const nan_val = math.nan(T); + + try std.testing.expect(math.isNan(generic_fmin(T, nan_val, nan_val))); + try std.testing.expectEqual(@as(T, 1.0), generic_fmin(T, nan_val, 1.0)); + try std.testing.expectEqual(@as(T, 1.0), generic_fmin(T, 1.0, nan_val)); + + try std.testing.expectEqual(@as(T, 1.0), generic_fmin(T, 1.0, 10.0)); + try std.testing.expectEqual(@as(T, -1.0), generic_fmin(T, 1.0, -1.0)); + } +} diff --git a/lib/std/special/compiler_rt/fmod.zig b/lib/std/special/compiler_rt/fmod.zig new file mode 100644 index 0000000000..b9a5710b9c --- /dev/null +++ b/lib/std/special/compiler_rt/fmod.zig @@ -0,0 +1,351 @@ +const builtin = @import("builtin"); +const std = @import("std"); +const math = std.math; +const assert = std.debug.assert; +const normalize = @import("divdf3.zig").normalize; + +pub fn __fmodh(x: f16, y: f16) callconv(.C) f16 { + // TODO: more efficient implementation + return @floatCast(f16, fmodf(x, y)); +} + +pub fn fmodf(x: f32, y: f32) callconv(.C) f32 { + return generic_fmod(f32, x, y); +} + +pub fn fmod(x: f64, y: f64) callconv(.C) f64 { + return generic_fmod(f64, x, y); +} + +/// fmodx - floating modulo large, returns the remainder of division for f80 types +/// Logic and flow heavily inspired by MUSL fmodl for 113 mantissa digits +pub fn __fmodx(a: f80, b: f80) callconv(.C) f80 { + @setRuntimeSafety(builtin.is_test); + + const T = f80; + const Z = std.meta.Int(.unsigned, @bitSizeOf(T)); + + const significandBits = math.floatMantissaBits(T); + const fractionalBits = math.floatFractionalBits(T); + const exponentBits = math.floatExponentBits(T); + + const signBit = (@as(Z, 1) << (significandBits + exponentBits)); + const maxExponent = ((1 << exponentBits) - 1); + + var aRep = @bitCast(Z, a); + var bRep = @bitCast(Z, b); + + const signA = aRep & signBit; + var expA = @intCast(i32, (@bitCast(Z, a) >> significandBits) & maxExponent); + var expB = @intCast(i32, (@bitCast(Z, b) >> significandBits) & maxExponent); + + // There are 3 cases where the answer is undefined, check for: + // - fmodx(val, 0) + // - fmodx(val, NaN) + // - fmodx(inf, val) + // The sign on checked values does not matter. + // Doing (a * b) / (a * b) procudes undefined results + // because the three cases always produce undefined calculations: + // - 0 / 0 + // - val * NaN + // - inf / inf + if (b == 0 or math.isNan(b) or expA == maxExponent) { + return (a * b) / (a * b); + } + + // Remove the sign from both + aRep &= ~signBit; + bRep &= ~signBit; + if (aRep <= bRep) { + if (aRep == bRep) { + return 0 * a; + } + return a; + } + + if (expA == 0) expA = normalize(f80, &aRep); + if (expB == 0) expB = normalize(f80, &bRep); + + var highA: u64 = 0; + var highB: u64 = 0; + var lowA: u64 = @truncate(u64, aRep); + var lowB: u64 = @truncate(u64, bRep); + + while (expA > expB) : (expA -= 1) { + var high = highA -% highB; + var low = lowA -% lowB; + if (lowA < lowB) { + high -%= 1; + } + if (high >> 63 == 0) { + if ((high | low) == 0) { + return 0 * a; + } + highA = 2 *% high + (low >> 63); + lowA = 2 *% low; + } else { + highA = 2 *% highA + (lowA >> 63); + lowA = 2 *% lowA; + } + } + + var high = highA -% highB; + var low = lowA -% lowB; + if (lowA < lowB) { + high -%= 1; + } + if (high >> 63 == 0) { + if ((high | low) == 0) { + return 0 * a; + } + highA = high; + lowA = low; + } + + while ((lowA >> fractionalBits) == 0) { + lowA = 2 *% lowA; + expA = expA - 1; + } + + // Combine the exponent with the sign and significand, normalize if happened to be denormalized + if (expA < -fractionalBits) { + return @bitCast(T, signA); + } else if (expA <= 0) { + return @bitCast(T, (lowA >> @intCast(math.Log2Int(u64), 1 - expA)) | signA); + } else { + return @bitCast(T, lowA | (@as(Z, @intCast(u16, expA)) << significandBits) | signA); + } +} + +/// fmodq - floating modulo large, returns the remainder of division for f128 types +/// Logic and flow heavily inspired by MUSL fmodl for 113 mantissa digits +pub fn fmodq(a: f128, b: f128) callconv(.C) f128 { + @setRuntimeSafety(builtin.is_test); + var amod = a; + var bmod = b; + const aPtr_u64 = @ptrCast([*]u64, &amod); + const bPtr_u64 = @ptrCast([*]u64, &bmod); + const aPtr_u16 = @ptrCast([*]u16, &amod); + const bPtr_u16 = @ptrCast([*]u16, &bmod); + + const exp_and_sign_index = comptime switch (builtin.target.cpu.arch.endian()) { + .Little => 7, + .Big => 0, + }; + const low_index = comptime switch (builtin.target.cpu.arch.endian()) { + .Little => 0, + .Big => 1, + }; + const high_index = comptime switch (builtin.target.cpu.arch.endian()) { + .Little => 1, + .Big => 0, + }; + + const signA = aPtr_u16[exp_and_sign_index] & 0x8000; + var expA = @intCast(i32, (aPtr_u16[exp_and_sign_index] & 0x7fff)); + var expB = @intCast(i32, (bPtr_u16[exp_and_sign_index] & 0x7fff)); + + // There are 3 cases where the answer is undefined, check for: + // - fmodq(val, 0) + // - fmodq(val, NaN) + // - fmodq(inf, val) + // The sign on checked values does not matter. + // Doing (a * b) / (a * b) procudes undefined results + // because the three cases always produce undefined calculations: + // - 0 / 0 + // - val * NaN + // - inf / inf + if (b == 0 or std.math.isNan(b) or expA == 0x7fff) { + return (a * b) / (a * b); + } + + // Remove the sign from both + aPtr_u16[exp_and_sign_index] = @bitCast(u16, @intCast(i16, expA)); + bPtr_u16[exp_and_sign_index] = @bitCast(u16, @intCast(i16, expB)); + if (amod <= bmod) { + if (amod == bmod) { + return 0 * a; + } + return a; + } + + if (expA == 0) { + amod *= 0x1p120; + expA = @as(i32, aPtr_u16[exp_and_sign_index]) - 120; + } + + if (expB == 0) { + bmod *= 0x1p120; + expB = @as(i32, bPtr_u16[exp_and_sign_index]) - 120; + } + + // OR in extra non-stored mantissa digit + var highA: u64 = (aPtr_u64[high_index] & (std.math.maxInt(u64) >> 16)) | 1 << 48; + var highB: u64 = (bPtr_u64[high_index] & (std.math.maxInt(u64) >> 16)) | 1 << 48; + var lowA: u64 = aPtr_u64[low_index]; + var lowB: u64 = bPtr_u64[low_index]; + + while (expA > expB) : (expA -= 1) { + var high = highA -% highB; + var low = lowA -% lowB; + if (lowA < lowB) { + high -%= 1; + } + if (high >> 63 == 0) { + if ((high | low) == 0) { + return 0 * a; + } + highA = 2 *% high + (low >> 63); + lowA = 2 *% low; + } else { + highA = 2 *% highA + (lowA >> 63); + lowA = 2 *% lowA; + } + } + + var high = highA -% highB; + var low = lowA -% lowB; + if (lowA < lowB) { + high -= 1; + } + if (high >> 63 == 0) { + if ((high | low) == 0) { + return 0 * a; + } + highA = high; + lowA = low; + } + + while (highA >> 48 == 0) { + highA = 2 *% highA + (lowA >> 63); + lowA = 2 *% lowA; + expA = expA - 1; + } + + // Overwrite the current amod with the values in highA and lowA + aPtr_u64[high_index] = highA; + aPtr_u64[low_index] = lowA; + + // Combine the exponent with the sign, normalize if happend to be denormalized + if (expA <= 0) { + aPtr_u16[exp_and_sign_index] = @truncate(u16, @bitCast(u32, (expA +% 120))) | signA; + amod *= 0x1p-120; + } else { + aPtr_u16[exp_and_sign_index] = @truncate(u16, @bitCast(u32, expA)) | signA; + } + + return amod; +} + +inline fn generic_fmod(comptime T: type, x: T, y: T) T { + @setRuntimeSafety(false); + + const bits = @typeInfo(T).Float.bits; + const uint = std.meta.Int(.unsigned, bits); + const log2uint = math.Log2Int(uint); + comptime assert(T == f32 or T == f64); + const digits = if (T == f32) 23 else 52; + const exp_bits = if (T == f32) 9 else 12; + const bits_minus_1 = bits - 1; + const mask = if (T == f32) 0xff else 0x7ff; + var ux = @bitCast(uint, x); + var uy = @bitCast(uint, y); + var ex = @intCast(i32, (ux >> digits) & mask); + var ey = @intCast(i32, (uy >> digits) & mask); + const sx = if (T == f32) @intCast(u32, ux & 0x80000000) else @intCast(i32, ux >> bits_minus_1); + var i: uint = undefined; + + if (uy << 1 == 0 or math.isNan(@bitCast(T, uy)) or ex == mask) + return (x * y) / (x * y); + + if (ux << 1 <= uy << 1) { + if (ux << 1 == uy << 1) + return 0 * x; + return x; + } + + // normalize x and y + if (ex == 0) { + i = ux << exp_bits; + while (i >> bits_minus_1 == 0) : ({ + ex -= 1; + i <<= 1; + }) {} + ux <<= @intCast(log2uint, @bitCast(u32, -ex + 1)); + } else { + ux &= math.maxInt(uint) >> exp_bits; + ux |= 1 << digits; + } + if (ey == 0) { + i = uy << exp_bits; + while (i >> bits_minus_1 == 0) : ({ + ey -= 1; + i <<= 1; + }) {} + uy <<= @intCast(log2uint, @bitCast(u32, -ey + 1)); + } else { + uy &= math.maxInt(uint) >> exp_bits; + uy |= 1 << digits; + } + + // x mod y + while (ex > ey) : (ex -= 1) { + i = ux -% uy; + if (i >> bits_minus_1 == 0) { + if (i == 0) + return 0 * x; + ux = i; + } + ux <<= 1; + } + i = ux -% uy; + if (i >> bits_minus_1 == 0) { + if (i == 0) + return 0 * x; + ux = i; + } + while (ux >> digits == 0) : ({ + ux <<= 1; + ex -= 1; + }) {} + + // scale result up + if (ex > 0) { + ux -%= 1 << digits; + ux |= @as(uint, @bitCast(u32, ex)) << digits; + } else { + ux >>= @intCast(log2uint, @bitCast(u32, -ex + 1)); + } + if (T == f32) { + ux |= sx; + } else { + ux |= @intCast(uint, sx) << bits_minus_1; + } + return @bitCast(T, ux); +} + +test "fmod, fmodf" { + inline for ([_]type{ f32, f64 }) |T| { + const nan_val = math.nan(T); + const inf_val = math.inf(T); + + try std.testing.expect(math.isNan(generic_fmod(T, nan_val, 1.0))); + try std.testing.expect(math.isNan(generic_fmod(T, 1.0, nan_val))); + try std.testing.expect(math.isNan(generic_fmod(T, inf_val, 1.0))); + try std.testing.expect(math.isNan(generic_fmod(T, 0.0, 0.0))); + try std.testing.expect(math.isNan(generic_fmod(T, 1.0, 0.0))); + + try std.testing.expectEqual(@as(T, 0.0), generic_fmod(T, 0.0, 2.0)); + try std.testing.expectEqual(@as(T, -0.0), generic_fmod(T, -0.0, 2.0)); + + try std.testing.expectEqual(@as(T, -2.0), generic_fmod(T, -32.0, 10.0)); + try std.testing.expectEqual(@as(T, -2.0), generic_fmod(T, -32.0, -10.0)); + try std.testing.expectEqual(@as(T, 2.0), generic_fmod(T, 32.0, 10.0)); + try std.testing.expectEqual(@as(T, 2.0), generic_fmod(T, 32.0, -10.0)); + } +} + +test { + _ = @import("fmodq_test.zig"); + _ = @import("fmodx_test.zig"); +} diff --git a/lib/std/special/compiler_rt/fmodq.zig b/lib/std/special/compiler_rt/fmodq.zig deleted file mode 100644 index 3f56c49796..0000000000 --- a/lib/std/special/compiler_rt/fmodq.zig +++ /dev/null @@ -1,126 +0,0 @@ -const builtin = @import("builtin"); -const std = @import("std"); - -// fmodq - floating modulo large, returns the remainder of division for f128 types -// Logic and flow heavily inspired by MUSL fmodl for 113 mantissa digits -pub fn fmodq(a: f128, b: f128) callconv(.C) f128 { - @setRuntimeSafety(builtin.is_test); - var amod = a; - var bmod = b; - const aPtr_u64 = @ptrCast([*]u64, &amod); - const bPtr_u64 = @ptrCast([*]u64, &bmod); - const aPtr_u16 = @ptrCast([*]u16, &amod); - const bPtr_u16 = @ptrCast([*]u16, &bmod); - - const exp_and_sign_index = comptime switch (builtin.target.cpu.arch.endian()) { - .Little => 7, - .Big => 0, - }; - const low_index = comptime switch (builtin.target.cpu.arch.endian()) { - .Little => 0, - .Big => 1, - }; - const high_index = comptime switch (builtin.target.cpu.arch.endian()) { - .Little => 1, - .Big => 0, - }; - - const signA = aPtr_u16[exp_and_sign_index] & 0x8000; - var expA = @intCast(i32, (aPtr_u16[exp_and_sign_index] & 0x7fff)); - var expB = @intCast(i32, (bPtr_u16[exp_and_sign_index] & 0x7fff)); - - // There are 3 cases where the answer is undefined, check for: - // - fmodq(val, 0) - // - fmodq(val, NaN) - // - fmodq(inf, val) - // The sign on checked values does not matter. - // Doing (a * b) / (a * b) procudes undefined results - // because the three cases always produce undefined calculations: - // - 0 / 0 - // - val * NaN - // - inf / inf - if (b == 0 or std.math.isNan(b) or expA == 0x7fff) { - return (a * b) / (a * b); - } - - // Remove the sign from both - aPtr_u16[exp_and_sign_index] = @bitCast(u16, @intCast(i16, expA)); - bPtr_u16[exp_and_sign_index] = @bitCast(u16, @intCast(i16, expB)); - if (amod <= bmod) { - if (amod == bmod) { - return 0 * a; - } - return a; - } - - if (expA == 0) { - amod *= 0x1p120; - expA = @as(i32, aPtr_u16[exp_and_sign_index]) - 120; - } - - if (expB == 0) { - bmod *= 0x1p120; - expB = @as(i32, bPtr_u16[exp_and_sign_index]) - 120; - } - - // OR in extra non-stored mantissa digit - var highA: u64 = (aPtr_u64[high_index] & (std.math.maxInt(u64) >> 16)) | 1 << 48; - var highB: u64 = (bPtr_u64[high_index] & (std.math.maxInt(u64) >> 16)) | 1 << 48; - var lowA: u64 = aPtr_u64[low_index]; - var lowB: u64 = bPtr_u64[low_index]; - - while (expA > expB) : (expA -= 1) { - var high = highA -% highB; - var low = lowA -% lowB; - if (lowA < lowB) { - high -%= 1; - } - if (high >> 63 == 0) { - if ((high | low) == 0) { - return 0 * a; - } - highA = 2 *% high + (low >> 63); - lowA = 2 *% low; - } else { - highA = 2 *% highA + (lowA >> 63); - lowA = 2 *% lowA; - } - } - - var high = highA -% highB; - var low = lowA -% lowB; - if (lowA < lowB) { - high -= 1; - } - if (high >> 63 == 0) { - if ((high | low) == 0) { - return 0 * a; - } - highA = high; - lowA = low; - } - - while (highA >> 48 == 0) { - highA = 2 *% highA + (lowA >> 63); - lowA = 2 *% lowA; - expA = expA - 1; - } - - // Overwrite the current amod with the values in highA and lowA - aPtr_u64[high_index] = highA; - aPtr_u64[low_index] = lowA; - - // Combine the exponent with the sign, normalize if happend to be denormalized - if (expA <= 0) { - aPtr_u16[exp_and_sign_index] = @truncate(u16, @bitCast(u32, (expA +% 120))) | signA; - amod *= 0x1p-120; - } else { - aPtr_u16[exp_and_sign_index] = @truncate(u16, @bitCast(u32, expA)) | signA; - } - - return amod; -} - -test { - _ = @import("fmodq_test.zig"); -} diff --git a/lib/std/special/compiler_rt/fmodq_test.zig b/lib/std/special/compiler_rt/fmodq_test.zig index b8baf8ae9b..07ddb8d182 100644 --- a/lib/std/special/compiler_rt/fmodq_test.zig +++ b/lib/std/special/compiler_rt/fmodq_test.zig @@ -1,24 +1,24 @@ const std = @import("std"); -const fmodq = @import("fmodq.zig"); +const fmod = @import("fmod.zig"); const testing = std.testing; fn test_fmodq(a: f128, b: f128, exp: f128) !void { - const res = fmodq.fmodq(a, b); + const res = fmod.fmodq(a, b); try testing.expect(exp == res); } fn test_fmodq_nans() !void { - try testing.expect(std.math.isNan(fmodq.fmodq(1.0, std.math.nan(f128)))); - try testing.expect(std.math.isNan(fmodq.fmodq(1.0, -std.math.nan(f128)))); - try testing.expect(std.math.isNan(fmodq.fmodq(std.math.nan(f128), 1.0))); - try testing.expect(std.math.isNan(fmodq.fmodq(-std.math.nan(f128), 1.0))); + try testing.expect(std.math.isNan(fmod.fmodq(1.0, std.math.nan(f128)))); + try testing.expect(std.math.isNan(fmod.fmodq(1.0, -std.math.nan(f128)))); + try testing.expect(std.math.isNan(fmod.fmodq(std.math.nan(f128), 1.0))); + try testing.expect(std.math.isNan(fmod.fmodq(-std.math.nan(f128), 1.0))); } fn test_fmodq_infs() !void { - try testing.expect(fmodq.fmodq(1.0, std.math.inf(f128)) == 1.0); - try testing.expect(fmodq.fmodq(1.0, -std.math.inf(f128)) == 1.0); - try testing.expect(std.math.isNan(fmodq.fmodq(std.math.inf(f128), 1.0))); - try testing.expect(std.math.isNan(fmodq.fmodq(-std.math.inf(f128), 1.0))); + try testing.expect(fmod.fmodq(1.0, std.math.inf(f128)) == 1.0); + try testing.expect(fmod.fmodq(1.0, -std.math.inf(f128)) == 1.0); + try testing.expect(std.math.isNan(fmod.fmodq(std.math.inf(f128), 1.0))); + try testing.expect(std.math.isNan(fmod.fmodq(-std.math.inf(f128), 1.0))); } test "fmodq" { diff --git a/lib/std/special/compiler_rt/fmodx.zig b/lib/std/special/compiler_rt/fmodx.zig deleted file mode 100644 index efe16f9f16..0000000000 --- a/lib/std/special/compiler_rt/fmodx.zig +++ /dev/null @@ -1,108 +0,0 @@ -const builtin = @import("builtin"); -const std = @import("std"); -const math = std.math; -const normalize = @import("divdf3.zig").normalize; - -// fmodx - floating modulo large, returns the remainder of division for f80 types -// Logic and flow heavily inspired by MUSL fmodl for 113 mantissa digits -pub fn fmodx(a: f80, b: f80) callconv(.C) f80 { - @setRuntimeSafety(builtin.is_test); - - const T = f80; - const Z = std.meta.Int(.unsigned, @bitSizeOf(T)); - - const significandBits = math.floatMantissaBits(T); - const fractionalBits = math.floatFractionalBits(T); - const exponentBits = math.floatExponentBits(T); - - const signBit = (@as(Z, 1) << (significandBits + exponentBits)); - const maxExponent = ((1 << exponentBits) - 1); - - var aRep = @bitCast(Z, a); - var bRep = @bitCast(Z, b); - - const signA = aRep & signBit; - var expA = @intCast(i32, (@bitCast(Z, a) >> significandBits) & maxExponent); - var expB = @intCast(i32, (@bitCast(Z, b) >> significandBits) & maxExponent); - - // There are 3 cases where the answer is undefined, check for: - // - fmodx(val, 0) - // - fmodx(val, NaN) - // - fmodx(inf, val) - // The sign on checked values does not matter. - // Doing (a * b) / (a * b) procudes undefined results - // because the three cases always produce undefined calculations: - // - 0 / 0 - // - val * NaN - // - inf / inf - if (b == 0 or math.isNan(b) or expA == maxExponent) { - return (a * b) / (a * b); - } - - // Remove the sign from both - aRep &= ~signBit; - bRep &= ~signBit; - if (aRep <= bRep) { - if (aRep == bRep) { - return 0 * a; - } - return a; - } - - if (expA == 0) expA = normalize(f80, &aRep); - if (expB == 0) expB = normalize(f80, &bRep); - - var highA: u64 = 0; - var highB: u64 = 0; - var lowA: u64 = @truncate(u64, aRep); - var lowB: u64 = @truncate(u64, bRep); - - while (expA > expB) : (expA -= 1) { - var high = highA -% highB; - var low = lowA -% lowB; - if (lowA < lowB) { - high -%= 1; - } - if (high >> 63 == 0) { - if ((high | low) == 0) { - return 0 * a; - } - highA = 2 *% high + (low >> 63); - lowA = 2 *% low; - } else { - highA = 2 *% highA + (lowA >> 63); - lowA = 2 *% lowA; - } - } - - var high = highA -% highB; - var low = lowA -% lowB; - if (lowA < lowB) { - high -%= 1; - } - if (high >> 63 == 0) { - if ((high | low) == 0) { - return 0 * a; - } - highA = high; - lowA = low; - } - - while ((lowA >> fractionalBits) == 0) { - lowA = 2 *% lowA; - expA = expA - 1; - } - - // Combine the exponent with the sign and significand, normalize if happened to be denormalized - if (expA < -fractionalBits) { - return @bitCast(T, signA); - } else if (expA <= 0) { - return @bitCast(T, (lowA >> @intCast(math.Log2Int(u64), 1 - expA)) | signA); - } else { - return @bitCast(T, lowA | (@as(Z, @intCast(u16, expA)) << significandBits) | signA); - } -} - -test { - _ = @import("fmodx_test.zig"); -} diff --git a/lib/std/special/compiler_rt/fmodx_test.zig b/lib/std/special/compiler_rt/fmodx_test.zig index a5d0887ea4..4bb1b5654a 100644 --- a/lib/std/special/compiler_rt/fmodx_test.zig +++ b/lib/std/special/compiler_rt/fmodx_test.zig @@ -1,24 +1,24 @@ const std = @import("std"); -const fmodx = @import("fmodx.zig"); +const fmod = @import("fmod.zig"); const testing = std.testing; fn test_fmodx(a: f80, b: f80, exp: f80) !void { - const res = fmodx.fmodx(a, b); + const res = fmod.__fmodx(a, b); try testing.expect(exp == res); } fn test_fmodx_nans() !void { - try testing.expect(std.math.isNan(fmodx.fmodx(1.0, std.math.nan(f80)))); - try testing.expect(std.math.isNan(fmodx.fmodx(1.0, -std.math.nan(f80)))); - try testing.expect(std.math.isNan(fmodx.fmodx(std.math.nan(f80), 1.0))); - try testing.expect(std.math.isNan(fmodx.fmodx(-std.math.nan(f80), 1.0))); + try testing.expect(std.math.isNan(fmod.__fmodx(1.0, std.math.nan(f80)))); + try testing.expect(std.math.isNan(fmod.__fmodx(1.0, -std.math.nan(f80)))); + try testing.expect(std.math.isNan(fmod.__fmodx(std.math.nan(f80), 1.0))); + try testing.expect(std.math.isNan(fmod.__fmodx(-std.math.nan(f80), 1.0))); } fn test_fmodx_infs() !void { - try testing.expect(fmodx.fmodx(1.0, std.math.inf(f80)) == 1.0); - try testing.expect(fmodx.fmodx(1.0, -std.math.inf(f80)) == 1.0); - try testing.expect(std.math.isNan(fmodx.fmodx(std.math.inf(f80), 1.0))); - try testing.expect(std.math.isNan(fmodx.fmodx(-std.math.inf(f80), 1.0))); + try testing.expect(fmod.__fmodx(1.0, std.math.inf(f80)) == 1.0); + try testing.expect(fmod.__fmodx(1.0, -std.math.inf(f80)) == 1.0); + try testing.expect(std.math.isNan(fmod.__fmodx(std.math.inf(f80), 1.0))); + try testing.expect(std.math.isNan(fmod.__fmodx(-std.math.inf(f80), 1.0))); } test "fmodx" { diff --git a/lib/std/special/compiler_rt/log.zig b/lib/std/special/compiler_rt/log.zig new file mode 100644 index 0000000000..8b09baac2e --- /dev/null +++ b/lib/std/special/compiler_rt/log.zig @@ -0,0 +1,168 @@ +// Ported from musl, which is licensed under the MIT license: +// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT +// +// https://git.musl-libc.org/cgit/musl/tree/src/math/lnf.c +// https://git.musl-libc.org/cgit/musl/tree/src/math/ln.c + +const std = @import("std"); +const math = std.math; +const testing = std.testing; + +pub fn __logh(a: f16) callconv(.C) f16 { + // TODO: more efficient implementation + return @floatCast(f16, logf(a)); +} + +pub fn logf(x_: f32) callconv(.C) f32 { + const ln2_hi: f32 = 6.9313812256e-01; + const ln2_lo: f32 = 9.0580006145e-06; + const Lg1: f32 = 0xaaaaaa.0p-24; + const Lg2: f32 = 0xccce13.0p-25; + const Lg3: f32 = 0x91e9ee.0p-25; + const Lg4: f32 = 0xf89e26.0p-26; + + var x = x_; + var ix = @bitCast(u32, x); + var k: i32 = 0; + + // x < 2^(-126) + if (ix < 0x00800000 or ix >> 31 != 0) { + // log(+-0) = -inf + if (ix << 1 == 0) { + return -math.inf(f32); + } + // log(-#) = nan + if (ix >> 31 != 0) { + return math.nan(f32); + } + + // subnormal, scale x + k -= 25; + x *= 0x1.0p25; + ix = @bitCast(u32, x); + } else if (ix >= 0x7F800000) { + return x; + } else if (ix == 0x3F800000) { + return 0; + } + + // x into [sqrt(2) / 2, sqrt(2)] + ix += 0x3F800000 - 0x3F3504F3; + k += @intCast(i32, ix >> 23) - 0x7F; + ix = (ix & 0x007FFFFF) + 0x3F3504F3; + x = @bitCast(f32, ix); + + const f = x - 1.0; + const s = f / (2.0 + f); + const z = s * s; + const w = z * z; + const t1 = w * (Lg2 + w * Lg4); + const t2 = z * (Lg1 + w * Lg3); + const R = t2 + t1; + const hfsq = 0.5 * f * f; + const dk = @intToFloat(f32, k); + + return s * (hfsq + R) + dk * ln2_lo - hfsq + f + dk * ln2_hi; +} + +pub fn log(x_: f64) callconv(.C) f64 { + const ln2_hi: f64 = 6.93147180369123816490e-01; + const ln2_lo: f64 = 1.90821492927058770002e-10; + const Lg1: f64 = 6.666666666666735130e-01; + const Lg2: f64 = 3.999999999940941908e-01; + const Lg3: f64 = 2.857142874366239149e-01; + const Lg4: f64 = 2.222219843214978396e-01; + const Lg5: f64 = 1.818357216161805012e-01; + const Lg6: f64 = 1.531383769920937332e-01; + const Lg7: f64 = 1.479819860511658591e-01; + + var x = x_; + var ix = @bitCast(u64, x); + var hx = @intCast(u32, ix >> 32); + var k: i32 = 0; + + if (hx < 0x00100000 or hx >> 31 != 0) { + // log(+-0) = -inf + if (ix << 1 == 0) { + return -math.inf(f64); + } + // log(-#) = nan + if (hx >> 31 != 0) { + return math.nan(f64); + } + + // subnormal, scale x + k -= 54; + x *= 0x1.0p54; + hx = @intCast(u32, @bitCast(u64, ix) >> 32); + } else if (hx >= 0x7FF00000) { + return x; + } else if (hx == 0x3FF00000 and ix << 32 == 0) { + return 0; + } + + // x into [sqrt(2) / 2, sqrt(2)] + hx += 0x3FF00000 - 0x3FE6A09E; + k += @intCast(i32, hx >> 20) - 0x3FF; + hx = (hx & 0x000FFFFF) + 0x3FE6A09E; + ix = (@as(u64, hx) << 32) | (ix & 0xFFFFFFFF); + x = @bitCast(f64, ix); + + const f = x - 1.0; + const hfsq = 0.5 * f * f; + const s = f / (2.0 + f); + const z = s * s; + const w = z * z; + const t1 = w * (Lg2 + w * (Lg4 + w * Lg6)); + const t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7))); + const R = t2 + t1; + const dk = @intToFloat(f64, k); + + return s * (hfsq + R) + dk * ln2_lo - hfsq + f + dk * ln2_hi; +} + +pub fn __logx(a: f80) callconv(.C) f80 { + // TODO: more efficient implementation + return @floatCast(f80, logq(a)); +} + +pub fn logq(a: f128) callconv(.C) f128 { + // TODO: more correct implementation + return log(@floatCast(f64, a)); +} + +test "ln32" { + const epsilon = 0.000001; + + try testing.expect(math.approxEqAbs(f32, logf(0.2), -1.609438, epsilon)); + try testing.expect(math.approxEqAbs(f32, logf(0.8923), -0.113953, epsilon)); + try testing.expect(math.approxEqAbs(f32, logf(1.5), 0.405465, epsilon)); + try testing.expect(math.approxEqAbs(f32, logf(37.45), 3.623007, epsilon)); + try testing.expect(math.approxEqAbs(f32, logf(89.123), 4.490017, epsilon)); + try testing.expect(math.approxEqAbs(f32, logf(123123.234375), 11.720941, epsilon)); +} + +test "ln64" { + const epsilon = 0.000001; + + try testing.expect(math.approxEqAbs(f64, log(0.2), -1.609438, epsilon)); + try testing.expect(math.approxEqAbs(f64, log(0.8923), -0.113953, epsilon)); + try testing.expect(math.approxEqAbs(f64, log(1.5), 0.405465, epsilon)); + try testing.expect(math.approxEqAbs(f64, log(37.45), 3.623007, epsilon)); + try testing.expect(math.approxEqAbs(f64, log(89.123), 4.490017, epsilon)); + try testing.expect(math.approxEqAbs(f64, log(123123.234375), 11.720941, epsilon)); +} + +test "ln32.special" { + try testing.expect(math.isPositiveInf(logf(math.inf(f32)))); + try testing.expect(math.isNegativeInf(logf(0.0))); + try testing.expect(math.isNan(logf(-1.0))); + try testing.expect(math.isNan(logf(math.nan(f32)))); +} + +test "ln64.special" { + try testing.expect(math.isPositiveInf(log(math.inf(f64)))); + try testing.expect(math.isNegativeInf(log(0.0))); + try testing.expect(math.isNan(log(-1.0))); + try testing.expect(math.isNan(log(math.nan(f64)))); +} diff --git a/lib/std/special/compiler_rt/log10.zig b/lib/std/special/compiler_rt/log10.zig new file mode 100644 index 0000000000..ce06d8c649 --- /dev/null +++ b/lib/std/special/compiler_rt/log10.zig @@ -0,0 +1,196 @@ +// Ported from musl, which is licensed under the MIT license: +// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT +// +// https://git.musl-libc.org/cgit/musl/tree/src/math/log10f.c +// https://git.musl-libc.org/cgit/musl/tree/src/math/log10.c + +const std = @import("std"); +const math = std.math; +const testing = std.testing; +const maxInt = std.math.maxInt; + +pub fn __log10h(a: f16) callconv(.C) f16 { + // TODO: more efficient implementation + return @floatCast(f16, log10f(a)); +} + +pub fn log10f(x_: f32) callconv(.C) f32 { + const ivln10hi: f32 = 4.3432617188e-01; + const ivln10lo: f32 = -3.1689971365e-05; + const log10_2hi: f32 = 3.0102920532e-01; + const log10_2lo: f32 = 7.9034151668e-07; + const Lg1: f32 = 0xaaaaaa.0p-24; + const Lg2: f32 = 0xccce13.0p-25; + const Lg3: f32 = 0x91e9ee.0p-25; + const Lg4: f32 = 0xf89e26.0p-26; + + var x = x_; + var u = @bitCast(u32, x); + var ix = u; + var k: i32 = 0; + + // x < 2^(-126) + if (ix < 0x00800000 or ix >> 31 != 0) { + // log(+-0) = -inf + if (ix << 1 == 0) { + return -math.inf(f32); + } + // log(-#) = nan + if (ix >> 31 != 0) { + return math.nan(f32); + } + + k -= 25; + x *= 0x1.0p25; + ix = @bitCast(u32, x); + } else if (ix >= 0x7F800000) { + return x; + } else if (ix == 0x3F800000) { + return 0; + } + + // x into [sqrt(2) / 2, sqrt(2)] + ix += 0x3F800000 - 0x3F3504F3; + k += @intCast(i32, ix >> 23) - 0x7F; + ix = (ix & 0x007FFFFF) + 0x3F3504F3; + x = @bitCast(f32, ix); + + const f = x - 1.0; + const s = f / (2.0 + f); + const z = s * s; + const w = z * z; + const t1 = w * (Lg2 + w * Lg4); + const t2 = z * (Lg1 + w * Lg3); + const R = t2 + t1; + const hfsq = 0.5 * f * f; + + var hi = f - hfsq; + u = @bitCast(u32, hi); + u &= 0xFFFFF000; + hi = @bitCast(f32, u); + const lo = f - hi - hfsq + s * (hfsq + R); + const dk = @intToFloat(f32, k); + + return dk * log10_2lo + (lo + hi) * ivln10lo + lo * ivln10hi + hi * ivln10hi + dk * log10_2hi; +} + +pub fn log10(x_: f64) callconv(.C) f64 { + const ivln10hi: f64 = 4.34294481878168880939e-01; + const ivln10lo: f64 = 2.50829467116452752298e-11; + const log10_2hi: f64 = 3.01029995663611771306e-01; + const log10_2lo: f64 = 3.69423907715893078616e-13; + const Lg1: f64 = 6.666666666666735130e-01; + const Lg2: f64 = 3.999999999940941908e-01; + const Lg3: f64 = 2.857142874366239149e-01; + const Lg4: f64 = 2.222219843214978396e-01; + const Lg5: f64 = 1.818357216161805012e-01; + const Lg6: f64 = 1.531383769920937332e-01; + const Lg7: f64 = 1.479819860511658591e-01; + + var x = x_; + var ix = @bitCast(u64, x); + var hx = @intCast(u32, ix >> 32); + var k: i32 = 0; + + if (hx < 0x00100000 or hx >> 31 != 0) { + // log(+-0) = -inf + if (ix << 1 == 0) { + return -math.inf(f32); + } + // log(-#) = nan + if (hx >> 31 != 0) { + return math.nan(f32); + } + + // subnormal, scale x + k -= 54; + x *= 0x1.0p54; + hx = @intCast(u32, @bitCast(u64, x) >> 32); + } else if (hx >= 0x7FF00000) { + return x; + } else if (hx == 0x3FF00000 and ix << 32 == 0) { + return 0; + } + + // x into [sqrt(2) / 2, sqrt(2)] + hx += 0x3FF00000 - 0x3FE6A09E; + k += @intCast(i32, hx >> 20) - 0x3FF; + hx = (hx & 0x000FFFFF) + 0x3FE6A09E; + ix = (@as(u64, hx) << 32) | (ix & 0xFFFFFFFF); + x = @bitCast(f64, ix); + + const f = x - 1.0; + const hfsq = 0.5 * f * f; + const s = f / (2.0 + f); + const z = s * s; + const w = z * z; + const t1 = w * (Lg2 + w * (Lg4 + w * Lg6)); + const t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7))); + const R = t2 + t1; + + // hi + lo = f - hfsq + s * (hfsq + R) ~ log(1 + f) + var hi = f - hfsq; + var hii = @bitCast(u64, hi); + hii &= @as(u64, maxInt(u64)) << 32; + hi = @bitCast(f64, hii); + const lo = f - hi - hfsq + s * (hfsq + R); + + // val_hi + val_lo ~ log10(1 + f) + k * log10(2) + var val_hi = hi * ivln10hi; + const dk = @intToFloat(f64, k); + const y = dk * log10_2hi; + var val_lo = dk * log10_2lo + (lo + hi) * ivln10lo + lo * ivln10hi; + + // Extra precision multiplication + const ww = y + val_hi; + val_lo += (y - ww) + val_hi; + val_hi = ww; + + return val_lo + val_hi; +} + +pub fn __log10x(a: f80) callconv(.C) f80 { + // TODO: more efficient implementation + return @floatCast(f80, log10q(a)); +} + +pub fn log10q(a: f128) callconv(.C) f128 { + // TODO: more correct implementation + return log10(@floatCast(f64, a)); +} + +test "log10_32" { + const epsilon = 0.000001; + + try testing.expect(math.approxEqAbs(f32, log10f(0.2), -0.698970, epsilon)); + try testing.expect(math.approxEqAbs(f32, log10f(0.8923), -0.049489, epsilon)); + try testing.expect(math.approxEqAbs(f32, log10f(1.5), 0.176091, epsilon)); + try testing.expect(math.approxEqAbs(f32, log10f(37.45), 1.573452, epsilon)); + try testing.expect(math.approxEqAbs(f32, log10f(89.123), 1.94999, epsilon)); + try testing.expect(math.approxEqAbs(f32, log10f(123123.234375), 5.09034, epsilon)); +} + +test "log10_64" { + const epsilon = 0.000001; + + try testing.expect(math.approxEqAbs(f64, log10(0.2), -0.698970, epsilon)); + try testing.expect(math.approxEqAbs(f64, log10(0.8923), -0.049489, epsilon)); + try testing.expect(math.approxEqAbs(f64, log10(1.5), 0.176091, epsilon)); + try testing.expect(math.approxEqAbs(f64, log10(37.45), 1.573452, epsilon)); + try testing.expect(math.approxEqAbs(f64, log10(89.123), 1.94999, epsilon)); + try testing.expect(math.approxEqAbs(f64, log10(123123.234375), 5.09034, epsilon)); +} + +test "log10_32.special" { + try testing.expect(math.isPositiveInf(log10f(math.inf(f32)))); + try testing.expect(math.isNegativeInf(log10f(0.0))); + try testing.expect(math.isNan(log10f(-1.0))); + try testing.expect(math.isNan(log10f(math.nan(f32)))); +} + +test "log10_64.special" { + try testing.expect(math.isPositiveInf(log10(math.inf(f64)))); + try testing.expect(math.isNegativeInf(log10(0.0))); + try testing.expect(math.isNan(log10(-1.0))); + try testing.expect(math.isNan(log10(math.nan(f64)))); +} diff --git a/lib/std/special/compiler_rt/log2.zig b/lib/std/special/compiler_rt/log2.zig new file mode 100644 index 0000000000..2c2d620c3d --- /dev/null +++ b/lib/std/special/compiler_rt/log2.zig @@ -0,0 +1,185 @@ +// Ported from musl, which is licensed under the MIT license: +// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT +// +// https://git.musl-libc.org/cgit/musl/tree/src/math/log2f.c +// https://git.musl-libc.org/cgit/musl/tree/src/math/log2.c + +const std = @import("std"); +const math = std.math; +const expect = std.testing.expect; +const maxInt = std.math.maxInt; + +pub fn __log2h(a: f16) callconv(.C) f16 { + // TODO: more efficient implementation + return @floatCast(f16, log2f(a)); +} + +pub fn log2f(x_: f32) callconv(.C) f32 { + const ivln2hi: f32 = 1.4428710938e+00; + const ivln2lo: f32 = -1.7605285393e-04; + const Lg1: f32 = 0xaaaaaa.0p-24; + const Lg2: f32 = 0xccce13.0p-25; + const Lg3: f32 = 0x91e9ee.0p-25; + const Lg4: f32 = 0xf89e26.0p-26; + + var x = x_; + var u = @bitCast(u32, x); + var ix = u; + var k: i32 = 0; + + // x < 2^(-126) + if (ix < 0x00800000 or ix >> 31 != 0) { + // log(+-0) = -inf + if (ix << 1 == 0) { + return -math.inf(f32); + } + // log(-#) = nan + if (ix >> 31 != 0) { + return math.nan(f32); + } + + k -= 25; + x *= 0x1.0p25; + ix = @bitCast(u32, x); + } else if (ix >= 0x7F800000) { + return x; + } else if (ix == 0x3F800000) { + return 0; + } + + // x into [sqrt(2) / 2, sqrt(2)] + ix += 0x3F800000 - 0x3F3504F3; + k += @intCast(i32, ix >> 23) - 0x7F; + ix = (ix & 0x007FFFFF) + 0x3F3504F3; + x = @bitCast(f32, ix); + + const f = x - 1.0; + const s = f / (2.0 + f); + const z = s * s; + const w = z * z; + const t1 = w * (Lg2 + w * Lg4); + const t2 = z * (Lg1 + w * Lg3); + const R = t2 + t1; + const hfsq = 0.5 * f * f; + + var hi = f - hfsq; + u = @bitCast(u32, hi); + u &= 0xFFFFF000; + hi = @bitCast(f32, u); + const lo = f - hi - hfsq + s * (hfsq + R); + return (lo + hi) * ivln2lo + lo * ivln2hi + hi * ivln2hi + @intToFloat(f32, k); +} + +pub fn log2(x_: f64) callconv(.C) f64 { + const ivln2hi: f64 = 1.44269504072144627571e+00; + const ivln2lo: f64 = 1.67517131648865118353e-10; + const Lg1: f64 = 6.666666666666735130e-01; + const Lg2: f64 = 3.999999999940941908e-01; + const Lg3: f64 = 2.857142874366239149e-01; + const Lg4: f64 = 2.222219843214978396e-01; + const Lg5: f64 = 1.818357216161805012e-01; + const Lg6: f64 = 1.531383769920937332e-01; + const Lg7: f64 = 1.479819860511658591e-01; + + var x = x_; + var ix = @bitCast(u64, x); + var hx = @intCast(u32, ix >> 32); + var k: i32 = 0; + + if (hx < 0x00100000 or hx >> 31 != 0) { + // log(+-0) = -inf + if (ix << 1 == 0) { + return -math.inf(f64); + } + // log(-#) = nan + if (hx >> 31 != 0) { + return math.nan(f64); + } + + // subnormal, scale x + k -= 54; + x *= 0x1.0p54; + hx = @intCast(u32, @bitCast(u64, x) >> 32); + } else if (hx >= 0x7FF00000) { + return x; + } else if (hx == 0x3FF00000 and ix << 32 == 0) { + return 0; + } + + // x into [sqrt(2) / 2, sqrt(2)] + hx += 0x3FF00000 - 0x3FE6A09E; + k += @intCast(i32, hx >> 20) - 0x3FF; + hx = (hx & 0x000FFFFF) + 0x3FE6A09E; + ix = (@as(u64, hx) << 32) | (ix & 0xFFFFFFFF); + x = @bitCast(f64, ix); + + const f = x - 1.0; + const hfsq = 0.5 * f * f; + const s = f / (2.0 + f); + const z = s * s; + const w = z * z; + const t1 = w * (Lg2 + w * (Lg4 + w * Lg6)); + const t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7))); + const R = t2 + t1; + + // hi + lo = f - hfsq + s * (hfsq + R) ~ log(1 + f) + var hi = f - hfsq; + var hii = @bitCast(u64, hi); + hii &= @as(u64, maxInt(u64)) << 32; + hi = @bitCast(f64, hii); + const lo = f - hi - hfsq + s * (hfsq + R); + + var val_hi = hi * ivln2hi; + var val_lo = (lo + hi) * ivln2lo + lo * ivln2hi; + + // spadd(val_hi, val_lo, y) + const y = @intToFloat(f64, k); + const ww = y + val_hi; + val_lo += (y - ww) + val_hi; + val_hi = ww; + + return val_lo + val_hi; +} + +pub fn __log2x(a: f80) callconv(.C) f80 { + // TODO: more efficient implementation + return @floatCast(f80, log2q(a)); +} + +pub fn log2q(a: f128) callconv(.C) f128 { + return math.log2(a); +} + +test "log2_32" { + const epsilon = 0.000001; + + try expect(math.approxEqAbs(f32, log2f(0.2), -2.321928, epsilon)); + try expect(math.approxEqAbs(f32, log2f(0.8923), -0.164399, epsilon)); + try expect(math.approxEqAbs(f32, log2f(1.5), 0.584962, epsilon)); + try expect(math.approxEqAbs(f32, log2f(37.45), 5.226894, epsilon)); + try expect(math.approxEqAbs(f32, log2f(123123.234375), 16.909744, epsilon)); +} + +test "log2_64" { + const epsilon = 0.000001; + + try expect(math.approxEqAbs(f64, log2(0.2), -2.321928, epsilon)); + try expect(math.approxEqAbs(f64, log2(0.8923), -0.164399, epsilon)); + try expect(math.approxEqAbs(f64, log2(1.5), 0.584962, epsilon)); + try expect(math.approxEqAbs(f64, log2(37.45), 5.226894, epsilon)); + try expect(math.approxEqAbs(f64, log2(123123.234375), 16.909744, epsilon)); +} + +test "log2_32.special" { + try expect(math.isPositiveInf(log2f(math.inf(f32)))); + try expect(math.isNegativeInf(log2f(0.0))); + try expect(math.isNan(log2f(-1.0))); + try expect(math.isNan(log2f(math.nan(f32)))); +} + +test "log2_64.special" { + try expect(math.isPositiveInf(log2(math.inf(f64)))); + try expect(math.isNegativeInf(log2(0.0))); + try expect(math.isNan(log2(-1.0))); + try expect(math.isNan(log2(math.nan(f64)))); +} diff --git a/lib/std/math/__rem_pio2.zig b/lib/std/special/compiler_rt/rem_pio2.zig index f01d8fe94a..73d477ee12 100644 --- a/lib/std/math/__rem_pio2.zig +++ b/lib/std/special/compiler_rt/rem_pio2.zig @@ -3,8 +3,8 @@ // // https://git.musl-libc.org/cgit/musl/tree/src/math/__rem_pio2.c -const std = @import("../std.zig"); -const __rem_pio2_large = @import("__rem_pio2_large.zig").__rem_pio2_large; +const std = @import("std"); +const rem_pio2_large = @import("rem_pio2_large.zig").rem_pio2_large; const math = std.math; const toint = 1.5 / math.floatEps(f64); @@ -82,10 +82,10 @@ fn medium(ix: u32, x: f64, y: *[2]f64) i32 { // Returns the remainder of x rem pi/2 in y[0]+y[1] // -// use __rem_pio2_large() for large x +// use rem_pio2_large() for large x // // caller must handle the case when reduction is not needed: |x| ~<= pi/4 */ -pub fn __rem_pio2(x: f64, y: *[2]f64) i32 { +pub fn rem_pio2(x: f64, y: *[2]f64) i32 { var z: f64 = undefined; var tx: [3]f64 = undefined; var ty: [2]f64 = undefined; @@ -186,7 +186,7 @@ pub fn __rem_pio2(x: f64, y: *[2]f64) i32 { while (tx[U(i)] == 0.0) { i -= 1; } - n = __rem_pio2_large(tx[0..], ty[0..], @intCast(i32, (ix >> 20)) - (0x3ff + 23), i + 1, 1); + n = rem_pio2_large(tx[0..], ty[0..], @intCast(i32, (ix >> 20)) - (0x3ff + 23), i + 1, 1); if (sign) { y[0] = -ty[0]; y[1] = -ty[1]; diff --git a/lib/std/math/__rem_pio2_large.zig b/lib/std/special/compiler_rt/rem_pio2_large.zig index 140e85f7f6..c8a53b741c 100644 --- a/lib/std/math/__rem_pio2_large.zig +++ b/lib/std/special/compiler_rt/rem_pio2_large.zig @@ -3,23 +3,22 @@ // // https://git.musl-libc.org/cgit/musl/tree/src/math/__rem_pio2_large.c -const std = @import("../std.zig"); +const std = @import("std"); const math = std.math; const init_jk = [_]i32{ 3, 4, 4, 6 }; // initial value for jk -// -// Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi -// -// integer array, contains the (24*i)-th to (24*i+23)-th -// bit of 2/pi after binary point. The corresponding -// floating value is -// -// ipio2[i] * 2^(-24(i+1)). -// -// NB: This table must have at least (e0-3)/24 + jk terms. -// For quad precision (e0 <= 16360, jk = 6), this is 686. /// +/// Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi +/// +/// integer array, contains the (24*i)-th to (24*i+23)-th +/// bit of 2/pi after binary point. The corresponding +/// floating value is +/// +/// ipio2[i] * 2^(-24(i+1)). +/// +/// NB: This table must have at least (e0-3)/24 + jk terms. +/// For quad precision (e0 <= 16360, jk = 6), this is 686. const ipio2 = [_]i32{ 0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, 0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, @@ -33,7 +32,6 @@ const ipio2 = [_]i32{ 0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, 0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B, - //#if LDBL_MAX_EXP > 1024 0x47C419, 0xC367CD, 0xDCE809, 0x2A8359, 0xC4768B, 0x961CA6, 0xDDAF44, 0xD15719, 0x053EA5, 0xFF0705, 0x3F7E33, 0xE832C2, 0xDE4F98, 0x327DBB, 0xC33D26, 0xEF6B1E, 0x5EF89F, 0x3A1F35, @@ -137,9 +135,7 @@ const ipio2 = [_]i32{ 0x237C7E, 0x32B90F, 0x8EF5A7, 0xE75614, 0x08F121, 0x2A9DB5, 0x4D7E6F, 0x5119A5, 0xABF9B5, 0xD6DF82, 0x61DD96, 0x023616, 0x9F3AC4, 0xA1A283, 0x6DED72, 0x7A8D39, 0xA9B882, 0x5C326B, - 0x5B2746, 0xED3400, 0x7700D2, 0x55F4FC, 0x4D5901, - 0x8071E0, - //#endif + 0x5B2746, 0xED3400, 0x7700D2, 0x55F4FC, 0x4D5901, 0x8071E0, }; const PIo2 = [_]f64{ @@ -157,109 +153,109 @@ fn U(x: anytype) usize { return @intCast(usize, x); } -// Returns the last three digits of N with y = x - N*pi/2 so that |y| < pi/2. -// -// The method is to compute the integer (mod 8) and fraction parts of -// (2/pi)*x without doing the full multiplication. In general we -// skip the part of the product that are known to be a huge integer ( -// more accurately, = 0 mod 8 ). Thus the number of operations are -// independent of the exponent of the input. -// -// (2/pi) is represented by an array of 24-bit integers in ipio2[]. -// -// Input parameters: -// x[] The input value (must be positive) is broken into nx -// pieces of 24-bit integers in double precision format. -// x[i] will be the i-th 24 bit of x. The scaled exponent -// of x[0] is given in input parameter e0 (i.e., x[0]*2^e0 -// match x's up to 24 bits. -// -// Example of breaking a double positive z into x[0]+x[1]+x[2]: -// e0 = ilogb(z)-23 -// z = scalbn(z,-e0) -// for i = 0,1,2 -// x[i] = floor(z) -// z = (z-x[i])*2**24 -// -// -// y[] ouput result in an array of double precision numbers. -// The dimension of y[] is: -// 24-bit precision 1 -// 53-bit precision 2 -// 64-bit precision 2 -// 113-bit precision 3 -// The actual value is the sum of them. Thus for 113-bit -// precison, one may have to do something like: -// -// long double t,w,r_head, r_tail; -// t = (long double)y[2] + (long double)y[1]; -// w = (long double)y[0]; -// r_head = t+w; -// r_tail = w - (r_head - t); -// -// e0 The exponent of x[0]. Must be <= 16360 or you need to -// expand the ipio2 table. -// -// nx dimension of x[] -// -// prec an integer indicating the precision: -// 0 24 bits (single) -// 1 53 bits (double) -// 2 64 bits (extended) -// 3 113 bits (quad) -// -// Here is the description of some local variables: -// -// jk jk+1 is the initial number of terms of ipio2[] needed -// in the computation. The minimum and recommended value -// for jk is 3,4,4,6 for single, double, extended, and quad. -// jk+1 must be 2 larger than you might expect so that our -// recomputation test works. (Up to 24 bits in the integer -// part (the 24 bits of it that we compute) and 23 bits in -// the fraction part may be lost to cancelation before we -// recompute.) -// -// jz local integer variable indicating the number of -// terms of ipio2[] used. -// -// jx nx - 1 -// -// jv index for pointing to the suitable ipio2[] for the -// computation. In general, we want -// ( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8 -// is an integer. Thus -// e0-3-24*jv >= 0 or (e0-3)/24 >= jv -// Hence jv = max(0,(e0-3)/24). -// -// jp jp+1 is the number of terms in PIo2[] needed, jp = jk. -// -// q[] double array with integral value, representing the -// 24-bits chunk of the product of x and 2/pi. -// -// q0 the corresponding exponent of q[0]. Note that the -// exponent for q[i] would be q0-24*i. -// -// PIo2[] double precision array, obtained by cutting pi/2 -// into 24 bits chunks. -// -// f[] ipio2[] in floating point -// -// iq[] integer array by breaking up q[] in 24-bits chunk. -// -// fq[] final product of x*(2/pi) in fq[0],..,fq[jk] -// -// ih integer. If >0 it indicates q[] is >= 0.5, hence -// it also indicates the *sign* of the result. -// +/// Returns the last three digits of N with y = x - N*pi/2 so that |y| < pi/2. /// -// -// Constants: -// The hexadecimal values are the intended ones for the following -// constants. The decimal values may be used, provided that the -// compiler will convert from decimal to binary accurately enough -// to produce the hexadecimal values shown. +/// The method is to compute the integer (mod 8) and fraction parts of +/// (2/pi)*x without doing the full multiplication. In general we +/// skip the part of the product that are known to be a huge integer ( +/// more accurately, = 0 mod 8 ). Thus the number of operations are +/// independent of the exponent of the input. +/// +/// (2/pi) is represented by an array of 24-bit integers in ipio2[]. +/// +/// Input parameters: +/// x[] The input value (must be positive) is broken into nx +/// pieces of 24-bit integers in double precision format. +/// x[i] will be the i-th 24 bit of x. The scaled exponent +/// of x[0] is given in input parameter e0 (i.e., x[0]*2^e0 +/// match x's up to 24 bits. +/// +/// Example of breaking a double positive z into x[0]+x[1]+x[2]: +/// e0 = ilogb(z)-23 +/// z = scalbn(z,-e0) +/// for i = 0,1,2 +/// x[i] = floor(z) +/// z = (z-x[i])*2**24 +/// +/// +/// y[] ouput result in an array of double precision numbers. +/// The dimension of y[] is: +/// 24-bit precision 1 +/// 53-bit precision 2 +/// 64-bit precision 2 +/// 113-bit precision 3 +/// The actual value is the sum of them. Thus for 113-bit +/// precison, one may have to do something like: +/// +/// long double t,w,r_head, r_tail; +/// t = (long double)y[2] + (long double)y[1]; +/// w = (long double)y[0]; +/// r_head = t+w; +/// r_tail = w - (r_head - t); +/// +/// e0 The exponent of x[0]. Must be <= 16360 or you need to +/// expand the ipio2 table. +/// +/// nx dimension of x[] +/// +/// prec an integer indicating the precision: +/// 0 24 bits (single) +/// 1 53 bits (double) +/// 2 64 bits (extended) +/// 3 113 bits (quad) +/// +/// Here is the description of some local variables: +/// +/// jk jk+1 is the initial number of terms of ipio2[] needed +/// in the computation. The minimum and recommended value +/// for jk is 3,4,4,6 for single, double, extended, and quad. +/// jk+1 must be 2 larger than you might expect so that our +/// recomputation test works. (Up to 24 bits in the integer +/// part (the 24 bits of it that we compute) and 23 bits in +/// the fraction part may be lost to cancelation before we +/// recompute.) +/// +/// jz local integer variable indicating the number of +/// terms of ipio2[] used. +/// +/// jx nx - 1 +/// +/// jv index for pointing to the suitable ipio2[] for the +/// computation. In general, we want +/// ( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8 +/// is an integer. Thus +/// e0-3-24*jv >= 0 or (e0-3)/24 >= jv +/// Hence jv = max(0,(e0-3)/24). +/// +/// jp jp+1 is the number of terms in PIo2[] needed, jp = jk. +/// +/// q[] double array with integral value, representing the +/// 24-bits chunk of the product of x and 2/pi. +/// +/// q0 the corresponding exponent of q[0]. Note that the +/// exponent for q[i] would be q0-24*i. +/// +/// PIo2[] double precision array, obtained by cutting pi/2 +/// into 24 bits chunks. +/// +/// f[] ipio2[] in floating point +/// +/// iq[] integer array by breaking up q[] in 24-bits chunk. +/// +/// fq[] final product of x*(2/pi) in fq[0],..,fq[jk] +/// +/// ih integer. If >0 it indicates q[] is >= 0.5, hence +/// it also indicates the *sign* of the result. +/// +/// +/// +/// Constants: +/// The hexadecimal values are the intended ones for the following +/// constants. The decimal values may be used, provided that the +/// compiler will convert from decimal to binary accurately enough +/// to produce the hexadecimal values shown. /// -pub fn __rem_pio2_large(x: []f64, y: []f64, e0: i32, nx: i32, prec: usize) i32 { +pub fn rem_pio2_large(x: []f64, y: []f64, e0: i32, nx: i32, prec: usize) i32 { var jz: i32 = undefined; var jx: i32 = undefined; var jv: i32 = undefined; @@ -333,7 +329,7 @@ pub fn __rem_pio2_large(x: []f64, y: []f64, e0: i32, nx: i32, prec: usize) i32 { // compute n z = math.scalbn(z, q0); // actual value of z - z -= 8.0 * math.floor(z * 0.125); // trim off integer >= 8 + z -= 8.0 * @floor(z * 0.125); // trim off integer >= 8 n = @floatToInt(i32, z); z -= @intToFloat(f64, n); ih = 0; diff --git a/lib/std/math/__rem_pio2f.zig b/lib/std/special/compiler_rt/rem_pio2f.zig index 5867fb30d9..34397dd734 100644 --- a/lib/std/math/__rem_pio2f.zig +++ b/lib/std/special/compiler_rt/rem_pio2f.zig @@ -3,8 +3,8 @@ // // https://git.musl-libc.org/cgit/musl/tree/src/math/__rem_pio2f.c -const std = @import("../std.zig"); -const __rem_pio2_large = @import("__rem_pio2_large.zig").__rem_pio2_large; +const std = @import("std"); +const rem_pio2_large = @import("rem_pio2_large.zig").rem_pio2_large; const math = std.math; const toint = 1.5 / math.floatEps(f64); @@ -19,8 +19,8 @@ const pio2_1t = 1.58932547735281966916e-08; // 0x3E5110b4, 0x611A6263 // Returns the remainder of x rem pi/2 in *y // use double precision for everything except passing x -// use __rem_pio2_large() for large x -pub fn __rem_pio2f(x: f32, y: *f64) i32 { +// use rem_pio2_large() for large x +pub fn rem_pio2f(x: f32, y: *f64) i32 { var tx: [1]f64 = undefined; var ty: [1]f64 = undefined; var @"fn": f64 = undefined; @@ -60,7 +60,7 @@ pub fn __rem_pio2f(x: f32, y: *f64) i32 { e0 = (ix >> 23) - (0x7f + 23); // e0 = ilogb(|x|)-23, positive ui = ix - (e0 << 23); tx[0] = @bitCast(f32, ui); - n = __rem_pio2_large(&tx, &ty, @intCast(i32, e0), 1, 0); + n = rem_pio2_large(&tx, &ty, @intCast(i32, e0), 1, 0); if (sign) { y.* = -ty[0]; return -n; diff --git a/lib/std/special/compiler_rt/round.zig b/lib/std/special/compiler_rt/round.zig new file mode 100644 index 0000000000..99201efcf8 --- /dev/null +++ b/lib/std/special/compiler_rt/round.zig @@ -0,0 +1,169 @@ +// Ported from musl, which is licensed under the MIT license: +// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT +// +// https://git.musl-libc.org/cgit/musl/tree/src/math/roundf.c +// https://git.musl-libc.org/cgit/musl/tree/src/math/round.c + +const std = @import("std"); +const math = std.math; +const expect = std.testing.expect; + +pub fn __roundh(x: f16) callconv(.C) f16 { + // TODO: more efficient implementation + return @floatCast(f16, roundf(x)); +} + +pub fn roundf(x_: f32) callconv(.C) f32 { + const f32_toint = 1.0 / math.floatEps(f32); + + var x = x_; + const u = @bitCast(u32, x); + const e = (u >> 23) & 0xFF; + var y: f32 = undefined; + + if (e >= 0x7F + 23) { + return x; + } + if (u >> 31 != 0) { + x = -x; + } + if (e < 0x7F - 1) { + math.doNotOptimizeAway(x + f32_toint); + return 0 * @bitCast(f32, u); + } + + y = x + f32_toint - f32_toint - x; + if (y > 0.5) { + y = y + x - 1; + } else if (y <= -0.5) { + y = y + x + 1; + } else { + y = y + x; + } + + if (u >> 31 != 0) { + return -y; + } else { + return y; + } +} + +pub fn round(x_: f64) callconv(.C) f64 { + const f64_toint = 1.0 / math.floatEps(f64); + + var x = x_; + const u = @bitCast(u64, x); + const e = (u >> 52) & 0x7FF; + var y: f64 = undefined; + + if (e >= 0x3FF + 52) { + return x; + } + if (u >> 63 != 0) { + x = -x; + } + if (e < 0x3ff - 1) { + math.doNotOptimizeAway(x + f64_toint); + return 0 * @bitCast(f64, u); + } + + y = x + f64_toint - f64_toint - x; + if (y > 0.5) { + y = y + x - 1; + } else if (y <= -0.5) { + y = y + x + 1; + } else { + y = y + x; + } + + if (u >> 63 != 0) { + return -y; + } else { + return y; + } +} + +pub fn __roundx(x: f80) callconv(.C) f80 { + // TODO: more efficient implementation + return @floatCast(f80, roundq(x)); +} + +pub fn roundq(x_: f128) callconv(.C) f128 { + const f128_toint = 1.0 / math.floatEps(f128); + + var x = x_; + const u = @bitCast(u128, x); + const e = (u >> 112) & 0x7FFF; + var y: f128 = undefined; + + if (e >= 0x3FFF + 112) { + return x; + } + if (u >> 127 != 0) { + x = -x; + } + if (e < 0x3FFF - 1) { + math.doNotOptimizeAway(x + f128_toint); + return 0 * @bitCast(f128, u); + } + + y = x + f128_toint - f128_toint - x; + if (y > 0.5) { + y = y + x - 1; + } else if (y <= -0.5) { + y = y + x + 1; + } else { + y = y + x; + } + + if (u >> 127 != 0) { + return -y; + } else { + return y; + } +} + +test "round32" { + try expect(roundf(1.3) == 1.0); + try expect(roundf(-1.3) == -1.0); + try expect(roundf(0.2) == 0.0); + try expect(roundf(1.8) == 2.0); +} + +test "round64" { + try expect(round(1.3) == 1.0); + try expect(round(-1.3) == -1.0); + try expect(round(0.2) == 0.0); + try expect(round(1.8) == 2.0); +} + +test "round128" { + try expect(roundq(1.3) == 1.0); + try expect(roundq(-1.3) == -1.0); + try expect(roundq(0.2) == 0.0); + try expect(roundq(1.8) == 2.0); +} + +test "round32.special" { + try expect(roundf(0.0) == 0.0); + try expect(roundf(-0.0) == -0.0); + try expect(math.isPositiveInf(roundf(math.inf(f32)))); + try expect(math.isNegativeInf(roundf(-math.inf(f32)))); + try expect(math.isNan(roundf(math.nan(f32)))); +} + +test "round64.special" { + try expect(round(0.0) == 0.0); + try expect(round(-0.0) == -0.0); + try expect(math.isPositiveInf(round(math.inf(f64)))); + try expect(math.isNegativeInf(round(-math.inf(f64)))); + try expect(math.isNan(round(math.nan(f64)))); +} + +test "round128.special" { + try expect(roundq(0.0) == 0.0); + try expect(roundq(-0.0) == -0.0); + try expect(math.isPositiveInf(roundq(math.inf(f128)))); + try expect(math.isNegativeInf(roundq(-math.inf(f128)))); + try expect(math.isNan(roundq(math.nan(f128)))); +} diff --git a/lib/std/math/sin.zig b/lib/std/special/compiler_rt/sin.zig index cf663b1d9e..aa77a961c7 100644 --- a/lib/std/math/sin.zig +++ b/lib/std/special/compiler_rt/sin.zig @@ -3,31 +3,21 @@ // // https://git.musl-libc.org/cgit/musl/tree/src/math/sinf.c // https://git.musl-libc.org/cgit/musl/tree/src/math/sin.c -// -const std = @import("../std.zig"); + +const std = @import("std"); const math = std.math; const expect = std.testing.expect; -const kernel = @import("__trig.zig"); -const __rem_pio2 = @import("__rem_pio2.zig").__rem_pio2; -const __rem_pio2f = @import("__rem_pio2f.zig").__rem_pio2f; - -/// Returns the sine of the radian value x. -/// -/// Special Cases: -/// - sin(+-0) = +-0 -/// - sin(+-inf) = nan -/// - sin(nan) = nan -pub fn sin(x: anytype) @TypeOf(x) { - const T = @TypeOf(x); - return switch (T) { - f32 => sin32(x), - f64 => sin64(x), - else => @compileError("sin not implemented for " ++ @typeName(T)), - }; +const kernel = @import("trig.zig"); +const rem_pio2 = @import("rem_pio2.zig").rem_pio2; +const rem_pio2f = @import("rem_pio2f.zig").rem_pio2f; + +pub fn __sinh(x: f16) callconv(.C) f16 { + // TODO: more efficient implementation + return @floatCast(f16, sinf(x)); } -fn sin32(x: f32) f32 { +pub fn sinf(x: f32) callconv(.C) f32 { // Small multiples of pi/2 rounded to double precision. const s1pio2: f64 = 1.0 * math.pi / 2.0; // 0x3FF921FB, 0x54442D18 const s2pio2: f64 = 2.0 * math.pi / 2.0; // 0x400921FB, 0x54442D18 @@ -73,7 +63,7 @@ fn sin32(x: f32) f32 { } var y: f64 = undefined; - const n = __rem_pio2f(x, &y); + const n = rem_pio2f(x, &y); return switch (n & 3) { 0 => kernel.__sindf(y), 1 => kernel.__cosdf(y), @@ -82,7 +72,7 @@ fn sin32(x: f32) f32 { }; } -fn sin64(x: f64) f64 { +pub fn sin(x: f64) callconv(.C) f64 { var ix = @bitCast(u64, x) >> 32; ix &= 0x7fffffff; @@ -102,7 +92,7 @@ fn sin64(x: f64) f64 { } var y: [2]f64 = undefined; - const n = __rem_pio2(x, &y); + const n = rem_pio2(x, &y); return switch (n & 3) { 0 => kernel.__sin(y[0], y[1], 1), 1 => kernel.__cos(y[0], y[1]), @@ -111,58 +101,68 @@ fn sin64(x: f64) f64 { }; } -test "math.sin" { - try expect(sin(@as(f32, 0.0)) == sin32(0.0)); - try expect(sin(@as(f64, 0.0)) == sin64(0.0)); +pub fn __sinx(x: f80) callconv(.C) f80 { + // TODO: more efficient implementation + return @floatCast(f80, sinq(x)); +} + +pub fn sinq(x: f128) callconv(.C) f128 { + // TODO: more correct implementation + return sin(@floatCast(f64, x)); +} + +test "sin" { + try expect(sin(@as(f32, 0.0)) == sinf(0.0)); + try expect(sin(@as(f64, 0.0)) == sin(0.0)); try expect(comptime (math.sin(@as(f64, 2))) == math.sin(@as(f64, 2))); } -test "math.sin32" { +test "sin32" { const epsilon = 0.00001; - try expect(math.approxEqAbs(f32, sin32(0.0), 0.0, epsilon)); - try expect(math.approxEqAbs(f32, sin32(0.2), 0.198669, epsilon)); - try expect(math.approxEqAbs(f32, sin32(0.8923), 0.778517, epsilon)); - try expect(math.approxEqAbs(f32, sin32(1.5), 0.997495, epsilon)); - try expect(math.approxEqAbs(f32, sin32(-1.5), -0.997495, epsilon)); - try expect(math.approxEqAbs(f32, sin32(37.45), -0.246544, epsilon)); - try expect(math.approxEqAbs(f32, sin32(89.123), 0.916166, epsilon)); + try expect(math.approxEqAbs(f32, sinf(0.0), 0.0, epsilon)); + try expect(math.approxEqAbs(f32, sinf(0.2), 0.198669, epsilon)); + try expect(math.approxEqAbs(f32, sinf(0.8923), 0.778517, epsilon)); + try expect(math.approxEqAbs(f32, sinf(1.5), 0.997495, epsilon)); + try expect(math.approxEqAbs(f32, sinf(-1.5), -0.997495, epsilon)); + try expect(math.approxEqAbs(f32, sinf(37.45), -0.246544, epsilon)); + try expect(math.approxEqAbs(f32, sinf(89.123), 0.916166, epsilon)); } -test "math.sin64" { +test "sin64" { const epsilon = 0.000001; - try expect(math.approxEqAbs(f64, sin64(0.0), 0.0, epsilon)); - try expect(math.approxEqAbs(f64, sin64(0.2), 0.198669, epsilon)); - try expect(math.approxEqAbs(f64, sin64(0.8923), 0.778517, epsilon)); - try expect(math.approxEqAbs(f64, sin64(1.5), 0.997495, epsilon)); - try expect(math.approxEqAbs(f64, sin64(-1.5), -0.997495, epsilon)); - try expect(math.approxEqAbs(f64, sin64(37.45), -0.246543, epsilon)); - try expect(math.approxEqAbs(f64, sin64(89.123), 0.916166, epsilon)); + try expect(math.approxEqAbs(f64, sin(0.0), 0.0, epsilon)); + try expect(math.approxEqAbs(f64, sin(0.2), 0.198669, epsilon)); + try expect(math.approxEqAbs(f64, sin(0.8923), 0.778517, epsilon)); + try expect(math.approxEqAbs(f64, sin(1.5), 0.997495, epsilon)); + try expect(math.approxEqAbs(f64, sin(-1.5), -0.997495, epsilon)); + try expect(math.approxEqAbs(f64, sin(37.45), -0.246543, epsilon)); + try expect(math.approxEqAbs(f64, sin(89.123), 0.916166, epsilon)); } -test "math.sin32.special" { - try expect(sin32(0.0) == 0.0); - try expect(sin32(-0.0) == -0.0); - try expect(math.isNan(sin32(math.inf(f32)))); - try expect(math.isNan(sin32(-math.inf(f32)))); - try expect(math.isNan(sin32(math.nan(f32)))); +test "sin32.special" { + try expect(sinf(0.0) == 0.0); + try expect(sinf(-0.0) == -0.0); + try expect(math.isNan(sinf(math.inf(f32)))); + try expect(math.isNan(sinf(-math.inf(f32)))); + try expect(math.isNan(sinf(math.nan(f32)))); } -test "math.sin64.special" { - try expect(sin64(0.0) == 0.0); - try expect(sin64(-0.0) == -0.0); - try expect(math.isNan(sin64(math.inf(f64)))); - try expect(math.isNan(sin64(-math.inf(f64)))); - try expect(math.isNan(sin64(math.nan(f64)))); +test "sin64.special" { + try expect(sin(0.0) == 0.0); + try expect(sin(-0.0) == -0.0); + try expect(math.isNan(sin(math.inf(f64)))); + try expect(math.isNan(sin(-math.inf(f64)))); + try expect(math.isNan(sin(math.nan(f64)))); } -test "math.sin32 #9901" { +test "sin32 #9901" { const float = @bitCast(f32, @as(u32, 0b11100011111111110000000000000000)); - _ = std.math.sin(float); + _ = sinf(float); } -test "math.sin64 #9901" { +test "sin64 #9901" { const float = @bitCast(f64, @as(u64, 0b1111111101000001000000001111110111111111100000000000000000000001)); - _ = std.math.sin(float); + _ = sin(float); } diff --git a/lib/std/special/compiler_rt/sincos.zig b/lib/std/special/compiler_rt/sincos.zig new file mode 100644 index 0000000000..fae326f182 --- /dev/null +++ b/lib/std/special/compiler_rt/sincos.zig @@ -0,0 +1,24 @@ +pub fn __sincosh(a: f16, r_sin: *f16, r_cos: *f16) callconv(.C) void { + r_sin.* = @sin(a); + r_cos.* = @cos(a); +} + +pub fn sincosf(a: f32, r_sin: *f32, r_cos: *f32) callconv(.C) void { + r_sin.* = @sin(a); + r_cos.* = @cos(a); +} + +pub fn sincos(a: f64, r_sin: *f64, r_cos: *f64) callconv(.C) void { + r_sin.* = @sin(a); + r_cos.* = @cos(a); +} + +pub fn __sincosx(a: f80, r_sin: *f80, r_cos: *f80) callconv(.C) void { + r_sin.* = @sin(a); + r_cos.* = @cos(a); +} + +pub fn sincosq(a: f128, r_sin: *f128, r_cos: *f128) callconv(.C) void { + r_sin.* = @sin(a); + r_cos.* = @cos(a); +} diff --git a/lib/std/special/compiler_rt/sqrt.zig b/lib/std/special/compiler_rt/sqrt.zig new file mode 100644 index 0000000000..ba07beb86e --- /dev/null +++ b/lib/std/special/compiler_rt/sqrt.zig @@ -0,0 +1,284 @@ +const std = @import("std"); +const math = std.math; + +pub fn __sqrth(x: f16) callconv(.C) f16 { + // TODO: more efficient implementation + return @floatCast(f16, sqrtf(x)); +} + +pub fn sqrtf(x: f32) callconv(.C) f32 { + const tiny: f32 = 1.0e-30; + const sign: i32 = @bitCast(i32, @as(u32, 0x80000000)); + var ix: i32 = @bitCast(i32, x); + + if ((ix & 0x7F800000) == 0x7F800000) { + return x * x + x; // sqrt(nan) = nan, sqrt(+inf) = +inf, sqrt(-inf) = snan + } + + // zero + if (ix <= 0) { + if (ix & ~sign == 0) { + return x; // sqrt (+-0) = +-0 + } + if (ix < 0) { + return math.snan(f32); + } + } + + // normalize + var m = ix >> 23; + if (m == 0) { + // subnormal + var i: i32 = 0; + while (ix & 0x00800000 == 0) : (i += 1) { + ix <<= 1; + } + m -= i - 1; + } + + m -= 127; // unbias exponent + ix = (ix & 0x007FFFFF) | 0x00800000; + + if (m & 1 != 0) { // odd m, double x to even + ix += ix; + } + + m >>= 1; // m = [m / 2] + + // sqrt(x) bit by bit + ix += ix; + var q: i32 = 0; // q = sqrt(x) + var s: i32 = 0; + var r: i32 = 0x01000000; // r = moving bit right -> left + + while (r != 0) { + const t = s + r; + if (t <= ix) { + s = t + r; + ix -= t; + q += r; + } + ix += ix; + r >>= 1; + } + + // floating add to find rounding direction + if (ix != 0) { + var z = 1.0 - tiny; // inexact + if (z >= 1.0) { + z = 1.0 + tiny; + if (z > 1.0) { + q += 2; + } else { + if (q & 1 != 0) { + q += 1; + } + } + } + } + + ix = (q >> 1) + 0x3f000000; + ix += m << 23; + return @bitCast(f32, ix); +} + +/// NOTE: The original code is full of implicit signed -> unsigned assumptions and u32 wraparound +/// behaviour. Most intermediate i32 values are changed to u32 where appropriate but there are +/// potentially some edge cases remaining that are not handled in the same way. +pub fn sqrt(x: f64) callconv(.C) f64 { + const tiny: f64 = 1.0e-300; + const sign: u32 = 0x80000000; + const u = @bitCast(u64, x); + + var ix0 = @intCast(u32, u >> 32); + var ix1 = @intCast(u32, u & 0xFFFFFFFF); + + // sqrt(nan) = nan, sqrt(+inf) = +inf, sqrt(-inf) = nan + if (ix0 & 0x7FF00000 == 0x7FF00000) { + return x * x + x; + } + + // sqrt(+-0) = +-0 + if (x == 0.0) { + return x; + } + // sqrt(-ve) = snan + if (ix0 & sign != 0) { + return math.snan(f64); + } + + // normalize x + var m = @intCast(i32, ix0 >> 20); + if (m == 0) { + // subnormal + while (ix0 == 0) { + m -= 21; + ix0 |= ix1 >> 11; + ix1 <<= 21; + } + + // subnormal + var i: u32 = 0; + while (ix0 & 0x00100000 == 0) : (i += 1) { + ix0 <<= 1; + } + m -= @intCast(i32, i) - 1; + ix0 |= ix1 >> @intCast(u5, 32 - i); + ix1 <<= @intCast(u5, i); + } + + // unbias exponent + m -= 1023; + ix0 = (ix0 & 0x000FFFFF) | 0x00100000; + if (m & 1 != 0) { + ix0 += ix0 + (ix1 >> 31); + ix1 = ix1 +% ix1; + } + m >>= 1; + + // sqrt(x) bit by bit + ix0 += ix0 + (ix1 >> 31); + ix1 = ix1 +% ix1; + + var q: u32 = 0; + var q1: u32 = 0; + var s0: u32 = 0; + var s1: u32 = 0; + var r: u32 = 0x00200000; + var t: u32 = undefined; + var t1: u32 = undefined; + + while (r != 0) { + t = s0 +% r; + if (t <= ix0) { + s0 = t + r; + ix0 -= t; + q += r; + } + ix0 = ix0 +% ix0 +% (ix1 >> 31); + ix1 = ix1 +% ix1; + r >>= 1; + } + + r = sign; + while (r != 0) { + t1 = s1 +% r; + t = s0; + if (t < ix0 or (t == ix0 and t1 <= ix1)) { + s1 = t1 +% r; + if (t1 & sign == sign and s1 & sign == 0) { + s0 += 1; + } + ix0 -= t; + if (ix1 < t1) { + ix0 -= 1; + } + ix1 = ix1 -% t1; + q1 += r; + } + ix0 = ix0 +% ix0 +% (ix1 >> 31); + ix1 = ix1 +% ix1; + r >>= 1; + } + + // rounding direction + if (ix0 | ix1 != 0) { + var z = 1.0 - tiny; // raise inexact + if (z >= 1.0) { + z = 1.0 + tiny; + if (q1 == 0xFFFFFFFF) { + q1 = 0; + q += 1; + } else if (z > 1.0) { + if (q1 == 0xFFFFFFFE) { + q += 1; + } + q1 += 2; + } else { + q1 += q1 & 1; + } + } + } + + ix0 = (q >> 1) + 0x3FE00000; + ix1 = q1 >> 1; + if (q & 1 != 0) { + ix1 |= 0x80000000; + } + + // NOTE: musl here appears to rely on signed twos-complement wraparound. +% has the same + // behaviour at least. + var iix0 = @intCast(i32, ix0); + iix0 = iix0 +% (m << 20); + + const uz = (@intCast(u64, iix0) << 32) | ix1; + return @bitCast(f64, uz); +} + +pub fn __sqrtx(x: f80) callconv(.C) f80 { + // TODO: more efficient implementation + return @floatCast(f80, sqrtq(x)); +} + +pub fn sqrtq(x: f128) callconv(.C) f128 { + // TODO: more correct implementation + return sqrt(@floatCast(f64, x)); +} + +test "sqrtf" { + const V = [_]f32{ + 0.0, + 4.089288054930154, + 7.538757127071935, + 8.97780793672623, + 5.304443821913729, + 5.682408965311888, + 0.5846878579110049, + 3.650338664297043, + 0.3178091951800732, + 7.1505232436382835, + 3.6589165881946464, + }; + + // Note that @sqrt will either generate the sqrt opcode (if supported by the + // target ISA) or a call to `sqrtf` otherwise. + for (V) |val| + try std.testing.expectEqual(@sqrt(val), sqrtf(val)); +} + +test "sqrtf special" { + try std.testing.expect(math.isPositiveInf(sqrtf(math.inf(f32)))); + try std.testing.expect(sqrtf(0.0) == 0.0); + try std.testing.expect(sqrtf(-0.0) == -0.0); + try std.testing.expect(math.isNan(sqrtf(-1.0))); + try std.testing.expect(math.isNan(sqrtf(math.nan(f32)))); +} + +test "sqrt" { + const V = [_]f64{ + 0.0, + 4.089288054930154, + 7.538757127071935, + 8.97780793672623, + 5.304443821913729, + 5.682408965311888, + 0.5846878579110049, + 3.650338664297043, + 0.3178091951800732, + 7.1505232436382835, + 3.6589165881946464, + }; + + // Note that @sqrt will either generate the sqrt opcode (if supported by the + // target ISA) or a call to `sqrtf` otherwise. + for (V) |val| + try std.testing.expectEqual(@sqrt(val), sqrt(val)); +} + +test "sqrt special" { + try std.testing.expect(math.isPositiveInf(sqrt(math.inf(f64)))); + try std.testing.expect(sqrt(0.0) == 0.0); + try std.testing.expect(sqrt(-0.0) == -0.0); + try std.testing.expect(math.isNan(sqrt(-1.0))); + try std.testing.expect(math.isNan(sqrt(math.nan(f64)))); +} diff --git a/lib/std/math/tan.zig b/lib/std/special/compiler_rt/tan.zig index fd5950df7c..d99f00b99e 100644 --- a/lib/std/math/tan.zig +++ b/lib/std/special/compiler_rt/tan.zig @@ -5,30 +5,20 @@ // https://git.musl-libc.org/cgit/musl/tree/src/math/tan.c // https://golang.org/src/math/tan.go -const std = @import("../std.zig"); +const std = @import("std"); const math = std.math; const expect = std.testing.expect; -const kernel = @import("__trig.zig"); -const __rem_pio2 = @import("__rem_pio2.zig").__rem_pio2; -const __rem_pio2f = @import("__rem_pio2f.zig").__rem_pio2f; - -/// Returns the tangent of the radian value x. -/// -/// Special Cases: -/// - tan(+-0) = +-0 -/// - tan(+-inf) = nan -/// - tan(nan) = nan -pub fn tan(x: anytype) @TypeOf(x) { - const T = @TypeOf(x); - return switch (T) { - f32 => tan32(x), - f64 => tan64(x), - else => @compileError("tan not implemented for " ++ @typeName(T)), - }; +const kernel = @import("trig.zig"); +const rem_pio2 = @import("rem_pio2.zig").rem_pio2; +const rem_pio2f = @import("rem_pio2f.zig").rem_pio2f; + +pub fn __tanh(x: f16) callconv(.C) f16 { + // TODO: more efficient implementation + return @floatCast(f16, tanf(x)); } -fn tan32(x: f32) f32 { +pub fn tanf(x: f32) callconv(.C) f32 { // Small multiples of pi/2 rounded to double precision. const t1pio2: f64 = 1.0 * math.pi / 2.0; // 0x3FF921FB, 0x54442D18 const t2pio2: f64 = 2.0 * math.pi / 2.0; // 0x400921FB, 0x54442D18 @@ -68,11 +58,11 @@ fn tan32(x: f32) f32 { } var y: f64 = undefined; - const n = __rem_pio2f(x, &y); + const n = rem_pio2f(x, &y); return kernel.__tandf(y, n & 1 != 0); } -fn tan64(x: f64) f64 { +pub fn tan(x: f64) callconv(.C) f64 { var ix = @bitCast(u64, x) >> 32; ix &= 0x7fffffff; @@ -92,49 +82,59 @@ fn tan64(x: f64) f64 { } var y: [2]f64 = undefined; - const n = __rem_pio2(x, &y); + const n = rem_pio2(x, &y); return kernel.__tan(y[0], y[1], n & 1 != 0); } -test "math.tan" { - try expect(tan(@as(f32, 0.0)) == tan32(0.0)); - try expect(tan(@as(f64, 0.0)) == tan64(0.0)); +pub fn __tanx(x: f80) callconv(.C) f80 { + // TODO: more efficient implementation + return @floatCast(f80, tanq(x)); +} + +pub fn tanq(x: f128) callconv(.C) f128 { + // TODO: more correct implementation + return tan(@floatCast(f64, x)); +} + +test "tan" { + try expect(tan(@as(f32, 0.0)) == tanf(0.0)); + try expect(tan(@as(f64, 0.0)) == tan(0.0)); } -test "math.tan32" { +test "tan32" { const epsilon = 0.00001; - try expect(math.approxEqAbs(f32, tan32(0.0), 0.0, epsilon)); - try expect(math.approxEqAbs(f32, tan32(0.2), 0.202710, epsilon)); - try expect(math.approxEqAbs(f32, tan32(0.8923), 1.240422, epsilon)); - try expect(math.approxEqAbs(f32, tan32(1.5), 14.101420, epsilon)); - try expect(math.approxEqAbs(f32, tan32(37.45), -0.254397, epsilon)); - try expect(math.approxEqAbs(f32, tan32(89.123), 2.285852, epsilon)); + try expect(math.approxEqAbs(f32, tanf(0.0), 0.0, epsilon)); + try expect(math.approxEqAbs(f32, tanf(0.2), 0.202710, epsilon)); + try expect(math.approxEqAbs(f32, tanf(0.8923), 1.240422, epsilon)); + try expect(math.approxEqAbs(f32, tanf(1.5), 14.101420, epsilon)); + try expect(math.approxEqAbs(f32, tanf(37.45), -0.254397, epsilon)); + try expect(math.approxEqAbs(f32, tanf(89.123), 2.285852, epsilon)); } -test "math.tan64" { +test "tan64" { const epsilon = 0.000001; - try expect(math.approxEqAbs(f64, tan64(0.0), 0.0, epsilon)); - try expect(math.approxEqAbs(f64, tan64(0.2), 0.202710, epsilon)); - try expect(math.approxEqAbs(f64, tan64(0.8923), 1.240422, epsilon)); - try expect(math.approxEqAbs(f64, tan64(1.5), 14.101420, epsilon)); - try expect(math.approxEqAbs(f64, tan64(37.45), -0.254397, epsilon)); - try expect(math.approxEqAbs(f64, tan64(89.123), 2.2858376, epsilon)); + try expect(math.approxEqAbs(f64, tan(0.0), 0.0, epsilon)); + try expect(math.approxEqAbs(f64, tan(0.2), 0.202710, epsilon)); + try expect(math.approxEqAbs(f64, tan(0.8923), 1.240422, epsilon)); + try expect(math.approxEqAbs(f64, tan(1.5), 14.101420, epsilon)); + try expect(math.approxEqAbs(f64, tan(37.45), -0.254397, epsilon)); + try expect(math.approxEqAbs(f64, tan(89.123), 2.2858376, epsilon)); } -test "math.tan32.special" { - try expect(tan32(0.0) == 0.0); - try expect(tan32(-0.0) == -0.0); - try expect(math.isNan(tan32(math.inf(f32)))); - try expect(math.isNan(tan32(-math.inf(f32)))); - try expect(math.isNan(tan32(math.nan(f32)))); +test "tan32.special" { + try expect(tanf(0.0) == 0.0); + try expect(tanf(-0.0) == -0.0); + try expect(math.isNan(tanf(math.inf(f32)))); + try expect(math.isNan(tanf(-math.inf(f32)))); + try expect(math.isNan(tanf(math.nan(f32)))); } -test "math.tan64.special" { - try expect(tan64(0.0) == 0.0); - try expect(tan64(-0.0) == -0.0); - try expect(math.isNan(tan64(math.inf(f64)))); - try expect(math.isNan(tan64(-math.inf(f64)))); - try expect(math.isNan(tan64(math.nan(f64)))); +test "tan64.special" { + try expect(tan(0.0) == 0.0); + try expect(tan(-0.0) == -0.0); + try expect(math.isNan(tan(math.inf(f64)))); + try expect(math.isNan(tan(-math.inf(f64)))); + try expect(math.isNan(tan(math.nan(f64)))); } diff --git a/lib/std/math/__trig.zig b/lib/std/special/compiler_rt/trig.zig index 0c08ed58bd..8ece83515e 100644 --- a/lib/std/math/__trig.zig +++ b/lib/std/special/compiler_rt/trig.zig @@ -8,41 +8,41 @@ // https://git.musl-libc.org/cgit/musl/tree/src/math/__tand.c // https://git.musl-libc.org/cgit/musl/tree/src/math/__tandf.c -// kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164 -// Input x is assumed to be bounded by ~pi/4 in magnitude. -// Input y is the tail of x. -// -// Algorithm -// 1. Since cos(-x) = cos(x), we need only to consider positive x. -// 2. if x < 2^-27 (hx<0x3e400000 0), return 1 with inexact if x!=0. -// 3. cos(x) is approximated by a polynomial of degree 14 on -// [0,pi/4] -// 4 14 -// cos(x) ~ 1 - x*x/2 + C1*x + ... + C6*x -// where the remez error is -// -// | 2 4 6 8 10 12 14 | -58 -// |cos(x)-(1-.5*x +C1*x +C2*x +C3*x +C4*x +C5*x +C6*x )| <= 2 -// | | -// -// 4 6 8 10 12 14 -// 4. let r = C1*x +C2*x +C3*x +C4*x +C5*x +C6*x , then -// cos(x) ~ 1 - x*x/2 + r -// since cos(x+y) ~ cos(x) - sin(x)*y -// ~ cos(x) - x*y, -// a correction term is necessary in cos(x) and hence -// cos(x+y) = 1 - (x*x/2 - (r - x*y)) -// For better accuracy, rearrange to -// cos(x+y) ~ w + (tmp + (r-x*y)) -// where w = 1 - x*x/2 and tmp is a tiny correction term -// (1 - x*x/2 == w + tmp exactly in infinite precision). -// The exactness of w + tmp in infinite precision depends on w -// and tmp having the same precision as x. If they have extra -// precision due to compiler bugs, then the extra precision is -// only good provided it is retained in all terms of the final -// expression for cos(). Retention happens in all cases tested -// under FreeBSD, so don't pessimize things by forcibly clipping -// any extra precision in w. +/// kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164 +/// Input x is assumed to be bounded by ~pi/4 in magnitude. +/// Input y is the tail of x. +/// +/// Algorithm +/// 1. Since cos(-x) = cos(x), we need only to consider positive x. +/// 2. if x < 2^-27 (hx<0x3e400000 0), return 1 with inexact if x!=0. +/// 3. cos(x) is approximated by a polynomial of degree 14 on +/// [0,pi/4] +/// 4 14 +/// cos(x) ~ 1 - x*x/2 + C1*x + ... + C6*x +/// where the remez error is +/// +/// | 2 4 6 8 10 12 14 | -58 +/// |cos(x)-(1-.5*x +C1*x +C2*x +C3*x +C4*x +C5*x +C6*x )| <= 2 +/// | | +/// +/// 4 6 8 10 12 14 +/// 4. let r = C1*x +C2*x +C3*x +C4*x +C5*x +C6*x , then +/// cos(x) ~ 1 - x*x/2 + r +/// since cos(x+y) ~ cos(x) - sin(x)*y +/// ~ cos(x) - x*y, +/// a correction term is necessary in cos(x) and hence +/// cos(x+y) = 1 - (x*x/2 - (r - x*y)) +/// For better accuracy, rearrange to +/// cos(x+y) ~ w + (tmp + (r-x*y)) +/// where w = 1 - x*x/2 and tmp is a tiny correction term +/// (1 - x*x/2 == w + tmp exactly in infinite precision). +/// The exactness of w + tmp in infinite precision depends on w +/// and tmp having the same precision as x. If they have extra +/// precision due to compiler bugs, then the extra precision is +/// only good provided it is retained in all terms of the final +/// expression for cos(). Retention happens in all cases tested +/// under FreeBSD, so don't pessimize things by forcibly clipping +/// any extra precision in w. pub fn __cos(x: f64, y: f64) f64 { const C1 = 4.16666666666666019037e-02; // 0x3FA55555, 0x5555554C const C2 = -1.38888888888741095749e-03; // 0xBF56C16C, 0x16C15177 @@ -73,33 +73,33 @@ pub fn __cosdf(x: f64) f32 { return @floatCast(f32, ((1.0 + z * C0) + w * C1) + (w * z) * r); } -// kernel sin function on ~[-pi/4, pi/4] (except on -0), pi/4 ~ 0.7854 -// Input x is assumed to be bounded by ~pi/4 in magnitude. -// Input y is the tail of x. -// Input iy indicates whether y is 0. (if iy=0, y assume to be 0). -// -// Algorithm -// 1. Since sin(-x) = -sin(x), we need only to consider positive x. -// 2. Callers must return sin(-0) = -0 without calling here since our -// odd polynomial is not evaluated in a way that preserves -0. -// Callers may do the optimization sin(x) ~ x for tiny x. -// 3. sin(x) is approximated by a polynomial of degree 13 on -// [0,pi/4] -// 3 13 -// sin(x) ~ x + S1*x + ... + S6*x -// where -// -// |sin(x) 2 4 6 8 10 12 | -58 -// |----- - (1+S1*x +S2*x +S3*x +S4*x +S5*x +S6*x )| <= 2 -// | x | -// -// 4. sin(x+y) = sin(x) + sin'(x')*y -// ~ sin(x) + (1-x*x/2)*y -// For better accuracy, let -// 3 2 2 2 2 -// r = x *(S2+x *(S3+x *(S4+x *(S5+x *S6)))) -// then 3 2 -// sin(x) = x + (S1*x + (x *(r-y/2)+y)) +/// kernel sin function on ~[-pi/4, pi/4] (except on -0), pi/4 ~ 0.7854 +/// Input x is assumed to be bounded by ~pi/4 in magnitude. +/// Input y is the tail of x. +/// Input iy indicates whether y is 0. (if iy=0, y assume to be 0). +/// +/// Algorithm +/// 1. Since sin(-x) = -sin(x), we need only to consider positive x. +/// 2. Callers must return sin(-0) = -0 without calling here since our +/// odd polynomial is not evaluated in a way that preserves -0. +/// Callers may do the optimization sin(x) ~ x for tiny x. +/// 3. sin(x) is approximated by a polynomial of degree 13 on +/// [0,pi/4] +/// 3 13 +/// sin(x) ~ x + S1*x + ... + S6*x +/// where +/// +/// |sin(x) 2 4 6 8 10 12 | -58 +/// |----- - (1+S1*x +S2*x +S3*x +S4*x +S5*x +S6*x )| <= 2 +/// | x | +/// +/// 4. sin(x+y) = sin(x) + sin'(x')*y +/// ~ sin(x) + (1-x*x/2)*y +/// For better accuracy, let +/// 3 2 2 2 2 +/// r = x *(S2+x *(S3+x *(S4+x *(S5+x *S6)))) +/// then 3 2 +/// sin(x) = x + (S1*x + (x *(r-y/2)+y)) pub fn __sin(x: f64, y: f64, iy: i32) f64 { const S1 = -1.66666666666666324348e-01; // 0xBFC55555, 0x55555549 const S2 = 8.33333333332248946124e-03; // 0x3F811111, 0x1110F8A6 @@ -134,38 +134,38 @@ pub fn __sindf(x: f64) f32 { return @floatCast(f32, (x + s * (S1 + z * S2)) + s * w * r); } -// kernel tan function on ~[-pi/4, pi/4] (except on -0), pi/4 ~ 0.7854 -// Input x is assumed to be bounded by ~pi/4 in magnitude. -// Input y is the tail of x. -// Input odd indicates whether tan (if odd = 0) or -1/tan (if odd = 1) is returned. -// -// Algorithm -// 1. Since tan(-x) = -tan(x), we need only to consider positive x. -// 2. Callers must return tan(-0) = -0 without calling here since our -// odd polynomial is not evaluated in a way that preserves -0. -// Callers may do the optimization tan(x) ~ x for tiny x. -// 3. tan(x) is approximated by a odd polynomial of degree 27 on -// [0,0.67434] -// 3 27 -// tan(x) ~ x + T1*x + ... + T13*x -// where -// -// |tan(x) 2 4 26 | -59.2 -// |----- - (1+T1*x +T2*x +.... +T13*x )| <= 2 -// | x | -// -// Note: tan(x+y) = tan(x) + tan'(x)*y -// ~ tan(x) + (1+x*x)*y -// Therefore, for better accuracy in computing tan(x+y), let -// 3 2 2 2 2 -// r = x *(T2+x *(T3+x *(...+x *(T12+x *T13)))) -// then -// 3 2 -// tan(x+y) = x + (T1*x + (x *(r+y)+y)) -// -// 4. For x in [0.67434,pi/4], let y = pi/4 - x, then -// tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y)) -// = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y))) +/// kernel tan function on ~[-pi/4, pi/4] (except on -0), pi/4 ~ 0.7854 +/// Input x is assumed to be bounded by ~pi/4 in magnitude. +/// Input y is the tail of x. +/// Input odd indicates whether tan (if odd = 0) or -1/tan (if odd = 1) is returned. +/// +/// Algorithm +/// 1. Since tan(-x) = -tan(x), we need only to consider positive x. +/// 2. Callers must return tan(-0) = -0 without calling here since our +/// odd polynomial is not evaluated in a way that preserves -0. +/// Callers may do the optimization tan(x) ~ x for tiny x. +/// 3. tan(x) is approximated by a odd polynomial of degree 27 on +/// [0,0.67434] +/// 3 27 +/// tan(x) ~ x + T1*x + ... + T13*x +/// where +/// +/// |tan(x) 2 4 26 | -59.2 +/// |----- - (1+T1*x +T2*x +.... +T13*x )| <= 2 +/// | x | +/// +/// Note: tan(x+y) = tan(x) + tan'(x)*y +/// ~ tan(x) + (1+x*x)*y +/// Therefore, for better accuracy in computing tan(x+y), let +/// 3 2 2 2 2 +/// r = x *(T2+x *(T3+x *(...+x *(T12+x *T13)))) +/// then +/// 3 2 +/// tan(x+y) = x + (T1*x + (x *(r+y)+y)) +/// +/// 4. For x in [0.67434,pi/4], let y = pi/4 - x, then +/// tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y)) +/// = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y))) pub fn __tan(x_: f64, y_: f64, odd: bool) f64 { var x = x_; var y = y_; diff --git a/lib/std/special/compiler_rt/trunc.zig b/lib/std/special/compiler_rt/trunc.zig new file mode 100644 index 0000000000..5406f9a02d --- /dev/null +++ b/lib/std/special/compiler_rt/trunc.zig @@ -0,0 +1,124 @@ +// Ported from musl, which is licensed under the MIT license: +// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT +// +// https://git.musl-libc.org/cgit/musl/tree/src/math/truncf.c +// https://git.musl-libc.org/cgit/musl/tree/src/math/trunc.c + +const std = @import("std"); +const math = std.math; +const expect = std.testing.expect; + +pub fn __trunch(x: f16) callconv(.C) f16 { + // TODO: more efficient implementation + return @floatCast(f16, truncf(x)); +} + +pub fn truncf(x: f32) callconv(.C) f32 { + const u = @bitCast(u32, x); + var e = @intCast(i32, ((u >> 23) & 0xFF)) - 0x7F + 9; + var m: u32 = undefined; + + if (e >= 23 + 9) { + return x; + } + if (e < 9) { + e = 1; + } + + m = @as(u32, math.maxInt(u32)) >> @intCast(u5, e); + if (u & m == 0) { + return x; + } else { + math.doNotOptimizeAway(x + 0x1p120); + return @bitCast(f32, u & ~m); + } +} + +pub fn trunc(x: f64) callconv(.C) f64 { + const u = @bitCast(u64, x); + var e = @intCast(i32, ((u >> 52) & 0x7FF)) - 0x3FF + 12; + var m: u64 = undefined; + + if (e >= 52 + 12) { + return x; + } + if (e < 12) { + e = 1; + } + + m = @as(u64, math.maxInt(u64)) >> @intCast(u6, e); + if (u & m == 0) { + return x; + } else { + math.doNotOptimizeAway(x + 0x1p120); + return @bitCast(f64, u & ~m); + } +} + +pub fn __truncx(x: f80) callconv(.C) f80 { + // TODO: more efficient implementation + return @floatCast(f80, truncq(x)); +} + +pub fn truncq(x: f128) callconv(.C) f128 { + const u = @bitCast(u128, x); + var e = @intCast(i32, ((u >> 112) & 0x7FFF)) - 0x3FFF + 16; + var m: u128 = undefined; + + if (e >= 112 + 16) { + return x; + } + if (e < 16) { + e = 1; + } + + m = @as(u128, math.maxInt(u128)) >> @intCast(u7, e); + if (u & m == 0) { + return x; + } else { + math.doNotOptimizeAway(x + 0x1p120); + return @bitCast(f128, u & ~m); + } +} + +test "trunc32" { + try expect(truncf(1.3) == 1.0); + try expect(truncf(-1.3) == -1.0); + try expect(truncf(0.2) == 0.0); +} + +test "trunc64" { + try expect(trunc(1.3) == 1.0); + try expect(trunc(-1.3) == -1.0); + try expect(trunc(0.2) == 0.0); +} + +test "trunc128" { + try expect(truncq(1.3) == 1.0); + try expect(truncq(-1.3) == -1.0); + try expect(truncq(0.2) == 0.0); +} + +test "trunc32.special" { + try expect(truncf(0.0) == 0.0); // 0x3F800000 + try expect(truncf(-0.0) == -0.0); + try expect(math.isPositiveInf(truncf(math.inf(f32)))); + try expect(math.isNegativeInf(truncf(-math.inf(f32)))); + try expect(math.isNan(truncf(math.nan(f32)))); +} + +test "trunc64.special" { + try expect(trunc(0.0) == 0.0); + try expect(trunc(-0.0) == -0.0); + try expect(math.isPositiveInf(trunc(math.inf(f64)))); + try expect(math.isNegativeInf(trunc(-math.inf(f64)))); + try expect(math.isNan(trunc(math.nan(f64)))); +} + +test "trunc128.special" { + try expect(truncq(0.0) == 0.0); + try expect(truncq(-0.0) == -0.0); + try expect(math.isPositiveInf(truncq(math.inf(f128)))); + try expect(math.isNegativeInf(truncq(-math.inf(f128)))); + try expect(math.isNan(truncq(math.nan(f128)))); +} diff --git a/lib/std/testing.zig b/lib/std/testing.zig index 004e2d0fa7..cfdf300c04 100644 --- a/lib/std/testing.zig +++ b/lib/std/testing.zig @@ -265,7 +265,7 @@ pub fn expectApproxEqRel(expected: anytype, actual: @TypeOf(expected), tolerance test "expectApproxEqRel" { inline for ([_]type{ f16, f32, f64, f128 }) |T| { const eps_value = comptime math.epsilon(T); - const sqrt_eps_value = comptime math.sqrt(eps_value); + const sqrt_eps_value = comptime @sqrt(eps_value); const pos_x: T = 12.0; const pos_y: T = pos_x + 2 * eps_value; |
