mulc3.zig (2275B) - Raw
1 const std = @import("std"); 2 const isNan = std.math.isNan; 3 const isInf = std.math.isInf; 4 const copysign = std.math.copysign; 5 6 pub 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) 14 pub 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 }