1const std = @import("std");
 2const isNan = std.math.isNan;
 3const isInf = std.math.isInf;
 4const scalbn = std.math.scalbn;
 5const ilogb = std.math.ilogb;
 6const maxInt = std.math.maxInt;
 7const minInt = std.math.minInt;
 8const isFinite = std.math.isFinite;
 9const copysign = std.math.copysign;
10const Complex = @import("mulc3.zig").Complex;
11
12/// Implementation based on Annex G of C17 Standard (N2176)
13pub inline fn divc3(comptime T: type, a: T, b: T, c_in: T, d_in: T) Complex(T) {
14    var c = c_in;
15    var d = d_in;
16
17    // logbw used to prevent under/over-flow
18    const logbw = ilogb(@max(@abs(c), @abs(d)));
19    const logbw_finite = logbw != maxInt(i32) and logbw != minInt(i32);
20    const ilogbw = if (logbw_finite) b: {
21        c = scalbn(c, -logbw);
22        d = scalbn(d, -logbw);
23        break :b logbw;
24    } else 0;
25    const denom = c * c + d * d;
26    const result = Complex(T){
27        .real = scalbn((a * c + b * d) / denom, -ilogbw),
28        .imag = scalbn((b * c - a * d) / denom, -ilogbw),
29    };
30
31    // Recover infinities and zeros that computed as NaN+iNaN;
32    // the only cases are non-zero/zero, infinite/finite, and finite/infinite, ...
33    if (isNan(result.real) and isNan(result.imag)) {
34        const zero: T = 0.0;
35        const one: T = 1.0;
36
37        if ((denom == 0.0) and (!isNan(a) or !isNan(b))) {
38            return .{
39                .real = copysign(std.math.inf(T), c) * a,
40                .imag = copysign(std.math.inf(T), c) * b,
41            };
42        } else if ((isInf(a) or isInf(b)) and isFinite(c) and isFinite(d)) {
43            const boxed_a = copysign(if (isInf(a)) one else zero, a);
44            const boxed_b = copysign(if (isInf(b)) one else zero, b);
45            return .{
46                .real = std.math.inf(T) * (boxed_a * c - boxed_b * d),
47                .imag = std.math.inf(T) * (boxed_b * c - boxed_a * d),
48            };
49        } else if (logbw == maxInt(i32) and isFinite(a) and isFinite(b)) {
50            const boxed_c = copysign(if (isInf(c)) one else zero, c);
51            const boxed_d = copysign(if (isInf(d)) one else zero, d);
52            return .{
53                .real = 0.0 * (a * boxed_c + b * boxed_d),
54                .imag = 0.0 * (b * boxed_c - a * boxed_d),
55            };
56        }
57    }
58
59    return result;
60}