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