diff options
Diffstat (limited to 'std/math')
28 files changed, 1461 insertions, 283 deletions
diff --git a/std/math/complex/abs.zig b/std/math/complex/abs.zig new file mode 100644 index 0000000000..4cd095c46b --- /dev/null +++ b/std/math/complex/abs.zig @@ -0,0 +1,18 @@ +const std = @import("../../index.zig"); +const debug = std.debug; +const math = std.math; +const cmath = math.complex; +const Complex = cmath.Complex; + +pub fn abs(z: var) @typeOf(z.re) { + const T = @typeOf(z.re); + return math.hypot(T, z.re, z.im); +} + +const epsilon = 0.0001; + +test "complex.cabs" { + const a = Complex(f32).new(5, 3); + const c = abs(a); + debug.assert(math.approxEq(f32, c, 5.83095, epsilon)); +} diff --git a/std/math/complex/acos.zig b/std/math/complex/acos.zig new file mode 100644 index 0000000000..a5760b4ace --- /dev/null +++ b/std/math/complex/acos.zig @@ -0,0 +1,21 @@ +const std = @import("../../index.zig"); +const debug = std.debug; +const math = std.math; +const cmath = math.complex; +const Complex = cmath.Complex; + +pub fn acos(z: var) Complex(@typeOf(z.re)) { + const T = @typeOf(z.re); + const q = cmath.asin(z); + return Complex(T).new(T(math.pi) / 2 - q.re, -q.im); +} + +const epsilon = 0.0001; + +test "complex.cacos" { + const a = Complex(f32).new(5, 3); + const c = acos(a); + + debug.assert(math.approxEq(f32, c.re, 0.546975, epsilon)); + debug.assert(math.approxEq(f32, c.im, -2.452914, epsilon)); +} diff --git a/std/math/complex/acosh.zig b/std/math/complex/acosh.zig new file mode 100644 index 0000000000..8dd91b2836 --- /dev/null +++ b/std/math/complex/acosh.zig @@ -0,0 +1,21 @@ +const std = @import("../../index.zig"); +const debug = std.debug; +const math = std.math; +const cmath = math.complex; +const Complex = cmath.Complex; + +pub fn acosh(z: var) Complex(@typeOf(z.re)) { + const T = @typeOf(z.re); + const q = cmath.acos(z); + return Complex(T).new(-q.im, q.re); +} + +const epsilon = 0.0001; + +test "complex.cacosh" { + const a = Complex(f32).new(5, 3); + const c = acosh(a); + + debug.assert(math.approxEq(f32, c.re, 2.452914, epsilon)); + debug.assert(math.approxEq(f32, c.im, 0.546975, epsilon)); +} diff --git a/std/math/complex/arg.zig b/std/math/complex/arg.zig new file mode 100644 index 0000000000..f24512ac73 --- /dev/null +++ b/std/math/complex/arg.zig @@ -0,0 +1,18 @@ +const std = @import("../../index.zig"); +const debug = std.debug; +const math = std.math; +const cmath = math.complex; +const Complex = cmath.Complex; + +pub fn arg(z: var) @typeOf(z.re) { + const T = @typeOf(z.re); + return math.atan2(T, z.im, z.re); +} + +const epsilon = 0.0001; + +test "complex.carg" { + const a = Complex(f32).new(5, 3); + const c = arg(a); + debug.assert(math.approxEq(f32, c, 0.540420, epsilon)); +} diff --git a/std/math/complex/asin.zig b/std/math/complex/asin.zig new file mode 100644 index 0000000000..584a3a1a9b --- /dev/null +++ b/std/math/complex/asin.zig @@ -0,0 +1,27 @@ +const std = @import("../../index.zig"); +const debug = std.debug; +const math = std.math; +const cmath = math.complex; +const Complex = cmath.Complex; + +pub fn asin(z: var) Complex(@typeOf(z.re)) { + const T = @typeOf(z.re); + const x = z.re; + const y = z.im; + + const p = Complex(T).new(1.0 - (x - y) * (x + y), -2.0 * x * y); + const q = Complex(T).new(-y, x); + const r = cmath.log(q.add(cmath.sqrt(p))); + + return Complex(T).new(r.im, -r.re); +} + +const epsilon = 0.0001; + +test "complex.casin" { + const a = Complex(f32).new(5, 3); + const c = asin(a); + + debug.assert(math.approxEq(f32, c.re, 1.023822, epsilon)); + debug.assert(math.approxEq(f32, c.im, 2.452914, epsilon)); +} diff --git a/std/math/complex/asinh.zig b/std/math/complex/asinh.zig new file mode 100644 index 0000000000..0c4dc2b6e4 --- /dev/null +++ b/std/math/complex/asinh.zig @@ -0,0 +1,22 @@ +const std = @import("../../index.zig"); +const debug = std.debug; +const math = std.math; +const cmath = math.complex; +const Complex = cmath.Complex; + +pub fn asinh(z: var) Complex(@typeOf(z.re)) { + const T = @typeOf(z.re); + const q = Complex(T).new(-z.im, z.re); + const r = cmath.asin(q); + return Complex(T).new(r.im, -r.re); +} + +const epsilon = 0.0001; + +test "complex.casinh" { + const a = Complex(f32).new(5, 3); + const c = asinh(a); + + debug.assert(math.approxEq(f32, c.re, 2.459831, epsilon)); + debug.assert(math.approxEq(f32, c.im, 0.533999, epsilon)); +} diff --git a/std/math/complex/atan.zig b/std/math/complex/atan.zig new file mode 100644 index 0000000000..b7bbf930eb --- /dev/null +++ b/std/math/complex/atan.zig @@ -0,0 +1,130 @@ +const std = @import("../../index.zig"); +const debug = std.debug; +const math = std.math; +const cmath = math.complex; +const Complex = cmath.Complex; + +pub fn atan(z: var) Complex(@typeOf(z.re)) { + const T = @typeOf(z.re); + return switch (T) { + f32 => atan32(z), + f64 => atan64(z), + else => @compileError("atan not implemented for " ++ @typeName(z)), + }; +} + +fn redupif32(x: f32) f32 { + const DP1 = 3.140625; + const DP2 = 9.67502593994140625e-4; + const DP3 = 1.509957990978376432e-7; + + var t = x / math.pi; + if (t >= 0.0) { + t += 0.5; + } else { + t -= 0.5; + } + + const u = f32(i32(t)); + return ((x - u * DP1) - u * DP2) - t * DP3; +} + +fn atan32(z: &const Complex(f32)) Complex(f32) { + const maxnum = 1.0e38; + + const x = z.re; + const y = z.im; + + if ((x == 0.0) and (y > 1.0)) { + // overflow + return Complex(f32).new(maxnum, maxnum); + } + + const x2 = x * x; + var a = 1.0 - x2 - (y * y); + if (a == 0.0) { + // overflow + return Complex(f32).new(maxnum, maxnum); + } + + var t = 0.5 * math.atan2(f32, 2.0 * x, a); + var w = redupif32(t); + + t = y - 1.0; + a = x2 + t * t; + if (a == 0.0) { + // overflow + return Complex(f32).new(maxnum, maxnum); + } + + t = y + 1.0; + a = (x2 + (t * t)) / a; + return Complex(f32).new(w, 0.25 * math.ln(a)); +} + +fn redupif64(x: f64) f64 { + const DP1 = 3.14159265160560607910; + const DP2 = 1.98418714791870343106e-9; + const DP3 = 1.14423774522196636802e-17; + + var t = x / math.pi; + if (t >= 0.0) { + t += 0.5; + } else { + t -= 0.5; + } + + const u = f64(i64(t)); + return ((x - u * DP1) - u * DP2) - t * DP3; +} + +fn atan64(z: &const Complex(f64)) Complex(f64) { + const maxnum = 1.0e308; + + const x = z.re; + const y = z.im; + + if ((x == 0.0) and (y > 1.0)) { + // overflow + return Complex(f64).new(maxnum, maxnum); + } + + const x2 = x * x; + var a = 1.0 - x2 - (y * y); + if (a == 0.0) { + // overflow + return Complex(f64).new(maxnum, maxnum); + } + + var t = 0.5 * math.atan2(f64, 2.0 * x, a); + var w = redupif64(t); + + t = y - 1.0; + a = x2 + t * t; + if (a == 0.0) { + // overflow + return Complex(f64).new(maxnum, maxnum); + } + + t = y + 1.0; + a = (x2 + (t * t)) / a; + return Complex(f64).new(w, 0.25 * math.ln(a)); +} + +const epsilon = 0.0001; + +test "complex.catan32" { + const a = Complex(f32).new(5, 3); + const c = atan(a); + + debug.assert(math.approxEq(f32, c.re, 1.423679, epsilon)); + debug.assert(math.approxEq(f32, c.im, 0.086569, epsilon)); +} + +test "complex.catan64" { + const a = Complex(f64).new(5, 3); + const c = atan(a); + + debug.assert(math.approxEq(f64, c.re, 1.423679, epsilon)); + debug.assert(math.approxEq(f64, c.im, 0.086569, epsilon)); +} diff --git a/std/math/complex/atanh.zig b/std/math/complex/atanh.zig new file mode 100644 index 0000000000..f70c741765 --- /dev/null +++ b/std/math/complex/atanh.zig @@ -0,0 +1,22 @@ +const std = @import("../../index.zig"); +const debug = std.debug; +const math = std.math; +const cmath = math.complex; +const Complex = cmath.Complex; + +pub fn atanh(z: var) Complex(@typeOf(z.re)) { + const T = @typeOf(z.re); + const q = Complex(T).new(-z.im, z.re); + const r = cmath.atan(q); + return Complex(T).new(r.im, -r.re); +} + +const epsilon = 0.0001; + +test "complex.catanh" { + const a = Complex(f32).new(5, 3); + const c = atanh(a); + + debug.assert(math.approxEq(f32, c.re, 0.146947, epsilon)); + debug.assert(math.approxEq(f32, c.im, 1.480870, epsilon)); +} diff --git a/std/math/complex/conj.zig b/std/math/complex/conj.zig new file mode 100644 index 0000000000..ad3e8b5036 --- /dev/null +++ b/std/math/complex/conj.zig @@ -0,0 +1,17 @@ +const std = @import("../../index.zig"); +const debug = std.debug; +const math = std.math; +const cmath = math.complex; +const Complex = cmath.Complex; + +pub fn conj(z: var) Complex(@typeOf(z.re)) { + const T = @typeOf(z.re); + return Complex(T).new(z.re, -z.im); +} + +test "complex.conj" { + const a = Complex(f32).new(5, 3); + const c = a.conjugate(); + + debug.assert(c.re == 5 and c.im == -3); +} diff --git a/std/math/complex/cos.zig b/std/math/complex/cos.zig new file mode 100644 index 0000000000..96e4ffcdb0 --- /dev/null +++ b/std/math/complex/cos.zig @@ -0,0 +1,21 @@ +const std = @import("../../index.zig"); +const debug = std.debug; +const math = std.math; +const cmath = math.complex; +const Complex = cmath.Complex; + +pub fn cos(z: var) Complex(@typeOf(z.re)) { + const T = @typeOf(z.re); + const p = Complex(T).new(-z.im, z.re); + return cmath.cosh(p); +} + +const epsilon = 0.0001; + +test "complex.ccos" { + const a = Complex(f32).new(5, 3); + const c = cos(a); + + debug.assert(math.approxEq(f32, c.re, 2.855815, epsilon)); + debug.assert(math.approxEq(f32, c.im, 9.606383, epsilon)); +} diff --git a/std/math/complex/cosh.zig b/std/math/complex/cosh.zig new file mode 100644 index 0000000000..96eac68556 --- /dev/null +++ b/std/math/complex/cosh.zig @@ -0,0 +1,165 @@ +const std = @import("../../index.zig"); +const debug = std.debug; +const math = std.math; +const cmath = math.complex; +const Complex = cmath.Complex; + +const ldexp_cexp = @import("ldexp.zig").ldexp_cexp; + +pub fn cosh(z: var) Complex(@typeOf(z.re)) { + const T = @typeOf(z.re); + return switch (T) { + f32 => cosh32(z), + f64 => cosh64(z), + else => @compileError("cosh not implemented for " ++ @typeName(z)), + }; +} + +fn cosh32(z: &const Complex(f32)) Complex(f32) { + const x = z.re; + const y = z.im; + + const hx = @bitCast(u32, x); + const ix = hx & 0x7fffffff; + + const hy = @bitCast(u32, y); + const iy = hy & 0x7fffffff; + + if (ix < 0x7f800000 and iy < 0x7f800000) { + if (iy == 0) { + return Complex(f32).new(math.cosh(x), y); + } + // small x: normal case + if (ix < 0x41100000) { + return Complex(f32).new(math.cosh(x) * math.cos(y), math.sinh(x) * math.sin(y)); + } + + // |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; + return Complex(f32).new(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).new(math.fabs(x), y); + const r = ldexp_cexp(v, -1); + return Complex(f32).new(x, y * math.copysign(f32, 1, x)); + } + // x >= 192.7: result always overflows + else { + const h = 0x1p127 * x; + return Complex(f32).new(h * h * math.cos(y), h * math.sin(y)); + } + } + + if (ix == 0 and iy >= 0x7f800000) { + return Complex(f32).new(y - y, math.copysign(f32, 0, x * (y - y))); + } + + if (iy == 0 and ix >= 0x7f800000) { + if (hx & 0x7fffff == 0) { + return Complex(f32).new(x * x, math.copysign(f32, 0, x) * y); + } + return Complex(f32).new(x, math.copysign(f32, 0, (x + x) * y)); + } + + if (ix < 0x7f800000 and iy >= 0x7f800000) { + return Complex(f32).new(y - y, x * (y - y)); + } + + if (ix >= 0x7f800000 and (hx & 0x7fffff) == 0) { + if (iy >= 0x7f800000) { + return Complex(f32).new(x * x, x * (y - y)); + } + return Complex(f32).new((x * x) * math.cos(y), x * math.sin(y)); + } + + return Complex(f32).new((x * x) * (y - y), (x + x) * (y - y)); +} + +fn cosh64(z: &const Complex(f64)) Complex(f64) { + const x = z.re; + const y = z.im; + + const fx = @bitCast(u64, x); + const hx = u32(fx >> 32); + const lx = @truncate(u32, fx); + const ix = hx & 0x7fffffff; + + const fy = @bitCast(u64, y); + const hy = u32(fy >> 32); + const ly = @truncate(u32, fy); + const iy = hy & 0x7fffffff; + + // nearly non-exceptional case where x, y are finite + if (ix < 0x7ff00000 and iy < 0x7ff00000) { + if (iy | ly == 0) { + return Complex(f64).new(math.cosh(x), x * y); + } + // small x: normal case + if (ix < 0x40360000) { + return Complex(f64).new(math.cosh(x) * math.cos(y), math.sinh(x) * math.sin(y)); + } + + // |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; + return Complex(f64).new(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).new(math.fabs(x), y); + const r = ldexp_cexp(v, -1); + return Complex(f64).new(x, y * math.copysign(f64, 1, x)); + } + // x >= 1455: result always overflows + else { + const h = 0x1p1023; + return Complex(f64).new(h * h * math.cos(y), h * math.sin(y)); + } + } + + if (ix | lx == 0 and iy >= 0x7ff00000) { + return Complex(f64).new(y - y, math.copysign(f64, 0, x * (y - y))); + } + + if (iy | ly == 0 and ix >= 0x7ff00000) { + if ((hx & 0xfffff) | lx == 0) { + return Complex(f64).new(x * x, math.copysign(f64, 0, x) * y); + } + return Complex(f64).new(x * x, math.copysign(f64, 0, (x + x) * y)); + } + + if (ix < 0x7ff00000 and iy >= 0x7ff00000) { + return Complex(f64).new(y - y, x * (y - y)); + } + + if (ix >= 0x7ff00000 and (hx & 0xfffff) | lx == 0) { + if (iy >= 0x7ff00000) { + return Complex(f64).new(x * x, x * (y - y)); + } + return Complex(f64).new(x * x * math.cos(y), x * math.sin(y)); + } + + return Complex(f64).new((x * x) * (y - y), (x + x) * (y - y)); +} + +const epsilon = 0.0001; + +test "complex.ccosh32" { + const a = Complex(f32).new(5, 3); + const c = cosh(a); + + debug.assert(math.approxEq(f32, c.re, -73.467300, epsilon)); + debug.assert(math.approxEq(f32, c.im, 10.471557, epsilon)); +} + +test "complex.ccosh64" { + const a = Complex(f64).new(5, 3); + const c = cosh(a); + + debug.assert(math.approxEq(f64, c.re, -73.467300, epsilon)); + debug.assert(math.approxEq(f64, c.im, 10.471557, epsilon)); +} diff --git a/std/math/complex/exp.zig b/std/math/complex/exp.zig new file mode 100644 index 0000000000..03f7f9e41b --- /dev/null +++ b/std/math/complex/exp.zig @@ -0,0 +1,140 @@ +const std = @import("../../index.zig"); +const debug = std.debug; +const math = std.math; +const cmath = math.complex; +const Complex = cmath.Complex; + +const ldexp_cexp = @import("ldexp.zig").ldexp_cexp; + +pub fn exp(z: var) Complex(@typeOf(z.re)) { + const T = @typeOf(z.re); + + return switch (T) { + f32 => exp32(z), + f64 => exp64(z), + else => @compileError("exp not implemented for " ++ @typeName(z)), + }; +} + +fn exp32(z: &const Complex(f32)) Complex(f32) { + @setFloatMode(this, @import("builtin").FloatMode.Strict); + + const exp_overflow = 0x42b17218; // max_exp * ln2 ~= 88.72283955 + const cexp_overflow = 0x43400074; // (max_exp - min_denom_exp) * ln2 + + const x = z.re; + const y = z.im; + + const hy = @bitCast(u32, y) & 0x7fffffff; + // cexp(x + i0) = exp(x) + i0 + if (hy == 0) { + return Complex(f32).new(math.exp(x), y); + } + + const hx = @bitCast(u32, x); + // cexp(0 + iy) = cos(y) + isin(y) + if ((hx & 0x7fffffff) == 0) { + return Complex(f32).new(math.cos(y), math.sin(y)); + } + + if (hy >= 0x7f800000) { + // cexp(finite|nan +- i inf|nan) = nan + i nan + if ((hx & 0x7fffffff) != 0x7f800000) { + return Complex(f32).new(y - y, y - y); + } + // cexp(-inf +- i inf|nan) = 0 + i0 + else if (hx & 0x80000000 != 0) { + return Complex(f32).new(0, 0); + } + // cexp(+inf +- i inf|nan) = inf + i nan + else { + return Complex(f32).new(x, y - y); + } + } + + // 88.7 <= x <= 192 so must scale + if (hx >= exp_overflow and hx <= cexp_overflow) { + return ldexp_cexp(z, 0); + } + // - x < exp_overflow => exp(x) won't overflow (common) + // - x > cexp_overflow, so exp(x) * s overflows for s > 0 + // - x = +-inf + // - x = nan + else { + const exp_x = math.exp(x); + return Complex(f32).new(exp_x * math.cos(y), exp_x * math.sin(y)); + } +} + +fn exp64(z: &const Complex(f64)) Complex(f64) { + const exp_overflow = 0x40862e42; // high bits of max_exp * ln2 ~= 710 + const cexp_overflow = 0x4096b8e4; // (max_exp - min_denorm_exp) * ln2 + + const x = z.re; + const y = z.im; + + const fy = @bitCast(u64, y); + const hy = u32(fy >> 32) & 0x7fffffff; + const ly = @truncate(u32, fy); + + // cexp(x + i0) = exp(x) + i0 + if (hy | ly == 0) { + return Complex(f64).new(math.exp(x), y); + } + + const fx = @bitCast(u64, x); + const hx = u32(fx >> 32); + const lx = @truncate(u32, fx); + + // cexp(0 + iy) = cos(y) + isin(y) + if ((hx & 0x7fffffff) | lx == 0) { + return Complex(f64).new(math.cos(y), math.sin(y)); + } + + if (hy >= 0x7ff00000) { + // cexp(finite|nan +- i inf|nan) = nan + i nan + if (lx != 0 or (hx & 0x7fffffff) != 0x7ff00000) { + return Complex(f64).new(y - y, y - y); + } + // cexp(-inf +- i inf|nan) = 0 + i0 + else if (hx & 0x80000000 != 0) { + return Complex(f64).new(0, 0); + } + // cexp(+inf +- i inf|nan) = inf + i nan + else { + return Complex(f64).new(x, y - y); + } + } + + // 709.7 <= x <= 1454.3 so must scale + if (hx >= exp_overflow and hx <= cexp_overflow) { + const r = ldexp_cexp(z, 0); + return *r; + } + // - x < exp_overflow => exp(x) won't overflow (common) + // - x > cexp_overflow, so exp(x) * s overflows for s > 0 + // - x = +-inf + // - x = nan + else { + const exp_x = math.exp(x); + return Complex(f64).new(exp_x * math.cos(y), exp_x * math.sin(y)); + } +} + +const epsilon = 0.0001; + +test "complex.cexp32" { + const a = Complex(f32).new(5, 3); + const c = exp(a); + + debug.assert(math.approxEq(f32, c.re, -146.927917, epsilon)); + debug.assert(math.approxEq(f32, c.im, 20.944065, epsilon)); +} + +test "complex.cexp64" { + const a = Complex(f32).new(5, 3); + const c = exp(a); + + debug.assert(math.approxEq(f64, c.re, -146.927917, epsilon)); + debug.assert(math.approxEq(f64, c.im, 20.944065, epsilon)); +} diff --git a/std/math/complex/index.zig b/std/math/complex/index.zig new file mode 100644 index 0000000000..a4d493307e --- /dev/null +++ b/std/math/complex/index.zig @@ -0,0 +1,171 @@ +const std = @import("../../index.zig"); +const debug = std.debug; +const math = std.math; + +pub const abs = @import("abs.zig").abs; +pub const acosh = @import("acosh.zig").acosh; +pub const acos = @import("acos.zig").acos; +pub const arg = @import("arg.zig").arg; +pub const asinh = @import("asinh.zig").asinh; +pub const asin = @import("asin.zig").asin; +pub const atanh = @import("atanh.zig").atanh; +pub const atan = @import("atan.zig").atan; +pub const conj = @import("conj.zig").conj; +pub const cosh = @import("cosh.zig").cosh; +pub const cos = @import("cos.zig").cos; +pub const exp = @import("exp.zig").exp; +pub const log = @import("log.zig").log; +pub const pow = @import("pow.zig").pow; +pub const proj = @import("proj.zig").proj; +pub const sinh = @import("sinh.zig").sinh; +pub const sin = @import("sin.zig").sin; +pub const sqrt = @import("sqrt.zig").sqrt; +pub const tanh = @import("tanh.zig").tanh; +pub const tan = @import("tan.zig").tan; + +pub fn Complex(comptime T: type) type { + return struct { + const Self = this; + + re: T, + im: T, + + pub fn new(re: T, im: T) Self { + return Self { + .re = re, + .im = im, + }; + } + + pub fn add(self: &const Self, other: &const Self) Self { + return Self { + .re = self.re + other.re, + .im = self.im + other.im, + }; + } + + pub fn sub(self: &const Self, other: &const Self) Self { + return Self { + .re = self.re - other.re, + .im = self.im - other.im, + }; + } + + pub fn mul(self: &const Self, other: &const Self) Self { + return Self { + .re = self.re * other.re - self.im * other.im, + .im = self.im * other.re + self.re * other.im, + }; + } + + pub fn div(self: &const Self, other: &const Self) Self { + const re_num = self.re * other.re + self.im * other.im; + const im_num = self.im * other.re - self.re * other.im; + const den = other.re * other.re + other.im * other.im; + + return Self { + .re = re_num / den, + .im = im_num / den, + }; + } + + pub fn conjugate(self: &const Self) Self { + return Self { + .re = self.re, + .im = -self.im, + }; + } + + pub fn reciprocal(self: &const Self) Self { + const m = self.re * self.re + self.im * self.im; + return Self { + .re = self.re / m, + .im = -self.im / m, + }; + } + + pub fn magnitude(self: &const Self) T { + return math.sqrt(self.re * self.re + self.im * self.im); + } + }; +} + +const epsilon = 0.0001; + +test "complex.add" { + const a = Complex(f32).new(5, 3); + const b = Complex(f32).new(2, 7); + const c = a.add(b); + + debug.assert(c.re == 7 and c.im == 10); +} + +test "complex.sub" { + const a = Complex(f32).new(5, 3); + const b = Complex(f32).new(2, 7); + const c = a.sub(b); + + debug.assert(c.re == 3 and c.im == -4); +} + +test "complex.mul" { + const a = Complex(f32).new(5, 3); + const b = Complex(f32).new(2, 7); + const c = a.mul(b); + + debug.assert(c.re == -11 and c.im == 41); +} + +test "complex.div" { + const a = Complex(f32).new(5, 3); + const b = Complex(f32).new(2, 7); + const c = a.div(b); + + debug.assert(math.approxEq(f32, c.re, f32(31)/53, epsilon) and + math.approxEq(f32, c.im, f32(-29)/53, epsilon)); +} + +test "complex.conjugate" { + const a = Complex(f32).new(5, 3); + const c = a.conjugate(); + + debug.assert(c.re == 5 and c.im == -3); +} + +test "complex.reciprocal" { + const a = Complex(f32).new(5, 3); + const c = a.reciprocal(); + + debug.assert(math.approxEq(f32, c.re, f32(5)/34, epsilon) and + math.approxEq(f32, c.im, f32(-3)/34, epsilon)); +} + +test "complex.magnitude" { + const a = Complex(f32).new(5, 3); + const c = a.magnitude(); + + debug.assert(math.approxEq(f32, c, 5.83095, epsilon)); +} + +test "complex.cmath" { + _ = @import("abs.zig"); + _ = @import("acosh.zig"); + _ = @import("acos.zig"); + _ = @import("arg.zig"); + _ = @import("asinh.zig"); + _ = @import("asin.zig"); + _ = @import("atanh.zig"); + _ = @import("atan.zig"); + _ = @import("conj.zig"); + _ = @import("cosh.zig"); + _ = @import("cos.zig"); + _ = @import("exp.zig"); + _ = @import("log.zig"); + _ = @import("pow.zig"); + _ = @import("proj.zig"); + _ = @import("sinh.zig"); + _ = @import("sin.zig"); + _ = @import("sqrt.zig"); + _ = @import("tanh.zig"); + _ = @import("tan.zig"); +} diff --git a/std/math/complex/ldexp.zig b/std/math/complex/ldexp.zig new file mode 100644 index 0000000000..4fb5a6815f --- /dev/null +++ b/std/math/complex/ldexp.zig @@ -0,0 +1,75 @@ +const std = @import("../../index.zig"); +const debug = std.debug; +const math = std.math; +const cmath = math.complex; +const Complex = cmath.Complex; + +pub fn ldexp_cexp(z: var, expt: i32) Complex(@typeOf(z.re)) { + const T = @typeOf(z.re); + + return switch (T) { + f32 => ldexp_cexp32(z, expt), + f64 => ldexp_cexp64(z, expt), + else => unreachable, + }; +} + +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 hx = @bitCast(u32, exp_x); + *expt = i32(hx >> 23) - (0x7f + 127) + k; + return @bitCast(f32, (hx & 0x7fffff) | ((0x7f + 127) << 23)); +} + +fn ldexp_cexp32(z: &const Complex(f32), expt: i32) Complex(f32) { + var ex_expt: i32 = undefined; + const exp_x = frexp_exp32(z.re, &ex_expt); + const exptf = expt + ex_expt; + + const half_expt1 = @divTrunc(exptf, 2); + const scale1 = @bitCast(f32, (0x7f + half_expt1) << 23); + + const half_expt2 = exptf - half_expt1; + const scale2 = @bitCast(f32, (0x7f + half_expt2) << 23); + + return Complex(f32).new( + math.cos(z.im) * exp_x * scale1 * scale2, + math.sin(z.im) * exp_x * scale1 * scale2, + ); +} + +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 fx = @bitCast(u64, x); + const hx = u32(fx >> 32); + const lx = @truncate(u32, fx); + + *expt = i32(hx >> 20) - (0x3ff + 1023) + k; + + const high_word = (hx & 0xfffff) | ((0x3ff + 1023) << 20); + return @bitCast(f64, (u64(high_word) << 32) | lx); +} + +fn ldexp_cexp64(z: &const Complex(f64), expt: i32) Complex(f64) { + var ex_expt: i32 = undefined; + const exp_x = frexp_exp64(z.re, &ex_expt); + const exptf = i64(expt + ex_expt); + + const half_expt1 = @divTrunc(exptf, 2); + const scale1 = @bitCast(f64, (0x3ff + half_expt1) << 20); + + const half_expt2 = exptf - half_expt1; + const scale2 = @bitCast(f64, (0x3ff + half_expt2) << 20); + + return Complex(f64).new( + math.cos(z.im) * exp_x * scale1 * scale2, + math.sin(z.im) * exp_x * scale1 * scale2, + ); +} diff --git a/std/math/complex/log.zig b/std/math/complex/log.zig new file mode 100644 index 0000000000..a4a1d1664f --- /dev/null +++ b/std/math/complex/log.zig @@ -0,0 +1,23 @@ +const std = @import("../../index.zig"); +const debug = std.debug; +const math = std.math; +const cmath = math.complex; +const Complex = cmath.Complex; + +pub fn log(z: var) Complex(@typeOf(z.re)) { + const T = @typeOf(z.re); + const r = cmath.abs(z); + const phi = cmath.arg(z); + + return Complex(T).new(math.ln(r), phi); +} + +const epsilon = 0.0001; + +test "complex.clog" { + const a = Complex(f32).new(5, 3); + const c = log(a); + + debug.assert(math.approxEq(f32, c.re, 1.763180, epsilon)); + debug.assert(math.approxEq(f32, c.im, 0.540419, epsilon)); +} diff --git a/std/math/complex/pow.zig b/std/math/complex/pow.zig new file mode 100644 index 0000000000..bef9fde542 --- /dev/null +++ b/std/math/complex/pow.zig @@ -0,0 +1,22 @@ +const std = @import("../../index.zig"); +const debug = std.debug; +const math = std.math; +const cmath = math.complex; +const Complex = cmath.Complex; + +pub fn pow(comptime T: type, z: &const T, c: &const T) T { + const p = cmath.log(z); + const q = c.mul(p); + return cmath.exp(q); +} + +const epsilon = 0.0001; + +test "complex.cpow" { + const a = Complex(f32).new(5, 3); + const b = Complex(f32).new(2.3, -1.3); + const c = pow(Complex(f32), a, b); + + debug.assert(math.approxEq(f32, c.re, 58.049110, epsilon)); + debug.assert(math.approxEq(f32, c.im, -101.003433, epsilon)); +} diff --git a/std/math/complex/proj.zig b/std/math/complex/proj.zig new file mode 100644 index 0000000000..b6c4cc046e --- /dev/null +++ b/std/math/complex/proj.zig @@ -0,0 +1,24 @@ +const std = @import("../../index.zig"); +const debug = std.debug; +const math = std.math; +const cmath = math.complex; +const Complex = cmath.Complex; + +pub fn proj(z: var) Complex(@typeOf(z.re)) { + const T = @typeOf(z.re); + + if (math.isInf(z.re) or math.isInf(z.im)) { + return Complex(T).new(math.inf(T), math.copysign(T, 0, z.re)); + } + + return Complex(T).new(z.re, z.im); +} + +const epsilon = 0.0001; + +test "complex.cproj" { + const a = Complex(f32).new(5, 3); + const c = proj(a); + + debug.assert(c.re == 5 and c.im == 3); +} diff --git a/std/math/complex/sin.zig b/std/math/complex/sin.zig new file mode 100644 index 0000000000..d32b771d3b --- /dev/null +++ b/std/math/complex/sin.zig @@ -0,0 +1,22 @@ +const std = @import("../../index.zig"); +const debug = std.debug; +const math = std.math; +const cmath = math.complex; +const Complex = cmath.Complex; + +pub fn sin(z: var) Complex(@typeOf(z.re)) { + const T = @typeOf(z.re); + const p = Complex(T).new(-z.im, z.re); + const q = cmath.sinh(p); + return Complex(T).new(q.im, -q.re); +} + +const epsilon = 0.0001; + +test "complex.csin" { + const a = Complex(f32).new(5, 3); + const c = sin(a); + + debug.assert(math.approxEq(f32, c.re, -9.654126, epsilon)); + debug.assert(math.approxEq(f32, c.im, 2.841692, epsilon)); +} diff --git a/std/math/complex/sinh.zig b/std/math/complex/sinh.zig new file mode 100644 index 0000000000..09a62ca058 --- /dev/null +++ b/std/math/complex/sinh.zig @@ -0,0 +1,164 @@ +const std = @import("../../index.zig"); +const debug = std.debug; +const math = std.math; +const cmath = math.complex; +const Complex = cmath.Complex; + +const ldexp_cexp = @import("ldexp.zig").ldexp_cexp; + +pub fn sinh(z: var) Complex(@typeOf(z.re)) { + const T = @typeOf(z.re); + return switch (T) { + f32 => sinh32(z), + f64 => sinh64(z), + else => @compileError("tan not implemented for " ++ @typeName(z)), + }; +} + +fn sinh32(z: &const Complex(f32)) Complex(f32) { + const x = z.re; + const y = z.im; + + const hx = @bitCast(u32, x); + const ix = hx & 0x7fffffff; + + const hy = @bitCast(u32, y); + const iy = hy & 0x7fffffff; + + if (ix < 0x7f800000 and iy < 0x7f800000) { + if (iy == 0) { + return Complex(f32).new(math.sinh(x), y); + } + // small x: normal case + if (ix < 0x41100000) { + return Complex(f32).new(math.sinh(x) * math.cos(y), math.cosh(x) * math.sin(y)); + } + + // |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; + return Complex(f32).new(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).new(math.fabs(x), y); + const r = ldexp_cexp(v, -1); + return Complex(f32).new(x * math.copysign(f32, 1, x), y); + } + // x >= 192.7: result always overflows + else { + const h = 0x1p127 * x; + return Complex(f32).new(h * math.cos(y), h * h * math.sin(y)); + } + } + + if (ix == 0 and iy >= 0x7f800000) { + return Complex(f32).new(math.copysign(f32, 0, x * (y - y)), y - y); + } + + if (iy == 0 and ix >= 0x7f800000) { + if (hx & 0x7fffff == 0) { + return Complex(f32).new(x, y); + } + return Complex(f32).new(x, math.copysign(f32, 0, y)); + } + + if (ix < 0x7f800000 and iy >= 0x7f800000) { + return Complex(f32).new(y - y, x * (y - y)); + } + + if (ix >= 0x7f800000 and (hx & 0x7fffff) == 0) { + if (iy >= 0x7f800000) { + return Complex(f32).new(x * x, x * (y - y)); + } + return Complex(f32).new(x * math.cos(y), math.inf_f32 * math.sin(y)); + } + + return Complex(f32).new((x * x) * (y - y), (x + x) * (y - y)); +} + +fn sinh64(z: &const Complex(f64)) Complex(f64) { + const x = z.re; + const y = z.im; + + const fx = @bitCast(u64, x); + const hx = u32(fx >> 32); + const lx = @truncate(u32, fx); + const ix = hx & 0x7fffffff; + + const fy = @bitCast(u64, y); + const hy = u32(fy >> 32); + const ly = @truncate(u32, fy); + const iy = hy & 0x7fffffff; + + if (ix < 0x7ff00000 and iy < 0x7ff00000) { + if (iy | ly == 0) { + return Complex(f64).new(math.sinh(x), y); + } + // small x: normal case + if (ix < 0x40360000) { + return Complex(f64).new(math.sinh(x) * math.cos(y), math.cosh(x) * math.sin(y)); + } + + // |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; + return Complex(f64).new(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).new(math.fabs(x), y); + const r = ldexp_cexp(v, -1); + return Complex(f64).new(x * math.copysign(f64, 1, x), y); + } + // x >= 1455: result always overflows + else { + const h = 0x1p1023 * x; + return Complex(f64).new(h * math.cos(y), h * h * math.sin(y)); + } + } + + if (ix | lx == 0 and iy >= 0x7ff00000) { + return Complex(f64).new(math.copysign(f64, 0, x * (y - y)), y - y); + } + + if (iy | ly == 0 and ix >= 0x7ff00000) { + if ((hx & 0xfffff) | lx == 0) { + return Complex(f64).new(x, y); + } + return Complex(f64).new(x, math.copysign(f64, 0, y)); + } + + if (ix < 0x7ff00000 and iy >= 0x7ff00000) { + return Complex(f64).new(y - y, x * (y - y)); + } + + if (ix >= 0x7ff00000 and (hx & 0xfffff) | lx == 0) { + if (iy >= 0x7ff00000) { + return Complex(f64).new(x * x, x * (y - y)); + } + return Complex(f64).new(x * math.cos(y), math.inf_f64 * math.sin(y)); + } + + return Complex(f64).new((x * x) * (y - y), (x + x) * (y - y)); +} + +const epsilon = 0.0001; + +test "complex.csinh32" { + const a = Complex(f32).new(5, 3); + const c = sinh(a); + + debug.assert(math.approxEq(f32, c.re, -73.460617, epsilon)); + debug.assert(math.approxEq(f32, c.im, 10.472508, epsilon)); +} + +test "complex.csinh64" { + const a = Complex(f64).new(5, 3); + const c = sinh(a); + + debug.assert(math.approxEq(f64, c.re, -73.460617, epsilon)); + debug.assert(math.approxEq(f64, c.im, 10.472508, epsilon)); +} diff --git a/std/math/complex/sqrt.zig b/std/math/complex/sqrt.zig new file mode 100644 index 0000000000..afda69f7c9 --- /dev/null +++ b/std/math/complex/sqrt.zig @@ -0,0 +1,133 @@ +const std = @import("../../index.zig"); +const debug = std.debug; +const math = std.math; +const cmath = math.complex; +const Complex = cmath.Complex; + +// TODO when #733 is solved this can be @typeOf(z) instead of Complex(@typeOf(z.re)) +pub fn sqrt(z: var) Complex(@typeOf(z.re)) { + const T = @typeOf(z.re); + + return switch (T) { + f32 => sqrt32(z), + f64 => sqrt64(z), + else => @compileError("sqrt not implemented for " ++ @typeName(z)), + }; +} + +fn sqrt32(z: &const Complex(f32)) Complex(f32) { + const x = z.re; + const y = z.im; + + if (x == 0 and y == 0) { + return Complex(f32).new(0, y); + } + if (math.isInf(y)) { + return Complex(f32).new(math.inf(f32), y); + } + if (math.isNan(x)) { + // raise invalid if y is not nan + const t = (y - y) / (y - y); + return Complex(f32).new(x, t); + } + if (math.isInf(x)) { + // sqrt(inf + i nan) = inf + nan i + // sqrt(inf + iy) = inf + i0 + // sqrt(-inf + i nan) = nan +- inf i + // sqrt(-inf + iy) = 0 + inf i + if (math.signbit(x)) { + return Complex(f32).new(math.fabs(x - y), math.copysign(f32, x, y)); + } else { + return Complex(f32).new(x, math.copysign(f32, y - y, y)); + } + } + + // y = nan special case is handled fine below + + // double-precision avoids overflow with correct rounding. + const dx = f64(x); + const dy = f64(y); + + if (dx >= 0) { + const t = math.sqrt((dx + math.hypot(f64, dx, dy)) * 0.5); + return Complex(f32).new(f32(t), f32(dy / (2.0 * t))); + } else { + const t = math.sqrt((-dx + math.hypot(f64, dx, dy)) * 0.5); + return Complex(f32).new(f32(math.fabs(y) / (2.0 * t)), f32(math.copysign(f64, t, y))); + } +} + +fn sqrt64(z: &const Complex(f64)) Complex(f64) { + // may encounter overflow for im,re >= DBL_MAX / (1 + sqrt(2)) + const threshold = 0x1.a827999fcef32p+1022; + + var x = z.re; + var y = z.im; + + if (x == 0 and y == 0) { + return Complex(f64).new(0, y); + } + if (math.isInf(y)) { + return Complex(f64).new(math.inf(f64), y); + } + if (math.isNan(x)) { + // raise invalid if y is not nan + const t = (y - y) / (y - y); + return Complex(f64).new(x, t); + } + if (math.isInf(x)) { + // sqrt(inf + i nan) = inf + nan i + // sqrt(inf + iy) = inf + i0 + // sqrt(-inf + i nan) = nan +- inf i + // sqrt(-inf + iy) = 0 + inf i + if (math.signbit(x)) { + return Complex(f64).new(math.fabs(x - y), math.copysign(f64, x, y)); + } else { + return Complex(f64).new(x, math.copysign(f64, y - y, y)); + } + } + + // y = nan special case is handled fine below + + // scale to avoid overflow + var scale = false; + if (math.fabs(x) >= threshold or math.fabs(y) >= threshold) { + x *= 0.25; + y *= 0.25; + scale = true; + } + + var result: Complex(f64) = undefined; + if (x >= 0) { + const t = math.sqrt((x + math.hypot(f64, x, y)) * 0.5); + result = Complex(f64).new(t, y / (2.0 * t)); + } else { + const t = math.sqrt((-x + math.hypot(f64, x, y)) * 0.5); + result = Complex(f64).new(math.fabs(y) / (2.0 * t), math.copysign(f64, t, y)); + } + + if (scale) { + result.re *= 2; + result.im *= 2; + } + + return result; +} + +const epsilon = 0.0001; + +test "complex.csqrt32" { + const a = Complex(f32).new(5, 3); + const c = sqrt(a); + + debug.assert(math.approxEq(f32, c.re, 2.327117, epsilon)); + debug.assert(math.approxEq(f32, c.im, 0.644574, epsilon)); +} + +test "complex.csqrt64" { + const a = Complex(f64).new(5, 3); + const c = sqrt(a); + + debug.assert(math.approxEq(f64, c.re, 2.3271175190399496, epsilon)); + debug.assert(math.approxEq(f64, c.im, 0.6445742373246469, epsilon)); +} diff --git a/std/math/complex/tan.zig b/std/math/complex/tan.zig new file mode 100644 index 0000000000..4ea5182fa7 --- /dev/null +++ b/std/math/complex/tan.zig @@ -0,0 +1,22 @@ +const std = @import("../../index.zig"); +const debug = std.debug; +const math = std.math; +const cmath = math.complex; +const Complex = cmath.Complex; + +pub fn tan(z: var) Complex(@typeOf(z.re)) { + const T = @typeOf(z.re); + const q = Complex(T).new(-z.im, z.re); + const r = cmath.tanh(q); + return Complex(T).new(r.im, -r.re); +} + +const epsilon = 0.0001; + +test "complex.ctan" { + const a = Complex(f32).new(5, 3); + const c = tan(a); + + debug.assert(math.approxEq(f32, c.re, -0.002708233, epsilon)); + debug.assert(math.approxEq(f32, c.im, 1.004165, epsilon)); +} diff --git a/std/math/complex/tanh.zig b/std/math/complex/tanh.zig new file mode 100644 index 0000000000..6af62f48ae --- /dev/null +++ b/std/math/complex/tanh.zig @@ -0,0 +1,111 @@ +const std = @import("../../index.zig"); +const debug = std.debug; +const math = std.math; +const cmath = math.complex; +const Complex = cmath.Complex; + +pub fn tanh(z: var) Complex(@typeOf(z.re)) { + const T = @typeOf(z.re); + return switch (T) { + f32 => tanh32(z), + f64 => tanh64(z), + else => @compileError("tan not implemented for " ++ @typeName(z)), + }; +} + +fn tanh32(z: &const Complex(f32)) Complex(f32) { + const x = z.re; + const y = z.im; + + const hx = @bitCast(u32, x); + const ix = hx & 0x7fffffff; + + if (ix >= 0x7f800000) { + if (ix & 0x7fffff != 0) { + const r = if (y == 0) y else x * y; + return Complex(f32).new(x, r); + } + const xx = @bitCast(f32, hx - 0x40000000); + const r = if (math.isInf(y)) y else math.sin(y) * math.cos(y); + return Complex(f32).new(xx, math.copysign(f32, 0, r)); + } + + if (!math.isFinite(y)) { + const r = if (ix != 0) y - y else x; + return Complex(f32).new(r, y - y); + } + + // x >= 11 + if (ix >= 0x41300000) { + const exp_mx = math.exp(-math.fabs(x)); + return Complex(f32).new(math.copysign(f32, 1, x), 4 * math.sin(y) * math.cos(y) * exp_mx * exp_mx); + } + + // Kahan's algorithm + const t = math.tan(y); + const beta = 1.0 + t * t; + const s = math.sinh(x); + const rho = math.sqrt(1 + s * s); + const den = 1 + beta * s * s; + + return Complex(f32).new((beta * rho * s) / den, t / den); +} + +fn tanh64(z: &const Complex(f64)) Complex(f64) { + const x = z.re; + const y = z.im; + + const fx = @bitCast(u64, x); + const hx = u32(fx >> 32); + const lx = @truncate(u32, fx); + const ix = hx & 0x7fffffff; + + if (ix >= 0x7ff00000) { + if ((ix & 0x7fffff) | lx != 0) { + const r = if (y == 0) y else x * y; + return Complex(f64).new(x, r); + } + + const xx = @bitCast(f64, (u64(hx - 0x40000000) << 32) | lx); + const r = if (math.isInf(y)) y else math.sin(y) * math.cos(y); + return Complex(f64).new(xx, math.copysign(f64, 0, r)); + } + + if (!math.isFinite(y)) { + const r = if (ix != 0) y - y else x; + return Complex(f64).new(r, y - y); + } + + // x >= 22 + if (ix >= 0x40360000) { + const exp_mx = math.exp(-math.fabs(x)); + return Complex(f64).new(math.copysign(f64, 1, x), 4 * math.sin(y) * math.cos(y) * exp_mx * exp_mx); + } + + // Kahan's algorithm + const t = math.tan(y); + const beta = 1.0 + t * t; + const s = math.sinh(x); + const rho = math.sqrt(1 + s * s); + const den = 1 + beta * s * s; + + return Complex(f64).new((beta * rho * s) / den, t / den); +} + +const epsilon = 0.0001; + +test "complex.ctanh32" { + const a = Complex(f32).new(5, 3); + const c = tanh(a); + + debug.assert(math.approxEq(f32, c.re, 0.999913, epsilon)); + debug.assert(math.approxEq(f32, c.im, -0.000025, epsilon)); +} + +test "complex.ctanh64" { + const a = Complex(f64).new(5, 3); + const c = tanh(a); + + debug.assert(math.approxEq(f64, c.re, 0.999913, epsilon)); + debug.assert(math.approxEq(f64, c.im, -0.000025, epsilon)); +} diff --git a/std/math/exp.zig b/std/math/exp.zig index 4032930a43..21aa558c57 100644 --- a/std/math/exp.zig +++ b/std/math/exp.zig @@ -6,6 +6,7 @@ const std = @import("../index.zig"); const math = std.math; const assert = std.debug.assert; +const builtin = @import("builtin"); pub fn exp(x: var) @typeOf(x) { const T = @typeOf(x); @@ -17,6 +18,8 @@ pub fn exp(x: var) @typeOf(x) { } fn exp32(x_: f32) f32 { + @setFloatMode(this, builtin.FloatMode.Strict); + const half = []f32 { 0.5, -0.5 }; const ln2hi = 6.9314575195e-1; const ln2lo = 1.4286067653e-6; @@ -94,6 +97,8 @@ fn exp32(x_: f32) f32 { } fn exp64(x_: f64) f64 { + @setFloatMode(this, builtin.FloatMode.Strict); + const half = []const f64 { 0.5, -0.5 }; const ln2hi: f64 = 6.93147180369123816490e-01; const ln2lo: f64 = 1.90821492927058770002e-10; diff --git a/std/math/index.zig b/std/math/index.zig index 477dafcbcc..a549a6bb61 100644 --- a/std/math/index.zig +++ b/std/math/index.zig @@ -129,6 +129,9 @@ pub const cos = @import("cos.zig").cos; pub const sin = @import("sin.zig").sin; pub const tan = @import("tan.zig").tan; +pub const complex = @import("complex/index.zig"); +pub const Complex = complex.Complex; + test "math" { _ = @import("nan.zig"); _ = @import("isnan.zig"); @@ -172,6 +175,8 @@ test "math" { _ = @import("sin.zig"); _ = @import("cos.zig"); _ = @import("tan.zig"); + + _ = @import("complex/index.zig"); } @@ -553,6 +558,32 @@ test "math.floorPowerOfTwo" { comptime testFloorPowerOfTwo(); } +pub fn log2_int(comptime T: type, x: T) Log2Int(T) { + assert(x != 0); + return Log2Int(T)(T.bit_count - 1 - @clz(x)); +} + +pub fn log2_int_ceil(comptime T: type, x: T) Log2Int(T) { + assert(x != 0); + const log2_val = log2_int(T, x); + if (T(1) << log2_val == x) + return log2_val; + return log2_val + 1; +} + +test "std.math.log2_int_ceil" { + assert(log2_int_ceil(u32, 1) == 0); + assert(log2_int_ceil(u32, 2) == 1); + assert(log2_int_ceil(u32, 3) == 2); + assert(log2_int_ceil(u32, 4) == 2); + assert(log2_int_ceil(u32, 5) == 3); + assert(log2_int_ceil(u32, 6) == 3); + assert(log2_int_ceil(u32, 7) == 3); + assert(log2_int_ceil(u32, 8) == 3); + assert(log2_int_ceil(u32, 9) == 4); + assert(log2_int_ceil(u32, 10) == 4); +} + fn testFloorPowerOfTwo() void { assert(floorPowerOfTwo(u32, 63) == 32); assert(floorPowerOfTwo(u32, 64) == 64); diff --git a/std/math/ln.zig b/std/math/ln.zig index c349ed7c6f..d09494b998 100644 --- a/std/math/ln.zig +++ b/std/math/ln.zig @@ -89,6 +89,8 @@ pub fn ln_32(x_: f32) f32 { } pub fn ln_64(x_: f64) f64 { + @setFloatMode(this, @import("builtin").FloatMode.Strict); + const ln2_hi: f64 = 6.93147180369123816490e-01; const ln2_lo: f64 = 1.90821492927058770002e-10; const Lg1: f64 = 6.666666666666735130e-01; diff --git a/std/math/log2.zig b/std/math/log2.zig index 998d6d6c5e..d5bbe385c2 100644 --- a/std/math/log2.zig +++ b/std/math/log2.zig @@ -31,17 +31,12 @@ pub fn log2(x: var) @typeOf(x) { return result; }, TypeId.Int => { - return log2_int(T, x); + return math.log2_int(T, x); }, else => @compileError("log2 not implemented for " ++ @typeName(T)), } } -pub fn log2_int(comptime T: type, x: T) T { - assert(x != 0); - return T.bit_count - 1 - T(@clz(x)); -} - pub fn log2_32(x_: f32) f32 { const ivln2hi: f32 = 1.4428710938e+00; const ivln2lo: f32 = -1.7605285393e-04; diff --git a/std/math/sqrt.zig b/std/math/sqrt.zig index 690f8b6901..982bd28b72 100644 --- a/std/math/sqrt.zig +++ b/std/math/sqrt.zig @@ -14,26 +14,8 @@ const TypeId = builtin.TypeId; pub fn sqrt(x: var) (if (@typeId(@typeOf(x)) == TypeId.Int) @IntType(false, @typeOf(x).bit_count / 2) else @typeOf(x)) { const T = @typeOf(x); switch (@typeId(T)) { - TypeId.FloatLiteral => { - return T(sqrt64(x)); - }, - TypeId.Float => { - switch (T) { - f32 => { - switch (builtin.arch) { - builtin.Arch.x86_64 => return @import("x86_64/sqrt.zig").sqrt32(x), - else => return sqrt32(x), - } - }, - f64 => { - switch (builtin.arch) { - builtin.Arch.x86_64 => return @import("x86_64/sqrt.zig").sqrt64(x), - else => return sqrt64(x), - } - }, - else => @compileError("sqrt not implemented for " ++ @typeName(T)), - } - }, + TypeId.FloatLiteral => return T(@sqrt(f64, x)), // TODO upgrade to f128 + TypeId.Float => return @sqrt(T, x), TypeId.IntLiteral => comptime { if (x > @maxValue(u128)) { @compileError("sqrt not implemented for comptime_int greater than 128 bits"); @@ -43,269 +25,58 @@ pub fn sqrt(x: var) (if (@typeId(@typeOf(x)) == TypeId.Int) @IntType(false, @typ } return T(sqrt_int(u128, x)); }, - TypeId.Int => { - return sqrt_int(T, x); - }, + TypeId.Int => return sqrt_int(T, x), else => @compileError("sqrt not implemented for " ++ @typeName(T)), } } -fn sqrt32(x: f32) f32 { - const tiny: f32 = 1.0e-30; - const sign: i32 = @bitCast(i32, u32(0x80000000)); - var ix: i32 = @bitCast(i32, x); - - if ((ix & 0x7F800000) == 0x7F800000) { - return x * x + x; // sqrt(nan) = nan, sqrt(+inf) = +inf, sqrt(-inf) = snan - } - - // zero - if (ix <= 0) { - if (ix & ~sign == 0) { - return x; // sqrt (+-0) = +-0 - } - if (ix < 0) { - return 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. -fn sqrt64(x: f64) f64 { - const tiny: f64 = 1.0e-300; - const sign: u32 = 0x80000000; - const u = @bitCast(u64, x); - - var ix0 = u32(u >> 32); - var ix1 = u32(u & 0xFFFFFFFF); - - // sqrt(nan) = nan, sqrt(+inf) = +inf, sqrt(-inf) = nan - if (ix0 & 0x7FF00000 == 0x7FF00000) { - return x * x + x; - } - - // sqrt(+-0) = +-0 - if (x == 0.0) { - return x; - } - // sqrt(-ve) = snan - if (ix0 & sign != 0) { - return math.snan(f64); - } - - // normalize x - var m = i32(ix0 >> 20); - if (m == 0) { - // subnormal - while (ix0 == 0) { - m -= 21; - ix0 |= ix1 >> 11; - ix1 <<= 21; - } - - // subnormal - var i: u32 = 0; - while (ix0 & 0x00100000 == 0) : (i += 1) { - ix0 <<= 1; - } - m -= i32(i) - 1; - ix0 |= ix1 >> u5(32 - i); - ix1 <<= 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) { - t = s1 +% r; - t = s0; - if (t < ix0 or (t == ix0 and t1 <= ix1)) { - s1 = t1 +% r; - if (t1 & sign == sign and s1 & sign == 0) { - s0 += 1; - } - ix0 -= t; - if (ix1 < t1) { - ix0 -= 1; - } - ix1 = ix1 -% t1; - q1 += r; - } - ix0 = ix0 +% ix0 +% (ix1 >> 31); - ix1 = ix1 +% ix1; - r >>= 1; - } - - // rounding direction - if (ix0 | ix1 != 0) { - var z = 1.0 - tiny; // raise inexact - if (z >= 1.0) { - z = 1.0 + tiny; - if (q1 == 0xFFFFFFFF) { - q1 = 0; - q += 1; - } else if (z > 1.0) { - if (q1 == 0xFFFFFFFE) { - q += 1; - } - q1 += 2; - } else { - q1 += q1 & 1; - } - } - } - - ix0 = (q >> 1) + 0x3FE00000; - ix1 = q1 >> 1; - if (q & 1 != 0) { - ix1 |= 0x80000000; - } - - // NOTE: musl here appears to rely on signed twos-complement wraparound. +% has the same - // behaviour at least. - var iix0 = i32(ix0); - iix0 = iix0 +% (m << 20); - - const uz = (u64(iix0) << 32) | ix1; - return @bitCast(f64, uz); -} - test "math.sqrt" { - assert(sqrt(f32(0.0)) == sqrt32(0.0)); - assert(sqrt(f64(0.0)) == sqrt64(0.0)); + assert(sqrt(f32(0.0)) == @sqrt(f32, 0.0)); + assert(sqrt(f64(0.0)) == @sqrt(f64, 0.0)); } test "math.sqrt32" { const epsilon = 0.000001; - assert(sqrt32(0.0) == 0.0); - assert(math.approxEq(f32, sqrt32(2.0), 1.414214, epsilon)); - assert(math.approxEq(f32, sqrt32(3.6), 1.897367, epsilon)); - assert(sqrt32(4.0) == 2.0); - assert(math.approxEq(f32, sqrt32(7.539840), 2.745877, epsilon)); - assert(math.approxEq(f32, sqrt32(19.230934), 4.385309, epsilon)); - assert(sqrt32(64.0) == 8.0); - assert(math.approxEq(f32, sqrt32(64.1), 8.006248, epsilon)); - assert(math.approxEq(f32, sqrt32(8942.230469), 94.563370, epsilon)); + assert(@sqrt(f32, 0.0) == 0.0); + assert(math.approxEq(f32, @sqrt(f32, 2.0), 1.414214, epsilon)); + assert(math.approxEq(f32, @sqrt(f32, 3.6), 1.897367, epsilon)); + assert(@sqrt(f32, 4.0) == 2.0); + assert(math.approxEq(f32, @sqrt(f32, 7.539840), 2.745877, epsilon)); + assert(math.approxEq(f32, @sqrt(f32, 19.230934), 4.385309, epsilon)); + assert(@sqrt(f32, 64.0) == 8.0); + assert(math.approxEq(f32, @sqrt(f32, 64.1), 8.006248, epsilon)); + assert(math.approxEq(f32, @sqrt(f32, 8942.230469), 94.563370, epsilon)); } test "math.sqrt64" { const epsilon = 0.000001; - assert(sqrt64(0.0) == 0.0); - assert(math.approxEq(f64, sqrt64(2.0), 1.414214, epsilon)); - assert(math.approxEq(f64, sqrt64(3.6), 1.897367, epsilon)); - assert(sqrt64(4.0) == 2.0); - assert(math.approxEq(f64, sqrt64(7.539840), 2.745877, epsilon)); - assert(math.approxEq(f64, sqrt64(19.230934), 4.385309, epsilon)); - assert(sqrt64(64.0) == 8.0); - assert(math.approxEq(f64, sqrt64(64.1), 8.006248, epsilon)); - assert(math.approxEq(f64, sqrt64(8942.230469), 94.563367, epsilon)); + assert(@sqrt(f64, 0.0) == 0.0); + assert(math.approxEq(f64, @sqrt(f64, 2.0), 1.414214, epsilon)); + assert(math.approxEq(f64, @sqrt(f64, 3.6), 1.897367, epsilon)); + assert(@sqrt(f64, 4.0) == 2.0); + assert(math.approxEq(f64, @sqrt(f64, 7.539840), 2.745877, epsilon)); + assert(math.approxEq(f64, @sqrt(f64, 19.230934), 4.385309, epsilon)); + assert(@sqrt(f64, 64.0) == 8.0); + assert(math.approxEq(f64, @sqrt(f64, 64.1), 8.006248, epsilon)); + assert(math.approxEq(f64, @sqrt(f64, 8942.230469), 94.563367, epsilon)); } test "math.sqrt32.special" { - assert(math.isPositiveInf(sqrt32(math.inf(f32)))); - assert(sqrt32(0.0) == 0.0); - assert(sqrt32(-0.0) == -0.0); - assert(math.isNan(sqrt32(-1.0))); - assert(math.isNan(sqrt32(math.nan(f32)))); + assert(math.isPositiveInf(@sqrt(f32, math.inf(f32)))); + assert(@sqrt(f32, 0.0) == 0.0); + assert(@sqrt(f32, -0.0) == -0.0); + assert(math.isNan(@sqrt(f32, -1.0))); + assert(math.isNan(@sqrt(f32, math.nan(f32)))); } test "math.sqrt64.special" { - assert(math.isPositiveInf(sqrt64(math.inf(f64)))); - assert(sqrt64(0.0) == 0.0); - assert(sqrt64(-0.0) == -0.0); - assert(math.isNan(sqrt64(-1.0))); - assert(math.isNan(sqrt64(math.nan(f64)))); + assert(math.isPositiveInf(@sqrt(f64, math.inf(f64)))); + assert(@sqrt(f64, 0.0) == 0.0); + assert(@sqrt(f64, -0.0) == -0.0); + assert(math.isNan(@sqrt(f64, -1.0))); + assert(math.isNan(@sqrt(f64, math.nan(f64)))); } fn sqrt_int(comptime T: type, value: T) @IntType(false, T.bit_count / 2) { diff --git a/std/math/x86_64/sqrt.zig b/std/math/x86_64/sqrt.zig deleted file mode 100644 index ad9ce0c96c..0000000000 --- a/std/math/x86_64/sqrt.zig +++ /dev/null @@ -1,15 +0,0 @@ -pub fn sqrt32(x: f32) f32 { - return asm ( - \\sqrtss %%xmm0, %%xmm0 - : [ret] "={xmm0}" (-> f32) - : [x] "{xmm0}" (x) - ); -} - -pub fn sqrt64(x: f64) f64 { - return asm ( - \\sqrtsd %%xmm0, %%xmm0 - : [ret] "={xmm0}" (-> f64) - : [x] "{xmm0}" (x) - ); -} |