1const std = @import("std");
 2const isNan = std.math.isNan;
 3const isInf = std.math.isInf;
 4const copysign = std.math.copysign;
 5
 6pub fn Complex(comptime T: type) type {
 7    return extern struct {
 8        real: T,
 9        imag: T,
10    };
11}
12
13/// Implementation based on Annex G of C17 Standard (N2176)
14pub inline fn mulc3(comptime T: type, a_in: T, b_in: T, c_in: T, d_in: T) Complex(T) {
15    var a = a_in;
16    var b = b_in;
17    var c = c_in;
18    var d = d_in;
19
20    const ac = a * c;
21    const bd = b * d;
22    const ad = a * d;
23    const bc = b * c;
24
25    const zero: T = 0.0;
26    const one: T = 1.0;
27
28    const z: Complex(T) = .{
29        .real = ac - bd,
30        .imag = ad + bc,
31    };
32    if (isNan(z.real) and isNan(z.imag)) {
33        var recalc: bool = false;
34
35        if (isInf(a) or isInf(b)) { // (a + ib) is infinite
36
37            // "Box" the infinity (+/-inf goes to +/-1, all finite values go to 0)
38            a = copysign(if (isInf(a)) one else zero, a);
39            b = copysign(if (isInf(b)) one else zero, b);
40
41            // Replace NaNs in the other factor with (signed) 0
42            if (isNan(c)) c = copysign(zero, c);
43            if (isNan(d)) d = copysign(zero, d);
44
45            recalc = true;
46        }
47
48        if (isInf(c) or isInf(d)) { // (c + id) is infinite
49
50            // "Box" the infinity (+/-inf goes to +/-1, all finite values go to 0)
51            c = copysign(if (isInf(c)) one else zero, c);
52            d = copysign(if (isInf(d)) one else zero, d);
53
54            // Replace NaNs in the other factor with (signed) 0
55            if (isNan(a)) a = copysign(zero, a);
56            if (isNan(b)) b = copysign(zero, b);
57
58            recalc = true;
59        }
60
61        if (!recalc and (isInf(ac) or isInf(bd) or isInf(ad) or isInf(bc))) {
62
63            // Recover infinities from overflow by changing NaNs to 0
64            if (isNan(a)) a = copysign(zero, a);
65            if (isNan(b)) b = copysign(zero, b);
66            if (isNan(c)) c = copysign(zero, c);
67            if (isNan(d)) d = copysign(zero, d);
68
69            recalc = true;
70        }
71        if (recalc) {
72            return .{
73                .real = std.math.inf(T) * (a * c - b * d),
74                .imag = std.math.inf(T) * (a * d + b * c),
75            };
76        }
77    }
78    return z;
79}