Add initial complex-number support
- Library type instead of builtin - All C complex functions implemented Partial WIP: Needs more tests for edge cases.master
Marc Tiehuis
parent
15bf0c1541
commit
d5e99cc05e
|
|
@ -498,6 +498,28 @@ set(ZIG_STD_FILES
|
|||
"math/tan.zig"
|
||||
"math/tanh.zig"
|
||||
"math/trunc.zig"
|
||||
"math/complex/abs.zig"
|
||||
"math/complex/acosh.zig"
|
||||
"math/complex/acos.zig"
|
||||
"math/complex/arg.zig"
|
||||
"math/complex/asinh.zig"
|
||||
"math/complex/asin.zig"
|
||||
"math/complex/atanh.zig"
|
||||
"math/complex/atan.zig"
|
||||
"math/complex/conj.zig"
|
||||
"math/complex/cosh.zig"
|
||||
"math/complex/cos.zig"
|
||||
"math/complex/exp.zig"
|
||||
"math/complex/index.zig"
|
||||
"math/complex/ldexp.zig"
|
||||
"math/complex/log.zig"
|
||||
"math/complex/pow.zig"
|
||||
"math/complex/proj.zig"
|
||||
"math/complex/sinh.zig"
|
||||
"math/complex/sin.zig"
|
||||
"math/complex/sqrt.zig"
|
||||
"math/complex/tanh.zig"
|
||||
"math/complex/tan.zig"
|
||||
"mem.zig"
|
||||
"net.zig"
|
||||
"os/child_process.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));
|
||||
}
|
||||
|
|
@ -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 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);
|
||||
|
||||
const re = c.re;
|
||||
const im = c.im;
|
||||
|
||||
debug.assert(math.approxEq(f32, re, 0.546975, epsilon));
|
||||
debug.assert(math.approxEq(f32, im, -2.452914, epsilon));
|
||||
}
|
||||
|
|
@ -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 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);
|
||||
|
||||
const re = c.re;
|
||||
const im = c.im;
|
||||
|
||||
debug.assert(math.approxEq(f32, re, 2.452914, epsilon));
|
||||
debug.assert(math.approxEq(f32, im, 0.546975, epsilon));
|
||||
}
|
||||
|
|
@ -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));
|
||||
}
|
||||
|
|
@ -0,0 +1,30 @@
|
|||
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);
|
||||
|
||||
const re = c.re;
|
||||
const im = c.im;
|
||||
|
||||
debug.assert(math.approxEq(f32, re, 1.023822, epsilon));
|
||||
debug.assert(math.approxEq(f32, im, 2.452914, epsilon));
|
||||
}
|
||||
|
|
@ -0,0 +1,25 @@
|
|||
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);
|
||||
|
||||
const re = c.re;
|
||||
const im = c.im;
|
||||
|
||||
debug.assert(math.approxEq(f32, re, 2.459831, epsilon));
|
||||
debug.assert(math.approxEq(f32, im, 0.533999, epsilon));
|
||||
}
|
||||
|
|
@ -0,0 +1,136 @@
|
|||
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.ctan32" {
|
||||
const a = Complex(f32).new(5, 3);
|
||||
const c = atan(a);
|
||||
|
||||
const re = c.re;
|
||||
const im = c.im;
|
||||
|
||||
debug.assert(math.approxEq(f32, re, 1.423679, epsilon));
|
||||
debug.assert(math.approxEq(f32, im, 0.086569, epsilon));
|
||||
}
|
||||
|
||||
test "complex.ctan64" {
|
||||
const a = Complex(f64).new(5, 3);
|
||||
const c = atan(a);
|
||||
|
||||
const re = c.re;
|
||||
const im = c.im;
|
||||
|
||||
debug.assert(math.approxEq(f64, re, 1.423679, epsilon));
|
||||
debug.assert(math.approxEq(f64, im, 0.086569, epsilon));
|
||||
}
|
||||
|
|
@ -0,0 +1,25 @@
|
|||
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);
|
||||
|
||||
const re = c.re;
|
||||
const im = c.im;
|
||||
|
||||
debug.assert(math.approxEq(f32, re, 0.146947, epsilon));
|
||||
debug.assert(math.approxEq(f32, im, 1.480870, epsilon));
|
||||
}
|
||||
|
|
@ -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);
|
||||
}
|
||||
|
|
@ -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 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);
|
||||
|
||||
const re = c.re;
|
||||
const im = c.im;
|
||||
|
||||
debug.assert(math.approxEq(f32, re, 2.855815, epsilon));
|
||||
debug.assert(math.approxEq(f32, im, 9.606383, epsilon));
|
||||
}
|
||||
|
|
@ -0,0 +1,171 @@
|
|||
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);
|
||||
|
||||
const re = c.re;
|
||||
const im = c.im;
|
||||
|
||||
debug.assert(math.approxEq(f32, re, -73.467300, epsilon));
|
||||
debug.assert(math.approxEq(f32, im, 10.471557, epsilon));
|
||||
}
|
||||
|
||||
test "complex.ccosh64" {
|
||||
const a = Complex(f64).new(5, 3);
|
||||
const c = cosh(a);
|
||||
|
||||
const re = c.re;
|
||||
const im = c.im;
|
||||
|
||||
debug.assert(math.approxEq(f64, re, -73.467300, epsilon));
|
||||
debug.assert(math.approxEq(f64, im, 10.471557, epsilon));
|
||||
}
|
||||
|
|
@ -0,0 +1,146 @@
|
|||
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);
|
||||
|
||||
const re = c.re;
|
||||
const im = c.im;
|
||||
|
||||
debug.assert(math.approxEq(f32, re, -146.927917, epsilon));
|
||||
debug.assert(math.approxEq(f32, im, 20.944065, epsilon));
|
||||
}
|
||||
|
||||
test "complex.cexp64" {
|
||||
const a = Complex(f32).new(5, 3);
|
||||
const c = exp(a);
|
||||
|
||||
const re = c.re;
|
||||
const im = c.im;
|
||||
|
||||
debug.assert(math.approxEq(f64, re, -146.927917, epsilon));
|
||||
debug.assert(math.approxEq(f64, im, 20.944065, epsilon));
|
||||
}
|
||||
|
|
@ -0,0 +1,172 @@
|
|||
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;
|
||||
|
||||
// NOTE: Make this a builtin at some point?
|
||||
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");
|
||||
}
|
||||
|
|
@ -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,
|
||||
);
|
||||
}
|
||||
|
|
@ -0,0 +1,26 @@
|
|||
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);
|
||||
|
||||
const re = c.re;
|
||||
const im = c.im;
|
||||
|
||||
debug.assert(math.approxEq(f32, re, 1.763180, epsilon));
|
||||
debug.assert(math.approxEq(f32, im, 0.540419, epsilon));
|
||||
}
|
||||
|
|
@ -0,0 +1,25 @@
|
|||
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);
|
||||
|
||||
const re = c.re;
|
||||
const im = c.im;
|
||||
|
||||
debug.assert(math.approxEq(f32, re, 58.049110, epsilon));
|
||||
debug.assert(math.approxEq(f32, im, -101.003433, epsilon));
|
||||
}
|
||||
|
|
@ -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);
|
||||
}
|
||||
|
|
@ -0,0 +1,25 @@
|
|||
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);
|
||||
|
||||
const re = c.re;
|
||||
const im = c.im;
|
||||
|
||||
debug.assert(math.approxEq(f32, re, -9.654126, epsilon));
|
||||
debug.assert(math.approxEq(f32, im, 2.841692, epsilon));
|
||||
}
|
||||
|
|
@ -0,0 +1,170 @@
|
|||
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);
|
||||
|
||||
const re = c.re;
|
||||
const im = c.im;
|
||||
|
||||
debug.assert(math.approxEq(f32, re, -73.460617, epsilon));
|
||||
debug.assert(math.approxEq(f32, im, 10.472508, epsilon));
|
||||
}
|
||||
|
||||
test "complex.csinh64" {
|
||||
const a = Complex(f64).new(5, 3);
|
||||
const c = sinh(a);
|
||||
|
||||
const re = c.re;
|
||||
const im = c.im;
|
||||
|
||||
debug.assert(math.approxEq(f64, re, -73.460617, epsilon));
|
||||
debug.assert(math.approxEq(f64, im, 10.472508, epsilon));
|
||||
}
|
||||
|
|
@ -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;
|
||||
|
||||
// NOTE: Returning @typeOf(z) here causes issues when trying to access the value. This is
|
||||
// why we currently assign re, im parts to a new value explicitly for all tests.
|
||||
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);
|
||||
|
||||
const re = c.re;
|
||||
const im = c.im;
|
||||
|
||||
debug.assert(math.approxEq(f32, re, 2.327117, epsilon));
|
||||
debug.assert(math.approxEq(f32, im, 0.644574, epsilon));
|
||||
}
|
||||
|
||||
test "complex.csqrt64" {
|
||||
const a = Complex(f64).new(5, 3);
|
||||
const c = sqrt(a);
|
||||
|
||||
const re = c.re;
|
||||
const im = c.im;
|
||||
|
||||
debug.assert(math.approxEq(f64, re, 2.3271175190399496, epsilon));
|
||||
debug.assert(math.approxEq(f64, im, 0.6445742373246469, epsilon));
|
||||
}
|
||||
|
|
@ -0,0 +1,25 @@
|
|||
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);
|
||||
|
||||
const re = c.re;
|
||||
const im = c.im;
|
||||
|
||||
debug.assert(math.approxEq(f32, re, -0.002708233, epsilon));
|
||||
debug.assert(math.approxEq(f32, im, 1.004165, epsilon));
|
||||
}
|
||||
|
|
@ -0,0 +1,117 @@
|
|||
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);
|
||||
|
||||
const re = c.re;
|
||||
const im = c.im;
|
||||
|
||||
debug.assert(math.approxEq(f32, re, 0.999913, epsilon));
|
||||
debug.assert(math.approxEq(f32, im, -0.000025, epsilon));
|
||||
}
|
||||
|
||||
test "complex.ctanh64" {
|
||||
const a = Complex(f64).new(5, 3);
|
||||
const c = tanh(a);
|
||||
|
||||
const re = c.re;
|
||||
const im = c.im;
|
||||
|
||||
debug.assert(math.approxEq(f64, re, 0.999913, epsilon));
|
||||
debug.assert(math.approxEq(f64, im, -0.000025, epsilon));
|
||||
}
|
||||
|
|
@ -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");
|
||||
}
|
||||
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue