Skip to content
215 changes: 143 additions & 72 deletions lib/std/math/complex.zig
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,6 @@ pub const asinh = @import("complex/asinh.zig").asinh;
pub const asin = @import("complex/asin.zig").asin;
pub const atanh = @import("complex/atanh.zig").atanh;
pub const atan = @import("complex/atan.zig").atan;
pub const conj = @import("complex/conj.zig").conj;
pub const cosh = @import("complex/cosh.zig").cosh;
pub const cos = @import("complex/cos.zig").cos;
pub const exp = @import("complex/exp.zig").exp;
Expand All @@ -23,7 +22,8 @@ pub const sqrt = @import("complex/sqrt.zig").sqrt;
pub const tanh = @import("complex/tanh.zig").tanh;
pub const tan = @import("complex/tan.zig").tan;

/// A complex number consisting of a real an imaginary part. T must be a floating-point value.
/// A complex number consisting of a real and imaginary part.
/// T must be a floating-point value.
pub fn Complex(comptime T: type) type {
return struct {
const Self = @This();
Expand All @@ -34,170 +34,242 @@ pub fn Complex(comptime T: type) type {
/// Imaginary part.
im: T,

/// Create a new Complex number from the given real and imaginary parts.
/// Imarinary unit that satisfies "i^2 = -1".
pub const i: Self = .init(0, 1);

/// Creates a new complex number from the given real and imaginary parts.
pub fn init(re: T, im: T) Self {
return Self{
return .{
.re = re,
.im = im,
};
}

/// Returns the sum of two complex numbers.
/// Calculates the sum of two complex numbers.
pub fn add(self: Self, other: Self) Self {
return Self{
return .{
.re = self.re + other.re,
.im = self.im + other.im,
};
}

/// Returns the subtraction of two complex numbers.
/// Calculates the subtraction of two complex numbers.
pub fn sub(self: Self, other: Self) Self {
return Self{
return .{
.re = self.re - other.re,
.im = self.im - other.im,
};
}

/// Returns the product of two complex numbers.
/// Calculates the product of two complex numbers.
pub fn mul(self: Self, other: Self) Self {
return Self{
return .{
.re = self.re * other.re - self.im * other.im,
.im = self.im * other.re + self.re * other.im,
};
}

/// Returns the quotient of two complex numbers.
/// Calculates the quotient of two complex numbers.
pub fn div(self: Self, other: 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{
return .{
.re = re_num / den,
.im = im_num / den,
};
}

/// Returns the complex conjugate of a number.
pub fn conjugate(self: Self) Self {
return Self{
.re = self.re,
/// Calculates the negation of a complex number.
pub fn neg(self: Self) Self {
return .{
.re = -self.re,
.im = -self.im,
};
}

/// Returns the negation of a complex number.
pub fn neg(self: Self) Self {
return Self{
.re = -self.re,
/// Calculates the complex conjugate of a complex number.
pub fn conj(self: Self) Self {
return .{
.re = self.re,
.im = -self.im,
};
}

/// Returns the product of complex number and i=sqrt(-1)
pub fn mulbyi(self: Self) Self {
return Self{
/// Calculates the product of a complex number and imaginary unit.
/// You should not manually does ".mul(.i, *)" instead of using this,
/// as its consumes more operations than this.
pub fn mulByI(self: Self) Self {
return .{
.re = -self.im,
.im = self.re,
};
}

/// Returns the reciprocal of a complex number.
pub fn reciprocal(self: Self) Self {
const m = self.re * self.re + self.im * self.im;
return Self{
.re = self.re / m,
.im = -self.im / m,
/// Calculates the product of a complex number and negation of imaginary unit,
/// thus this rotates 90 degrees clockwise on the complex plane.
/// You should not manually does "*.mul(.i).neg()" (or "*.neg().mul(.i)") instead of using this,
/// as its consumes more operations than this.
pub fn mulByMinusI(self: Self) Self {
return .{
.re = self.im,
.im = -self.re,
};
}

/// Returns the magnitude of a complex number.
/// Calculates the reciprocal of a complex number.
pub fn recip(self: Self) Self {
const magnitude_sq = self.squaredMagnitude();

return .{
.re = self.re / magnitude_sq,
.im = -self.im / magnitude_sq,
};
}

/// Calculates the magnitude of a complex number.
pub fn magnitude(self: Self) T {
return @sqrt(self.re * self.re + self.im * self.im);
return @sqrt(self.squaredMagnitude());
}

/// Calculates the squared magnitude of a complex number.
pub fn squaredMagnitude(self: Self) T {
return self.re * self.re + self.im * self.im;
}
};
}

const epsilon = 0.0001;
const TestingComplex = Complex(f32);

const testing_epsilon = 0.0001;

test "add" {
const a = Complex(f32).init(5, 3);
const b = Complex(f32).init(2, 7);
const c = a.add(b);
const a: TestingComplex = .init(5, 3);
const b: TestingComplex = .init(2, 7);

const a_add_b = a.add(b);

try testing.expect(c.re == 7 and c.im == 10);
try testing.expectEqual(7, a_add_b.re);
try testing.expectEqual(10, a_add_b.im);
}

test "sub" {
const a = Complex(f32).init(5, 3);
const b = Complex(f32).init(2, 7);
const c = a.sub(b);
const a: TestingComplex = .init(5, 3);
const b: TestingComplex = .init(2, 7);

try testing.expect(c.re == 3 and c.im == -4);
const a_sub_b = a.sub(b);

try testing.expectEqual(3, a_sub_b.re);
try testing.expectEqual(-4, a_sub_b.im);
}

test "mul" {
const a = Complex(f32).init(5, 3);
const b = Complex(f32).init(2, 7);
const c = a.mul(b);
const a: TestingComplex = .init(5, 3);
const b: TestingComplex = .init(2, 7);

const a_mul_b = a.mul(b);

try testing.expect(c.re == -11 and c.im == 41);
try testing.expectEqual(-11, a_mul_b.re);
try testing.expectEqual(41, a_mul_b.im);
}

test "div" {
const a = Complex(f32).init(5, 3);
const b = Complex(f32).init(2, 7);
const c = a.div(b);
const a: TestingComplex = .init(5, 3);
const b: TestingComplex = .init(2, 7);

try testing.expect(math.approxEqAbs(f32, c.re, @as(f32, 31) / 53, epsilon) and
math.approxEqAbs(f32, c.im, @as(f32, -29) / 53, epsilon));
const a_div_b = a.div(b);

try testing.expectApproxEqAbs(@as(f32, 31) / 53, a_div_b.re, testing_epsilon);
try testing.expectApproxEqAbs(@as(f32, -29) / 53, a_div_b.im, testing_epsilon);
}

test "conjugate" {
const a = Complex(f32).init(5, 3);
const c = a.conjugate();
test "conj" {
const a: TestingComplex = .init(5, 3);
const a_conj = a.conj();

try testing.expect(c.re == 5 and c.im == -3);
try testing.expectEqual(5, a_conj.re);
try testing.expectEqual(-3, a_conj.im);
}

test "neg" {
const a = Complex(f32).init(5, 3);
const c = a.neg();
const a: TestingComplex = .init(5, 3);
const neg_a = a.neg();

try testing.expectEqual(-5, neg_a.re);
try testing.expectEqual(-3, neg_a.im);
}

test "mulByI" {
const a: TestingComplex = .init(5, 3);
const i_a = a.mulByI();

try testing.expectEqual(-3, i_a.re);
try testing.expectEqual(5, i_a.im);
}

test "multiplication by i yields same result as mulByI" {
const a: TestingComplex = .init(5, 3);

const i_a_natural = a.mulByI();
const i_a_unnatural: TestingComplex = .mul(.i, a);

try testing.expectEqual(i_a_unnatural.re, i_a_natural.re);
try testing.expectEqual(i_a_unnatural.im, i_a_natural.im);
}

test "mulByMinusI" {
const a: TestingComplex = .init(5, 3);
const minus_i_a = a.mulByMinusI();

try testing.expectEqual(3, minus_i_a.re);
try testing.expectEqual(-5, minus_i_a.im);
}

test "multiplication by negation of i yields same result as mulByMinusI" {
const a: TestingComplex = .init(5, 3);

const minus_i_a_natural = a.mulByMinusI();
const minus_i_a_unnatural: TestingComplex = a.mul(.i).neg(); // x.mul(.i).neg() -> -ix

try testing.expectEqual(minus_i_a_unnatural.re, minus_i_a_natural.re);
try testing.expectEqual(minus_i_a_unnatural.im, minus_i_a_natural.im);
}

test "i^2 equals to -1" {
const a: TestingComplex = .mul(.i, .i);

try testing.expect(c.re == -5 and c.im == -3);
try testing.expectEqual(-1, a.re);
try testing.expectEqual(0, a.im);
}

test "mulbyi" {
const a = Complex(f32).init(5, 3);
const c = a.mulbyi();
test "(-i)^2 equals to -1" {
const a: TestingComplex = .mul(.neg(.i), .neg(.i));

try testing.expect(c.re == -3 and c.im == 5);
try testing.expectEqual(-1, a.re);
try testing.expectEqual(0, a.im);
}

test "reciprocal" {
const a = Complex(f32).init(5, 3);
const c = a.reciprocal();
test "recip" {
const a: TestingComplex = .init(5, 3);
const a_recip = a.recip();

try testing.expect(math.approxEqAbs(f32, c.re, @as(f32, 5) / 34, epsilon) and
math.approxEqAbs(f32, c.im, @as(f32, -3) / 34, epsilon));
try testing.expectApproxEqAbs(@as(f32, 5) / 34, a_recip.re, testing_epsilon);
try testing.expectApproxEqAbs(@as(f32, -3) / 34, a_recip.im, testing_epsilon);
}

test "magnitude" {
const a = Complex(f32).init(5, 3);
const c = a.magnitude();
const a: TestingComplex = .init(5, 3);
const a_magnitude = a.magnitude();

try testing.expect(math.approxEqAbs(f32, c, 5.83095, epsilon));
try testing.expectApproxEqAbs(5.83095, a_magnitude, testing_epsilon);
}

test "squaredMagnitude" {
const a = Complex(f32).init(5, 3);
const c = a.squaredMagnitude();
const a: TestingComplex = .init(5, 3);
const a_magnitude_sq = a.squaredMagnitude();

try testing.expect(math.approxEqAbs(f32, c, math.pow(f32, a.magnitude(), 2), epsilon));
try testing.expectApproxEqAbs(math.pow(f32, a.magnitude(), 2), a_magnitude_sq, testing_epsilon);
}

test {
Expand All @@ -209,7 +281,6 @@ test {
_ = @import("complex/asin.zig");
_ = @import("complex/atanh.zig");
_ = @import("complex/atan.zig");
_ = @import("complex/conj.zig");
_ = @import("complex/cosh.zig");
_ = @import("complex/cos.zig");
_ = @import("complex/exp.zig");
Expand Down
13 changes: 7 additions & 6 deletions lib/std/math/complex/abs.zig
Original file line number Diff line number Diff line change
@@ -1,17 +1,18 @@
const std = @import("../../std.zig");
const testing = std.testing;
const math = std.math;
const cmath = math.complex;
const Complex = cmath.Complex;
const Complex = math.Complex;

/// Returns the absolute value (modulus) of z.
/// Calculates the absolute value (modulus) of a complex number.
pub fn abs(z: anytype) @TypeOf(z.re, z.im) {
return math.hypot(z.re, z.im);
}

test abs {
const epsilon = math.floatEps(f32);
const a = Complex(f32).init(5, 3);
const c = abs(a);
try testing.expectApproxEqAbs(5.8309517, c, epsilon);

const a: Complex(f32) = .init(5, 3);
const a_abs = abs(a);

try testing.expectApproxEqAbs(5.8309517, a_abs, epsilon);
}
22 changes: 13 additions & 9 deletions lib/std/math/complex/acos.zig
Original file line number Diff line number Diff line change
@@ -1,21 +1,25 @@
const std = @import("../../std.zig");
const testing = std.testing;
const math = std.math;
const cmath = math.complex;
const Complex = cmath.Complex;
const Complex = math.Complex;

/// Returns the arc-cosine of z.
const asin = @import("asin.zig").asin;

/// Calculates the arc-cosine of a complex number.
pub fn acos(z: anytype) Complex(@TypeOf(z.re, z.im)) {
const T = @TypeOf(z.re, z.im);
const q = cmath.asin(z);
return Complex(T).init(@as(T, math.pi) / 2 - q.re, -q.im);

const q = asin(z);

return .init(@as(T, math.pi) / 2 - q.re, -q.im);
}

test acos {
const epsilon = math.floatEps(f32);
const a = Complex(f32).init(5, 3);
const c = acos(a);

try testing.expectApproxEqAbs(0.5469737, c.re, epsilon);
try testing.expectApproxEqAbs(-2.4529128, c.im, epsilon);
const a: Complex(f32) = .init(5, 3);
const acos_a = acos(a);

try testing.expectApproxEqAbs(0.5469737, acos_a.re, epsilon);
try testing.expectApproxEqAbs(-2.4529128, acos_a.im, epsilon);
}
Loading