zig

divc3.zig (2280B) - Raw

---

 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;
 }