1const std = @import("std");
2const isNan = std.math.isNan;
3const isInf = std.math.isInf;
4const copysign = std.math.copysign;
5
6const compiler_rt = @import("../compiler_rt.zig");
7const symbol = compiler_rt.symbol;
8const Complex = compiler_rt.Complex;
9
10comptime {
11 if (@import("builtin").zig_backend != .stage2_c) {
12 symbol(&__mulhc3, "__mulhc3");
13 symbol(&__mulsc3, "__mulsc3");
14 symbol(&__muldc3, "__muldc3");
15 symbol(&__mulxc3, "__mulxc3");
16 if (compiler_rt.want_ppc_abi) {
17 symbol(&__multc3, "__mulkc3");
18 } else {
19 symbol(&__multc3, "__multc3");
20 }
21 }
22}
23
24fn __mulhc3(lhs_real: compiler_rt.f16.Abi, lhs_imag: compiler_rt.f16.Abi, rhs_real: compiler_rt.f16.Abi, rhs_imag: compiler_rt.f16.Abi) callconv(.c) compiler_rt.f16.complex.Abi {
25 return compiler_rt.f16.complex.toAbi(mul_cf16(
26 compiler_rt.f16.complex.fromAbi(.{ .real = lhs_real, .imag = lhs_imag }),
27 compiler_rt.f16.complex.fromAbi(.{ .real = rhs_real, .imag = rhs_imag }),
28 ));
29}
30pub fn mul_cf16(a: Complex(f16), b: Complex(f16)) Complex(f16) {
31 return mulc3(f16, a, b);
32}
33
34fn __mulsc3(lhs_real: compiler_rt.f32.Abi, lhs_imag: compiler_rt.f32.Abi, rhs_real: compiler_rt.f32.Abi, rhs_imag: compiler_rt.f32.Abi) callconv(.c) compiler_rt.f32.complex.Abi {
35 return compiler_rt.f32.complex.toAbi(mul_cf32(
36 compiler_rt.f32.complex.fromAbi(.{ .real = lhs_real, .imag = lhs_imag }),
37 compiler_rt.f32.complex.fromAbi(.{ .real = rhs_real, .imag = rhs_imag }),
38 ));
39}
40pub fn mul_cf32(a: Complex(f32), b: Complex(f32)) Complex(f32) {
41 return mulc3(f32, a, b);
42}
43
44fn __muldc3(lhs_real: compiler_rt.f64.Abi, lhs_imag: compiler_rt.f64.Abi, rhs_real: compiler_rt.f64.Abi, rhs_imag: compiler_rt.f64.Abi) callconv(.c) compiler_rt.f64.complex.Abi {
45 return compiler_rt.f64.complex.toAbi(mul_cf64(
46 compiler_rt.f64.complex.fromAbi(.{ .real = lhs_real, .imag = lhs_imag }),
47 compiler_rt.f64.complex.fromAbi(.{ .real = rhs_real, .imag = rhs_imag }),
48 ));
49}
50pub fn mul_cf64(a: Complex(f64), b: Complex(f64)) Complex(f64) {
51 return mulc3(f64, a, b);
52}
53
54fn __mulxc3(lhs_real: compiler_rt.f80.Abi, lhs_imag: compiler_rt.f80.Abi, rhs_real: compiler_rt.f80.Abi, rhs_imag: compiler_rt.f80.Abi) callconv(.c) compiler_rt.f80.complex.Abi {
55 return compiler_rt.f80.complex.toAbi(mul_cf80(
56 compiler_rt.f80.complex.fromAbi(.{ .real = lhs_real, .imag = lhs_imag }),
57 compiler_rt.f80.complex.fromAbi(.{ .real = rhs_real, .imag = rhs_imag }),
58 ));
59}
60pub fn mul_cf80(a: Complex(f80), b: Complex(f80)) Complex(f80) {
61 return mulc3(f80, a, b);
62}
63
64fn __multc3(lhs_real: compiler_rt.f128.Abi, lhs_imag: compiler_rt.f128.Abi, rhs_real: compiler_rt.f128.Abi, rhs_imag: compiler_rt.f128.Abi) callconv(.c) compiler_rt.f128.complex.Abi {
65 return compiler_rt.f128.complex.toAbi(mul_cf128(
66 compiler_rt.f128.complex.fromAbi(.{ .real = lhs_real, .imag = lhs_imag }),
67 compiler_rt.f128.complex.fromAbi(.{ .real = rhs_real, .imag = rhs_imag }),
68 ));
69}
70pub fn mul_cf128(a: Complex(f128), b: Complex(f128)) Complex(f128) {
71 return mulc3(f128, a, b);
72}
73
74/// Implementation based on Annex G of C17 Standard (N2176)
75inline fn mulc3(comptime T: type, lhs: Complex(T), rhs: Complex(T)) Complex(T) {
76 var a = lhs.real;
77 var b = lhs.imag;
78 var c = rhs.real;
79 var d = rhs.imag;
80
81 const ac = a * c;
82 const bd = b * d;
83 const ad = a * d;
84 const bc = b * c;
85
86 const zero: T = 0.0;
87 const one: T = 1.0;
88
89 const z: Complex(T) = .{
90 .real = ac - bd,
91 .imag = ad + bc,
92 };
93 if (isNan(z.real) and isNan(z.imag)) {
94 var recalc: bool = false;
95
96 if (isInf(a) or isInf(b)) { // (a + ib) is infinite
97
98 // "Box" the infinity (+/-inf goes to +/-1, all finite values go to 0)
99 a = copysign(if (isInf(a)) one else zero, a);
100 b = copysign(if (isInf(b)) one else zero, b);
101
102 // Replace NaNs in the other factor with (signed) 0
103 if (isNan(c)) c = copysign(zero, c);
104 if (isNan(d)) d = copysign(zero, d);
105
106 recalc = true;
107 }
108
109 if (isInf(c) or isInf(d)) { // (c + id) is infinite
110
111 // "Box" the infinity (+/-inf goes to +/-1, all finite values go to 0)
112 c = copysign(if (isInf(c)) one else zero, c);
113 d = copysign(if (isInf(d)) one else zero, d);
114
115 // Replace NaNs in the other factor with (signed) 0
116 if (isNan(a)) a = copysign(zero, a);
117 if (isNan(b)) b = copysign(zero, b);
118
119 recalc = true;
120 }
121
122 if (!recalc and (isInf(ac) or isInf(bd) or isInf(ad) or isInf(bc))) {
123
124 // Recover infinities from overflow by changing NaNs to 0
125 if (isNan(a)) a = copysign(zero, a);
126 if (isNan(b)) b = copysign(zero, b);
127 if (isNan(c)) c = copysign(zero, c);
128 if (isNan(d)) d = copysign(zero, d);
129
130 recalc = true;
131 }
132 if (recalc) {
133 return .{
134 .real = std.math.inf(T) * (a * c - b * d),
135 .imag = std.math.inf(T) * (a * d + b * c),
136 };
137 }
138 }
139 return z;
140}
141
142test {
143 _ = @import("mulc3_test.zig");
144}