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;
10
11const compiler_rt = @import("../compiler_rt.zig");
12const symbol = compiler_rt.symbol;
13const Complex = compiler_rt.Complex;
14
15comptime {
16 if (@import("builtin").zig_backend != .stage2_c) {
17 symbol(&__divhc3, "__divhc3");
18 symbol(&__divsc3, "__divsc3");
19 symbol(&__divdc3, "__divdc3");
20 symbol(&__divxc3, "__divxc3");
21 if (compiler_rt.want_ppc_abi) {
22 symbol(&__divtc3, "__divkc3");
23 } else {
24 symbol(&__divtc3, "__divtc3");
25 }
26 }
27}
28
29fn __divhc3(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 {
30 return compiler_rt.f16.complex.toAbi(div_cf16(
31 compiler_rt.f16.complex.fromAbi(.{ .real = lhs_real, .imag = lhs_imag }),
32 compiler_rt.f16.complex.fromAbi(.{ .real = rhs_real, .imag = rhs_imag }),
33 ));
34}
35pub fn div_cf16(a: Complex(f16), b: Complex(f16)) Complex(f16) {
36 return divc3(f16, a, b);
37}
38
39fn __divsc3(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 {
40 return compiler_rt.f32.complex.toAbi(div_cf32(
41 compiler_rt.f32.complex.fromAbi(.{ .real = lhs_real, .imag = lhs_imag }),
42 compiler_rt.f32.complex.fromAbi(.{ .real = rhs_real, .imag = rhs_imag }),
43 ));
44}
45pub fn div_cf32(a: Complex(f32), b: Complex(f32)) Complex(f32) {
46 return divc3(f32, a, b);
47}
48
49fn __divdc3(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 {
50 return compiler_rt.f64.complex.toAbi(div_cf64(
51 compiler_rt.f64.complex.fromAbi(.{ .real = lhs_real, .imag = lhs_imag }),
52 compiler_rt.f64.complex.fromAbi(.{ .real = rhs_real, .imag = rhs_imag }),
53 ));
54}
55pub fn div_cf64(a: Complex(f64), b: Complex(f64)) Complex(f64) {
56 return divc3(f64, a, b);
57}
58
59fn __divxc3(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 {
60 return compiler_rt.f80.complex.toAbi(div_cf80(
61 compiler_rt.f80.complex.fromAbi(.{ .real = lhs_real, .imag = lhs_imag }),
62 compiler_rt.f80.complex.fromAbi(.{ .real = rhs_real, .imag = rhs_imag }),
63 ));
64}
65pub fn div_cf80(a: Complex(f80), b: Complex(f80)) Complex(f80) {
66 return divc3(f80, a, b);
67}
68
69fn __divtc3(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 {
70 return compiler_rt.f128.complex.toAbi(div_cf128(
71 compiler_rt.f128.complex.fromAbi(.{ .real = lhs_real, .imag = lhs_imag }),
72 compiler_rt.f128.complex.fromAbi(.{ .real = rhs_real, .imag = rhs_imag }),
73 ));
74}
75pub fn div_cf128(a: Complex(f128), b: Complex(f128)) Complex(f128) {
76 return divc3(f128, a, b);
77}
78
79/// Implementation based on Annex G of C17 Standard (N2176)
80inline fn divc3(comptime T: type, lhs: Complex(T), rhs: Complex(T)) Complex(T) {
81 const a = lhs.real;
82 const b = lhs.imag;
83 var c = rhs.real;
84 var d = rhs.imag;
85
86 // logbw used to prevent under/over-flow
87 const logbw = ilogb(@max(@abs(c), @abs(d)));
88 const logbw_finite = logbw != maxInt(i32) and logbw != minInt(i32);
89 const ilogbw = if (logbw_finite) b: {
90 c = scalbn(c, -logbw);
91 d = scalbn(d, -logbw);
92 break :b logbw;
93 } else 0;
94 const denom = c * c + d * d;
95 const result: Complex(T) = .{
96 .real = scalbn((a * c + b * d) / denom, -ilogbw),
97 .imag = scalbn((b * c - a * d) / denom, -ilogbw),
98 };
99
100 // Recover infinities and zeros that computed as NaN+iNaN;
101 // the only cases are non-zero/zero, infinite/finite, and finite/infinite, ...
102 if (isNan(result.real) and isNan(result.imag)) {
103 const zero: T = 0.0;
104 const one: T = 1.0;
105
106 if ((denom == 0.0) and (!isNan(a) or !isNan(b))) {
107 return .{
108 .real = copysign(std.math.inf(T), c) * a,
109 .imag = copysign(std.math.inf(T), c) * b,
110 };
111 } else if ((isInf(a) or isInf(b)) and isFinite(c) and isFinite(d)) {
112 const boxed_a = copysign(if (isInf(a)) one else zero, a);
113 const boxed_b = copysign(if (isInf(b)) one else zero, b);
114 return .{
115 .real = std.math.inf(T) * (boxed_a * c - boxed_b * d),
116 .imag = std.math.inf(T) * (boxed_b * c - boxed_a * d),
117 };
118 } else if (logbw == maxInt(i32) and isFinite(a) and isFinite(b)) {
119 const boxed_c = copysign(if (isInf(c)) one else zero, c);
120 const boxed_d = copysign(if (isInf(d)) one else zero, d);
121 return .{
122 .real = 0.0 * (a * boxed_c + b * boxed_d),
123 .imag = 0.0 * (b * boxed_c - a * boxed_d),
124 };
125 }
126 }
127
128 return result;
129}
130
131test {
132 _ = @import("divc3_test.zig");
133}