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