diff options
| author | Marc Tiehuis <marctiehuis@gmail.com> | 2017-06-16 20:26:10 +1200 |
|---|---|---|
| committer | Marc Tiehuis <marctiehuis@gmail.com> | 2017-06-16 20:32:31 +1200 |
| commit | 4c16f9a3c35b23b9917f2a27b91ba8cd20e6fd82 (patch) | |
| tree | 778f0f06734f7dc17e9269ee1cf5b513f7b504c0 | |
| parent | 865b53f2860405a718262abf9a794d2bf9529dbc (diff) | |
| download | zig-4c16f9a3c35b23b9917f2a27b91ba8cd20e6fd82.tar.gz zig-4c16f9a3c35b23b9917f2a27b91ba8cd20e6fd82.zip | |
Add math library
This covers the majority of the functions as covered by the C99
specification for a math library.
Code is adapted primarily from musl libc, with the pow and standard
trigonometric functions adapted from the Go stdlib.
Changes:
- Remove assert expose in index and import as needed.
- Add float log function and merge with existing base 2 integer
implementation.
See https://github.com/tiehuis/zig-fmath.
See #374.
47 files changed, 6150 insertions, 174 deletions
diff --git a/std/math/_expo2.zig b/std/math/_expo2.zig new file mode 100644 index 0000000000..e0d8a8130d --- /dev/null +++ b/std/math/_expo2.zig @@ -0,0 +1,28 @@ +const math = @import("index.zig"); + +pub fn expo2(x: var) -> @typeOf(x) { + const T = @typeOf(x); + switch (T) { + f32 => expo2f(x), + f64 => expo2d(x), + else => @compileError("expo2 not implemented for " ++ @typeName(T)), + } +} + +fn expo2f(x: f32) -> f32 { + const k: u32 = 235; + const kln2 = 0x1.45C778p+7; + + const u = (0x7F + k / 2) << 23; + const scale = @bitCast(f32, u); + math.exp(x - kln2) * scale * scale +} + +fn expo2d(x: f64) -> f64 { + const k: u32 = 2043; + const kln2 = 0x1.62066151ADD8BP+10; + + const u = (0x3FF + k / 2) << 20; + const scale = @bitCast(f64, u64(u) << 32); + math.exp(x - kln2) * scale * scale +} diff --git a/std/math/acos.zig b/std/math/acos.zig new file mode 100644 index 0000000000..c6661a272d --- /dev/null +++ b/std/math/acos.zig @@ -0,0 +1,165 @@ +const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; + +pub fn acos(x: var) -> @typeOf(x) { + const T = @typeOf(x); + switch (T) { + f32 => @inlineCall(acos32, x), + f64 => @inlineCall(acos64, x), + else => @compileError("acos not implemented for " ++ @typeName(T)), + } +} + +fn r32(z: f32) -> f32 { + const pS0 = 1.6666586697e-01; + const pS1 = -4.2743422091e-02; + const pS2 = -8.6563630030e-03; + const qS1 = -7.0662963390e-01; + + const p = z * (pS0 + z * (pS1 + z * pS2)); + const q = 1.0 + z * qS1; + p / q +} + +fn acos32(x: f32) -> f32 { + const pio2_hi = 1.5707962513e+00; + const pio2_lo = 7.5497894159e-08; + + const hx: u32 = @bitCast(u32, x); + const ix: u32 = hx & 0x7FFFFFFF; + + // |x| >= 1 or nan + if (ix >= 0x3F800000) { + if (ix == 0x3F800000) { + if (hx >> 31 != 0) { + return 2.0 * pio2_hi + 0x1.0p-120; + } else { + return 0; + } + } else { + return 0 / (x - x); + } + } + + // |x| < 0.5 + if (ix < 0x3F000000) { + if (ix <= 0x32800000) { // |x| < 2^(-26) + return pio2_hi + 0x1.0p-120; + } else { + return pio2_hi - (x - (pio2_lo - x * r32(x * x))); + } + } + + // x < -0.5 + if (hx >> 31 != 0) { + const z = (1 + x) * 0.5; + const s = math.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 jx = @bitCast(u32, s); + const df = @bitCast(f32, jx & 0xFFFFF000); + const c = (z - df * df) / (s + df); + const w = r32(z) * s + c; + 2 * (df + w) +} + +fn r64(z: f64) -> f64 { + const pS0: f64 = 1.66666666666666657415e-01; + const pS1: f64 = -3.25565818622400915405e-01; + const pS2: f64 = 2.01212532134862925881e-01; + const pS3: f64 = -4.00555345006794114027e-02; + const pS4: f64 = 7.91534994289814532176e-04; + const pS5: f64 = 3.47933107596021167570e-05; + const qS1: f64 = -2.40339491173441421878e+00; + const qS2: f64 = 2.02094576023350569471e+00; + const qS3: f64 = -6.88283971605453293030e-01; + const qS4: f64 = 7.70381505559019352791e-02; + + const p = z * (pS0 + z * (pS1 + z * (pS2 + z * (pS3 + z * (pS4 + z * pS5))))); + const q = 1.0 + z * (qS1 + z * (qS2 + z * (qS3 + z * qS4))); + p / q +} + +fn acos64(x: f64) -> f64 { + const pio2_hi: f64 = 1.57079632679489655800e+00; + const pio2_lo: f64 = 6.12323399573676603587e-17; + + const ux = @bitCast(u64, x); + const hx = u32(ux >> 32); + const ix = hx & 0x7FFFFFFF; + + // |x| >= 1 or nan + if (ix >= 0x3FF00000) { + const lx = u32(ux & 0xFFFFFFFF); + + // acos(1) = 0, acos(-1) = pi + if ((ix - 0x3FF00000) | lx == 0) { + if (hx >> 31 != 0) { + return 2 * pio2_hi + 0x1.0p-120; + } else { + return 0; + } + } + + return 0 / (x - x); + } + + // |x| < 0.5 + if (ix < 0x3FE00000) { + // |x| < 2^(-57) + if (ix <= 0x3C600000) { + return pio2_hi + 0x1.0p-120; + } else { + return pio2_hi - (x - (pio2_lo - x * r64(x * x))); + } + } + + // x < -0.5 + if (hx >> 31 != 0) { + const z = (1.0 + x) * 0.5; + const s = math.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 jx = @bitCast(u64, s); + const df = @bitCast(f64, jx & 0xFFFFFFFF00000000); + const c = (z - df * df) / (s + df); + const w = r64(z) * s + c; + 2 * (df + w) +} + +test "acos" { + assert(acos(f32(0.0)) == acos32(0.0)); + assert(acos(f64(0.0)) == acos64(0.0)); +} + +test "acos32" { + const epsilon = 0.000001; + + assert(math.approxEq(f32, acos32(0.0), 1.570796, epsilon)); + assert(math.approxEq(f32, acos32(0.2), 1.369438, epsilon)); + assert(math.approxEq(f32, acos32(0.3434), 1.220262, epsilon)); + assert(math.approxEq(f32, acos32(0.5), 1.047198, epsilon)); + assert(math.approxEq(f32, acos32(0.8923), 0.468382, epsilon)); + assert(math.approxEq(f32, acos32(-0.2), 1.772154, epsilon)); +} + +test "acos64" { + const epsilon = 0.000001; + + assert(math.approxEq(f64, acos64(0.0), 1.570796, epsilon)); + assert(math.approxEq(f64, acos64(0.2), 1.369438, epsilon)); + assert(math.approxEq(f64, acos64(0.3434), 1.220262, epsilon)); + assert(math.approxEq(f64, acos64(0.5), 1.047198, epsilon)); + assert(math.approxEq(f64, acos64(0.8923), 0.468382, epsilon)); + assert(math.approxEq(f64, acos64(-0.2), 1.772154, epsilon)); +} diff --git a/std/math/acosh.zig b/std/math/acosh.zig new file mode 100644 index 0000000000..7da81360a9 --- /dev/null +++ b/std/math/acosh.zig @@ -0,0 +1,71 @@ +const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; + +pub fn acosh(x: var) -> @typeOf(x) { + const T = @typeOf(x); + switch (T) { + f32 => @inlineCall(acoshf, x), + f64 => @inlineCall(acoshd, x), + else => @compileError("acosh not implemented for " ++ @typeName(T)), + } +} + +// acosh(x) = log(x + sqrt(x * x - 1)) +fn acoshf(x: f32) -> f32 { + const u = @bitCast(u32, x); + const i = u & 0x7FFFFFFF; + + // |x| < 2, invalid if x < 1 or nan + if (i < 0x3F800000 + (1 << 23)) { + math.log1p(x - 1 + math.sqrt((x - 1) * (x - 1) + 2 * (x - 1))) + } + // |x| < 0x1p12 + else if (i < 0x3F800000 + (12 << 23)) { + math.ln(2 * x - 1 / (x + math.sqrt(x * x - 1))) + } + // |x| >= 0x1p12 + else { + math.ln(x) + 0.693147180559945309417232121458176568 + } +} + +fn acoshd(x: f64) -> f64 { + const u = @bitCast(u64, x); + const e = (u >> 52) & 0x7FF; + + // |x| < 2, invalid if x < 1 or nan + if (e < 0x3FF + 1) { + math.log1p(x - 1 + math.sqrt((x - 1) * (x - 1) + 2 * (x - 1))) + } + // |x| < 0x1p26 + else if (e < 0x3FF + 26) { + math.ln(2 * x - 1 / (x + math.sqrt(x * x - 1))) + } + // |x| >= 0x1p26 or nan + else { + math.ln(x) + 0.693147180559945309417232121458176568 + } +} + +test "acosh" { + assert(acosh(f32(1.5)) == acoshf(1.5)); + assert(acosh(f64(1.5)) == acoshd(1.5)); +} + +test "acoshf" { + const epsilon = 0.000001; + + assert(math.approxEq(f32, acoshf(1.5), 0.962424, epsilon)); + assert(math.approxEq(f32, acoshf(37.45), 4.315976, epsilon)); + assert(math.approxEq(f32, acoshf(89.123), 5.183133, epsilon)); + assert(math.approxEq(f32, acoshf(123123.234375), 12.414088, epsilon)); +} + +test "acoshd" { + const epsilon = 0.000001; + + assert(math.approxEq(f64, acoshd(1.5), 0.962424, epsilon)); + assert(math.approxEq(f64, acoshd(37.45), 4.315976, epsilon)); + assert(math.approxEq(f64, acoshd(89.123), 5.183133, epsilon)); + assert(math.approxEq(f64, acoshd(123123.234375), 12.414088, epsilon)); +} diff --git a/std/math/asin.zig b/std/math/asin.zig new file mode 100644 index 0000000000..e671477b3c --- /dev/null +++ b/std/math/asin.zig @@ -0,0 +1,157 @@ +const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; + +pub fn asin(x: var) -> @typeOf(x) { + const T = @typeOf(x); + switch (T) { + f32 => @inlineCall(asin32, x), + f64 => @inlineCall(asin64, x), + else => @compileError("asin not implemented for " ++ @typeName(T)), + } +} + +fn r32(z: f32) -> f32 { + const pS0 = 1.6666586697e-01; + const pS1 = -4.2743422091e-02; + const pS2 = -8.6563630030e-03; + const qS1 = -7.0662963390e-01; + + const p = z * (pS0 + z * (pS1 + z * pS2)); + const q = 1.0 + z * qS1; + p / q +} + +fn asin32(x: f32) -> f32 { + const pio2 = 1.570796326794896558e+00; + + const hx: u32 = @bitCast(u32, x); + const ix: u32 = hx & 0x7FFFFFFF; + + // |x| >= 1 + if (ix >= 0x3F800000) { + // |x| >= 1 + if (ix == 0x3F800000) { + return x * pio2 + 0x1.0p-120; // asin(+-1) = +-pi/2 with inexact + } else { + return 0 / (x - x); // asin(|x| > 1) is nan + } + } + + // |x| < 0.5 + if (ix < 0x3F000000) { + // 0x1p-126 <= |x| < 0x1p-12 + if (ix < 0x39800000 and ix >= 0x00800000) { + return x; + } else { + return x + x * r32(x * x); + } + } + + // 1 > |x| >= 0.5 + const z = (1 - math.fabs(x)) * 0.5; + const s = math.sqrt(z); + const fx = pio2 - 2 * (s + s * r32(z)); + + if (hx >> 31 != 0) { + -fx + } else { + fx + } +} + +fn r64(z: f64) -> f64 { + const pS0: f64 = 1.66666666666666657415e-01; + const pS1: f64 = -3.25565818622400915405e-01; + const pS2: f64 = 2.01212532134862925881e-01; + const pS3: f64 = -4.00555345006794114027e-02; + const pS4: f64 = 7.91534994289814532176e-04; + const pS5: f64 = 3.47933107596021167570e-05; + const qS1: f64 = -2.40339491173441421878e+00; + const qS2: f64 = 2.02094576023350569471e+00; + const qS3: f64 = -6.88283971605453293030e-01; + const qS4: f64 = 7.70381505559019352791e-02; + + const p = z * (pS0 + z * (pS1 + z * (pS2 + z * (pS3 + z * (pS4 + z * pS5))))); + const q = 1.0 + z * (qS1 + z * (qS2 + z * (qS3 + z * qS4))); + p / q +} + +fn asin64(x: f64) -> f64 { + const pio2_hi: f64 = 1.57079632679489655800e+00; + const pio2_lo: f64 = 6.12323399573676603587e-17; + + const ux = @bitCast(u64, x); + const hx = u32(ux >> 32); + const ix = hx & 0x7FFFFFFF; + + // |x| >= 1 or nan + if (ix >= 0x3FF00000) { + const lx = u32(ux & 0xFFFFFFFF); + + // asin(1) = +-pi/2 with inexact + if ((ix - 0x3FF00000) | lx == 0) { + return x * pio2_hi + 0x1.0p-120; + } else { + return 0/ (x - x); + } + } + + // |x| < 0.5 + if (ix < 0x3FE00000) { + // if 0x1p-1022 <= |x| < 0x1p-26 avoid raising overflow + if (ix < 0x3E500000 and ix >= 0x00100000) { + return x; + } else { + return x + x * r64(x * x); + } + } + + // 1 > |x| >= 0.5 + const z = (1 - math.fabs(x)) * 0.5; + const s = math.sqrt(z); + const r = r64(z); + var fx: f64 = undefined; + + // |x| > 0.975 + if (ix >= 0x3FEF3333) { + fx = pio2_hi - 2 * (s + s * r) + } else { + const jx = @bitCast(u64, s); + const df = @bitCast(f64, jx & 0xFFFFFFFF00000000); + const c = (z - df * df) / (s + df); + fx = 0.5 * pio2_hi - (2 * s * r - (pio2_lo - 2 * c) - (0.5 * pio2_hi - 2 * df)); + } + + if (hx >> 31 != 0) { + -fx + } else { + fx + } +} + +test "asin" { + assert(asin(f32(0.0)) == asin32(0.0)); + assert(asin(f64(0.0)) == asin64(0.0)); +} + +test "asin32" { + const epsilon = 0.000001; + + assert(math.approxEq(f32, asin32(0.0), 0.0, epsilon)); + assert(math.approxEq(f32, asin32(0.2), 0.201358, epsilon)); + assert(math.approxEq(f32, asin32(-0.2), -0.201358, epsilon)); + assert(math.approxEq(f32, asin32(0.3434), 0.350535, epsilon)); + assert(math.approxEq(f32, asin32(0.5), 0.523599, epsilon)); + assert(math.approxEq(f32, asin32(0.8923), 1.102415, epsilon)); +} + +test "asin64" { + const epsilon = 0.000001; + + assert(math.approxEq(f64, asin64(0.0), 0.0, epsilon)); + assert(math.approxEq(f64, asin64(0.2), 0.201358, epsilon)); + assert(math.approxEq(f64, asin64(-0.2), -0.201358, epsilon)); + assert(math.approxEq(f64, asin64(0.3434), 0.350535, epsilon)); + assert(math.approxEq(f64, asin64(0.5), 0.523599, epsilon)); + assert(math.approxEq(f64, asin64(0.8923), 1.102415, epsilon)); +} diff --git a/std/math/asinh.zig b/std/math/asinh.zig new file mode 100644 index 0000000000..2fa75c9d76 --- /dev/null +++ b/std/math/asinh.zig @@ -0,0 +1,95 @@ +const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; + +pub fn asinh(x: var) -> @typeOf(x) { + const T = @typeOf(x); + switch (T) { + f32 => @inlineCall(asinhf, x), + f64 => @inlineCall(asinhd, x), + else => @compileError("asinh not implemented for " ++ @typeName(T)), + } +} + +// asinh(x) = sign(x) * log(|x| + sqrt(x * x + 1)) ~= x - x^3/6 + o(x^5) +fn asinhf(x: f32) -> f32 { + const u = @bitCast(u32, x); + const i = u & 0x7FFFFFFF; + const s = i >> 31; + + var rx = @bitCast(f32, i); // |x| + + // |x| >= 0x1p12 or inf or nan + if (i >= 0x3F800000 + (12 << 23)) { + rx = math.ln(rx) + 0.69314718055994530941723212145817656; + } + // |x| >= 2 + else if (i >= 0x3F800000 + (1 << 23)) { + rx = math.ln(2 * x + 1 / (math.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)); + } + // |x| < 0x1p-12, inexact if x != 0 + else { + math.forceEval(x + 0x1.0p120); + } + + if (s != 0) -rx else rx +} + +fn asinhd(x: f64) -> f64 { + const u = @bitCast(u64, x); + const e = (u >> 52) & 0x7FF; + const s = u >> 63; + + var rx = @bitCast(f64, u & (@maxValue(u64) >> 1)); // |x| + + // |x| >= 0x1p26 or inf or nan + if (e >= 0x3FF + 26) { + rx = math.ln(rx) + 0.693147180559945309417232121458176568; + } + // |x| >= 2 + else if (e >= 0x3FF + 1) { + rx = math.ln(2 * x + 1 / (math.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)); + } + // |x| < 0x1p-12, inexact if x != 0 + else { + math.forceEval(x + 0x1.0p120); + } + + if (s != 0) -rx else rx +} + +test "asinh" { + assert(asinh(f32(0.0)) == asinhf(0.0)); + assert(asinh(f64(0.0)) == asinhd(0.0)); +} + +test "asinhf" { + const epsilon = 0.000001; + + assert(math.approxEq(f32, asinhf(0.0), 0.0, epsilon)); + assert(math.approxEq(f32, asinhf(0.2), 0.198690, epsilon)); + assert(math.approxEq(f32, asinhf(0.8923), 0.803133, epsilon)); + assert(math.approxEq(f32, asinhf(1.5), 1.194763, epsilon)); + assert(math.approxEq(f32, asinhf(37.45), 4.316332, epsilon)); + assert(math.approxEq(f32, asinhf(89.123), 5.183196, epsilon)); + assert(math.approxEq(f32, asinhf(123123.234375), 12.414088, epsilon)); +} + +test "asinhd" { + const epsilon = 0.000001; + + assert(math.approxEq(f64, asinhd(0.0), 0.0, epsilon)); + assert(math.approxEq(f64, asinhd(0.2), 0.198690, epsilon)); + assert(math.approxEq(f64, asinhd(0.8923), 0.803133, epsilon)); + assert(math.approxEq(f64, asinhd(1.5), 1.194763, epsilon)); + assert(math.approxEq(f64, asinhd(37.45), 4.316332, epsilon)); + assert(math.approxEq(f64, asinhd(89.123), 5.183196, epsilon)); + assert(math.approxEq(f64, asinhd(123123.234375), 12.414088, epsilon)); +} diff --git a/std/math/atan.zig b/std/math/atan.zig new file mode 100644 index 0000000000..094827e000 --- /dev/null +++ b/std/math/atan.zig @@ -0,0 +1,227 @@ +const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; + +pub fn atan(x: var) -> @typeOf(x) { + const T = @typeOf(x); + switch (T) { + f32 => @inlineCall(atan32, x), + f64 => @inlineCall(atan64, x), + else => @compileError("atan not implemented for " ++ @typeName(T)), + } +} + +fn atan32(x_: f32) -> f32 { + const atanhi = []const f32 { + 4.6364760399e-01, // atan(0.5)hi + 7.8539812565e-01, // atan(1.0)hi + 9.8279368877e-01, // atan(1.5)hi + 1.5707962513e+00, // atan(inf)hi + }; + + const atanlo = []const f32 { + 5.0121582440e-09, // atan(0.5)lo + 3.7748947079e-08, // atan(1.0)lo + 3.4473217170e-08, // atan(1.5)lo + 7.5497894159e-08, // atan(inf)lo + }; + + const aT = []const f32 { + 3.3333328366e-01, + -1.9999158382e-01, + 1.4253635705e-01, + -1.0648017377e-01, + 6.1687607318e-02, + }; + + var x = x_; + var ix: u32 = @bitCast(u32, x); + const sign = ix >> 31; + ix &= 0x7FFFFFFF; + + // |x| >= 2^26 + if (ix >= 0x4C800000) { + if (math.isNan(x)) { + return x; + } else { + const z = atanhi[3] + 0x1.0p-120; + return if (sign != 0) -z else z; + } + } + + var id: ?usize = undefined; + + // |x| < 0.4375 + if (ix < 0x3EE00000) { + // |x| < 2^(-12) + if (ix < 0x39800000) { + if (ix < 0x00800000) { + math.forceEval(x * x); + } + return x; + } + id = null; + } else { + x = math.fabs(x); + // |x| < 1.1875 + if (ix < 0x3F980000) { + // 7/16 <= |x| < 11/16 + if (ix < 0x3F300000) { + id = 0; + x = (2.0 * x - 1.0) / (2.0 + x); + } + // 11/16 <= |x| < 19/16 + else { + id = 1; + x = (x - 1.0) / (x + 1.0); + } + } + else { + // |x| < 2.4375 + if (ix < 0x401C0000) { + id = 2; + x = (x - 1.5) / (1.0 + 1.5 * x); + } + // 2.4375 <= |x| < 2^26 + else { + id = 3; + x = -1.0 / x; + } + } + } + + const z = x * x; + const w = z * z; + const s1 = z * (aT[0] + w * (aT[2] + w * aT[4])); + const s2 = w * (aT[1] + w * aT[3]); + + if (id == null) { + x - x * (s1 + s2) + } else { + const zz = atanhi[??id] - ((x * (s1 + s2) - atanlo[??id]) - x); + if (sign != 0) -zz else zz + } +} + +fn atan64(x_: f64) -> f64 { + const atanhi = []const f64 { + 4.63647609000806093515e-01, // atan(0.5)hi + 7.85398163397448278999e-01, // atan(1.0)hi + 9.82793723247329054082e-01, // atan(1.5)hi + 1.57079632679489655800e+00, // atan(inf)hi + }; + + const atanlo = []const f64 { + 2.26987774529616870924e-17, // atan(0.5)lo + 3.06161699786838301793e-17, // atan(1.0)lo + 1.39033110312309984516e-17, // atan(1.5)lo + 6.12323399573676603587e-17, // atan(inf)lo + }; + + const aT = []const f64 { + 3.33333333333329318027e-01, + -1.99999999998764832476e-01, + 1.42857142725034663711e-01, + -1.11111104054623557880e-01, + 9.09088713343650656196e-02, + -7.69187620504482999495e-02, + 6.66107313738753120669e-02, + -5.83357013379057348645e-02, + 4.97687799461593236017e-02, + -3.65315727442169155270e-02, + 1.62858201153657823623e-02, + }; + + var x = x_; + var ux = @bitCast(u64, x); + var ix = u32(ux >> 32); + const sign = ix >> 31; + ix &= 0x7FFFFFFF; + + // |x| >= 2^66 + if (ix >= 0x44100000) { + if (math.isNan(x)) { + return x; + } else { + const z = atanhi[3] + 0x1.0p-120; + return if (sign != 0) -z else z; + } + } + + var id: ?usize = undefined; + + // |x| < 0.4375 + if (ix < 0x3DFC0000) { + // |x| < 2^(-27) + if (ix < 0x3E400000) { + if (ix < 0x00100000) { + math.forceEval(f32(x)); + } + return x; + } + id = null; + } else { + x = math.fabs(x); + // |x| < 1.1875 + if (ix < 0x3FF30000) { + // 7/16 <= |x| < 11/16 + if (ix < 0x3FE60000) { + id = 0; + x = (2.0 * x - 1.0) / (2.0 + x); + } + // 11/16 <= |x| < 19/16 + else { + id = 1; + x = (x - 1.0) / (x + 1.0); + } + } + else { + // |x| < 2.4375 + if (ix < 0x40038000) { + id = 2; + x = (x - 1.5) / (1.0 + 1.5 * x); + } + // 2.4375 <= |x| < 2^66 + else { + id = 3; + x = -1.0 / x; + } + } + } + + const z = x * x; + const w = z * z; + const s1 = z * (aT[0] + w * (aT[2] + w * (aT[4] + w * (aT[6] + w * (aT[8] + w * aT[10]))))); + const s2 = w * (aT[1] + w * (aT[3] + w * (aT[5] + w * (aT[7] + w * aT[9])))); + + if (id == null) { + x - x * (s1 + s2) + } else { + const zz = atanhi[??id] - ((x * (s1 + s2) - atanlo[??id]) - x); + if (sign != 0) -zz else zz + } +} + +test "atan" { + assert(atan(f32(0.2)) == atan32(0.2)); + assert(atan(f64(0.2)) == atan64(0.2)); +} + +test "atan32" { + const epsilon = 0.000001; + + assert(math.approxEq(f32, atan32(0.2), 0.197396, epsilon)); + assert(math.approxEq(f32, atan32(-0.2), -0.197396, epsilon)); + assert(math.approxEq(f32, atan32(0.3434), 0.330783, epsilon)); + assert(math.approxEq(f32, atan32(0.8923), 0.728545, epsilon)); + assert(math.approxEq(f32, atan32(1.5), 0.982794, epsilon)); +} + +test "atan64" { + const epsilon = 0.000001; + + assert(math.approxEq(f64, atan64(0.2), 0.197396, epsilon)); + assert(math.approxEq(f64, atan64(-0.2), -0.197396, epsilon)); + assert(math.approxEq(f64, atan64(0.3434), 0.330783, epsilon)); + assert(math.approxEq(f64, atan64(0.8923), 0.728545, epsilon)); + assert(math.approxEq(f64, atan64(1.5), 0.982794, epsilon)); +} diff --git a/std/math/atan2.zig b/std/math/atan2.zig new file mode 100644 index 0000000000..14bcc60ad2 --- /dev/null +++ b/std/math/atan2.zig @@ -0,0 +1,214 @@ +const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; + +pub fn atan2(comptime T: type, x: T, y: T) -> T { + switch (T) { + f32 => @inlineCall(atan2f, x, y), + f64 => @inlineCall(atan2d, x, y), + else => @compileError("atan2 not implemented for " ++ @typeName(T)), + } +} + +fn atan2f(y: f32, x: f32) -> f32 { + const pi: f32 = 3.1415927410e+00; + const pi_lo: f32 = -8.7422776573e-08; + + if (math.isNan(x) or math.isNan(y)) { + return x + y; + } + + var ix = @bitCast(u32, x); + var iy = @bitCast(u32, y); + + // x = 1.0 + if (ix == 0x3F800000) { + return math.atan(y); + } + + // 2 * sign(x) + sign(y) + const m = ((iy >> 31) & 1) | ((ix >> 30) & 2); + ix &= 0x7FFFFFFF; + iy &= 0x7FFFFFFF; + + if (iy == 0) { + switch (m) { + 0, 1 => return y, // atan(+-0, +...) + 2 => return pi, // atan(+0, -...) + 3 => return -pi, // atan(-0, -...) + else => unreachable, + } + } + + if (ix == 0) { + if (m & 1 != 0) { + return -pi / 2; + } else { + return pi / 2; + } + } + + if (ix == 0x7F800000) { + if (iy == 0x7F800000) { + switch (m) { + 0 => return pi / 4, // atan(+inf, +inf) + 1 => return -pi / 4, // atan(-inf, +inf) + 2 => return 3*pi / 4, // atan(+inf, -inf) + 3 => return -3*pi / 4, // atan(-inf, -inf) + else => unreachable, + } + } else { + switch (m) { + 0 => return 0.0, // atan(+..., +inf) + 1 => return -0.0, // atan(-..., +inf) + 2 => return pi, // atan(+..., -inf) + 3 => return -pi, // atan(-...f, -inf) + else => unreachable, + } + } + } + + // |y / x| > 0x1p26 + if (ix + (26 << 23) < iy or iy == 0x7F800000) { + if (m & 1 != 0) { + return -pi / 2; + } else { + return pi / 2; + } + } + + // z = atan(|y / x|) with correct underflow + var z = { + if ((m & 2) != 0 and iy + (26 << 23) < ix) { + 0.0 + } else { + math.atan(math.fabs(y / x)) + } + }; + + switch (m) { + 0 => return z, // atan(+, +) + 1 => return -z, // atan(-, +) + 2 => return pi - (z - pi_lo), // atan(+, -) + 3 => return (z - pi_lo) - pi, // atan(-, -) + else => unreachable, + } +} + +fn atan2d(y: f64, x: f64) -> f64 { + const pi: f64 = 3.1415926535897931160E+00; + const pi_lo: f64 = 1.2246467991473531772E-16; + + if (math.isNan(x) or math.isNan(y)) { + return x + y; + } + + var ux = @bitCast(u64, x); + var ix = u32(ux >> 32); + var lx = u32(ux & 0xFFFFFFFF); + + var uy = @bitCast(u64, y); + var iy = u32(uy >> 32); + var ly = u32(uy & 0xFFFFFFFF); + + // x = 1.0 + if ((ix -% 0x3FF00000) | lx == 0) { + return math.atan(y); + } + + // 2 * sign(x) + sign(y) + const m = ((iy >> 31) & 1) | ((ix >> 30) & 2); + ix &= 0x7FFFFFFF; + iy &= 0x7FFFFFFF; + + if (iy | ly == 0) { + switch (m) { + 0, 1 => return y, // atan(+-0, +...) + 2 => return pi, // atan(+0, -...) + 3 => return -pi, // atan(-0, -...) + else => unreachable, + } + } + + if (ix | lx == 0) { + if (m & 1 != 0) { + return -pi / 2; + } else { + return pi / 2; + } + } + + if (ix == 0x7FF00000) { + if (iy == 0x7FF00000) { + switch (m) { + 0 => return pi / 4, // atan(+inf, +inf) + 1 => return -pi / 4, // atan(-inf, +inf) + 2 => return 3*pi / 4, // atan(+inf, -inf) + 3 => return -3*pi / 4, // atan(-inf, -inf) + else => unreachable, + } + } else { + switch (m) { + 0 => return 0.0, // atan(+..., +inf) + 1 => return -0.0, // atan(-..., +inf) + 2 => return pi, // atan(+..., -inf) + 3 => return -pi, // atan(-...f, -inf) + else => unreachable, + } + } + } + + // |y / x| > 0x1p64 + if (ix +% (64 << 20) < iy or iy == 0x7FF00000) { + if (m & 1 != 0) { + return -pi / 2; + } else { + return pi / 2; + } + } + + // z = atan(|y / x|) with correct underflow + var z = { + if ((m & 2) != 0 and iy +% (64 << 20) < ix) { + 0.0 + } else { + math.atan(math.fabs(y / x)) + } + }; + + switch (m) { + 0 => return z, // atan(+, +) + 1 => return -z, // atan(-, +) + 2 => return pi - (z - pi_lo), // atan(+, -) + 3 => return (z - pi_lo) - pi, // atan(-, -) + else => unreachable, + } +} + +test "atan2" { + assert(atan2(f32, 0.2, 0.21) == atan2f(0.2, 0.21)); + assert(atan2(f64, 0.2, 0.21) == atan2d(0.2, 0.21)); +} + +test "atan2f" { + const epsilon = 0.000001; + + assert(math.approxEq(f32, atan2f(0.0, 0.0), 0.0, epsilon)); + assert(math.approxEq(f32, atan2f(0.2, 0.2), 0.785398, epsilon)); + assert(math.approxEq(f32, atan2f(-0.2, 0.2), -0.785398, epsilon)); + assert(math.approxEq(f32, atan2f(0.2, -0.2), 2.356194, epsilon)); + assert(math.approxEq(f32, atan2f(-0.2, -0.2), -2.356194, epsilon)); + assert(math.approxEq(f32, atan2f(0.34, -0.4), 2.437099, epsilon)); + assert(math.approxEq(f32, atan2f(0.34, 1.243), 0.267001, epsilon)); +} + +test "atan2d" { + const epsilon = 0.000001; + + assert(math.approxEq(f64, atan2d(0.0, 0.0), 0.0, epsilon)); + assert(math.approxEq(f64, atan2d(0.2, 0.2), 0.785398, epsilon)); + assert(math.approxEq(f64, atan2d(-0.2, 0.2), -0.785398, epsilon)); + assert(math.approxEq(f64, atan2d(0.2, -0.2), 2.356194, epsilon)); + assert(math.approxEq(f64, atan2d(-0.2, -0.2), -2.356194, epsilon)); + assert(math.approxEq(f64, atan2d(0.34, -0.4), 2.437099, epsilon)); + assert(math.approxEq(f64, atan2d(0.34, 1.243), 0.267001, epsilon)); +} diff --git a/std/math/atanh.zig b/std/math/atanh.zig new file mode 100644 index 0000000000..9098616c53 --- /dev/null +++ b/std/math/atanh.zig @@ -0,0 +1,85 @@ +const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; + +pub fn atanh(x: var) -> @typeOf(x) { + const T = @typeOf(x); + switch (T) { + f32 => @inlineCall(atanhf, x), + f64 => @inlineCall(atanhd, x), + else => @compileError("atanh not implemented for " ++ @typeName(T)), + } +} + +// atanh(x) = log((1 + x) / (1 - x)) / 2 = log1p(2x / (1 - x)) / 2 ~= x + x^3 / 3 + o(x^5) +fn atanhf(x: f32) -> f32 { + const u = @bitCast(u32, x); + const i = u & 0x7FFFFFFF; + const s = u >> 31; + + var y = @bitCast(f32, i); // |x| + + if (u < 0x3F800000 - (1 << 23)) { + if (u < 0x3F800000 - (32 << 23)) { + // underflow + if (u < (1 << 23)) { + math.forceEval(y * y) + } + } + // |x| < 0.5 + else { + y = 0.5 * math.log1p(2 * y + 2 * y * y / (1 - y)); + } + } else { + // avoid overflow + y = 0.5 * math.log1p(2 * (y / (1 - y))); + } + + if (s != 0) -y else y +} + +fn atanhd(x: f64) -> f64 { + const u = @bitCast(u64, x); + const e = (u >> 52) & 0x7FF; + const s = u >> 63; + + var y = @bitCast(f64, u & (@maxValue(u64) >> 1)); // |x| + + if (e < 0x3FF - 1) { + if (e < 0x3FF - 32) { + // underflow + if (e == 0) { + math.forceEval(f32(y)); + } + } + // |x| < 0.5 + else { + y = 0.5 * math.log1p(2 * y + 2 * y * y / (1 - y)); + } + } else { + // avoid overflow + y = 0.5 * math.log1p(2 * (y / (1 - y))); + } + + if (s != 0) -y else y +} + +test "atanh" { + assert(atanh(f32(0.0)) == atanhf(0.0)); + assert(atanh(f64(0.0)) == atanhd(0.0)); +} + +test "atanhf" { + const epsilon = 0.000001; + + assert(math.approxEq(f32, atanhf(0.0), 0.0, epsilon)); + assert(math.approxEq(f32, atanhf(0.2), 0.202733, epsilon)); + assert(math.approxEq(f32, atanhf(0.8923), 1.433099, epsilon)); +} + +test "atanhd" { + const epsilon = 0.000001; + + assert(math.approxEq(f64, atanhd(0.0), 0.0, epsilon)); + assert(math.approxEq(f64, atanhd(0.2), 0.202733, epsilon)); + assert(math.approxEq(f64, atanhd(0.8923), 1.433099, epsilon)); +} diff --git a/std/math/cbrt.zig b/std/math/cbrt.zig new file mode 100644 index 0000000000..e286f94ab4 --- /dev/null +++ b/std/math/cbrt.zig @@ -0,0 +1,134 @@ +const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; + +pub fn cbrt(x: var) -> @typeOf(x) { + const T = @typeOf(x); + switch (T) { + f32 => @inlineCall(cbrt32, x), + f64 => @inlineCall(cbrt64, x), + else => @compileError("cbrt not implemented for " ++ @typeName(T)), + } +} + +fn cbrt32(x: f32) -> f32 { + const B1: u32 = 709958130; // (127 - 127.0 / 3 - 0.03306235651) * 2^23 + const B2: u32 = 642849266; // (127 - 127.0 / 3 - 24 / 3 - 0.03306235651) * 2^23 + + var u = @bitCast(u32, x); + var hx = u & 0x7FFFFFFF; + + // cbrt(nan, inf) = itself + if (hx >= 0x7F800000) { + return x + x; + } + + // cbrt to ~5bits + if (hx < 0x00800000) { + // cbrt(+-0) = itself + if (hx == 0) { + return x; + } + u = @bitCast(u32, x * 0x1.0p24); + hx = u & 0x7FFFFFFF; + hx = hx / 3 + B2; + } else { + hx = hx / 3 + B1; + } + + u &= 0x80000000; + u |= hx; + + // first step newton to 16 bits + var t: f64 = @bitCast(f32, u); + var r: f64 = t * t * t; + t = t * (f64(x) + x + r) / (x + r + r); + + // second step newton to 47 bits + r = t * t * t; + t = t * (f64(x) + x + r) / (x + r + r); + + f32(t) +} + +fn cbrt64(x: f64) -> f64 { + const B1: u32 = 715094163; // (1023 - 1023 / 3 - 0.03306235651 * 2^20 + const B2: u32 = 696219795; // (1023 - 1023 / 3 - 54 / 3 - 0.03306235651 * 2^20 + + // |1 / cbrt(x) - p(x)| < 2^(23.5) + const P0: f64 = 1.87595182427177009643; + const P1: f64 = -1.88497979543377169875; + const P2: f64 = 1.621429720105354466140; + const P3: f64 = -0.758397934778766047437; + const P4: f64 = 0.145996192886612446982; + + var u = @bitCast(u64, x); + var hx = u32(u >> 32) & 0x7FFFFFFF; + + // cbrt(nan, inf) = itself + if (hx >= 0x7FF00000) { + return x + x; + } + + // cbrt to ~5bits + if (hx < 0x00100000) { + u = @bitCast(u64, x * 0x1.0p54); + hx = u32(u >> 32) & 0x7FFFFFFF; + + // cbrt(0) is itself + if (hx == 0) { + return 0; + } + hx = hx / 3 + B2; + } else { + hx = hx / 3 + B1; + } + + u &= 1 << 63; + u |= u64(hx) << 32; + var t = @bitCast(f64, u); + + // cbrt to 23 bits + // cbrt(x) = t * cbrt(x / t^3) ~= t * P(t^3 / x) + var r = (t * t) * (t / x); + t = t * ((P0 + r * (P1 + r * P2)) + ((r * r) * r) * (P3 + r * P4)); + + // Round t away from 0 to 23 bits + u = @bitCast(u64, t); + u = (u + 0x80000000) & 0xFFFFFFFFC0000000; + t = @bitCast(f64, u); + + // one step newton to 53 bits + const s = t * t; + var q = x / s; + var w = t + t; + q = (q - t) / (w + q); + + t + t * q +} + +test "cbrt" { + assert(cbrt(f32(0.0)) == cbrt32(0.0)); + assert(cbrt(f64(0.0)) == cbrt64(0.0)); +} + +test "cbrt32" { + const epsilon = 0.000001; + + assert(cbrt32(0.0) == 0.0); + assert(math.approxEq(f32, cbrt32(0.2), 0.584804, epsilon)); + assert(math.approxEq(f32, cbrt32(0.8923), 0.962728, epsilon)); + assert(math.approxEq(f32, cbrt32(1.5), 1.144714, epsilon)); + assert(math.approxEq(f32, cbrt32(37.45), 3.345676, epsilon)); + assert(math.approxEq(f32, cbrt32(123123.234375), 49.748501, epsilon)); +} + +test "cbrt64" { + const epsilon = 0.000001; + + assert(cbrt64(0.0) == 0.0); + assert(math.approxEq(f64, cbrt64(0.2), 0.584804, epsilon)); + assert(math.approxEq(f64, cbrt64(0.8923), 0.962728, epsilon)); + assert(math.approxEq(f64, cbrt64(1.5), 1.144714, epsilon)); + assert(math.approxEq(f64, cbrt64(37.45), 3.345676, epsilon)); + assert(math.approxEq(f64, cbrt64(123123.234375), 49.748501, epsilon)); +} diff --git a/std/math/ceil.zig b/std/math/ceil.zig new file mode 100644 index 0000000000..f80e493f1f --- /dev/null +++ b/std/math/ceil.zig @@ -0,0 +1,89 @@ +const builtin = @import("builtin"); +const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; + +pub fn ceil(x: var) -> @typeOf(x) { + const T = @typeOf(x); + switch (T) { + f32 => @inlineCall(ceil32, x), + f64 => @inlineCall(ceil64, x), + else => @compileError("ceil not implemented for " ++ @typeName(T)), + } +} + +fn ceil32(x: f32) -> f32 { + var u = @bitCast(u32, x); + var e = i32((u >> 23) & 0xFF) - 0x7F; + var m: u32 = undefined; + + if (e >= 23) { + return x; + } + else if (e >= 0) { + m = 0x007FFFFF >> u32(e); + if (u & m == 0) { + return x; + } + math.forceEval(x + 0x1.0p120); + if (u >> 31 == 0) { + u += m; + } + u &= ~m; + @bitCast(f32, u) + } else { + math.forceEval(x + 0x1.0p120); + if (u >> 31 != 0) { + return -0.0; + } else { + 1.0 + } + } +} + +fn ceil64(x: f64) -> f64 { + const u = @bitCast(u64, x); + const e = (u >> 52) & 0x7FF; + var y: f64 = undefined; + + if (e >= 0x3FF+52 or x == 0) { + return x; + } + + if (u >> 63 != 0) { + @setFloatMode(this, builtin.FloatMode.Strict); + y = x - math.f64_toint + math.f64_toint - x; + } else { + @setFloatMode(this, builtin.FloatMode.Strict); + y = x + math.f64_toint - math.f64_toint - x; + } + + if (e <= 0x3FF-1) { + math.forceEval(y); + if (u >> 63 != 0) { + return -0.0; // Compiler requires return. + } else { + 1.0 + } + } else if (y < 0) { + x + y + 1 + } else { + x + y + } +} + +test "ceil" { + assert(ceil(f32(0.0)) == ceil32(0.0)); + assert(ceil(f64(0.0)) == ceil64(0.0)); +} + +test "ceil32" { + assert(ceil32(1.3) == 2.0); + assert(ceil32(-1.3) == -1.0); + assert(ceil32(0.2) == 1.0); +} + +test "ceil64" { + assert(ceil64(1.3) == 2.0); + assert(ceil64(-1.3) == -1.0); + assert(ceil64(0.2) == 1.0); +} diff --git a/std/math/copysign.zig b/std/math/copysign.zig new file mode 100644 index 0000000000..7043cc4852 --- /dev/null +++ b/std/math/copysign.zig @@ -0,0 +1,47 @@ +const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; + +pub fn copysign(comptime T: type, x: T, y: T) -> T { + switch (T) { + f32 => @inlineCall(copysign32, x, y), + f64 => @inlineCall(copysign64, x, y), + else => @compileError("copysign not implemented for " ++ @typeName(T)), + } +} + +fn copysign32(x: f32, y: f32) -> f32 { + const ux = @bitCast(u32, x); + const uy = @bitCast(u32, y); + + const h1 = ux & (@maxValue(u32) / 2); + const h2 = uy & (u32(1) << 31); + @bitCast(f32, h1 | h2) +} + +fn copysign64(x: f64, y: f64) -> f64 { + const ux = @bitCast(u64, x); + const uy = @bitCast(u64, y); + + const h1 = ux & (@maxValue(u64) / 2); + const h2 = uy & (u64(1) << 63); + @bitCast(f64, h1 | h2) +} + +test "copysign" { + assert(copysign(f32, 1.0, 1.0) == copysign32(1.0, 1.0)); + assert(copysign(f64, 1.0, 1.0) == copysign64(1.0, 1.0)); +} + +test "copysign32" { + assert(copysign32(5.0, 1.0) == 5.0); + assert(copysign32(5.0, -1.0) == -5.0); + assert(copysign32(-5.0, -1.0) == -5.0); + assert(copysign32(-5.0, 1.0) == 5.0); +} + +test "copysign64" { + assert(copysign64(5.0, 1.0) == 5.0); + assert(copysign64(5.0, -1.0) == -5.0); + assert(copysign64(-5.0, -1.0) == -5.0); + assert(copysign64(-5.0, 1.0) == 5.0); +} diff --git a/std/math/cos.zig b/std/math/cos.zig new file mode 100644 index 0000000000..5c90853297 --- /dev/null +++ b/std/math/cos.zig @@ -0,0 +1,159 @@ +const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; + +pub fn cos(x: var) -> @typeOf(x) { + const T = @typeOf(x); + switch (T) { + f32 => @inlineCall(cos32, x), + f64 => @inlineCall(cos64, x), + else => @compileError("cos not implemented for " ++ @typeName(T)), + } +} + +// sin polynomial coefficients +const S0 = 1.58962301576546568060E-10; +const S1 = -2.50507477628578072866E-8; +const S2 = 2.75573136213857245213E-6; +const S3 = -1.98412698295895385996E-4; +const S4 = 8.33333333332211858878E-3; +const S5 = -1.66666666666666307295E-1; + +// cos polynomial coeffiecients +const C0 = -1.13585365213876817300E-11; +const C1 = 2.08757008419747316778E-9; +const C2 = -2.75573141792967388112E-7; +const C3 = 2.48015872888517045348E-5; +const C4 = -1.38888888888730564116E-3; +const C5 = 4.16666666666665929218E-2; + +// NOTE: This is taken from the go stdlib. The musl implementation is much more complex. +// +// This may have slight differences on some edge cases and may need to replaced if so. +fn cos32(x_: f32) -> f32 { + const pi4a = 7.85398125648498535156e-1; + const pi4b = 3.77489470793079817668E-8; + const pi4c = 2.69515142907905952645E-15; + const m4pi = 1.273239544735162542821171882678754627704620361328125; + + var x = x_; + if (math.isNan(x) or math.isInf(x)) { + return math.nan(f32); + } + + var sign = false; + if (x < 0) { + x = -x; + } + + var y = math.floor(x * m4pi); + var j = i64(y); + + if (j & 1 == 1) { + j += 1; + y += 1; + } + + j &= 7; + if (j > 3) { + j -= 4; + sign = !sign; + } + if (j > 1) { + sign = !sign; + } + + const z = ((x - y * pi4a) - y * pi4b) - y * pi4c; + const w = z * z; + + const r = { + if (j == 1 or j == 2) { + z + z * w * (S5 + w * (S4 + w * (S3 + w * (S2 + w * (S1 + w * S0))))) + } else { + 1.0 - 0.5 * w + w * w * (C5 + w * (C4 + w * (C3 + w * (C2 + w * (C1 + w * C0))))) + } + }; + + if (sign) { + -r + } else { + r + } +} + +fn cos64(x_: f64) -> f64 { + const pi4a = 7.85398125648498535156e-1; + const pi4b = 3.77489470793079817668E-8; + const pi4c = 2.69515142907905952645E-15; + const m4pi = 1.273239544735162542821171882678754627704620361328125; + + var x = x_; + if (math.isNan(x) or math.isInf(x)) { + return math.nan(f64); + } + + var sign = false; + if (x < 0) { + x = -x; + } + + var y = math.floor(x * m4pi); + var j = i64(y); + + if (j & 1 == 1) { + j += 1; + y += 1; + } + + j &= 7; + if (j > 3) { + j -= 4; + sign = !sign; + } + if (j > 1) { + sign = !sign; + } + + const z = ((x - y * pi4a) - y * pi4b) - y * pi4c; + const w = z * z; + + const r = { + if (j == 1 or j == 2) { + z + z * w * (S5 + w * (S4 + w * (S3 + w * (S2 + w * (S1 + w * S0))))) + } else { + 1.0 - 0.5 * w + w * w * (C5 + w * (C4 + w * (C3 + w * (C2 + w * (C1 + w * C0))))) + } + }; + + if (sign) { + -r + } else { + r + } +} + +test "cos" { + assert(cos(f32(0.0)) == cos32(0.0)); + assert(cos(f64(0.0)) == cos64(0.0)); +} + +test "cos32" { + const epsilon = 0.000001; + + assert(math.approxEq(f32, cos32(0.0), 1.0, epsilon)); + assert(math.approxEq(f32, cos32(0.2), 0.980067, epsilon)); + assert(math.approxEq(f32, cos32(0.8923), 0.627623, epsilon)); + assert(math.approxEq(f32, cos32(1.5), 0.070737, epsilon)); + assert(math.approxEq(f32, cos32(37.45), 0.969132, epsilon)); + assert(math.approxEq(f32, cos32(89.123), 0.400798, epsilon)); +} + +test "cos64" { + const epsilon = 0.000001; + + assert(math.approxEq(f64, cos64(0.0), 1.0, epsilon)); + assert(math.approxEq(f64, cos64(0.2), 0.980067, epsilon)); + assert(math.approxEq(f64, cos64(0.8923), 0.627623, epsilon)); + assert(math.approxEq(f64, cos64(1.5), 0.070737, epsilon)); + assert(math.approxEq(f64, cos64(37.45), 0.969132, epsilon)); + assert(math.approxEq(f64, cos64(89.123), 0.40080, epsilon)); +} diff --git a/std/math/cosh.zig b/std/math/cosh.zig new file mode 100644 index 0000000000..60bff79b2b --- /dev/null +++ b/std/math/cosh.zig @@ -0,0 +1,91 @@ +const math = @import("index.zig"); +const expo2 = @import("_expo2.zig").expo2; +const assert = @import("../debug.zig").assert; + +pub fn cosh(x: var) -> @typeOf(x) { + const T = @typeOf(x); + switch (T) { + f32 => @inlineCall(coshf, x), + f64 => @inlineCall(coshd, x), + else => @compileError("cosh not implemented for " ++ @typeName(T)), + } +} + +// cosh(x) = (exp(x) + 1 / exp(x)) / 2 +// = 1 + 0.5 * (exp(x) - 1) * (exp(x) - 1) / exp(x) +// = 1 + (x * x) / 2 + o(x^4) +fn coshf(x: f32) -> f32 { + const u = @bitCast(u32, x); + const ux = u & 0x7FFFFFFF; + const ax = @bitCast(f32, ux); + + // |x| < log(2) + if (ux < 0x3F317217) { + if (ux < 0x3F800000 - (12 << 23)) { + math.raiseOverflow(); + return 1.0; + } + const t = math.expm1(ax); + return 1 + t * t / (2 * (1 + t)); + } + + // |x| < log(FLT_MAX) + if (ux < 0x42B17217) { + const t = math.exp(ax); + return 0.5 * (t + 1 / t); + } + + // |x| > log(FLT_MAX) or nan + expo2(ax) +} + +fn coshd(x: f64) -> f64 { + const u = @bitCast(u64, x); + const w = u32(u >> 32); + const ax = @bitCast(f64, u & (@maxValue(u64) >> 1)); + + // |x| < log(2) + if (w < 0x3FE62E42) { + if (w < 0x3FF00000 - (26 << 20)) { + if (x != 0) { + math.raiseInexact(); + } + return 1.0; + } + const t = math.expm1(ax); + return 1 + t * t / (2 * (1 + t)); + } + + // |x| < log(DBL_MAX) + if (w < 0x40862E42) { + const t = math.exp(ax); + // NOTE: If x > log(0x1p26) then 1/t is not required. + return 0.5 * (t + 1 / t); + } + + // |x| > log(CBL_MAX) or nan + expo2(ax) +} + +test "cosh" { + assert(cosh(f32(1.5)) == coshf(1.5)); + assert(cosh(f64(1.5)) == coshd(1.5)); +} + +test "coshf" { + const epsilon = 0.000001; + + assert(math.approxEq(f32, coshf(0.0), 1.0, epsilon)); + assert(math.approxEq(f32, coshf(0.2), 1.020067, epsilon)); + assert(math.approxEq(f32, coshf(0.8923), 1.425225, epsilon)); + assert(math.approxEq(f32, coshf(1.5), 2.352410, epsilon)); +} + +test "coshd" { + const epsilon = 0.000001; + + assert(math.approxEq(f64, coshd(0.0), 1.0, epsilon)); + assert(math.approxEq(f64, coshd(0.2), 1.020067, epsilon)); + assert(math.approxEq(f64, coshd(0.8923), 1.425225, epsilon)); + assert(math.approxEq(f64, coshd(1.5), 2.352410, epsilon)); +} diff --git a/std/math/exp.zig b/std/math/exp.zig new file mode 100644 index 0000000000..53cac3facd --- /dev/null +++ b/std/math/exp.zig @@ -0,0 +1,191 @@ +const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; + +pub fn exp(x: var) -> @typeOf(x) { + const T = @typeOf(x); + switch (T) { + f32 => @inlineCall(exp32, x), + f64 => @inlineCall(exp64, x), + else => @compileError("exp not implemented for " ++ @typeName(T)), + } +} + +fn exp32(x_: f32) -> f32 { + const half = []const f32 { 0.5, -0.5 }; + const ln2hi = 6.9314575195e-1; + const ln2lo = 1.4286067653e-6; + const invln2 = 1.4426950216e+0; + const P1 = 1.6666625440e-1; + const P2 = -2.7667332906e-3; + + var x = x_; + var hx = @bitCast(u32, x); + const sign = i32(hx >> 31); + hx &= 0x7FFFFFFF; + + // |x| >= -87.33655 or nan + if (hx >= 0x42AEAC50) { + // nan + if (hx > 0x7F800000) { + return x; + } + // x >= 88.722839 + if (hx >= 0x42b17218 and sign == 0) { + return x * 0x1.0p127; + } + if (sign != 0) { + math.forceEval(-0x1.0p-149 / x); // overflow + // x <= -103.972084 + if (hx >= 0x42CFF1B5) { + return 0; + } + } + } + + var k: i32 = undefined; + var hi: f32 = undefined; + var lo: f32 = undefined; + + // |x| > 0.5 * ln2 + if (hx > 0x3EB17218) { + // |x| > 1.5 * ln2 + if (hx > 0x3F851592) { + k = i32(invln2 * x + half[usize(sign)]); + } + else { + k = 1 - sign - sign; + } + + const fk = f32(k); + hi = x - fk * ln2hi; + lo = fk * ln2lo; + x = hi - lo; + } + // |x| > 2^(-14) + else if (hx > 0x39000000) { + k = 0; + hi = x; + lo = 0; + } + else { + math.forceEval(0x1.0p127 + x); // inexact + return 1 + x; + } + + const xx = x * x; + const c = x - xx * (P1 + xx * P2); + const y = 1 + (x * c / (2 - c) - lo + hi); + + if (k == 0) { + y + } else { + math.scalbn(y, k) + } +} + +fn exp64(x_: f64) -> f64 { + const half = []const f64 { 0.5, -0.5 }; + const ln2hi: f64 = 6.93147180369123816490e-01; + const ln2lo: f64 = 1.90821492927058770002e-10; + const invln2: f64 = 1.44269504088896338700e+00; + const P1: f64 = 1.66666666666666019037e-01; + const P2: f64 = -2.77777777770155933842e-03; + const P3: f64 = 6.61375632143793436117e-05; + const P4: f64 = -1.65339022054652515390e-06; + const P5: f64 = 4.13813679705723846039e-08; + + var x = x_; + var ux = @bitCast(u64, x); + var hx = ux >> 32; + const sign = i32(hx >> 31); + hx &= 0x7FFFFFFF; + + // |x| >= 708.39 or nan + if (hx >= 0x4086232B) { + // nan + if (hx > 0x7FF00000) { + return x; + } + if (x > 709.782712893383973096) { + // overflow if x != inf + if (!math.isInf(x)) { + math.raiseOverflow(); + } + return math.inf(f64); + } + if (x < -708.39641853226410622) { + // underflow if x != -inf + // math.forceEval(f32(-0x1.0p-149 / x)); + if (x < -745.13321910194110842) { + return 0; + } + } + } + + // argument reduction + var k: i32 = undefined; + var hi: f64 = undefined; + var lo: f64 = undefined; + + // |x| > 0.5 * ln2 + if (hx > 0x3EB17218) { + // |x| >= 1.5 * ln2 + if (hx > 0x3FF0A2B2) { + k = i32(invln2 * x + half[usize(sign)]); + } + else { + k = 1 - sign - sign; + } + + const dk = f64(k); + hi = x - dk * ln2hi; + lo = dk * ln2lo; + x = hi - lo; + } + // |x| > 2^(-28) + else if (hx > 0x3E300000) { + k = 0; + hi = x; + lo = 0; + } + else { + // inexact if x != 0 + // math.forceEval(0x1.0p1023 + x); + return 1 + x; + } + + const xx = x * x; + const c = x - xx * (P1 + xx * (P2 + xx * (P3 + xx * (P4 + xx * P5)))); + const y = 1 + (x * c / (2 - c) - lo + hi); + + if (k == 0) { + y + } else { + math.scalbn(y, k) + } +} + +test "exp" { + assert(exp(f32(0.0)) == exp32(0.0)); + assert(exp(f64(0.0)) == exp64(0.0)); +} + +test "exp32" { + const epsilon = 0.000001; + + assert(exp32(0.0) == 1.0); + assert(math.approxEq(f32, exp32(0.0), 1.0, epsilon)); + assert(math.approxEq(f32, exp32(0.2), 1.221403, epsilon)); + assert(math.approxEq(f32, exp32(0.8923), 2.440737, epsilon)); + assert(math.approxEq(f32, exp32(1.5), 4.481689, epsilon)); +} + +test "exp64" { + const epsilon = 0.000001; + + assert(exp64(0.0) == 1.0); + assert(math.approxEq(f64, exp64(0.0), 1.0, epsilon)); + assert(math.approxEq(f64, exp64(0.2), 1.221403, epsilon)); + assert(math.approxEq(f64, exp64(0.8923), 2.440737, epsilon)); + assert(math.approxEq(f64, exp64(1.5), 4.481689, epsilon)); +} diff --git a/std/math/exp2.zig b/std/math/exp2.zig new file mode 100644 index 0000000000..3968357251 --- /dev/null +++ b/std/math/exp2.zig @@ -0,0 +1,429 @@ +const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; + +pub fn exp2(x: var) -> @typeOf(x) { + const T = @typeOf(x); + switch (T) { + f32 => @inlineCall(exp2f, x), + f64 => @inlineCall(exp2d, x), + else => @compileError("exp2 not implemented for " ++ @typeName(T)), + } +} + +const exp2ft = []const 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 exp2f(x: f32) -> f32 { + const tblsiz = u32(exp2ft.len); + const redux: f32 = 0x1.8p23 / f32(tblsiz); + const P1: f32 = 0x1.62e430p-1; + const P2: f32 = 0x1.ebfbe0p-3; + const P3: f32 = 0x1.c6b348p-5; + const P4: f32 = 0x1.3b2c9cp-7; + + var u = @bitCast(u32, x); + const ix = u & 0x7FFFFFFF; + + // |x| > 126 + if (ix > 0x42FC0000) { + // nan + if (ix > 0x7F800000) { + return x; + } + // x >= 128 + if (u >= 0x43000000 and u < 0x80000000) { + return x * 0x1.0p127; + } + // x < -126 + if (u >= 0x80000000) { + if (u >= 0xC3160000 or u & 0x000FFFF != 0) { + math.forceEval(-0x1.0p-149 / x); + } + // x <= -150 + if (u >= 0x3160000) { + return 0; + } + } + } + // |x| <= 0x1p-25 + else if (ix <= 0x33000000) { + return 1.0 + x; + } + + var uf = x + redux; + var i0 = @bitCast(u32, uf); + i0 += tblsiz / 2; + + const k = i0 / tblsiz; + // NOTE: musl relies on undefined overflow shift behaviour. Appears that this produces the + // intended result but should confirm how GCC/Clang handle this to ensure. + const uk = @bitCast(f64, u64(0x3FF + k) <<% 52); + i0 &= tblsiz - 1; + uf -= redux; + + const z: f64 = x - uf; + var r: f64 = exp2ft[i0]; + const t: f64 = r * z; + r = r + t * (P1 + z * P2) + t * (z * z) * (P3 + z * P4); + f32(r * uk) +} + +const exp2dt = []f64 { + // exp2(z + eps) eps + 0x1.6a09e667f3d5dp-1, 0x1.9880p-44, + 0x1.6b052fa751744p-1, 0x1.8000p-50, + 0x1.6c012750bd9fep-1, -0x1.8780p-45, + 0x1.6cfdcddd476bfp-1, 0x1.ec00p-46, + 0x1.6dfb23c651a29p-1, -0x1.8000p-50, + 0x1.6ef9298593ae3p-1, -0x1.c000p-52, + 0x1.6ff7df9519386p-1, -0x1.fd80p-45, + 0x1.70f7466f42da3p-1, -0x1.c880p-45, + 0x1.71f75e8ec5fc3p-1, 0x1.3c00p-46, + 0x1.72f8286eacf05p-1, -0x1.8300p-44, + 0x1.73f9a48a58152p-1, -0x1.0c00p-47, + 0x1.74fbd35d7ccfcp-1, 0x1.f880p-45, + 0x1.75feb564267f1p-1, 0x1.3e00p-47, + 0x1.77024b1ab6d48p-1, -0x1.7d00p-45, + 0x1.780694fde5d38p-1, -0x1.d000p-50, + 0x1.790b938ac1d00p-1, 0x1.3000p-49, + 0x1.7a11473eb0178p-1, -0x1.d000p-49, + 0x1.7b17b0976d060p-1, 0x1.0400p-45, + 0x1.7c1ed0130c133p-1, 0x1.0000p-53, + 0x1.7d26a62ff8636p-1, -0x1.6900p-45, + 0x1.7e2f336cf4e3bp-1, -0x1.2e00p-47, + 0x1.7f3878491c3e8p-1, -0x1.4580p-45, + 0x1.80427543e1b4ep-1, 0x1.3000p-44, + 0x1.814d2add1071ap-1, 0x1.f000p-47, + 0x1.82589994ccd7ep-1, -0x1.1c00p-45, + 0x1.8364c1eb942d0p-1, 0x1.9d00p-45, + 0x1.8471a4623cab5p-1, 0x1.7100p-43, + 0x1.857f4179f5bbcp-1, 0x1.2600p-45, + 0x1.868d99b4491afp-1, -0x1.2c40p-44, + 0x1.879cad931a395p-1, -0x1.3000p-45, + 0x1.88ac7d98a65b8p-1, -0x1.a800p-45, + 0x1.89bd0a4785800p-1, -0x1.d000p-49, + 0x1.8ace5422aa223p-1, 0x1.3280p-44, + 0x1.8be05bad619fap-1, 0x1.2b40p-43, + 0x1.8cf3216b54383p-1, -0x1.ed00p-45, + 0x1.8e06a5e08664cp-1, -0x1.0500p-45, + 0x1.8f1ae99157807p-1, 0x1.8280p-45, + 0x1.902fed0282c0ep-1, -0x1.cb00p-46, + 0x1.9145b0b91ff96p-1, -0x1.5e00p-47, + 0x1.925c353aa2ff9p-1, 0x1.5400p-48, + 0x1.93737b0cdc64ap-1, 0x1.7200p-46, + 0x1.948b82b5f98aep-1, -0x1.9000p-47, + 0x1.95a44cbc852cbp-1, 0x1.5680p-45, + 0x1.96bdd9a766f21p-1, -0x1.6d00p-44, + 0x1.97d829fde4e2ap-1, -0x1.1000p-47, + 0x1.98f33e47a23a3p-1, 0x1.d000p-45, + 0x1.9a0f170ca0604p-1, -0x1.8a40p-44, + 0x1.9b2bb4d53ff89p-1, 0x1.55c0p-44, + 0x1.9c49182a3f15bp-1, 0x1.6b80p-45, + 0x1.9d674194bb8c5p-1, -0x1.c000p-49, + 0x1.9e86319e3238ep-1, 0x1.7d00p-46, + 0x1.9fa5e8d07f302p-1, 0x1.6400p-46, + 0x1.a0c667b5de54dp-1, -0x1.5000p-48, + 0x1.a1e7aed8eb8f6p-1, 0x1.9e00p-47, + 0x1.a309bec4a2e27p-1, 0x1.ad80p-45, + 0x1.a42c980460a5dp-1, -0x1.af00p-46, + 0x1.a5503b23e259bp-1, 0x1.b600p-47, + 0x1.a674a8af46213p-1, 0x1.8880p-44, + 0x1.a799e1330b3a7p-1, 0x1.1200p-46, + 0x1.a8bfe53c12e8dp-1, 0x1.6c00p-47, + 0x1.a9e6b5579fcd2p-1, -0x1.9b80p-45, + 0x1.ab0e521356fb8p-1, 0x1.b700p-45, + 0x1.ac36bbfd3f381p-1, 0x1.9000p-50, + 0x1.ad5ff3a3c2780p-1, 0x1.4000p-49, + 0x1.ae89f995ad2a3p-1, -0x1.c900p-45, + 0x1.afb4ce622f367p-1, 0x1.6500p-46, + 0x1.b0e07298db790p-1, 0x1.fd40p-45, + 0x1.b20ce6c9a89a9p-1, 0x1.2700p-46, + 0x1.b33a2b84f1a4bp-1, 0x1.d470p-43, + 0x1.b468415b747e7p-1, -0x1.8380p-44, + 0x1.b59728de5593ap-1, 0x1.8000p-54, + 0x1.b6c6e29f1c56ap-1, 0x1.ad00p-47, + 0x1.b7f76f2fb5e50p-1, 0x1.e800p-50, + 0x1.b928cf22749b2p-1, -0x1.4c00p-47, + 0x1.ba5b030a10603p-1, -0x1.d700p-47, + 0x1.bb8e0b79a6f66p-1, 0x1.d900p-47, + 0x1.bcc1e904bc1ffp-1, 0x1.2a00p-47, + 0x1.bdf69c3f3a16fp-1, -0x1.f780p-46, + 0x1.bf2c25bd71db8p-1, -0x1.0a00p-46, + 0x1.c06286141b2e9p-1, -0x1.1400p-46, + 0x1.c199bdd8552e0p-1, 0x1.be00p-47, + 0x1.c2d1cd9fa64eep-1, -0x1.9400p-47, + 0x1.c40ab5fffd02fp-1, -0x1.ed00p-47, + 0x1.c544778fafd15p-1, 0x1.9660p-44, + 0x1.c67f12e57d0cbp-1, -0x1.a100p-46, + 0x1.c7ba88988c1b6p-1, -0x1.8458p-42, + 0x1.c8f6d9406e733p-1, -0x1.a480p-46, + 0x1.ca3405751c4dfp-1, 0x1.b000p-51, + 0x1.cb720dcef9094p-1, 0x1.1400p-47, + 0x1.ccb0f2e6d1689p-1, 0x1.0200p-48, + 0x1.cdf0b555dc412p-1, 0x1.3600p-48, + 0x1.cf3155b5bab3bp-1, -0x1.6900p-47, + 0x1.d072d4a0789bcp-1, 0x1.9a00p-47, + 0x1.d1b532b08c8fap-1, -0x1.5e00p-46, + 0x1.d2f87080d8a85p-1, 0x1.d280p-46, + 0x1.d43c8eacaa203p-1, 0x1.1a00p-47, + 0x1.d5818dcfba491p-1, 0x1.f000p-50, + 0x1.d6c76e862e6a1p-1, -0x1.3a00p-47, + 0x1.d80e316c9834ep-1, -0x1.cd80p-47, + 0x1.d955d71ff6090p-1, 0x1.4c00p-48, + 0x1.da9e603db32aep-1, 0x1.f900p-48, + 0x1.dbe7cd63a8325p-1, 0x1.9800p-49, + 0x1.dd321f301b445p-1, -0x1.5200p-48, + 0x1.de7d5641c05bfp-1, -0x1.d700p-46, + 0x1.dfc97337b9aecp-1, -0x1.6140p-46, + 0x1.e11676b197d5ep-1, 0x1.b480p-47, + 0x1.e264614f5a3e7p-1, 0x1.0ce0p-43, + 0x1.e3b333b16ee5cp-1, 0x1.c680p-47, + 0x1.e502ee78b3fb4p-1, -0x1.9300p-47, + 0x1.e653924676d68p-1, -0x1.5000p-49, + 0x1.e7a51fbc74c44p-1, -0x1.7f80p-47, + 0x1.e8f7977cdb726p-1, -0x1.3700p-48, + 0x1.ea4afa2a490e8p-1, 0x1.5d00p-49, + 0x1.eb9f4867ccae4p-1, 0x1.61a0p-46, + 0x1.ecf482d8e680dp-1, 0x1.5500p-48, + 0x1.ee4aaa2188514p-1, 0x1.6400p-51, + 0x1.efa1bee615a13p-1, -0x1.e800p-49, + 0x1.f0f9c1cb64106p-1, -0x1.a880p-48, + 0x1.f252b376bb963p-1, -0x1.c900p-45, + 0x1.f3ac948dd7275p-1, 0x1.a000p-53, + 0x1.f50765b6e4524p-1, -0x1.4f00p-48, + 0x1.f6632798844fdp-1, 0x1.a800p-51, + 0x1.f7bfdad9cbe38p-1, 0x1.abc0p-48, + 0x1.f91d802243c82p-1, -0x1.4600p-50, + 0x1.fa7c1819e908ep-1, -0x1.b0c0p-47, + 0x1.fbdba3692d511p-1, -0x1.0e00p-51, + 0x1.fd3c22b8f7194p-1, -0x1.0de8p-46, + 0x1.fe9d96b2a23eep-1, 0x1.e430p-49, + 0x1.0000000000000p+0, 0x0.0000p+0, + 0x1.00b1afa5abcbep+0, -0x1.3400p-52, + 0x1.0163da9fb3303p+0, -0x1.2170p-46, + 0x1.02168143b0282p+0, 0x1.a400p-52, + 0x1.02c9a3e77806cp+0, 0x1.f980p-49, + 0x1.037d42e11bbcap+0, -0x1.7400p-51, + 0x1.04315e86e7f89p+0, 0x1.8300p-50, + 0x1.04e5f72f65467p+0, -0x1.a3f0p-46, + 0x1.059b0d315855ap+0, -0x1.2840p-47, + 0x1.0650a0e3c1f95p+0, 0x1.1600p-48, + 0x1.0706b29ddf71ap+0, 0x1.5240p-46, + 0x1.07bd42b72a82dp+0, -0x1.9a00p-49, + 0x1.0874518759bd0p+0, 0x1.6400p-49, + 0x1.092bdf66607c8p+0, -0x1.0780p-47, + 0x1.09e3ecac6f383p+0, -0x1.8000p-54, + 0x1.0a9c79b1f3930p+0, 0x1.fa00p-48, + 0x1.0b5586cf988fcp+0, -0x1.ac80p-48, + 0x1.0c0f145e46c8ap+0, 0x1.9c00p-50, + 0x1.0cc922b724816p+0, 0x1.5200p-47, + 0x1.0d83b23395dd8p+0, -0x1.ad00p-48, + 0x1.0e3ec32d3d1f3p+0, 0x1.bac0p-46, + 0x1.0efa55fdfa9a6p+0, -0x1.4e80p-47, + 0x1.0fb66affed2f0p+0, -0x1.d300p-47, + 0x1.1073028d7234bp+0, 0x1.1500p-48, + 0x1.11301d0125b5bp+0, 0x1.c000p-49, + 0x1.11edbab5e2af9p+0, 0x1.6bc0p-46, + 0x1.12abdc06c31d5p+0, 0x1.8400p-49, + 0x1.136a814f2047dp+0, -0x1.ed00p-47, + 0x1.1429aaea92de9p+0, 0x1.8e00p-49, + 0x1.14e95934f3138p+0, 0x1.b400p-49, + 0x1.15a98c8a58e71p+0, 0x1.5300p-47, + 0x1.166a45471c3dfp+0, 0x1.3380p-47, + 0x1.172b83c7d5211p+0, 0x1.8d40p-45, + 0x1.17ed48695bb9fp+0, -0x1.5d00p-47, + 0x1.18af9388c8d93p+0, -0x1.c880p-46, + 0x1.1972658375d66p+0, 0x1.1f00p-46, + 0x1.1a35beb6fcba7p+0, 0x1.0480p-46, + 0x1.1af99f81387e3p+0, -0x1.7390p-43, + 0x1.1bbe084045d54p+0, 0x1.4e40p-45, + 0x1.1c82f95281c43p+0, -0x1.a200p-47, + 0x1.1d4873168b9b2p+0, 0x1.3800p-49, + 0x1.1e0e75eb44031p+0, 0x1.ac00p-49, + 0x1.1ed5022fcd938p+0, 0x1.1900p-47, + 0x1.1f9c18438cdf7p+0, -0x1.b780p-46, + 0x1.2063b88628d8fp+0, 0x1.d940p-45, + 0x1.212be3578a81ep+0, 0x1.8000p-50, + 0x1.21f49917ddd41p+0, 0x1.b340p-45, + 0x1.22bdda2791323p+0, 0x1.9f80p-46, + 0x1.2387a6e7561e7p+0, -0x1.9c80p-46, + 0x1.2451ffb821427p+0, 0x1.2300p-47, + 0x1.251ce4fb2a602p+0, -0x1.3480p-46, + 0x1.25e85711eceb0p+0, 0x1.2700p-46, + 0x1.26b4565e27d16p+0, 0x1.1d00p-46, + 0x1.2780e341de00fp+0, 0x1.1ee0p-44, + 0x1.284dfe1f5633ep+0, -0x1.4c00p-46, + 0x1.291ba7591bb30p+0, -0x1.3d80p-46, + 0x1.29e9df51fdf09p+0, 0x1.8b00p-47, + 0x1.2ab8a66d10e9bp+0, -0x1.27c0p-45, + 0x1.2b87fd0dada3ap+0, 0x1.a340p-45, + 0x1.2c57e39771af9p+0, -0x1.0800p-46, + 0x1.2d285a6e402d9p+0, -0x1.ed00p-47, + 0x1.2df961f641579p+0, -0x1.4200p-48, + 0x1.2ecafa93e2ecfp+0, -0x1.4980p-45, + 0x1.2f9d24abd8822p+0, -0x1.6300p-46, + 0x1.306fe0a31b625p+0, -0x1.2360p-44, + 0x1.31432edeea50bp+0, -0x1.0df8p-40, + 0x1.32170fc4cd7b8p+0, -0x1.2480p-45, + 0x1.32eb83ba8e9a2p+0, -0x1.5980p-45, + 0x1.33c08b2641766p+0, 0x1.ed00p-46, + 0x1.3496266e3fa27p+0, -0x1.c000p-50, + 0x1.356c55f929f0fp+0, -0x1.0d80p-44, + 0x1.36431a2de88b9p+0, 0x1.2c80p-45, + 0x1.371a7373aaa39p+0, 0x1.0600p-45, + 0x1.37f26231e74fep+0, -0x1.6600p-46, + 0x1.38cae6d05d838p+0, -0x1.ae00p-47, + 0x1.39a401b713ec3p+0, -0x1.4720p-43, + 0x1.3a7db34e5a020p+0, 0x1.8200p-47, + 0x1.3b57fbfec6e95p+0, 0x1.e800p-44, + 0x1.3c32dc313a8f2p+0, 0x1.f800p-49, + 0x1.3d0e544ede122p+0, -0x1.7a00p-46, + 0x1.3dea64c1234bbp+0, 0x1.6300p-45, + 0x1.3ec70df1c4eccp+0, -0x1.8a60p-43, + 0x1.3fa4504ac7e8cp+0, -0x1.cdc0p-44, + 0x1.40822c367a0bbp+0, 0x1.5b80p-45, + 0x1.4160a21f72e95p+0, 0x1.ec00p-46, + 0x1.423fb27094646p+0, -0x1.3600p-46, + 0x1.431f5d950a920p+0, 0x1.3980p-45, + 0x1.43ffa3f84b9ebp+0, 0x1.a000p-48, + 0x1.44e0860618919p+0, -0x1.6c00p-48, + 0x1.45c2042a7d201p+0, -0x1.bc00p-47, + 0x1.46a41ed1d0016p+0, -0x1.2800p-46, + 0x1.4786d668b3326p+0, 0x1.0e00p-44, + 0x1.486a2b5c13c00p+0, -0x1.d400p-45, + 0x1.494e1e192af04p+0, 0x1.c200p-47, + 0x1.4a32af0d7d372p+0, -0x1.e500p-46, + 0x1.4b17dea6db801p+0, 0x1.7800p-47, + 0x1.4bfdad53629e1p+0, -0x1.3800p-46, + 0x1.4ce41b817c132p+0, 0x1.0800p-47, + 0x1.4dcb299fddddbp+0, 0x1.c700p-45, + 0x1.4eb2d81d8ab96p+0, -0x1.ce00p-46, + 0x1.4f9b2769d2d02p+0, 0x1.9200p-46, + 0x1.508417f4531c1p+0, -0x1.8c00p-47, + 0x1.516daa2cf662ap+0, -0x1.a000p-48, + 0x1.5257de83f51eap+0, 0x1.a080p-43, + 0x1.5342b569d4edap+0, -0x1.6d80p-45, + 0x1.542e2f4f6ac1ap+0, -0x1.2440p-44, + 0x1.551a4ca5d94dbp+0, 0x1.83c0p-43, + 0x1.56070dde9116bp+0, 0x1.4b00p-45, + 0x1.56f4736b529dep+0, 0x1.15a0p-43, + 0x1.57e27dbe2c40ep+0, -0x1.9e00p-45, + 0x1.58d12d497c76fp+0, -0x1.3080p-45, + 0x1.59c0827ff0b4cp+0, 0x1.dec0p-43, + 0x1.5ab07dd485427p+0, -0x1.4000p-51, + 0x1.5ba11fba87af4p+0, 0x1.0080p-44, + 0x1.5c9268a59460bp+0, -0x1.6c80p-45, + 0x1.5d84590998e3fp+0, 0x1.69a0p-43, + 0x1.5e76f15ad20e1p+0, -0x1.b400p-46, + 0x1.5f6a320dcebcap+0, 0x1.7700p-46, + 0x1.605e1b976dcb8p+0, 0x1.6f80p-45, + 0x1.6152ae6cdf715p+0, 0x1.1000p-47, + 0x1.6247eb03a5531p+0, -0x1.5d00p-46, + 0x1.633dd1d1929b5p+0, -0x1.2d00p-46, + 0x1.6434634ccc313p+0, -0x1.a800p-49, + 0x1.652b9febc8efap+0, -0x1.8600p-45, + 0x1.6623882553397p+0, 0x1.1fe0p-40, + 0x1.671c1c708328ep+0, -0x1.7200p-44, + 0x1.68155d44ca97ep+0, 0x1.6800p-49, + 0x1.690f4b19e9471p+0, -0x1.9780p-45, +}; + +fn exp2d(x: f64) -> f64 { + const tblsiz = u32(exp2dt.len / 2); + const redux: f64 = 0x1.8p52 / 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 = u32(ux >> 32) & 0x7FFFFFFF; + + // |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.forceEval(f32(-0x1.0p-149 / x)); + } + if (x <= -1075) { + return 0; + } + } + } + // |x| < 0x1p-54 + else if (ix < 0x3C900000) { + return 1.0 + x; + } + + // reduce x + var uf = x + redux; + // NOTE: musl performs an implicit 64-bit to 32-bit u32 truncation here + var i0 = @truncate(u32, @bitCast(u64, uf)); + i0 += tblsiz / 2; + + const k: u32 = i0 / tblsiz * tblsiz; + const ik = @bitCast(i32, k / tblsiz); + i0 %= tblsiz; + uf -= redux; + + // r = exp2(y) = exp2t[i0] * p(z - eps[i]) + var z = x - uf; + const t = exp2dt[2 * i0]; + z -= exp2dt[2 * i0 + 1]; + const r = t + t * z * (P1 + z * (P2 + z * (P3 + z * (P4 + z * P5)))); + + math.scalbn(r, ik) +} + +test "exp2" { + assert(exp2(f32(0.8923)) == exp2f(0.8923)); + assert(exp2(f64(0.8923)) == exp2d(0.8923)); +} + +test "exp2f" { + const epsilon = 0.000001; + + assert(exp2f(0.0) == 1.0); + assert(math.approxEq(f32, exp2f(0.2), 1.148698, epsilon)); + assert(math.approxEq(f32, exp2f(0.8923), 1.856133, epsilon)); + assert(math.approxEq(f32, exp2f(1.5), 2.828427, epsilon)); + assert(math.approxEq(f32, exp2f(37.45), 187747237888, epsilon)); +} + +test "exp2d" { + const epsilon = 0.000001; + + assert(exp2d(0.0) == 1.0); + assert(math.approxEq(f64, exp2d(0.2), 1.148698, epsilon)); + assert(math.approxEq(f64, exp2d(0.8923), 1.856133, epsilon)); + assert(math.approxEq(f64, exp2d(1.5), 2.828427, epsilon)); + // assert(math.approxEq(f64, exp2d(37.45), 18379273786760560.000000, epsilon)); +} diff --git a/std/math/expm1.zig b/std/math/expm1.zig new file mode 100644 index 0000000000..7b089e023b --- /dev/null +++ b/std/math/expm1.zig @@ -0,0 +1,282 @@ +const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; + +pub fn expm1(x: var) -> @typeOf(x) { + const T = @typeOf(x); + switch (T) { + f32 => @inlineCall(expm1f, x), + f64 => @inlineCall(expm1d, x), + else => @compileError("exp1m not implemented for " ++ @typeName(T)), + } +} + +fn expm1f(x_: f32) -> f32 { + const o_threshold: f32 = 8.8721679688e+01; + const ln2_hi: f32 = 6.9313812256e-01; + const ln2_lo: f32 = 9.0580006145e-06; + const invln2: f32 = 1.4426950216e+00; + const Q1: f32 = -3.3333212137e-2; + const Q2: f32 = 1.5807170421e-3; + + var x = x_; + const ux = @bitCast(u32, x); + const hx = ux & 0x7FFFFFFF; + const sign = hx >> 31; + + // |x| >= 27 * ln2 + if (hx >= 0x4195B844) { + // nan + if (hx > 0x7F800000) { + return x; + } + if (sign != 0) { + return -1; + } + if (x > o_threshold) { + x *= 0x1.0p127; + return x; + } + } + + var hi: f32 = undefined; + var lo: f32 = undefined; + var c: f32 = undefined; + var k: i32 = undefined; + + // |x| > 0.5 * ln2 + if (hx > 0x3EB17218) { + // |x| < 1.5 * ln2 + if (hx < 0x3F851592) { + if (sign == 0) { + hi = x - ln2_hi; + lo = ln2_lo; + k = 1; + } else { + hi = x + ln2_hi; + lo = -ln2_lo; + k = -1; + } + } else { + var kf = invln2 * x; + if (sign != 0) { + kf -= 0.5; + } else { + kf += 0.5; + } + + k = i32(kf); + const t = f32(k); + hi = x - t * ln2_hi; + lo = t * ln2_lo; + } + + x = hi - lo; + c = (hi - x) - lo; + } + // |x| < 2^(-25) + else if (hx < 0x33000000) { + if (hx < 0x00800000) { + math.forceEval(x * x); + } + return x; + } + else { + k = 0; + } + + const hfx = 0.5 * x; + const hxs = x * hfx; + const r1 = 1.0 + hxs * (Q1 + hxs * Q2); + const t = 3.0 - r1 * hfx; + var e = hxs * ((r1 - t) / (6.0 - x * t)); + + // c is 0 + if (k == 0) { + return x - (x * e - hxs); + } + + e = x * (e - c) - c; + e -= hxs; + + // exp(x) ~ 2^k (x_reduced - e + 1) + if (k == -1) { + return 0.5 * (x - e) - 0.5; + } + if (k == 1) { + if (x < -0.25) { + return -2.0 * (e - (x + 0.5)); + } else { + return 1.0 + 2.0 * (x - e); + } + } + + const twopk = @bitCast(f32, u32(0x7F + k) << 23); + + if (k < 0 or k > 56) { + var y = x - e + 1.0; + if (k == 128) { + y = y * 2.0 * 0x1.0p127; + } else { + y = y * twopk; + } + + return y - 1.0; + } + + const uf = @bitCast(f32, u32(0x7F - k) << 23); + if (k < 23) { + return (x - e + (1 - uf)) * twopk; + } else { + return (x - (e + uf) + 1) * twopk; + } +} + +fn expm1d(x_: f64) -> f64 { + const o_threshold: f64 = 7.09782712893383973096e+02; + const ln2_hi: f64 = 6.93147180369123816490e-01; + const ln2_lo: f64 = 1.90821492927058770002e-10; + const invln2: f64 = 1.44269504088896338700e+00; + const Q1: f64 = -3.33333333333331316428e-02; + const Q2: f64 = 1.58730158725481460165e-03; + const Q3: f64 = -7.93650757867487942473e-05; + const Q4: f64 = 4.00821782732936239552e-06; + const Q5: f64 = -2.01099218183624371326e-07; + + var x = x_; + const ux = @bitCast(u64, x); + const hx = u32(ux >> 32) & 0x7FFFFFFF; + const sign = hx >> 63; + + // |x| >= 56 * ln2 + if (hx >= 0x4043687A) { + // exp1md(nan) = nan + if (hx > 0x7FF00000) { + return x; + } + // exp1md(-ve) = -1 + if (sign != 0) { + return -1; + } + if (x > o_threshold) { + math.raiseOverflow(); + return math.nan(f64); + } + } + + var hi: f64 = undefined; + var lo: f64 = undefined; + var c: f64 = undefined; + var k: i32 = undefined; + + // |x| > 0.5 * ln2 + if (hx > 0x3FD62E42) { + // |x| < 1.5 * ln2 + if (hx < 0x3FF0A2B2) { + if (sign == 0) { + hi = x - ln2_hi; + lo = ln2_lo; + k = 1; + } else { + hi = x + ln2_hi; + lo = -ln2_lo; + k = -1; + } + } else { + var kf = invln2 * x; + if (sign != 0) { + kf -= 0.5; + } else { + kf += 0.5; + } + + k = i32(kf); + const t = f64(k); + hi = x - t * ln2_hi; + lo = t * ln2_lo; + } + + x = hi - lo; + c = (hi - x) - lo; + } + // |x| < 2^(-54) + else if (hx < 0x3C900000) { + if (hx < 0x00100000) { + math.forceEval(f32(x)); + } + return x; + } + else { + k = 0; + } + + const hfx = 0.5 * x; + const hxs = x * hfx; + const r1 = 1.0 + hxs * (Q1 + hxs * (Q2 + hxs * (Q3 + hxs * (Q4 + hxs * Q5)))); + const t = 3.0 - r1 * hfx; + var e = hxs * ((r1 - t) / (6.0 - x * t)); + + // c is 0 + if (k == 0) { + return x - (x * e - hxs); + } + + e = x * (e - c) - c; + e -= hxs; + + // exp(x) ~ 2^k (x_reduced - e + 1) + if (k == -1) { + return 0.5 * (x - e) - 0.5; + } + if (k == 1) { + if (x < -0.25) { + return -2.0 * (e - (x + 0.5)); + } else { + return 1.0 + 2.0 * (x - e); + } + } + + const twopk = @bitCast(f64, u64(0x3FF + k) << 52); + + if (k < 0 or k > 56) { + var y = x - e + 1.0; + if (k == 1024) { + y = y * 2.0; // TODO: * 0x1.0p1023; + } else { + y = y * twopk; + } + + return y - 1.0; + } + + const uf = @bitCast(f64, u64(0x3FF - k) << 52); + if (k < 20) { + return (x - e + (1 - uf)) * twopk; + } else { + return (x - (e + uf) + 1) * twopk; + } +} + +test "exp1m" { + assert(expm1(f32(0.0)) == expm1f(0.0)); + assert(expm1(f64(0.0)) == expm1d(0.0)); +} + +test "expm1f" { + const epsilon = 0.000001; + + assert(expm1f(0.0) == 0.0); + assert(math.approxEq(f32, expm1f(0.0), 0.0, epsilon)); + assert(math.approxEq(f32, expm1f(0.2), 0.221403, epsilon)); + assert(math.approxEq(f32, expm1f(0.8923), 1.440737, epsilon)); + assert(math.approxEq(f32, expm1f(1.5), 3.481689, epsilon)); +} + +test "expm1d" { + const epsilon = 0.000001; + + assert(expm1d(0.0) == 0.0); + assert(math.approxEq(f64, expm1d(0.0), 0.0, epsilon)); + assert(math.approxEq(f64, expm1d(0.2), 0.221403, epsilon)); + assert(math.approxEq(f64, expm1d(0.8923), 1.440737, epsilon)); + assert(math.approxEq(f64, expm1d(1.5), 3.481689, epsilon)); +} diff --git a/std/math/fabs.zig b/std/math/fabs.zig index 268a339b18..128d4eae76 100644 --- a/std/math/fabs.zig +++ b/std/math/fabs.zig @@ -1,10 +1,11 @@ +const math = @import("index.zig"); const assert = @import("../debug.zig").assert; pub fn fabs(x: var) -> @typeOf(x) { const T = @typeOf(x); switch (T) { - f32 => fabs32(x), - f64 => fabs64(x), + f32 => @inlineCall(fabs32, x), + f64 => @inlineCall(fabs64, x), else => @compileError("fabs not implemented for " ++ @typeName(T)), } } @@ -24,26 +25,14 @@ fn fabs64(x: f64) -> f64 { test "fabs" { assert(fabs(f32(1.0)) == fabs32(1.0)); assert(fabs(f64(1.0)) == fabs64(1.0)); - comptime { - assert(fabs(f32(1.0)) == fabs32(1.0)); - assert(fabs(f64(1.0)) == fabs64(1.0)); - } } test "fabs32" { assert(fabs64(1.0) == 1.0); assert(fabs64(-1.0) == 1.0); - comptime { - assert(fabs64(1.0) == 1.0); - assert(fabs64(-1.0) == 1.0); - } } test "fabs64" { assert(fabs64(1.0) == 1.0); assert(fabs64(-1.0) == 1.0); - comptime { - assert(fabs64(1.0) == 1.0); - assert(fabs64(-1.0) == 1.0); - } } diff --git a/std/math/floor.zig b/std/math/floor.zig new file mode 100644 index 0000000000..af1ee98bed --- /dev/null +++ b/std/math/floor.zig @@ -0,0 +1,89 @@ +const builtin = @import("builtin"); +const assert = @import("../debug.zig").assert; +const math = @import("index.zig"); + +pub fn floor(x: var) -> @typeOf(x) { + const T = @typeOf(x); + switch (T) { + f32 => @inlineCall(floor32, x), + f64 => @inlineCall(floor64, x), + else => @compileError("floor not implemented for " ++ @typeName(T)), + } +} + +fn floor32(x: f32) -> f32 { + var u = @bitCast(u32, x); + const e = i32((u >> 23) & 0xFF) - 0x7F; + var m: u32 = undefined; + + if (e >= 23) { + return x; + } + + if (e >= 0) { + m = 0x007FFFFF >> u32(e); + if (u & m == 0) { + return x; + } + math.forceEval(x + 0x1.0p120); + if (u >> 31 != 0) { + u += m; + } + @bitCast(f32, u & ~m) + } else { + math.forceEval(x + 0x1.0p120); + if (u >> 31 == 0) { + return 0.0; // Compiler requires return + } else { + -1.0 + } + } +} + +fn floor64(x: f64) -> f64 { + const u = @bitCast(u64, x); + const e = (u >> 52) & 0x7FF; + var y: f64 = undefined; + + if (e >= 0x3FF+52 or x == 0) { + return x; + } + + if (u >> 63 != 0) { + @setFloatMode(this, builtin.FloatMode.Strict); + y = x - math.f64_toint + math.f64_toint - x; + } else { + @setFloatMode(this, builtin.FloatMode.Strict); + y = x + math.f64_toint - math.f64_toint - x; + } + + if (e <= 0x3FF-1) { + math.forceEval(y); + if (u >> 63 != 0) { + return -1.0; // Compiler requires return. + } else { + 0.0 + } + } else if (y > 0) { + x + y - 1 + } else { + x + y + } +} + +test "floor" { + assert(floor(f32(1.3)) == floor32(1.3)); + assert(floor(f64(1.3)) == floor64(1.3)); +} + +test "floor32" { + assert(floor32(1.3) == 1.0); + assert(floor32(-1.3) == -2.0); + assert(floor32(0.2) == 0.0); +} + +test "floor64" { + assert(floor64(1.3) == 1.0); + assert(floor64(-1.3) == -2.0); + assert(floor64(0.2) == 0.0); +} diff --git a/std/math/fma.zig b/std/math/fma.zig new file mode 100644 index 0000000000..2ff57b968d --- /dev/null +++ b/std/math/fma.zig @@ -0,0 +1,160 @@ +const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; + +pub fn fma(comptime T: type, x: T, y: T, z: T) -> T { + switch (T) { + f32 => @inlineCall(fma32, x, y, z), + f64 => @inlineCall(fma64, x, y ,z), + else => @compileError("fma not implemented for " ++ @typeName(T)), + } +} + +fn fma32(x: f32, y: f32, z: f32) -> f32 { + const xy = f64(x) * y; + const xy_z = xy + z; + const u = @bitCast(u64, xy_z); + const e = (u >> 52) & 0x7FF; + + if ((u & 0x1FFFFFFF) != 0x10000000 or e == 0x7FF or xy_z - xy == z) { + f32(xy_z) + } else { + // TODO: Handle inexact case with double-rounding + f32(xy_z) + } +} + +fn fma64(x: f64, y: f64, z: f64) -> f64 { + 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 <= 53 * 2) { + zs = math.scalbn(zs, -spread); + } else { + zs = math.copysign(f64, math.f64_min, zs); + } + + const xy = dd_mul(xs, ys); + const r = dd_add(xy.hi, zs); + spread = ex + ey; + + if (r.hi == 0.0) { + return xy.hi + zs + math.scalbn(xy.lo, spread); + } + + const adj = add_adjusted(r.lo, xy.lo); + if (spread + math.ilogb(r.hi) > -1023) { + math.scalbn(r.hi + adj, spread) + } else { + add_and_denorm(r.hi, adj, spread) + } +} + +const dd = struct { hi: f64, lo: f64, }; + +fn dd_add(a: f64, b: f64) -> dd { + var ret: dd = undefined; + ret.hi = a + b; + const s = ret.hi - a; + ret.lo = (a - (ret.hi - s)) + (b - s); + ret +} + +fn dd_mul(a: f64, b: f64) -> dd { + var ret: dd = undefined; + const split: f64 = 0x1.0p27 + 1.0; + + var p = a * split; + var ha = a - p; + ha += p; + var la = a - ha; + + p = b * split; + var hb = b - p; + hb += p; + var lb = b - hb; + + p = ha * hb; + var q = ha * lb + la * hb; + + ret.hi = p + q; + ret.lo = p - ret.hi + q + la * lb; + ret +} + +fn add_adjusted(a: f64, b: f64) -> f64 { + var sum = dd_add(a, b); + if (sum.lo != 0) { + var uhii = @bitCast(u64, sum.hi); + if (uhii & 1 == 0) { + // hibits += copysign(1.0, sum.hi, sum.lo) + const uloi = @bitCast(u64, sum.lo); + uhii += 1 - ((uhii ^ uloi) >> 62); + sum.hi = @bitCast(f64, uhii); + } + } + sum.hi +} + +fn add_and_denorm(a: f64, b: f64, scale: i32) -> f64 { + var sum = dd_add(a, b); + if (sum.lo != 0) { + var uhii = @bitCast(u64, sum.hi); + const bits_lost = -i32((uhii >> 52) & 0x7FF) - scale + 1; + if ((bits_lost != 1) == (uhii & 1 != 0)) { + const uloi = @bitCast(u64, sum.lo); + uhii += 1 - (((uhii ^ uloi) >> 62) & 2); + sum.hi = @bitCast(f64, uhii); + } + } + math.scalbn(sum.hi, scale) +} + +test "fma" { + assert(fma(f32, 0.0, 1.0, 1.0) == fma32(0.0, 1.0, 1.0)); + assert(fma(f64, 0.0, 1.0, 1.0) == fma64(0.0, 1.0, 1.0)); +} + +test "fma32" { + const epsilon = 0.000001; + + assert(math.approxEq(f32, fma32(0.0, 5.0, 9.124), 9.124, epsilon)); + assert(math.approxEq(f32, fma32(0.2, 5.0, 9.124), 10.124, epsilon)); + assert(math.approxEq(f32, fma32(0.8923, 5.0, 9.124), 13.5855, epsilon)); + assert(math.approxEq(f32, fma32(1.5, 5.0, 9.124), 16.624, epsilon)); + assert(math.approxEq(f32, fma32(37.45, 5.0, 9.124), 196.374004, epsilon)); + assert(math.approxEq(f32, fma32(89.123, 5.0, 9.124), 454.739005, epsilon)); + assert(math.approxEq(f32, fma32(123123.234375, 5.0, 9.124), 615625.295875, epsilon)); +} + +test "fma64" { + const epsilon = 0.000001; + + assert(math.approxEq(f64, fma64(0.0, 5.0, 9.124), 9.124, epsilon)); + assert(math.approxEq(f64, fma64(0.2, 5.0, 9.124), 10.124, epsilon)); + assert(math.approxEq(f64, fma64(0.8923, 5.0, 9.124), 13.5855, epsilon)); + assert(math.approxEq(f64, fma64(1.5, 5.0, 9.124), 16.624, epsilon)); + assert(math.approxEq(f64, fma64(37.45, 5.0, 9.124), 196.374, epsilon)); + assert(math.approxEq(f64, fma64(89.123, 5.0, 9.124), 454.739, epsilon)); + assert(math.approxEq(f64, fma64(123123.234375, 5.0, 9.124), 615625.295875, epsilon)); +} diff --git a/std/math/fmod.zig b/std/math/fmod.zig new file mode 100644 index 0000000000..8cc9586bc0 --- /dev/null +++ b/std/math/fmod.zig @@ -0,0 +1,190 @@ +const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; + +pub fn fmod(comptime T: type, x: T, y: T) -> T { + switch (T) { + f32 => @inlineCall(fmod32, x, y), + f64 => @inlineCall(fmod64, x, y), + else => @compileError("fmod not implemented for " ++ @typeName(T)), + } +} + +fn fmod32(x: f32, y: f32) -> f32 { + var ux = @bitCast(u32, x); + var uy = @bitCast(u32, y); + var ex = i32(ux >> 23) & 0xFF; + var ey = i32(ux >> 23) & 0xFF; + const sx = ux & 0x80000000; + + if (uy << 1 == 0 or math.isNan(y) or ex == 0xFF) { + return (x * y) / (x * y); + } + if (ux << 1 <= uy << 1) { + if (ux << 1 == uy << 1) { + return 0 * x; + } else { + return x; + } + } + + // normalize x and y + if (ex == 0) { + var i = ux << 9; + while (i >> 31 == 0) : (i <<= 1) { + ex -= 1; + } + ux <<= u32(-ex + 1); + } else { + ux &= @maxValue(u32) >> 9; + ux |= 1 << 23; + } + + if (ey == 0) { + var i = uy << 9; + while (i >> 31 == 0) : (i <<= 1) { + ey -= 1; + } + uy <<= u32(-ey + 1); + } else { + uy &= @maxValue(u32) >> 9; + uy |= 1 << 23; + } + + // x mod y + while (ex > ey) : (ex -= 1) { + const i = ux - uy; + if (i >> 31 == 0) { + if (i == 0) { + return 0 * x; + } + ux = i; + } + ux <<= 1; + } + { + const i = ux - uy; + if (i >> 31 == 0) { + if (i == 0) { + return 0 * x; + } + ux = i; + } + } + + while (ux >> 23 == 0) : (ux <<= 1) { + ex -= 1; + } + + // scale result up + if (ex > 0) { + ux -= 1 << 23; + ux |= u32(ex) << 23; + } else { + ux >>= u32(-ex + 1); + } + + ux |= sx; + @bitCast(f32, ux) +} + +fn fmod64(x: f64, y: f64) -> f64 { + var ux = @bitCast(u64, x); + var uy = @bitCast(u64, y); + var ex = i32(ux >> 52) & 0x7FF; + var ey = i32(ux >> 52) & 0x7FF; + const sx = ux >> 63; + + if (uy << 1 == 0 or math.isNan(y) or ex == 0x7FF) { + return (x * y) / (x * y); + } + if (ux << 1 <= uy << 1) { + if (ux << 1 == uy << 1) { + return 0 * x; + } else { + return x; + } + } + + // normalize x and y + if (ex == 0) { + var i = ux << 12; + while (i >> 63 == 0) : (i <<= 1) { + ex -= 1; + } + ux <<= u64(-ex + 1); + } else { + ux &= @maxValue(u64) >> 12; + ux |= 1 << 52; + } + + if (ey == 0) { + var i = uy << 12; + while (i >> 63 == 0) : (i <<= 1) { + ey -= 1; + } + uy <<= u64(-ey + 1); + } else { + uy &= @maxValue(u64) >> 12; + uy |= 1 << 52; + } + + // x mod y + while (ex > ey) : (ex -= 1) { + const i = ux - uy; + if (i >> 63 == 0) { + if (i == 0) { + return 0 * x; + } + ux = i; + } + ux <<= 1; + } + { + const i = ux - uy; + if (i >> 63 == 0) { + if (i == 0) { + return 0 * x; + } + ux = i; + } + } + + while (ux >> 52 == 0) : (ux <<= 1) { + ex -= 1; + } + + // scale result up + if (ex > 0) { + ux -= 1 << 52; + ux |= u64(ex) << 52; + } else { + ux >>= u64(-ex + 1); + } + + ux |= sx << 63; + @bitCast(f64, ux) +} + +// duplicate symbol clash with `fmod` test name +test "fmod_" { + assert(fmod(f32, 1.3, 2.5) == fmod32(1.3, 2.5)); + assert(fmod(f64, 1.3, 2.5) == fmod64(1.3, 2.5)); +} + +test "fmod32" { + const epsilon = 0.000001; + + assert(math.approxEq(f32, fmod32(5.2, 2.0), 1.2, epsilon)); + assert(math.approxEq(f32, fmod32(18.5, 4.2), 1.7, epsilon)); + assert(math.approxEq(f32, fmod32(23, 48.34), 23.0, epsilon)); + assert(math.approxEq(f32, fmod32(123.340890, 2398.2314), 123.340889, epsilon)); +} + +test "fmod64" { + const epsilon = 0.000001; + + assert(math.approxEq(f64, fmod64(5.2, 2.0), 1.2, epsilon)); + assert(math.approxEq(f64, fmod64(18.5, 4.2), 1.7, epsilon)); + assert(math.approxEq(f64, fmod64(23, 48.34), 23.0, epsilon)); + assert(math.approxEq(f64, fmod64(123.340890, 2398.2314), 123.340889, epsilon)); +} diff --git a/std/math/frexp.zig b/std/math/frexp.zig index 981d19d3c5..a96ae8f37b 100644 --- a/std/math/frexp.zig +++ b/std/math/frexp.zig @@ -1,89 +1,113 @@ -const assert = @import("../debug.zig").assert; const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; -pub fn frexp(x: var, e: &i32) -> @typeOf(x) { +fn frexp_result(comptime T: type) -> type { + struct { + significand: T, + exponent: i32, + } +} +pub const frexp32_result = frexp_result(f32); +pub const frexp64_result = frexp_result(f64); + +pub fn frexp(x: var) -> frexp_result(@typeOf(x)) { const T = @typeOf(x); switch (T) { - f32 => frexp32(x, e), - f64 => frexp64(x, e), + f32 => @inlineCall(frexp32, x), + f64 => @inlineCall(frexp64, x), else => @compileError("frexp not implemented for " ++ @typeName(T)), } } -fn frexp32(x_: f32, e: &i32) -> f32 { - var x = x_; +fn frexp32(x: f32) -> frexp32_result { + var result: frexp32_result = undefined; + var y = @bitCast(u32, x); - const ee = i32(y >> 23) & 0xFF; + const e = i32(y >> 23) & 0xFF; - if (ee == 0) { + if (e == 0) { if (x != 0) { - x = frexp32(x * 0x1.0p64, e); - *e -= 64; + // subnormal + result = frexp32(x * 0x1.0p64); + result.exponent -= 64; } else { - *e = 0; + // frexp(+-0) = (+-0, 0) + result.significand = x; + result.exponent = 0; } - return x; - } else if (ee == 0xFF) { - return x; + return result; + } else if (e == 0xFF) { + // frexp(nan) = (nan, 0) + result.significand = x; + result.exponent = 0; + return result; } - *e = ee - 0x7E; + result.exponent = e - 0x7E; y &= 0x807FFFFF; y |= 0x3F000000; - @bitCast(f32, y) + result.significand = @bitCast(f32, y); + result } -fn frexp64(x_: f64, e: &i32) -> f64 { - var x = x_; +fn frexp64(x: f64) -> frexp64_result { + var result: frexp64_result = undefined; + var y = @bitCast(u64, x); - const ee = i32(y >> 52) & 0x7FF; + const e = i32(y >> 52) & 0x7FF; - if (ee == 0) { + if (e == 0) { if (x != 0) { - x = frexp64(x * 0x1.0p64, e); - *e -= 64; + // subnormal + result = frexp64(x * 0x1.0p64); + result.exponent -= 64; } else { - *e = 0; + // frexp(+-0) = (+-0, 0) + result.significand = x; + result.exponent = 0; } - return x; - } else if (ee == 0x7FF) { - return x; + return result; + } else if (e == 0x7FF) { + // frexp(nan) = (nan, 0) + result.significand = x; + return result; } - *e = ee - 0x3FE; + result.exponent = e - 0x3FE; y &= 0x800FFFFFFFFFFFFF; y |= 0x3FE0000000000000; - @bitCast(f64, y) + result.significand = @bitCast(f64, y); + result } test "frexp" { - var i0: i32 = undefined; - var i1: i32 = undefined; + const a = frexp(f32(1.3)); + const b = frexp32(1.3); + assert(a.significand == b.significand and a.exponent == b.exponent); - assert(frexp(f32(1.3), &i0) == frexp32(1.3, &i1)); - assert(frexp(f64(1.3), &i0) == frexp64(1.3, &i1)); + const c = frexp(f64(1.3)); + const d = frexp64(1.3); + assert(c.significand == d.significand and c.exponent == d.exponent); } test "frexp32" { const epsilon = 0.000001; - var i: i32 = undefined; - var d: f32 = undefined; + var r: frexp32_result = undefined; - d = frexp32(1.3, &i); - assert(math.approxEq(f32, d, 0.65, epsilon) and i == 1); + r = frexp32(1.3); + assert(math.approxEq(f32, r.significand, 0.65, epsilon) and r.exponent == 1); - d = frexp32(78.0234, &i); - assert(math.approxEq(f32, d, 0.609558, epsilon) and i == 7); + r = frexp32(78.0234); + assert(math.approxEq(f32, r.significand, 0.609558, epsilon) and r.exponent == 7); } test "frexp64" { const epsilon = 0.000001; - var i: i32 = undefined; - var d: f64 = undefined; + var r: frexp64_result = undefined; - d = frexp64(1.3, &i); - assert(math.approxEq(f64, d, 0.65, epsilon) and i == 1); + r = frexp64(1.3); + assert(math.approxEq(f64, r.significand, 0.65, epsilon) and r.exponent == 1); - d = frexp64(78.0234, &i); - assert(math.approxEq(f64, d, 0.609558, epsilon) and i == 7); + r = frexp64(78.0234); + assert(math.approxEq(f64, r.significand, 0.609558, epsilon) and r.exponent == 7); } diff --git a/std/math/hypot.zig b/std/math/hypot.zig new file mode 100644 index 0000000000..61a37f1403 --- /dev/null +++ b/std/math/hypot.zig @@ -0,0 +1,135 @@ +const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; + +pub fn hypot(comptime T: type, x: T, y: T) -> T { + switch (T) { + f32 => @inlineCall(hypot32, x, y), + f64 => @inlineCall(hypot64, x, y), + else => @compileError("hypot not implemented for " ++ @typeName(T)), + } +} + +fn hypot32(x: f32, y: f32) -> f32 { + var ux = @bitCast(u32, x); + var uy = @bitCast(u32, y); + + ux &= @maxValue(u32) >> 1; + uy &= @maxValue(u32) >> 1; + if (ux < uy) { + const tmp = ux; + ux = uy; + uy = tmp; + } + + var xx = @bitCast(f32, ux); + var yy = @bitCast(f32, uy); + if (uy == 0xFF << 23) { + return yy; + } + if (ux >= 0xFF << 23 or uy == 0 or ux - uy >= (25 << 23)) { + return xx + yy; + } + + var z: f32 = 1.0; + if (ux >= (0x7F+60) << 23) { + z = 0x1.0p90; + xx *= 0x1.0p-90; + yy *= 0x1.0p-90; + } else if (uy < (0x7F-60) << 23) { + z = 0x1.0p-90; + xx *= 0x1.0p-90; + yy *= 0x1.0p-90; + } + + z * math.sqrt(f32(f64(x) * x + f64(y) * y)) +} + +fn sq(hi: &f64, lo: &f64, x: f64) { + const split: f64 = 0x1.0p27 + 1.0; + const xc = x * split; + const xh = x - xc + xc; + const xl = x - xh; + *hi = x * x; + *lo = xh * xh - *hi + 2 * xh * xl + xl * xl; +} + +fn hypot64(x: f64, y: f64) -> f64 { + var ux = @bitCast(u64, x); + var uy = @bitCast(u64, y); + + ux &= @maxValue(u64) >> 1; + uy &= @maxValue(u64) >> 1; + if (ux < uy) { + const tmp = ux; + ux = uy; + uy = tmp; + } + + const ex = ux >> 52; + const ey = uy >> 52; + var xx = @bitCast(f64, ux); + var yy = @bitCast(f64, uy); + + // hypot(inf, nan) == inf + if (ey == 0x7FF) { + return yy; + } + if (ex == 0x7FF or uy == 0) { + return xx; + } + + // hypot(x, y) ~= x + y * y / x / 2 with inexact for small y/x + if (ex - ey > 64) { + return xx + yy; + } + + var z: f64 = 1; + if (ex > 0x3FF + 510) { + z = 0x1.0p700; + xx *= 0x1.0p-700; + yy *= 0x1.0p-700; + } else if (ey < 0x3FF - 450) { + z = 0x1.0p-700; + xx *= 0x1.0p700; + yy *= 0x1.0p700; + } + + var hx: f64 = undefined; + var lx: f64 = undefined; + var hy: f64 = undefined; + var ly: f64 = undefined; + + sq(&hx, &lx, x); + sq(&hy, &ly, y); + + z * math.sqrt(ly + lx + hy + hx) +} + +test "hypot" { + assert(hypot(f32, 0.0, -1.2) == hypot32(0.0, -1.2)); + assert(hypot(f64, 0.0, -1.2) == hypot64(0.0, -1.2)); +} + +test "hypot32" { + const epsilon = 0.000001; + + assert(math.approxEq(f32, hypot32(0.0, -1.2), 1.2, epsilon)); + assert(math.approxEq(f32, hypot32(0.2, -0.34), 0.394462, epsilon)); + assert(math.approxEq(f32, hypot32(0.8923, 2.636890), 2.783772, epsilon)); + assert(math.approxEq(f32, hypot32(1.5, 5.25), 5.460083, epsilon)); + assert(math.approxEq(f32, hypot32(37.45, 159.835), 164.163742, epsilon)); + assert(math.approxEq(f32, hypot32(89.123, 382.028905), 392.286865, epsilon)); + assert(math.approxEq(f32, hypot32(123123.234375, 529428.707813), 543556.875, epsilon)); +} + +test "hypot64" { + const epsilon = 0.000001; + + assert(math.approxEq(f64, hypot64(0.0, -1.2), 1.2, epsilon)); + assert(math.approxEq(f64, hypot64(0.2, -0.34), 0.394462, epsilon)); + assert(math.approxEq(f64, hypot64(0.8923, 2.636890), 2.783772, epsilon)); + assert(math.approxEq(f64, hypot64(1.5, 5.25), 5.460082, epsilon)); + assert(math.approxEq(f64, hypot64(37.45, 159.835), 164.163728, epsilon)); + assert(math.approxEq(f64, hypot64(89.123, 382.028905), 392.286876, epsilon)); + assert(math.approxEq(f64, hypot64(123123.234375, 529428.707813), 543556.885247, epsilon)); +} diff --git a/std/math/ilogb.zig b/std/math/ilogb.zig new file mode 100644 index 0000000000..8528a19e22 --- /dev/null +++ b/std/math/ilogb.zig @@ -0,0 +1,100 @@ +const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; + +pub fn ilogb(x: var) -> i32 { + const T = @typeOf(x); + switch (T) { + f32 => @inlineCall(ilogb32, x), + f64 => @inlineCall(ilogb64, x), + else => @compileError("ilogb not implemented for " ++ @typeName(T)), + } +} + +// NOTE: Should these be exposed publically? +const fp_ilogbnan = -1 - i32(@maxValue(u32) >> 1); +const fp_ilogb0 = fp_ilogbnan; + +fn ilogb32(x: f32) -> i32 { + var u = @bitCast(u32, x); + var e = i32((u >> 23) & 0xFF); + + if (e == 0) { + u <<= 9; + if (u == 0) { + math.raiseInvalid(); + return fp_ilogb0; + } + + // subnormal + e = -0x7F; + while (u >> 31 == 0) : (u <<= 1) { + e -= 1; + } + return e; + } + + if (e == 0xFF) { + math.raiseInvalid(); + if (u << 9 != 0) { + return fp_ilogbnan; + } else { + return @maxValue(i32); + } + } + + e - 0x7F +} + +fn ilogb64(x: f64) -> i32 { + var u = @bitCast(u64, x); + var e = i32((u >> 52) & 0x7FF); + + if (e == 0) { + u <<= 12; + if (u == 0) { + math.raiseInvalid(); + return fp_ilogb0; + } + + // subnormal + e = -0x3FF; + while (u >> 63 == 0) : (u <<= 1) { + e -= 1; + } + return e; + } + + if (e == 0x7FF) { + math.raiseInvalid(); + if (u << 12 != 0) { + return fp_ilogbnan; + } else { + return @maxValue(i32); + } + } + + e - 0x3FF +} + +test "ilogb" { + assert(ilogb(f32(0.2)) == ilogb32(0.2)); + assert(ilogb(f64(0.2)) == ilogb64(0.2)); +} + +test "ilogb32" { + assert(ilogb32(0.0) == fp_ilogb0); + assert(ilogb32(0.5) == -1); + assert(ilogb32(0.8923) == -1); + assert(ilogb32(10.0) == 3); + assert(ilogb32(-123984) == 16); + assert(ilogb32(2398.23) == 11); +} + +test "ilogb64" { + assert(ilogb64(0.0) == fp_ilogb0); + assert(ilogb64(0.5) == -1); + assert(ilogb64(0.8923) == -1); + assert(ilogb64(10.0) == 3); + assert(ilogb64(-123984) == 16); + assert(ilogb64(2398.23) == 11); +} diff --git a/std/math/index.zig b/std/math/index.zig index 82e5af82ea..b031bba25d 100644 --- a/std/math/index.zig +++ b/std/math/index.zig @@ -1,7 +1,189 @@ -const assert = @import("../debug.zig").assert; const builtin = @import("builtin"); +const TypeId = builtin.TypeId; +const assert = @import("../debug.zig").assert; + +pub const e = 2.7182818284590452354; // e +pub const log2_e = 1.4426950408889634074; // log_2(e) +pub const log10_e = 0.43429448190325182765; // log_10(e) +pub const ln_2 = 0.69314718055994530942; // log_e(2) +pub const ln_10 = 2.30258509299404568402; // log_e(10) +pub const pi = 3.14159265358979323846; // pi +pub const pi_2 = 1.57079632679489661923; // pi/2 +pub const pi_4 = 0.78539816339744830962; // pi/4 +pub const r1_pi = 0.31830988618379067154; // 1/pi +pub const r2_pi = 0.63661977236758134308; // 2/pi +pub const r2_sqrtpi = 1.12837916709551257390; // 2/sqrt(pi) +pub const sqrt2 = 1.41421356237309504880; // sqrt(2) +pub const r1_sqrt2 = 0.70710678118654752440; // 1/sqrt(2) + +// float.h details +pub const f64_true_min = 4.94065645841246544177e-324; +pub const f64_min = 2.22507385850720138309e-308; +pub const f64_max = 1.79769313486231570815e+308; +pub const f64_epsilon = 2.22044604925031308085e-16; +pub const f64_toint = 1.0 / f64_epsilon; + +pub const f32_true_min = 1.40129846432481707092e-45; +pub const f32_min = 1.17549435082228750797e-38; +pub const f32_max = 3.40282346638528859812e+38; +pub const f32_epsilon = 1.1920928955078125e-07; +pub const f32_toint = 1.0 / f32_epsilon; + +pub const nan_u32 = u32(0x7F800001); +pub const nan_f32 = @bitCast(f32, nan_u32); + +pub const inf_u32 = u32(0x7F800000); +pub const inf_f32 = @bitCast(f32, inf_u32); + +pub const nan_u64 = u64(0x7FF << 52) | 1; +pub const nan_f64 = @bitCast(f64, nan_u64); + +pub const inf_u64 = u64(0x7FF << 52); +pub const inf_f64 = @bitCast(f64, inf_u64); + +pub const nan = @import("nan.zig").nan; +pub const inf = @import("inf.zig").inf; + +pub fn approxEq(comptime T: type, x: T, y: T, epsilon: T) -> bool { + assert(@typeId(T) == TypeId.Float); + fabs(x - y) < epsilon +} +// TODO: Hide the following in an internal module. +pub fn forceEval(value: var) { + const T = @typeOf(value); + switch (T) { + f32 => { + var x: f32 = undefined; + const p = @ptrCast(&volatile f32, &x); + *p = x; + }, + f64 => { + var x: f64 = undefined; + const p = @ptrCast(&volatile f64, &x); + *p = x; + }, + else => { + @compileError("forceEval not implemented for " ++ @typeName(T)); + }, + } +} + +pub fn raiseInvalid() { + // Raise INVALID fpu exception +} + +pub fn raiseUnderflow() { + // Raise UNDERFLOW fpu exception +} + +pub fn raiseOverflow() { + // Raise OVERFLOW fpu exception +} + +pub fn raiseInexact() { + // Raise INEXACT fpu exception +} + +pub fn raiseDivByZero() { + // Raise INEXACT fpu exception +} + +pub const isNan = @import("isnan.zig").isNan; +pub const fabs = @import("fabs.zig").fabs; +pub const ceil = @import("ceil.zig").ceil; +pub const floor = @import("floor.zig").floor; +pub const trunc = @import("floor.zig").trunc; +pub const round = @import("round.zig").round; pub const frexp = @import("frexp.zig").frexp; +pub const frexp32_result = @import("frexp.zig").frexp32_result; +pub const frexp64_result = @import("frexp.zig").frexp64_result; +pub const fmod = @import("fmod.zig").fmod; +pub const modf = @import("modf.zig").modf; +pub const modf32_result = @import("modf.zig").modf32_result; +pub const modf64_result = @import("modf.zig").modf64_result; +pub const copysign = @import("copysign.zig").copysign; +pub const isFinite = @import("isfinite.zig").isFinite; +pub const isInf = @import("isinf.zig").isInf; +pub const isPositiveInf = @import("isinf.zig").isPositiveInf; +pub const isNegativeInf = @import("isinf.zig").isNegativeInf; +pub const isNormal = @import("isnormal.zig").isNormal; +pub const signbit = @import("signbit.zig").signbit; +pub const scalbn = @import("scalbn.zig").scalbn; +pub const pow = @import("pow.zig").pow; +pub const sqrt = @import("sqrt.zig").sqrt; +pub const cbrt = @import("cbrt.zig").cbrt; +pub const acos = @import("acos.zig").acos; +pub const asin = @import("asin.zig").asin; +pub const atan = @import("atan.zig").atan; +pub const atan2 = @import("atan2.zig").atan2; +pub const hypot = @import("hypot.zig").hypot; +pub const exp = @import("exp.zig").exp; +pub const exp2 = @import("exp2.zig").exp2; +pub const expm1 = @import("expm1.zig").expm1; +pub const ilogb = @import("ilogb.zig").ilogb; +pub const ln = @import("ln.zig").ln; +pub const log = @import("log.zig").log; +pub const log2 = @import("log2.zig").log2; +pub const log10 = @import("log10.zig").log10; +pub const log1p = @import("log1p.zig").log1p; +pub const fma = @import("fma.zig").fma; +pub const asinh = @import("asinh.zig").asinh; +pub const acosh = @import("acosh.zig").acosh; +pub const atanh = @import("atanh.zig").atanh; +pub const sinh = @import("sinh.zig").sinh; +pub const cosh = @import("cosh.zig").cosh; +pub const tanh = @import("tanh.zig").tanh; +pub const cos = @import("cos.zig").cos; +pub const sin = @import("sin.zig").sin; +pub const tan = @import("tan.zig").tan; + +test "math" { + _ = @import("nan.zig"); + _ = @import("isnan.zig"); + _ = @import("fabs.zig"); + _ = @import("ceil.zig"); + _ = @import("floor.zig"); + _ = @import("trunc.zig"); + _ = @import("round.zig"); + _ = @import("frexp.zig"); + _ = @import("fmod.zig"); + _ = @import("modf.zig"); + _ = @import("copysign.zig"); + _ = @import("isfinite.zig"); + _ = @import("isinf.zig"); + _ = @import("isnormal.zig"); + _ = @import("signbit.zig"); + _ = @import("scalbn.zig"); + _ = @import("pow.zig"); + _ = @import("sqrt.zig"); + _ = @import("cbrt.zig"); + _ = @import("acos.zig"); + _ = @import("asin.zig"); + _ = @import("atan.zig"); + _ = @import("atan2.zig"); + _ = @import("hypot.zig"); + _ = @import("exp.zig"); + _ = @import("exp2.zig"); + _ = @import("expm1.zig"); + _ = @import("ilogb.zig"); + _ = @import("ln.zig"); + _ = @import("log.zig"); + _ = @import("log2.zig"); + _ = @import("log10.zig"); + _ = @import("log1p.zig"); + _ = @import("fma.zig"); + _ = @import("asinh.zig"); + _ = @import("acosh.zig"); + _ = @import("atanh.zig"); + _ = @import("sinh.zig"); + _ = @import("cosh.zig"); + _ = @import("tanh.zig"); + _ = @import("sin.zig"); + _ = @import("cos.zig"); + _ = @import("tan.zig"); +} + pub const Cmp = enum { Less, @@ -66,25 +248,6 @@ fn testOverflow() { } -pub fn log(comptime base: usize, value: var) -> @typeOf(value) { - const T = @typeOf(value); - switch (@typeId(T)) { - builtin.TypeId.Int => { - if (base == 2) { - return T.bit_count - 1 - @clz(value); - } else { - @compileError("TODO implement log for non base 2 integers"); - } - }, - builtin.TypeId.Float => { - @compileError("TODO implement log for floats"); - }, - else => { - @compileError("log expects integer or float, found '" ++ @typeName(T) ++ "'"); - }, - } -} - error Overflow; pub fn absInt(x: var) -> %@typeOf(x) { const T = @typeOf(x); @@ -244,92 +407,6 @@ fn testRem() { if (rem(f32, 10, 0)) |_| unreachable else |err| assert(err == error.DivisionByZero); } -fn isNan(comptime T: type, x: T) -> bool { - assert(@typeId(T) == builtin.TypeId.Float); - if (T == f32) { - const bits = @bitCast(u32, x); - return (bits & 0x7fffffff) > 0x7f800000; - } else if (T == f64) { - const bits = @bitCast(u64, x); - return (bits & (@maxValue(u64) >> 1)) > (u64(0x7ff) << 52); - } else if (T == c_longdouble) { - @compileError("TODO support isNan for c_longdouble"); - } else { - unreachable; - } -} - -pub fn floor(x: var) -> @typeOf(x) { - switch (@typeOf(x)) { - f32 => floor_f32(x), - f64 => floor_f64(x), - c_longdouble => @compileError("TODO support floor for c_longdouble"), - else => @compileError("Invalid type for floor: " ++ @typeName(@typeOf(x))), - } -} - -fn floor_f32(x: f32) -> f32 { - var i = @bitCast(u32, x); - const e = i32((i >> 23) & 0xff) -% 0x7f; - if (e >= 23) - return x; - if (e >= 0) { - const m = @bitCast(u32, 0x007fffff >> e); - if ((i & m) == 0) - return x; - if (i >> 31 != 0) - i +%= m; - i &= ~m; - } else { - if (i >> 31 == 0) - return 0; - if (i <<% 1 != 0) - return -1.0; - } - return @bitCast(f32, i); -} - -fn floor_f64(x: f64) -> f64 { - const DBL_EPSILON = 2.22044604925031308085e-16; - const toint = 1.0 / DBL_EPSILON; - - var i = @bitCast(u64, x); - const e = (i >> 52) & 0x7ff; - - if (e >= 0x3ff +% 52 or x == 0) - return x; - // y = int(x) - x, where int(x) is an integer neighbor of x - const y = { - @setFloatMode(this, builtin.FloatMode.Strict); - if (i >> 63 != 0) { - x - toint + toint - x - } else { - x + toint - toint - x - } - }; - // special case because of non-nearest rounding modes - if (e <= 0x3ff - 1) { - if (i >> 63 != 0) - return -1.0; - return 0.0; - } - if (y > 0) - return x + y - 1; - return x + y; -} - -test "math.floor" { - assert(floor(f32(1.234)) == 1.0); - assert(floor(f32(-1.234)) == -2.0); - assert(floor(f32(999.0)) == 999.0); - assert(floor(f32(-999.0)) == -999.0); - - assert(floor(f64(1.234)) == 1.0); - assert(floor(f64(-1.234)) == -2.0); - assert(floor(f64(999.0)) == 999.0); - assert(floor(f64(-999.0)) == -999.0); -} - /// Returns the absolute value of the integer parameter. /// Result is an unsigned integer. pub fn absCast(x: var) -> @IntType(false, @typeOf(x).bit_count) { @@ -377,13 +454,3 @@ test "math.negateCast" { if (negateCast(u32(@maxValue(i32) + 10))) |_| unreachable else |err| assert(err == error.Overflow); } - -test "math" { - _ = @import("frexp.zig"); -} - - -pub fn approxEq(comptime T: type, x: T, y: T, epsilon: T) -> bool { - comptime assert(@typeId(T) == builtin.TypeId.Float); - absFloat(x - y) < epsilon -} diff --git a/std/math/inf.zig b/std/math/inf.zig new file mode 100644 index 0000000000..ed559c313e --- /dev/null +++ b/std/math/inf.zig @@ -0,0 +1,10 @@ +const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; + +pub fn inf(comptime T: type) -> T { + switch (T) { + f32 => @bitCast(f32, math.inf_u32), + f64 => @bitCast(f64, math.inf_u64), + else => @compileError("inf not implemented for " ++ @typeName(T)), + } +} diff --git a/std/math/isfinite.zig b/std/math/isfinite.zig new file mode 100644 index 0000000000..a820cb2f52 --- /dev/null +++ b/std/math/isfinite.zig @@ -0,0 +1,30 @@ +const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; + +pub fn isFinite(x: var) -> bool { + const T = @typeOf(x); + switch (T) { + f32 => { + const bits = @bitCast(u32, x); + bits & 0x7FFFFFFF < 0x7F800000 + }, + f64 => { + const bits = @bitCast(u64, x); + bits & (@maxValue(u64) >> 1) < (0x7FF << 52) + }, + else => { + @compileError("isFinite not implemented for " ++ @typeName(T)); + }, + } +} + +test "isFinite" { + assert(isFinite(f32(0.0))); + assert(isFinite(f32(-0.0))); + assert(isFinite(f64(0.0))); + assert(isFinite(f64(-0.0))); + assert(!isFinite(math.inf(f32))); + assert(!isFinite(-math.inf(f32))); + assert(!isFinite(math.inf(f64))); + assert(!isFinite(-math.inf(f64))); +} diff --git a/std/math/isinf.zig b/std/math/isinf.zig new file mode 100644 index 0000000000..cfb8aae55a --- /dev/null +++ b/std/math/isinf.zig @@ -0,0 +1,82 @@ +const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; + +pub fn isInf(x: var) -> bool { + const T = @typeOf(x); + switch (T) { + f32 => { + const bits = @bitCast(u32, x); + bits & 0x7FFFFFFF == 0x7F800000 + }, + f64 => { + const bits = @bitCast(u64, x); + bits & (@maxValue(u64) >> 1) == (0x7FF << 52) + }, + else => { + @compileError("isInf not implemented for " ++ @typeName(T)); + }, + } +} + +pub fn isPositiveInf(x: var) -> bool { + const T = @typeOf(x); + switch (T) { + f32 => { + @bitCast(u32, x) == 0x7F800000 + }, + f64 => { + @bitCast(u64, x) == 0x7FF << 52 + }, + else => { + @compileError("isPositiveInf not implemented for " ++ @typeName(T)); + }, + } +} + +pub fn isNegativeInf(x: var) -> bool { + const T = @typeOf(x); + switch (T) { + f32 => { + @bitCast(u32, x) == 0xFF800000 + }, + f64 => { + @bitCast(u64, x) == 0xFFF << 52 + }, + else => { + @compileError("isNegativeInf not implemented for " ++ @typeName(T)); + }, + } +} + +test "isInf" { + assert(!isInf(f32(0.0))); + assert(!isInf(f32(-0.0))); + assert(!isInf(f64(0.0))); + assert(!isInf(f64(-0.0))); + assert(isInf(math.inf(f32))); + assert(isInf(-math.inf(f32))); + assert(isInf(math.inf(f64))); + assert(isInf(-math.inf(f64))); +} + +test "isPositiveInf" { + assert(!isPositiveInf(f32(0.0))); + assert(!isPositiveInf(f32(-0.0))); + assert(!isPositiveInf(f64(0.0))); + assert(!isPositiveInf(f64(-0.0))); + assert(isPositiveInf(math.inf(f32))); + assert(!isPositiveInf(-math.inf(f32))); + assert(isPositiveInf(math.inf(f64))); + assert(!isPositiveInf(-math.inf(f64))); +} + +test "isNegativeInf" { + assert(!isNegativeInf(f32(0.0))); + assert(!isNegativeInf(f32(-0.0))); + assert(!isNegativeInf(f64(0.0))); + assert(!isNegativeInf(f64(-0.0))); + assert(!isNegativeInf(math.inf(f32))); + assert(isNegativeInf(-math.inf(f32))); + assert(!isNegativeInf(math.inf(f64))); + assert(isNegativeInf(-math.inf(f64))); +} diff --git a/std/math/isnan.zig b/std/math/isnan.zig new file mode 100644 index 0000000000..6d3951a3f9 --- /dev/null +++ b/std/math/isnan.zig @@ -0,0 +1,26 @@ +const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; + +pub fn isNan(x: var) -> bool { + const T = @typeOf(x); + switch (T) { + f32 => { + const bits = @bitCast(u32, x); + bits & 0x7FFFFFFF > 0x7F800000 + }, + f64 => { + const bits = @bitCast(u64, x); + (bits & (@maxValue(u64) >> 1)) > (u64(0x7FF) << 52) + }, + else => { + @compileError("isNan not implemented for " ++ @typeName(T)); + }, + } +} + +test "isNan" { + assert(isNan(math.nan(f32))); + assert(isNan(math.nan(f64))); + assert(!isNan(f32(1.0))); + assert(!isNan(f64(1.0))); +} diff --git a/std/math/isnormal.zig b/std/math/isnormal.zig new file mode 100644 index 0000000000..597a7d8849 --- /dev/null +++ b/std/math/isnormal.zig @@ -0,0 +1,26 @@ +const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; + +pub fn isNormal(x: var) -> bool { + const T = @typeOf(x); + switch (T) { + f32 => { + const bits = @bitCast(u32, x); + (bits + 0x00800000) & 0x7FFFFFFF >= 0x01000000 + }, + f64 => { + const bits = @bitCast(u64, x); + (bits + (1 << 52)) & (@maxValue(u64) >> 1) >= (1 << 53) + }, + else => { + @compileError("isNormal not implemented for " ++ @typeName(T)); + }, + } +} + +test "isNormal" { + assert(!isNormal(math.nan(f32))); + assert(!isNormal(math.nan(f64))); + assert(isNormal(f32(1.0))); + assert(isNormal(f64(1.0))); +} diff --git a/std/math/ln.zig b/std/math/ln.zig new file mode 100644 index 0000000000..aff84303b7 --- /dev/null +++ b/std/math/ln.zig @@ -0,0 +1,148 @@ +const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; + +pub fn ln(x: var) -> @typeOf(x) { + const T = @typeOf(x); + switch (T) { + f32 => @inlineCall(lnf, x), + f64 => @inlineCall(lnd, x), + else => @compileError("ln not implemented for " ++ @typeName(T)), + } +} + +fn lnf(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 -1 / (x * x); + } + // log(-#) = nan + if (ix >> 31 != 0) { + return (x - x) / 0.0 + } + + // 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 += 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 = f32(k); + + s * (hfsq + R) + dk * ln2_lo - hfsq + f + dk * ln2_hi +} + +fn lnd(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 = u32(ix >> 32); + var k: i32 = 0; + + if (hx < 0x00100000 or hx >> 31 != 0) { + // log(+-0) = -inf + if (ix << 1 == 0) { + return -1 / (x * x); + } + // log(-#) = nan + if (hx >> 31 != 0) { + return (x - x) / 0.0; + } + + // subnormal, scale x + k -= 54; + x *= 0x1.0p54; + hx = 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 += i32(hx >> 20) - 0x3FF; + hx = (hx & 0x000FFFFF) + 0x3FE6A09E; + ix = (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 = f64(k); + + s * (hfsq + R) + dk * ln2_lo - hfsq + f + dk * ln2_hi +} + +test "log" { + assert(ln(f32(0.2)) == lnf(0.2)); + assert(ln(f64(0.2)) == lnd(0.2)); +} + +test "logf" { + const epsilon = 0.000001; + + assert(math.approxEq(f32, lnf(0.2), -1.609438, epsilon)); + assert(math.approxEq(f32, lnf(0.8923), -0.113953, epsilon)); + assert(math.approxEq(f32, lnf(1.5), 0.405465, epsilon)); + assert(math.approxEq(f32, lnf(37.45), 3.623007, epsilon)); + assert(math.approxEq(f32, lnf(89.123), 4.490017, epsilon)); + assert(math.approxEq(f32, lnf(123123.234375), 11.720941, epsilon)); +} + +test "logd" { + const epsilon = 0.000001; + + assert(math.approxEq(f64, lnd(0.2), -1.609438, epsilon)); + assert(math.approxEq(f64, lnd(0.8923), -0.113953, epsilon)); + assert(math.approxEq(f64, lnd(1.5), 0.405465, epsilon)); + assert(math.approxEq(f64, lnd(37.45), 3.623007, epsilon)); + assert(math.approxEq(f64, lnd(89.123), 4.490017, epsilon)); + assert(math.approxEq(f64, lnd(123123.234375), 11.720941, epsilon)); +} diff --git a/std/math/log.zig b/std/math/log.zig new file mode 100644 index 0000000000..abd480b6a6 --- /dev/null +++ b/std/math/log.zig @@ -0,0 +1,72 @@ +const math = @import("index.zig"); +const builtin = @import("builtin"); +const assert = @import("../debug.zig").assert; + +pub fn log(comptime base: usize, x: var) -> @typeOf(x) { + const T = @typeOf(x); + switch (@typeId(T)) { + builtin.TypeId.Int => { + if (base == 2) { + return T.bit_count - 1 - @clz(x); + } else { + @compileError("TODO implement log for non base 2 integers"); + } + }, + + builtin.TypeId.Float => { + return logf(base, x); + }, + + else => { + @compileError("log expects integer or float, found '" ++ @typeName(T) ++ "'"); + }, + } +} + +fn logf(comptime base: usize, x: var) -> @typeOf(x) { + const T = @typeOf(x); + switch (T) { + f32 => { + switch (base) { + 2 => return math.log2(x), + 10 => return math.log10(x), + else => return f32(math.ln(f64(x)) / math.ln(f64(base))), + } + }, + + f64 => { + switch (base) { + 2 => return math.log2(x), + 10 => return math.log10(x), + // NOTE: This likely is computed with reduced accuracy. + else => return math.ln(x) / math.ln(f64(base)), + } + }, + + else => @compileError("log not implemented for " ++ @typeName(T)), + } +} + +test "log_integer" { + assert(log(2, u8(0x1)) == 0); + assert(log(2, u8(0x2)) == 1); + assert(log(2, i16(0x72)) == 6); + assert(log(2, u32(0xFFFFFF)) == 23); + assert(log(2, u64(0x7FF0123456789ABC)) == 62); +} + +test "log_float" { + const epsilon = 0.000001; + + assert(math.approxEq(f32, log(6, f32(0.23947)), -0.797723, epsilon)); + assert(math.approxEq(f32, log(89, f32(0.23947)), -0.318432, epsilon)); + assert(math.approxEq(f64, log(123897, f64(12389216414)), 1.981724596, epsilon)); +} + +test "log_float_special" { + assert(log(2, f32(0.2301974)) == math.log2(f32(0.2301974))); + assert(log(10, f32(0.2301974)) == math.log10(f32(0.2301974))); + + assert(log(2, f64(213.23019799993)) == math.log2(f64(213.23019799993))); + assert(log(10, f64(213.23019799993)) == math.log10(f64(213.23019799993))); +} diff --git a/std/math/log10.zig b/std/math/log10.zig new file mode 100644 index 0000000000..4a4d2a02e7 --- /dev/null +++ b/std/math/log10.zig @@ -0,0 +1,175 @@ +const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; + +pub fn log10(x: var) -> @typeOf(x) { + const T = @typeOf(x); + switch (T) { + f32 => @inlineCall(log10f, x), + f64 => @inlineCall(log10d, x), + else => @compileError("log10 not implemented for " ++ @typeName(T)), + } +} + +fn log10f(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 -1 / (x * x); + } + // log(-#) = nan + if (ix >> 31 != 0) { + return (x - x) / 0.0 + } + + 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 += 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 = f32(k); + + dk * log10_2lo + (lo + hi) * ivln10lo + lo * ivln10hi + hi * ivln10hi + dk * log10_2hi +} + +fn log10d(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 = u32(ix >> 32); + var k: i32 = 0; + + if (hx < 0x00100000 or hx >> 31 != 0) { + // log(+-0) = -inf + if (ix << 1 == 0) { + return -1 / (x * x); + } + // log(-#) = nan + if (hx >> 31 != 0) { + return (x - x) / 0.0; + } + + // subnormal, scale x + k -= 54; + x *= 0x1.0p54; + hx = 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 += i32(hx >> 20) - 0x3FF; + hx = (hx & 0x000FFFFF) + 0x3FE6A09E; + ix = (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 &= @maxValue(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 = 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; + + val_lo + val_hi +} + +test "log10" { + assert(log10(f32(0.2)) == log10f(0.2)); + assert(log10(f64(0.2)) == log10d(0.2)); +} + +test "log10f" { + const epsilon = 0.000001; + + assert(math.approxEq(f32, log10f(0.2), -0.698970, epsilon)); + assert(math.approxEq(f32, log10f(0.8923), -0.049489, epsilon)); + assert(math.approxEq(f32, log10f(1.5), 0.176091, epsilon)); + assert(math.approxEq(f32, log10f(37.45), 1.573452, epsilon)); + assert(math.approxEq(f32, log10f(89.123), 1.94999, epsilon)); + assert(math.approxEq(f32, log10f(123123.234375), 5.09034, epsilon)); +} + +test "log10d" { + const epsilon = 0.000001; + + assert(math.approxEq(f64, log10d(0.2), -0.698970, epsilon)); + assert(math.approxEq(f64, log10d(0.8923), -0.049489, epsilon)); + assert(math.approxEq(f64, log10d(1.5), 0.176091, epsilon)); + assert(math.approxEq(f64, log10d(37.45), 1.573452, epsilon)); + assert(math.approxEq(f64, log10d(89.123), 1.94999, epsilon)); + assert(math.approxEq(f64, log10d(123123.234375), 5.09034, epsilon)); +} diff --git a/std/math/log1p.zig b/std/math/log1p.zig new file mode 100644 index 0000000000..811bd9e728 --- /dev/null +++ b/std/math/log1p.zig @@ -0,0 +1,197 @@ +const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; + +pub fn log1p(x: var) -> @typeOf(x) { + const T = @typeOf(x); + switch (T) { + f32 => @inlineCall(log1pf, x), + f64 => @inlineCall(log1pd, x), + else => @compileError("log1p not implemented for " ++ @typeName(T)), + } +} + +fn log1pf(x: f32) -> f32 { + const ln2_hi = 6.9313812256e-01; + const ln2_lo = 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; + + const u = @bitCast(u32, x); + var ix = u; + var k: i32 = 1; + var f: f32 = undefined; + var c: f32 = undefined; + + // 1 + x < sqrt(2)+ + if (ix < 0x3ED413D0 or ix >> 31 != 0) { + // x <= -1.0 + if (ix >= 0xBF800000) { + // log1p(-1) = +inf + if (x == -1) { + return x / 0.0; + } + // log1p(x < -1) = nan + else { + return (x - x) / 0.0; + } + } + // |x| < 2^(-24) + if ((ix << 1) < (0x33800000 << 1)) { + // underflow if subnormal + if (ix & 0x7F800000 == 0) { + math.forceEval(x * x); + } + return x; + } + // sqrt(2) / 2- <= 1 + x < sqrt(2)+ + if (ix <= 0xBE95F619) { + k = 0; + c = 0; + f = x; + } + } else if (ix >= 0x7F800000) { + return x; + } + + if (k != 0) { + const uf = 1 + x; + var iu = @bitCast(u32, uf); + iu += 0x3F800000 - 0x3F3504F3; + k = i32(iu >> 23) - 0x7F; + + // correction to avoid underflow in c / u + if (k < 25) { + c = if (k >= 2) 1 - (uf - x) else x - (uf - 1); + c /= uf; + } else { + c = 0; + } + + // u into [sqrt(2)/2, sqrt(2)] + iu = (iu & 0x007FFFFF) + 0x3F3504F3; + f = @bitCast(f32, iu) - 1; + } + + 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 = f32(k); + + s * (hfsq + R) + (dk * ln2_lo + c) - hfsq + f + dk * ln2_hi +} + +fn log1pd(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 ix = @bitCast(u64, x); + var hx = u32(ix >> 32); + var k: i32 = 1; + var c: f64 = undefined; + var f: f64 = undefined; + + // 1 + x < sqrt(2) + if (hx < 0x3FDA827A or hx >> 31 != 0) { + // x <= -1.0 + if (hx >= 0xBFF00000) { + // log1p(-1) = -inf + if (x == 1) { + return x / 0.0; + } + // log1p(x < -1) = nan + else { + return (x - x) / 0.0; + } + } + // |x| < 2^(-53) + if ((hx << 1) < (0x3CA00000 << 1)) { + if ((hx & 0x7FF00000) == 0) { + math.raiseUnderflow(); + } + return x; + } + // sqrt(2) / 2- <= 1 + x < sqrt(2)+ + if (hx <= 0xBFD2BEC4) { + k = 0; + c = 0; + f = x; + } + } + else if (hx >= 0x7FF00000) { + return x; + } + + if (k != 0) { + const uf = 1 + x; + const hu = @bitCast(u64, uf); + var iu = u32(hu >> 32); + iu += 0x3FF00000 - 0x3FE6A09E; + k = i32(iu >> 20) - 0x3FF; + + // correction to avoid underflow in c / u + if (k < 54) { + c = if (k >= 2) 1 - (uf - x) else x - (uf - 1); + c /= uf; + } else { + c = 0; + } + + // u into [sqrt(2)/2, sqrt(2)] + iu = (iu & 0x000FFFFF) + 0x3FE6A09E; + const iq = (u64(iu) << 32) | (hu & 0xFFFFFFFF); + f = @bitCast(f64, iq) - 1; + } + + 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 = f64(k); + + s * (hfsq + R) + (dk * ln2_lo + c) - hfsq + f + dk * ln2_hi +} + +test "log1p" { + assert(log1p(f32(0.0)) == log1pf(0.0)); + assert(log1p(f64(0.0)) == log1pd(0.0)); +} + +test "log1pf" { + const epsilon = 0.000001; + + assert(math.approxEq(f32, log1pf(0.0), 0.0, epsilon)); + assert(math.approxEq(f32, log1pf(0.2), 0.182322, epsilon)); + assert(math.approxEq(f32, log1pf(0.8923), 0.637793, epsilon)); + assert(math.approxEq(f32, log1pf(1.5), 0.916291, epsilon)); + assert(math.approxEq(f32, log1pf(37.45), 3.649359, epsilon)); + assert(math.approxEq(f32, log1pf(89.123), 4.501175, epsilon)); + assert(math.approxEq(f32, log1pf(123123.234375), 11.720949, epsilon)); +} + +test "log1pd" { + const epsilon = 0.000001; + + assert(math.approxEq(f64, log1pd(0.0), 0.0, epsilon)); + assert(math.approxEq(f64, log1pd(0.2), 0.182322, epsilon)); + assert(math.approxEq(f64, log1pd(0.8923), 0.637793, epsilon)); + assert(math.approxEq(f64, log1pd(1.5), 0.916291, epsilon)); + assert(math.approxEq(f64, log1pd(37.45), 3.649359, epsilon)); + assert(math.approxEq(f64, log1pd(89.123), 4.501175, epsilon)); + assert(math.approxEq(f64, log1pd(123123.234375), 11.720949, epsilon)); +} diff --git a/std/math/log2.zig b/std/math/log2.zig new file mode 100644 index 0000000000..67eb2c5e08 --- /dev/null +++ b/std/math/log2.zig @@ -0,0 +1,165 @@ +const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; + +pub fn log2(x: var) -> @typeOf(x) { + const T = @typeOf(x); + switch (T) { + f32 => @inlineCall(log2f, x), + f64 => @inlineCall(log2d, x), + else => @compileError("log2 not implemented for " ++ @typeName(T)), + } +} + +fn log2f(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 -1 / (x * x); + } + // log(-#) = nan + if (ix >> 31 != 0) { + return (x - x) / 0.0 + } + + 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 += 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); + (lo + hi) * ivln2lo + lo * ivln2hi + hi * ivln2hi + f32(k) +} + +fn log2d(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 = u32(ix >> 32); + var k: i32 = 0; + + if (hx < 0x00100000 or hx >> 31 != 0) { + // log(+-0) = -inf + if (ix << 1 == 0) { + return -1 / (x * x); + } + // log(-#) = nan + if (hx >> 31 != 0) { + return (x - x) / 0.0; + } + + // subnormal, scale x + k -= 54; + x *= 0x1.0p54; + hx = 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 += i32(hx >> 20) - 0x3FF; + hx = (hx & 0x000FFFFF) + 0x3FE6A09E; + ix = (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 &= @maxValue(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 = f64(k); + const ww = y + val_hi; + val_lo += (y - ww) + val_hi; + val_hi = ww; + + val_lo + val_hi +} + +test "log2" { + assert(log2(f32(0.2)) == log2f(0.2)); + assert(log2(f64(0.2)) == log2d(0.2)); +} + +test "log2f" { + const epsilon = 0.000001; + + assert(math.approxEq(f32, log2f(0.2), -2.321928, epsilon)); + assert(math.approxEq(f32, log2f(0.8923), -0.164399, epsilon)); + assert(math.approxEq(f32, log2f(1.5), 0.584962, epsilon)); + assert(math.approxEq(f32, log2f(37.45), 5.226894, epsilon)); + assert(math.approxEq(f32, log2f(123123.234375), 16.909744, epsilon)); +} + +test "log2d" { + const epsilon = 0.000001; + + assert(math.approxEq(f64, log2d(0.2), -2.321928, epsilon)); + assert(math.approxEq(f64, log2d(0.8923), -0.164399, epsilon)); + assert(math.approxEq(f64, log2d(1.5), 0.584962, epsilon)); + assert(math.approxEq(f64, log2d(37.45), 5.226894, epsilon)); + assert(math.approxEq(f64, log2d(123123.234375), 16.909744, epsilon)); +} diff --git a/std/math/modf.zig b/std/math/modf.zig new file mode 100644 index 0000000000..666fa04c81 --- /dev/null +++ b/std/math/modf.zig @@ -0,0 +1,157 @@ +const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; + +fn modf_result(comptime T: type) -> type { + struct { + fpart: T, + ipart: T, + } +} +pub const modf32_result = modf_result(f32); +pub const modf64_result = modf_result(f64); + +pub fn modf(x: var) -> modf_result(@typeOf(x)) { + const T = @typeOf(x); + switch (T) { + f32 => @inlineCall(modf32, x), + f64 => @inlineCall(modf64, x), + else => @compileError("modf not implemented for " ++ @typeName(T)), + } +} + +fn modf32(x: f32) -> modf32_result { + var result: modf32_result = undefined; + + const u = @bitCast(u32, x); + const e = i32((u >> 23) & 0xFF) - 0x7F; + const us = u & 0x80000000; + + // no fractional part + if (e >= 23) { + result.ipart = x; + if (e == 0x80 and u << 9 != 0) { // nan + result.fpart = x; + } else { + result.fpart = @bitCast(f32, us); + } + return result; + } + + // no integral part + if (e < 0) { + result.ipart = @bitCast(f32, us); + result.fpart = x; + return result; + } + + const mask = 0x007FFFFF >> u32(e); + if (u & mask == 0) { + result.ipart = x; + result.fpart = @bitCast(f32, us); + return result; + } + + const uf = @bitCast(f32, u & ~mask); + result.ipart = uf; + result.fpart = x - uf; + result +} + +fn modf64(x: f64) -> modf64_result { + var result: modf64_result = undefined; + + const u = @bitCast(u64, x); + const e = i32((u >> 52) & 0x7FF) - 0x3FF; + const us = u & (1 << 63); + + // no fractional part + if (e >= 52) { + result.ipart = x; + if (e == 0x400 and u << 12 != 0) { // nan + result.fpart = x; + } else { + result.fpart = @bitCast(f64, us); + } + return result; + } + + // no integral part + if (e < 0) { + result.ipart = @bitCast(f64, us); + result.fpart = x; + return result; + } + + const mask = @maxValue(u64) >> 12 >> u64(e); + if (u & mask == 0) { + result.ipart = x; + result.fpart = @bitCast(f64, us); + return result; + } + + const uf = @bitCast(f64, u & ~mask); + result.ipart = uf; + result.fpart = x - uf; + result +} + +test "modf" { + const a = modf(f32(1.0)); + const b = modf32(1.0); + // NOTE: No struct comparison on generic return type function? non-named, makes sense, but still. + assert(a.ipart == b.ipart and a.fpart == b.fpart); + + const c = modf(f64(1.0)); + const d = modf64(1.0); + assert(a.ipart == b.ipart and a.fpart == b.fpart); +} + +test "modf32" { + const epsilon = 0.000001; + var r: modf32_result = undefined; + + r = modf32(1.0); + assert(math.approxEq(f32, r.ipart, 1.0, epsilon)); + assert(math.approxEq(f32, r.fpart, 0.0, epsilon)); + + r = modf32(2.545); + assert(math.approxEq(f32, r.ipart, 2.0, epsilon)); + assert(math.approxEq(f32, r.fpart, 0.545, epsilon)); + + r = modf32(3.978123); + assert(math.approxEq(f32, r.ipart, 3.0, epsilon)); + assert(math.approxEq(f32, r.fpart, 0.978123, epsilon)); + + r = modf32(43874.3); + assert(math.approxEq(f32, r.ipart, 43874, epsilon)); + assert(math.approxEq(f32, r.fpart, 0.300781, epsilon)); + + r = modf32(1234.340780); + assert(math.approxEq(f32, r.ipart, 1234, epsilon)); + assert(math.approxEq(f32, r.fpart, 0.340820, epsilon)); +} + +test "modf64" { + const epsilon = 0.000001; + var r: modf64_result = undefined; + + r = modf64(1.0); + assert(math.approxEq(f64, r.ipart, 1.0, epsilon)); + assert(math.approxEq(f64, r.fpart, 0.0, epsilon)); + + r = modf64(2.545); + assert(math.approxEq(f64, r.ipart, 2.0, epsilon)); + assert(math.approxEq(f64, r.fpart, 0.545, epsilon)); + + r = modf64(3.978123); + assert(math.approxEq(f64, r.ipart, 3.0, epsilon)); + assert(math.approxEq(f64, r.fpart, 0.978123, epsilon)); + + r = modf64(43874.3); + assert(math.approxEq(f64, r.ipart, 43874, epsilon)); + assert(math.approxEq(f64, r.fpart, 0.3, epsilon)); + + r = modf64(1234.340780); + assert(math.approxEq(f64, r.ipart, 1234, epsilon)); + assert(math.approxEq(f64, r.fpart, 0.340780, epsilon)); +} diff --git a/std/math/nan.zig b/std/math/nan.zig new file mode 100644 index 0000000000..bd079026a5 --- /dev/null +++ b/std/math/nan.zig @@ -0,0 +1,9 @@ +const math = @import("index.zig"); + +pub fn nan(comptime T: type) -> T { + switch (T) { + f32 => @bitCast(f32, math.nan_u32), + f64 => @bitCast(f64, math.nan_u64), + else => @compileError("nan not implemented for " ++ @typeName(T)), + } +} diff --git a/std/math/oindex.zig b/std/math/oindex.zig new file mode 100644 index 0000000000..312833249c --- /dev/null +++ b/std/math/oindex.zig @@ -0,0 +1,274 @@ +const assert = @import("../debug.zig").assert; +const builtin = @import("builtin"); + +pub const frexp = @import("frexp.zig").frexp; + +pub const Cmp = enum { + Less, + Equal, + Greater, +}; + +pub fn min(x: var, y: var) -> @typeOf(x + y) { + if (x < y) x else y +} + +test "math.min" { + assert(min(i32(-1), i32(2)) == -1); +} + +pub fn max(x: var, y: var) -> @typeOf(x + y) { + if (x > y) x else y +} + +test "math.max" { + assert(max(i32(-1), i32(2)) == 2); +} + +error Overflow; +pub fn mul(comptime T: type, a: T, b: T) -> %T { + var answer: T = undefined; + if (@mulWithOverflow(T, a, b, &answer)) error.Overflow else answer +} + +error Overflow; +pub fn add(comptime T: type, a: T, b: T) -> %T { + var answer: T = undefined; + if (@addWithOverflow(T, a, b, &answer)) error.Overflow else answer +} + +error Overflow; +pub fn sub(comptime T: type, a: T, b: T) -> %T { + var answer: T = undefined; + if (@subWithOverflow(T, a, b, &answer)) error.Overflow else answer +} + +pub fn negate(x: var) -> %@typeOf(x) { + return sub(@typeOf(x), 0, x); +} + +error Overflow; +pub fn shl(comptime T: type, a: T, b: T) -> %T { + var answer: T = undefined; + if (@shlWithOverflow(T, a, b, &answer)) error.Overflow else answer +} + +test "math overflow functions" { + testOverflow(); + comptime testOverflow(); +} + +fn testOverflow() { + assert(%%mul(i32, 3, 4) == 12); + assert(%%add(i32, 3, 4) == 7); + assert(%%sub(i32, 3, 4) == -1); + assert(%%shl(i32, 0b11, 4) == 0b110000); +} + + +error Overflow; +pub fn absInt(x: var) -> %@typeOf(x) { + const T = @typeOf(x); + comptime assert(@typeId(T) == builtin.TypeId.Int); // must pass an integer to absInt + comptime assert(T.is_signed); // must pass a signed integer to absInt + if (x == @minValue(@typeOf(x))) + return error.Overflow; + { + @setDebugSafety(this, false); + return if (x < 0) -x else x; + } +} + +test "math.absInt" { + testAbsInt(); + comptime testAbsInt(); +} +fn testAbsInt() { + assert(%%absInt(i32(-10)) == 10); + assert(%%absInt(i32(10)) == 10); +} + +pub const absFloat = @import("fabs.zig").fabs; + +error DivisionByZero; +error Overflow; +pub fn divTrunc(comptime T: type, numerator: T, denominator: T) -> %T { + @setDebugSafety(this, false); + if (denominator == 0) + return error.DivisionByZero; + if (@typeId(T) == builtin.TypeId.Int and T.is_signed and numerator == @minValue(T) and denominator == -1) + return error.Overflow; + return @divTrunc(numerator, denominator); +} + +test "math.divTrunc" { + testDivTrunc(); + comptime testDivTrunc(); +} +fn testDivTrunc() { + assert(%%divTrunc(i32, 5, 3) == 1); + assert(%%divTrunc(i32, -5, 3) == -1); + if (divTrunc(i8, -5, 0)) |_| unreachable else |err| assert(err == error.DivisionByZero); + if (divTrunc(i8, -128, -1)) |_| unreachable else |err| assert(err == error.Overflow); + + assert(%%divTrunc(f32, 5.0, 3.0) == 1.0); + assert(%%divTrunc(f32, -5.0, 3.0) == -1.0); +} + +error DivisionByZero; +error Overflow; +pub fn divFloor(comptime T: type, numerator: T, denominator: T) -> %T { + @setDebugSafety(this, false); + if (denominator == 0) + return error.DivisionByZero; + if (@typeId(T) == builtin.TypeId.Int and T.is_signed and numerator == @minValue(T) and denominator == -1) + return error.Overflow; + return @divFloor(numerator, denominator); +} + +test "math.divFloor" { + testDivFloor(); + comptime testDivFloor(); +} +fn testDivFloor() { + assert(%%divFloor(i32, 5, 3) == 1); + assert(%%divFloor(i32, -5, 3) == -2); + if (divFloor(i8, -5, 0)) |_| unreachable else |err| assert(err == error.DivisionByZero); + if (divFloor(i8, -128, -1)) |_| unreachable else |err| assert(err == error.Overflow); + + assert(%%divFloor(f32, 5.0, 3.0) == 1.0); + assert(%%divFloor(f32, -5.0, 3.0) == -2.0); +} + +error DivisionByZero; +error Overflow; +error UnexpectedRemainder; +pub fn divExact(comptime T: type, numerator: T, denominator: T) -> %T { + @setDebugSafety(this, false); + if (denominator == 0) + return error.DivisionByZero; + if (@typeId(T) == builtin.TypeId.Int and T.is_signed and numerator == @minValue(T) and denominator == -1) + return error.Overflow; + const result = @divTrunc(numerator, denominator); + if (result * denominator != numerator) + return error.UnexpectedRemainder; + return result; +} + +test "math.divExact" { + testDivExact(); + comptime testDivExact(); +} +fn testDivExact() { + assert(%%divExact(i32, 10, 5) == 2); + assert(%%divExact(i32, -10, 5) == -2); + if (divExact(i8, -5, 0)) |_| unreachable else |err| assert(err == error.DivisionByZero); + if (divExact(i8, -128, -1)) |_| unreachable else |err| assert(err == error.Overflow); + if (divExact(i32, 5, 2)) |_| unreachable else |err| assert(err == error.UnexpectedRemainder); + + assert(%%divExact(f32, 10.0, 5.0) == 2.0); + assert(%%divExact(f32, -10.0, 5.0) == -2.0); + if (divExact(f32, 5.0, 2.0)) |_| unreachable else |err| assert(err == error.UnexpectedRemainder); +} + +error DivisionByZero; +error NegativeDenominator; +pub fn mod(comptime T: type, numerator: T, denominator: T) -> %T { + @setDebugSafety(this, false); + if (denominator == 0) + return error.DivisionByZero; + if (denominator < 0) + return error.NegativeDenominator; + return @mod(numerator, denominator); +} + +test "math.mod" { + testMod(); + comptime testMod(); +} +fn testMod() { + assert(%%mod(i32, -5, 3) == 1); + assert(%%mod(i32, 5, 3) == 2); + if (mod(i32, 10, -1)) |_| unreachable else |err| assert(err == error.NegativeDenominator); + if (mod(i32, 10, 0)) |_| unreachable else |err| assert(err == error.DivisionByZero); + + assert(%%mod(f32, -5, 3) == 1); + assert(%%mod(f32, 5, 3) == 2); + if (mod(f32, 10, -1)) |_| unreachable else |err| assert(err == error.NegativeDenominator); + if (mod(f32, 10, 0)) |_| unreachable else |err| assert(err == error.DivisionByZero); +} + +error DivisionByZero; +error NegativeDenominator; +pub fn rem(comptime T: type, numerator: T, denominator: T) -> %T { + @setDebugSafety(this, false); + if (denominator == 0) + return error.DivisionByZero; + if (denominator < 0) + return error.NegativeDenominator; + return @rem(numerator, denominator); +} + +test "math.rem" { + testRem(); + comptime testRem(); +} +fn testRem() { + assert(%%rem(i32, -5, 3) == -2); + assert(%%rem(i32, 5, 3) == 2); + if (rem(i32, 10, -1)) |_| unreachable else |err| assert(err == error.NegativeDenominator); + if (rem(i32, 10, 0)) |_| unreachable else |err| assert(err == error.DivisionByZero); + + assert(%%rem(f32, -5, 3) == -2); + assert(%%rem(f32, 5, 3) == 2); + if (rem(f32, 10, -1)) |_| unreachable else |err| assert(err == error.NegativeDenominator); + if (rem(f32, 10, 0)) |_| unreachable else |err| assert(err == error.DivisionByZero); +} + +/// Returns the absolute value of the integer parameter. +/// Result is an unsigned integer. +pub fn absCast(x: var) -> @IntType(false, @typeOf(x).bit_count) { + const uint = @IntType(false, @typeOf(x).bit_count); + if (x >= 0) + return uint(x); + + return uint(-(x + 1)) + 1; +} + +test "math.absCast" { + assert(absCast(i32(-999)) == 999); + assert(@typeOf(absCast(i32(-999))) == u32); + + assert(absCast(i32(999)) == 999); + assert(@typeOf(absCast(i32(999))) == u32); + + assert(absCast(i32(@minValue(i32))) == -@minValue(i32)); + assert(@typeOf(absCast(i32(@minValue(i32)))) == u32); +} + +/// Returns the negation of the integer parameter. +/// Result is a signed integer. +error Overflow; +pub fn negateCast(x: var) -> %@IntType(true, @typeOf(x).bit_count) { + if (@typeOf(x).is_signed) + return negate(x); + + const int = @IntType(true, @typeOf(x).bit_count); + if (x > -@minValue(int)) + return error.Overflow; + + if (x == -@minValue(int)) + return @minValue(int); + + return -int(x); +} + +test "math.negateCast" { + assert(%%negateCast(u32(999)) == -999); + assert(@typeOf(%%negateCast(u32(999))) == i32); + + assert(%%negateCast(u32(-@minValue(i32))) == @minValue(i32)); + assert(@typeOf(%%negateCast(u32(-@minValue(i32)))) == i32); + + if (negateCast(u32(@maxValue(i32) + 10))) |_| unreachable else |err| assert(err == error.Overflow); +} diff --git a/std/math/pow.zig b/std/math/pow.zig new file mode 100644 index 0000000000..ab281b67a0 --- /dev/null +++ b/std/math/pow.zig @@ -0,0 +1,316 @@ +const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; + +pub fn pow(comptime T: type, x: T, y: T) -> T { + switch (T) { + f32 => @inlineCall(pow32, x, y), + f64 => @inlineCall(pow64, x, y), + else => @compileError("pow not implemented for " ++ @typeName(T)), + } +} + +fn isOddInteger(x: f64) -> bool { + const r = math.modf(x); + r.fpart == 0.0 and i64(r.ipart) & 1 == 1 +} + +// This implementation is taken from the go stlib, musl is a bit more complex. +fn pow32(x: f32, y: f32) -> f32 { + // pow(x, +-0) = 1 for all x + // pow(1, y) = 1 for all y + if (y == 0 or x == 1) { + return 1; + } + + // pow(nan, y) = nan for all y + // pow(x, nan) = nan for all x + if (math.isNan(x) or math.isNan(y)) { + return math.nan(f32); + } + + // pow(x, 1) = x for all x + if (y == 1) { + return x; + } + + // special case sqrt + if (y == 0.5) { + return math.sqrt(x); + } + + if (y == -0.5) { + return 1 / math.sqrt(x); + } + + if (x == 0) { + if (y < 0) { + // pow(+-0, y) = +- 0 for y an odd integer + if (isOddInteger(y)) { + return math.copysign(f32, math.inf(f32), x); + } + // pow(+-0, y) = +inf for y an even integer + else { + return math.inf(f32); + } + } else { + if (isOddInteger(y)) { + return x; + } else { + return 0; + } + } + } + + if (math.isInf(y)) { + // pow(-1, inf) = -1 for all x + if (x == -1) { + return -1; + } + // pow(x, +inf) = +0 for |x| < 1 + // pow(x, -inf) = +0 for |x| > 1 + else if ((math.fabs(x) < 1) == math.isPositiveInf(y)) { + return 0; + } + // pow(x, -inf) = +inf for |x| < 1 + // pow(x, +inf) = +inf for |x| > 1 + else { + return math.inf(f32); + } + } + + if (math.isInf(x)) { + if (math.isNegativeInf(x)) { + return pow32(1 / x, -y); + } + // pow(+inf, y) = +0 for y < 0 + else if (y < 0) { + return 0; + } + // pow(+inf, y) = +0 for y > 0 + else if (y > 0) { + return math.inf(f32); + } + } + + var ay = y; + var flip = false; + if (ay < 0) { + ay = -ay; + flip = true; + } + + const r1 = math.modf(ay); + var yi = r1.ipart; + var yf = r1.fpart; + + if (yf != 0 and x < 0) { + return math.nan(f32); + } + if (yi >= 1 << 31) { + return math.exp(y * math.ln(x)); + } + + // a = a1 * 2^ae + var a1: f32 = 1.0; + var ae: i32 = 0; + + // a *= x^yf + if (yf != 0) { + if (yf > 0.5) { + yf -= 1; + yi += 1; + } + a1 = math.exp(yf * math.ln(x)); + } + + // a *= x^yi + const r2 = math.frexp(x); + var xe = r2.exponent; + var x1 = r2.significand; + + var i = i32(yi); + while (i != 0) : (i >>= 1) { + if (i & 1 == 1) { + a1 *= x1; + ae += xe; + } + x1 *= x1; + xe <<= 1; + if (x1 < 0.5) { + x1 += x1; + xe -= 1; + } + } + + // a *= a1 * 2^ae + if (flip) { + a1 = 1 / a1; + ae = -ae; + } + + math.scalbn(a1, ae) +} + +// This implementation is taken from the go stlib, musl is a bit more complex. +fn pow64(x: f64, y: f64) -> f64 { + // pow(x, +-0) = 1 for all x + // pow(1, y) = 1 for all y + if (y == 0 or x == 1) { + return 1; + } + + // pow(nan, y) = nan for all y + // pow(x, nan) = nan for all x + if (math.isNan(x) or math.isNan(y)) { + return math.nan(f64); + } + + // pow(x, 1) = x for all x + if (y == 1) { + return x; + } + + // special case sqrt + if (y == 0.5) { + return math.sqrt(x); + } + + if (y == -0.5) { + return 1 / math.sqrt(x); + } + + if (x == 0) { + if (y < 0) { + // pow(+-0, y) = +- 0 for y an odd integer + if (isOddInteger(y)) { + return math.copysign(f64, math.inf(f64), x); + } + // pow(+-0, y) = +inf for y an even integer + else { + return math.inf(f64); + } + } else { + if (isOddInteger(y)) { + return x; + } else { + return 0; + } + } + } + + if (math.isInf(y)) { + // pow(-1, inf) = -1 for all x + if (x == -1) { + return -1; + } + // pow(x, +inf) = +0 for |x| < 1 + // pow(x, -inf) = +0 for |x| > 1 + else if ((math.fabs(x) < 1) == math.isInf(y)) { + return 0; + } + // pow(x, -inf) = +inf for |x| < 1 + // pow(x, +inf) = +inf for |x| > 1 + else { + return math.inf(f64); + } + } + + if (math.isInf(x)) { + if (math.isInf(x)) { + return pow64(1 / x, -y); + } + // pow(+inf, y) = +0 for y < 0 + else if (y < 0) { + return 0; + } + // pow(+inf, y) = +0 for y > 0 + else if (y > 0) { + return math.inf(f64); + } + } + + var ay = y; + var flip = false; + if (ay < 0) { + ay = -ay; + flip = true; + } + + const r1 = math.modf(ay); + var yi = r1.ipart; + var yf = r1.fpart; + + if (yf != 0 and x < 0) { + return math.nan(f64); + } + if (yi >= 1 << 63) { + return math.exp(y * math.ln(x)); + } + + // a = a1 * 2^ae + var a1: f64 = 1.0; + var ae: i32 = 0; + + // a *= x^yf + if (yf != 0) { + if (yf > 0.5) { + yf -= 1; + yi += 1; + } + a1 = math.exp(yf * math.ln(x)); + } + + // a *= x^yi + const r2 = math.frexp(x); + var xe = r2.exponent; + var x1 = r2.significand; + + var i = i64(yi); + while (i != 0) : (i >>= 1) { + if (i & 1 == 1) { + a1 *= x1; + ae += xe; + } + x1 *= x1; + xe <<= 1; + if (x1 < 0.5) { + x1 += x1; + xe -= 1; + } + } + + // a *= a1 * 2^ae + if (flip) { + a1 = 1 / a1; + ae = -ae; + } + + math.scalbn(a1, ae) +} + +test "pow" { + assert(pow(f32, 0.2, 3.3) == pow32(0.2, 3.3)); + assert(pow(f64, 0.2, 3.3) == pow64(0.2, 3.3)); +} + +test "pow32" { + const epsilon = 0.000001; + + // assert(math.approxEq(f32, pow32(0.0, 3.3), 0.0, epsilon)); // TODO: Handle div zero + assert(math.approxEq(f32, pow32(0.8923, 3.3), 0.686572, epsilon)); + assert(math.approxEq(f32, pow32(0.2, 3.3), 0.004936, epsilon)); + assert(math.approxEq(f32, pow32(1.5, 3.3), 3.811546, epsilon)); + assert(math.approxEq(f32, pow32(37.45, 3.3), 155736.703125, epsilon)); + assert(math.approxEq(f32, pow32(89.123, 3.3), 2722489.5, epsilon)); +} + +test "pow64" { + const epsilon = 0.000001; + + // assert(math.approxEq(f32, pow32(0.0, 3.3), 0.0, epsilon)); // TODO: Handle div zero + assert(math.approxEq(f64, pow64(0.8923, 3.3), 0.686572, epsilon)); + assert(math.approxEq(f64, pow64(0.2, 3.3), 0.004936, epsilon)); + assert(math.approxEq(f64, pow64(1.5, 3.3), 3.811546, epsilon)); + assert(math.approxEq(f64, pow64(37.45, 3.3), 155736.7160616, epsilon)); + assert(math.approxEq(f64, pow64(89.123, 3.3), 2722490.231436, epsilon)); +} diff --git a/std/math/round.zig b/std/math/round.zig new file mode 100644 index 0000000000..41768b7ba0 --- /dev/null +++ b/std/math/round.zig @@ -0,0 +1,105 @@ +const builtin = @import("builtin"); +const assert = @import("../debug.zig").assert; +const math = @import("index.zig"); + +pub fn round(x: var) -> @typeOf(x) { + const T = @typeOf(x); + switch (T) { + f32 => @inlineCall(round32, x), + f64 => @inlineCall(round64, x), + else => @compileError("round not implemented for " ++ @typeName(T)), + } +} + +fn round32(x_: f32) -> 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.forceEval(x + math.f32_toint); + return 0 * @bitCast(f32, u); + } + + { + @setFloatMode(this, builtin.FloatMode.Strict); + y = x + math.f32_toint - math.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) { + -y + } else { + y + } +} + +fn round64(x_: f64) -> 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.forceEval(x + math.f64_toint); + return 0 * @bitCast(f64, u); + } + + { + @setFloatMode(this, builtin.FloatMode.Strict); + y = x + math.f64_toint - math.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) { + -y + } else { + y + } +} + +test "round" { + assert(round(f32(1.3)) == round32(1.3)); + assert(round(f64(1.3)) == round64(1.3)); +} + +test "round32" { + assert(round32(1.3) == 1.0); + assert(round32(-1.3) == -1.0); + assert(round32(0.2) == 0.0); + assert(round32(1.8) == 2.0); +} + +test "round64" { + assert(round64(1.3) == 1.0); + assert(round64(-1.3) == -1.0); + assert(round64(0.2) == 0.0); + assert(round64(1.8) == 2.0); +} diff --git a/std/math/scalbn.zig b/std/math/scalbn.zig new file mode 100644 index 0000000000..b0b207e9d9 --- /dev/null +++ b/std/math/scalbn.zig @@ -0,0 +1,85 @@ +const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; + +pub fn scalbn(x: var, n: i32) -> @typeOf(x) { + const T = @typeOf(x); + switch (T) { + f32 => @inlineCall(scalbn32, x, n), + f64 => @inlineCall(scalbn64, x, n), + else => @compileError("scalbn not implemented for " ++ @typeName(T)), + } +} + +fn scalbn32(x: f32, n_: i32) -> f32 { + var y = x; + var n = n_; + + if (n > 127) { + y *= 0x1.0p127; + n -= 127; + if (n > 1023) { + y *= 0x1.0p127; + n -= 127; + if (n > 127) { + n = 127; + } + } + } else if (n < -126) { + y *= 0x1.0p-126 * 0x1.0p24; + n += 126 - 24; + if (n < -126) { + y *= 0x1.0p-126 * 0x1.0p24; + n += 126 - 24; + if (n < -126) { + n = -126; + } + } + } + + const u = u32(n +% 0x7F) << 23; + y * @bitCast(f32, u) +} + +fn scalbn64(x: f64, n_: i32) -> f64 { + var y = x; + var n = n_; + + if (n > 1023) { + // TODO: Determine how to do the following. + // y *= 0x1.0p1023; + n -= 1023; + if (n > 1023) { + // y *= 0x1.0p1023; + n -= 1023; + if (n > 1023) { + n = 1023; + } + } + } else if (n < -1022) { + y *= 0x1.0p-1022 * 0x1.0p53; + n += 1022 - 53; + if (n < -1022) { + y *= 0x1.0p-1022 * 0x1.0p53; + n += 1022 - 53; + if (n < -1022) { + n = -1022; + } + } + } + + const u = u64(n +% 0x3FF) << 52; + y * @bitCast(f64, u) +} + +test "scalbn" { + assert(scalbn(f32(1.5), 4) == scalbn32(1.5, 4)); + assert(scalbn(f64(1.5), 4) == scalbn64(1.5, 4)); +} + +test "scalbn32" { + assert(scalbn32(1.5, 4) == 24.0); +} + +test "scalbn64" { + assert(scalbn64(1.5, 4) == 24.0); +} diff --git a/std/math/signbit.zig b/std/math/signbit.zig new file mode 100644 index 0000000000..fcc7224a80 --- /dev/null +++ b/std/math/signbit.zig @@ -0,0 +1,36 @@ +const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; + +pub fn signbit(x: var) -> bool { + const T = @typeOf(x); + switch (T) { + f32 => @inlineCall(signbit32, x), + f64 => @inlineCall(signbit64, x), + else => @compileError("signbit not implemented for " ++ @typeName(T)), + } +} + +fn signbit32(x: f32) -> bool { + const bits = @bitCast(u32, x); + bits >> 31 != 0 +} + +fn signbit64(x: f64) -> bool { + const bits = @bitCast(u64, x); + bits >> 63 != 0 +} + +test "signbit" { + assert(signbit(f32(4.0)) == signbit32(4.0)); + assert(signbit(f64(4.0)) == signbit64(4.0)); +} + +test "signbit32" { + assert(!signbit32(4.0)); + assert(signbit32(-3.0)); +} + +test "signbit64" { + assert(!signbit64(4.0)); + assert(signbit64(-3.0)); +} diff --git a/std/math/sin.zig b/std/math/sin.zig new file mode 100644 index 0000000000..3c8a30c0b5 --- /dev/null +++ b/std/math/sin.zig @@ -0,0 +1,161 @@ +const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; + +pub fn sin(x: var) -> @typeOf(x) { + const T = @typeOf(x); + switch (T) { + f32 => @inlineCall(sin32, x), + f64 => @inlineCall(sin64, x), + else => @compileError("sin not implemented for " ++ @typeName(T)), + } +} + +// sin polynomial coefficients +const S0 = 1.58962301576546568060E-10; +const S1 = -2.50507477628578072866E-8; +const S2 = 2.75573136213857245213E-6; +const S3 = -1.98412698295895385996E-4; +const S4 = 8.33333333332211858878E-3; +const S5 = -1.66666666666666307295E-1; + +// cos polynomial coeffiecients +const C0 = -1.13585365213876817300E-11; +const C1 = 2.08757008419747316778E-9; +const C2 = -2.75573141792967388112E-7; +const C3 = 2.48015872888517045348E-5; +const C4 = -1.38888888888730564116E-3; +const C5 = 4.16666666666665929218E-2; + +// NOTE: This is taken from the go stdlib. The musl implementation is much more complex. +// +// This may have slight differences on some edge cases and may need to replaced if so. +fn sin32(x_: f32) -> f32 { + const pi4a = 7.85398125648498535156e-1; + const pi4b = 3.77489470793079817668E-8; + const pi4c = 2.69515142907905952645E-15; + const m4pi = 1.273239544735162542821171882678754627704620361328125; + + var x = x_; + if (x == 0 or math.isNan(x)) { + return x; + } + if (math.isInf(x)) { + return math.nan(f32); + } + + var sign = false; + if (x < 0) { + x = -x; + sign = true; + } + + var y = math.floor(x * m4pi); + var j = i64(y); + + if (j & 1 == 1) { + j += 1; + y += 1; + } + + j &= 7; + if (j > 3) { + j -= 4; + sign = !sign; + } + + const z = ((x - y * pi4a) - y * pi4b) - y * pi4c; + const w = z * z; + + const r = { + if (j == 1 or j == 2) { + 1.0 - 0.5 * w + w * w * (C5 + w * (C4 + w * (C3 + w * (C2 + w * (C1 + w * C0))))) + } else { + z + z * w * (S5 + w * (S4 + w * (S3 + w * (S2 + w * (S1 + w * S0))))) + } + }; + + if (sign) { + -r + } else { + r + } +} + +fn sin64(x_: f64) -> f64 { + const pi4a = 7.85398125648498535156e-1; + const pi4b = 3.77489470793079817668E-8; + const pi4c = 2.69515142907905952645E-15; + const m4pi = 1.273239544735162542821171882678754627704620361328125; + + var x = x_; + if (x == 0 or math.isNan(x)) { + return x; + } + if (math.isInf(x)) { + return math.nan(f64); + } + + var sign = false; + if (x < 0) { + x = -x; + sign = true; + } + + var y = math.floor(x * m4pi); + var j = i64(y); + + if (j & 1 == 1) { + j += 1; + y += 1; + } + + j &= 7; + if (j > 3) { + j -= 4; + sign = !sign; + } + + const z = ((x - y * pi4a) - y * pi4b) - y * pi4c; + const w = z * z; + + const r = { + if (j == 1 or j == 2) { + 1.0 - 0.5 * w + w * w * (C5 + w * (C4 + w * (C3 + w * (C2 + w * (C1 + w * C0))))) + } else { + z + z * w * (S5 + w * (S4 + w * (S3 + w * (S2 + w * (S1 + w * S0))))) + } + }; + + if (sign) { + -r + } else { + r + } +} + +test "sin" { + assert(sin(f32(0.0)) == sin32(0.0)); + assert(sin(f64(0.0)) == sin64(0.0)); +} + +test "sin32" { + const epsilon = 0.000001; + + assert(math.approxEq(f32, sin32(0.0), 0.0, epsilon)); + assert(math.approxEq(f32, sin32(0.2), 0.198669, epsilon)); + assert(math.approxEq(f32, sin32(0.8923), 0.778517, epsilon)); + assert(math.approxEq(f32, sin32(1.5), 0.997495, epsilon)); + assert(math.approxEq(f32, sin32(37.45), -0.246544, epsilon)); + assert(math.approxEq(f32, sin32(89.123), 0.916166, epsilon)); +} + +test "sin64" { + const epsilon = 0.000001; + + assert(math.approxEq(f64, sin64(0.0), 0.0, epsilon)); + assert(math.approxEq(f64, sin64(0.2), 0.198669, epsilon)); + assert(math.approxEq(f64, sin64(0.8923), 0.778517, epsilon)); + assert(math.approxEq(f64, sin64(1.5), 0.997495, epsilon)); + assert(math.approxEq(f64, sin64(37.45), -0.246543, epsilon)); + assert(math.approxEq(f64, sin64(89.123), 0.916166, epsilon)); +} diff --git a/std/math/sinh.zig b/std/math/sinh.zig new file mode 100644 index 0000000000..7535bc2b5a --- /dev/null +++ b/std/math/sinh.zig @@ -0,0 +1,93 @@ +const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; +const expo2 = @import("_expo2.zig").expo2; + +pub fn sinh(x: var) -> @typeOf(x) { + const T = @typeOf(x); + switch (T) { + f32 => @inlineCall(sinhf, x), + f64 => @inlineCall(sinhd, x), + else => @compileError("sinh not implemented for " ++ @typeName(T)), + } +} + +// sinh(x) = (exp(x) - 1 / exp(x)) / 2 +// = (exp(x) - 1 + (exp(x) - 1) / exp(x)) / 2 +// = x + x^3 / 6 + o(x^5) +fn sinhf(x: f32) -> f32 { + const u = @bitCast(u32, x); + const ux = u & 0x7FFFFFFF; + const ax = @bitCast(f32, ux); + + var h: f32 = 0.5; + if (u >> 31 != 0) { + h = -h; + } + + // |x| < log(FLT_MAX) + if (ux < 0x42B17217) { + const t = math.expm1(ax); + if (ux < 0x3F800000) { + if (ux < 0x3F800000 - (12 << 23)) { + return x; + } else { + return h * (2 * t - t * t / (t + 1)); + } + } + return h * (t + t / (t + 1)); + } + + // |x| > log(FLT_MAX) or nan + 2 * h * expo2(ax) +} + +fn sinhd(x: f64) -> f64 { + const u = @bitCast(u64, x); + const w = u32(u >> 32); + const ax = @bitCast(f64, u & (@maxValue(u64) >> 1)); + + var h: f32 = 0.5; + if (u >> 63 != 0) { + h = -h; + } + + // |x| < log(FLT_MAX) + if (w < 0x40862E42) { + const t = math.expm1(ax); + if (w < 0x3FF00000) { + if (w < 0x3FF00000 - (26 << 20)) { + return x; + } else { + return h * (2 * t - t * t / (t + 1)); + } + } + // NOTE: |x| > log(0x1p26) + eps could be h * exp(x) + return h * (t + t / (t + 1)); + } + + // |x| > log(DBL_MAX) or nan + 2 * h * expo2(ax) +} + +test "sinh" { + assert(sinh(f32(1.5)) == sinhf(1.5)); + assert(sinh(f64(1.5)) == sinhd(1.5)); +} + +test "sinhf" { + const epsilon = 0.000001; + + assert(math.approxEq(f32, sinhf(0.0), 0.0, epsilon)); + assert(math.approxEq(f32, sinhf(0.2), 0.201336, epsilon)); + assert(math.approxEq(f32, sinhf(0.8923), 1.015512, epsilon)); + assert(math.approxEq(f32, sinhf(1.5), 2.129279, epsilon)); +} + +test "sinhd" { + const epsilon = 0.000001; + + assert(math.approxEq(f64, sinhd(0.0), 0.0, epsilon)); + assert(math.approxEq(f64, sinhd(0.2), 0.201336, epsilon)); + assert(math.approxEq(f64, sinhd(0.8923), 1.015512, epsilon)); + assert(math.approxEq(f64, sinhd(1.5), 2.129279, epsilon)); +} diff --git a/std/math/sqrt.zig b/std/math/sqrt.zig new file mode 100644 index 0000000000..f47a5972c3 --- /dev/null +++ b/std/math/sqrt.zig @@ -0,0 +1,253 @@ +const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; + +pub fn sqrt(x: var) -> @typeOf(x) { + const T = @typeOf(x); + switch (T) { + f32 => @inlineCall(sqrt32, x), + f64 => @inlineCall(sqrt64, x), + else => @compileError("sqrt not implemented for " ++ @typeName(T)), + } +} + +fn sqrt32(x: f32) -> f32 { + const tiny: f32 = 1.0e-30; + const sign: i32 = @bitCast(i32, 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 (x - x) / (x - x); // sqrt(-ve) = snan + } + } + + // 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; + @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. +fn sqrt64(x: f64) -> f64 { + const tiny: f64 = 1.0e-300; + const sign: u32 = 0x80000000; + const u = @bitCast(u64, x); + + var ix0 = u32(u >> 32); + var ix1 = u32(u & 0xFFFFFFFF); + + // sqrt(nan) = nan, sqrt(+inf) = +inf, sqrt(-inf) = nan + if (ix0 & 0x7FF00000 == 0x7FF00000) { + return x * x + x; + } + + // sqrt(+-0) = +-0 + if ((ix0 & ~sign) | ix0 == 0) { + return x; + } + // sqrt(-ve) = snan + if (ix0 & sign != 0) { + return (x - x) / (x - x); + } + + // normalize x + var m = 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 -= i32(i) - 1; + ix0 |= ix1 >> (32 - i); + ix1 <<= 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) { + t = 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 = i32(ix0); + iix0 = iix0 +% (m << 20); + + const uz = (u64(iix0) << 32) | ix1; + @bitCast(f64, uz) +} + +test "sqrt" { + assert(sqrt(f32(0.0)) == sqrt32(0.0)); + assert(sqrt(f64(0.0)) == sqrt64(0.0)); +} + +test "sqrt32" { + const epsilon = 0.000001; + + assert(sqrt32(0.0) == 0.0); + assert(math.approxEq(f32, sqrt32(2.0), 1.414214, epsilon)); + assert(math.approxEq(f32, sqrt32(3.6), 1.897367, epsilon)); + assert(sqrt32(4.0) == 2.0); + assert(math.approxEq(f32, sqrt32(7.539840), 2.745877, epsilon)); + assert(math.approxEq(f32, sqrt32(19.230934), 4.385309, epsilon)); + assert(sqrt32(64.0) == 8.0); + assert(math.approxEq(f32, sqrt32(64.1), 8.006248, epsilon)); + assert(math.approxEq(f32, sqrt32(8942.230469), 94.563370, epsilon)); +} + +test "sqrt64" { + const epsilon = 0.000001; + + assert(sqrt64(0.0) == 0.0); + assert(math.approxEq(f64, sqrt64(2.0), 1.414214, epsilon)); + assert(math.approxEq(f64, sqrt64(3.6), 1.897367, epsilon)); + assert(sqrt64(4.0) == 2.0); + assert(math.approxEq(f64, sqrt64(7.539840), 2.745877, epsilon)); + assert(math.approxEq(f64, sqrt64(19.230934), 4.385309, epsilon)); + assert(sqrt64(64.0) == 8.0); + assert(math.approxEq(f64, sqrt64(64.1), 8.006248, epsilon)); + assert(math.approxEq(f64, sqrt64(8942.230469), 94.563367, epsilon)); +} diff --git a/std/math/tan.zig b/std/math/tan.zig new file mode 100644 index 0000000000..a8e907a59d --- /dev/null +++ b/std/math/tan.zig @@ -0,0 +1,148 @@ +const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; + +pub fn tan(x: var) -> @typeOf(x) { + const T = @typeOf(x); + switch (T) { + f32 => @inlineCall(tan32, x), + f64 => @inlineCall(tan64, x), + else => @compileError("tan not implemented for " ++ @typeName(T)), + } +} + +const Tp0 = -1.30936939181383777646E4; +const Tp1 = 1.15351664838587416140E6; +const Tp2 = -1.79565251976484877988E7; + +const Tq1 = 1.36812963470692954678E4; +const Tq2 = -1.32089234440210967447E6; +const Tq3 = 2.50083801823357915839E7; +const Tq4 = -5.38695755929454629881E7; + +// NOTE: This is taken from the go stdlib. The musl implementation is much more complex. +// +// This may have slight differences on some edge cases and may need to replaced if so. +fn tan32(x_: f32) -> f32 { + const pi4a = 7.85398125648498535156e-1; + const pi4b = 3.77489470793079817668E-8; + const pi4c = 2.69515142907905952645E-15; + const m4pi = 1.273239544735162542821171882678754627704620361328125; + + var x = x_; + if (x == 0 or math.isNan(x)) { + return x; + } + if (math.isInf(x)) { + return math.nan(f32); + } + + var sign = false; + if (x < 0) { + x = -x; + sign = true; + } + + var y = math.floor(x * m4pi); + var j = i64(y); + + if (j & 1 == 1) { + j += 1; + y += 1; + } + + const z = ((x - y * pi4a) - y * pi4b) - y * pi4c; + const w = z * z; + + var r = { + if (w > 1e-14) { + z + z * (w * ((Tp0 * w + Tp1) * w + Tp2) / ((((w + Tq1) * w + Tq2) * w + Tq3) * w + Tq4)) + } else { + z + } + }; + + if (j & 2 == 2) { + r = -1 / r; + } + if (sign) { + r = -r; + } + + r +} + +fn tan64(x_: f64) -> f64 { + const pi4a = 7.85398125648498535156e-1; + const pi4b = 3.77489470793079817668E-8; + const pi4c = 2.69515142907905952645E-15; + const m4pi = 1.273239544735162542821171882678754627704620361328125; + + var x = x_; + if (x == 0 or math.isNan(x)) { + return x; + } + if (math.isInf(x)) { + return math.nan(f64); + } + + var sign = false; + if (x < 0) { + x = -x; + sign = true; + } + + var y = math.floor(x * m4pi); + var j = i64(y); + + if (j & 1 == 1) { + j += 1; + y += 1; + } + + const z = ((x - y * pi4a) - y * pi4b) - y * pi4c; + const w = z * z; + + var r = { + if (w > 1e-14) { + z + z * (w * ((Tp0 * w + Tp1) * w + Tp2) / ((((w + Tq1) * w + Tq2) * w + Tq3) * w + Tq4)) + } else { + z + } + }; + + if (j & 2 == 2) { + r = -1 / r; + } + if (sign) { + r = -r; + } + + r +} + +test "tan" { + assert(tan(f32(0.0)) == tan32(0.0)); + assert(tan(f64(0.0)) == tan64(0.0)); +} + +test "tan32" { + const epsilon = 0.000001; + + assert(math.approxEq(f32, tan32(0.0), 0.0, epsilon)); + assert(math.approxEq(f32, tan32(0.2), 0.202710, epsilon)); + assert(math.approxEq(f32, tan32(0.8923), 1.240422, epsilon)); + assert(math.approxEq(f32, tan32(1.5), 14.101420, epsilon)); + assert(math.approxEq(f32, tan32(37.45), -0.254397, epsilon)); + assert(math.approxEq(f32, tan32(89.123), 2.285852, epsilon)); +} + +test "tan64" { + const epsilon = 0.000001; + + assert(math.approxEq(f64, tan64(0.0), 0.0, epsilon)); + assert(math.approxEq(f64, tan64(0.2), 0.202710, epsilon)); + assert(math.approxEq(f64, tan64(0.8923), 1.240422, epsilon)); + assert(math.approxEq(f64, tan64(1.5), 14.101420, epsilon)); + assert(math.approxEq(f64, tan64(37.45), -0.254397, epsilon)); + assert(math.approxEq(f64, tan64(89.123), 2.2858376, epsilon)); +} diff --git a/std/math/tanh.zig b/std/math/tanh.zig new file mode 100644 index 0000000000..6f64a7ef87 --- /dev/null +++ b/std/math/tanh.zig @@ -0,0 +1,120 @@ +const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; +const expo2 = @import("_expo2.zig").expo2; + +pub fn tanh(x: var) -> @typeOf(x) { + const T = @typeOf(x); + switch (T) { + f32 => @inlineCall(tanhf, x), + f64 => @inlineCall(tanhd, x), + else => @compileError("tanh not implemented for " ++ @typeName(T)), + } +} + +// tanh(x) = (exp(x) - exp(-x)) / (exp(x) + exp(-x)) +// = (exp(2x) - 1) / (exp(2x) - 1 + 2) +// = (1 - exp(-2x)) / (exp(-2x) - 1 + 2) +fn tanhf(x: f32) -> f32 { + const u = @bitCast(u32, x); + const ux = u & 0x7FFFFFFF; + const ax = @bitCast(f32, ux); + + var t: f32 = undefined; + + // |x| < log(3) / 2 ~= 0.5493 or nan + if (ux > 0x3F0C9F54) { + // |x| > 10 + if (ux > 0x41200000) { + t = 1.0 + 0 / x; + } else { + t = math.expm1(2 * x); + t = 1 - 2 / (t + 2); + } + } + // |x| > log(5 / 3) / 2 ~= 0.2554 + else if (ux > 0x3E82C578) { + t = math.expm1(2 * x); + t = t / (t + 2); + } + // |x| >= 0x1.0p-126 + else if (ux >= 0x00800000) { + t = math.expm1(-2 * x); + t = -t / (t + 2); + } + // |x| is subnormal + else { + math.forceEval(x * x); + t = x; + } + + if (u >> 31 != 0) { + -t + } else { + t + } +} + +fn tanhd(x: f64) -> f64 { + const u = @bitCast(u64, x); + const w = u32(u >> 32); + const ax = @bitCast(f64, u & (@maxValue(u64) >> 1)); + + var t: f64 = undefined; + + // |x| < log(3) / 2 ~= 0.5493 or nan + if (w > 0x3Fe193EA) { + // |x| > 20 or nan + if (w > 0x40340000) { + t = 1.0 + 0 / x; + } else { + t = math.expm1(2 * x); + t = 1 - 2 / (t + 2); + } + } + // |x| > log(5 / 3) / 2 ~= 0.2554 + else if (w > 0x3FD058AE) { + t = math.expm1(2 * x); + t = t / (t + 2); + } + // |x| >= 0x1.0p-1022 + else if (w >= 0x00100000) { + t = math.expm1(-2 * x); + t = -t / (t + 2); + } + // |x| is subnormal + else { + math.forceEval(f32(x)); + t = x; + } + + if (u >> 63 != 0) { + -t + } else { + t + } +} + +test "tanh" { + assert(tanh(f32(1.5)) == tanhf(1.5)); + assert(tanh(f64(1.5)) == tanhd(1.5)); +} + +test "tanhf" { + const epsilon = 0.000001; + + assert(math.approxEq(f32, tanhf(0.0), 0.0, epsilon)); + assert(math.approxEq(f32, tanhf(0.2), 0.197375, epsilon)); + assert(math.approxEq(f32, tanhf(0.8923), 0.712528, epsilon)); + assert(math.approxEq(f32, tanhf(1.5), 0.905148, epsilon)); + assert(math.approxEq(f32, tanhf(37.45), 1.0, epsilon)); +} + +test "tanhd" { + const epsilon = 0.000001; + + assert(math.approxEq(f64, tanhd(0.0), 0.0, epsilon)); + assert(math.approxEq(f64, tanhd(0.2), 0.197375, epsilon)); + assert(math.approxEq(f64, tanhd(0.8923), 0.712528, epsilon)); + assert(math.approxEq(f64, tanhd(1.5), 0.905148, epsilon)); + assert(math.approxEq(f64, tanhd(37.45), 1.0, epsilon)); +} diff --git a/std/math/trunc.zig b/std/math/trunc.zig new file mode 100644 index 0000000000..8c2c1b0496 --- /dev/null +++ b/std/math/trunc.zig @@ -0,0 +1,70 @@ +const math = @import("index.zig"); +const assert = @import("../debug.zig").assert; + +pub fn trunc(x: var) -> @typeOf(x) { + const T = @typeOf(x); + switch (T) { + f32 => @inlineCall(trunc32, x), + f64 => @inlineCall(trunc64, x), + else => @compileError("trunc not implemented for " ++ @typeName(T)), + } +} + +fn trunc32(x: f32) -> f32 { + const u = @bitCast(u32, x); + var e = i32(((u >> 23) & 0xFF)) - 0x7F + 9; + var m: u32 = undefined; + + if (e >= 23 + 9) { + return x; + } + if (e < 9) { + e = 1; + } + + m = @maxValue(u32) >> u32(e); + if (u & m == 0) { + x + } else { + math.forceEval(x + 0x1p120); + @bitCast(f32, u & ~m) + } +} + +fn trunc64(x: f64) -> f64 { + const u = @bitCast(u64, x); + var e = i32(((u >> 52) & 0x7FF)) - 0x3FF + 12; + var m: u64 = undefined; + + if (e >= 52 + 12) { + return x; + } + if (e < 12) { + e = 1; + } + + m = @maxValue(u64) >> u64(e); + if (u & m == 0) { + x + } else { + math.forceEval(x + 0x1p120); + @bitCast(f64, u & ~m) + } +} + +test "trunc" { + assert(trunc(f32(1.3)) == trunc32(1.3)); + assert(trunc(f64(1.3)) == trunc64(1.3)); +} + +test "trunc32" { + assert(trunc32(1.3) == 1.0); + assert(trunc32(-1.3) == -1.0); + assert(trunc32(0.2) == 0.0); +} + +test "trunc64" { + assert(trunc64(1.3) == 1.0); + assert(trunc64(-1.3) == -1.0); + assert(trunc64(0.2) == 0.0); +} |