compiler-rt: Implement mulXc3 and divXc3 functions

These are the standard complex multiplication/division functions
required by the C standard (Annex G).

Don't get me started on the standard's handling of complex-infinity...

lib/compiler_rt.zig

@@ -15,11 +15,25 @@ comptime {
     _ = @import("compiler_rt/subxf3.zig");
 
     _ = @import("compiler_rt/mulf3.zig");
-    _ = @import("compiler_rt/muldf3.zig");
     _ = @import("compiler_rt/mulsf3.zig");
+    _ = @import("compiler_rt/muldf3.zig");
     _ = @import("compiler_rt/multf3.zig");
     _ = @import("compiler_rt/mulxf3.zig");
 
+    _ = @import("compiler_rt/mulc3.zig");
+    _ = @import("compiler_rt/mulhc3.zig");
+    _ = @import("compiler_rt/mulsc3.zig");
+    _ = @import("compiler_rt/muldc3.zig");
+    _ = @import("compiler_rt/mulxc3.zig");
+    _ = @import("compiler_rt/multc3.zig");
+
+    _ = @import("compiler_rt/divc3.zig");
+    _ = @import("compiler_rt/divhc3.zig");
+    _ = @import("compiler_rt/divsc3.zig");
+    _ = @import("compiler_rt/divdc3.zig");
+    _ = @import("compiler_rt/divxc3.zig");
+    _ = @import("compiler_rt/divtc3.zig");
+
     _ = @import("compiler_rt/negsf2.zig");
     _ = @import("compiler_rt/negdf2.zig");
     _ = @import("compiler_rt/negtf2.zig");

lib/compiler_rt/divc3.zig

@@ -0,0 +1,62 @@
+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 max = std.math.max;
+const fabs = std.math.fabs;
+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(fabs(c), fabs(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;
+}

lib/compiler_rt/divc3_test.zig

@@ -0,0 +1,77 @@
+const std = @import("std");
+const math = std.math;
+const expect = std.testing.expect;
+
+const Complex = @import("./mulc3.zig").Complex;
+const __divhc3 = @import("./divhc3.zig").__divhc3;
+const __divsc3 = @import("./divsc3.zig").__divsc3;
+const __divdc3 = @import("./divdc3.zig").__divdc3;
+const __divxc3 = @import("./divxc3.zig").__divxc3;
+const __divtc3 = @import("./divtc3.zig").__divtc3;
+
+test {
+    try testDiv(f16, __divhc3);
+    try testDiv(f32, __divsc3);
+    try testDiv(f64, __divdc3);
+    try testDiv(f80, __divxc3);
+    try testDiv(f128, __divtc3);
+}
+
+fn testDiv(comptime T: type, comptime f: fn (T, T, T, T) callconv(.C) Complex(T)) !void {
+    {
+        var a: T = 1.0;
+        var b: T = 0.0;
+        var c: T = -1.0;
+        var d: T = 0.0;
+
+        const result = f(a, b, c, d);
+        try expect(result.real == -1.0);
+        try expect(result.imag == 0.0);
+    }
+    {
+        var a: T = 1.0;
+        var b: T = 0.0;
+        var c: T = -4.0;
+        var d: T = 0.0;
+
+        const result = f(a, b, c, d);
+        try expect(result.real == -0.25);
+        try expect(result.imag == 0.0);
+    }
+    {
+        // if the first operand is an infinity and the second operand is a finite number, then the
+        // result of the / operator is an infinity;
+        var a: T = -math.inf(T);
+        var b: T = 0.0;
+        var c: T = -4.0;
+        var d: T = 1.0;
+
+        const result = f(a, b, c, d);
+        try expect(result.real == math.inf(T));
+        try expect(result.imag == math.inf(T));
+    }
+    {
+        // if the first operand is a finite number and the second operand is an infinity, then the
+        // result of the / operator is a zero;
+        var a: T = 17.2;
+        var b: T = 0.0;
+        var c: T = -math.inf(T);
+        var d: T = 0.0;
+
+        const result = f(a, b, c, d);
+        try expect(result.real == -0.0);
+        try expect(result.imag == 0.0);
+    }
+    {
+        // if the first operand is a nonzero finite number or an infinity and the second operand is
+        // a zero, then the result of the / operator is an infinity
+        var a: T = 1.1;
+        var b: T = 0.1;
+        var c: T = 0.0;
+        var d: T = 0.0;
+
+        const result = f(a, b, c, d);
+        try expect(result.real == math.inf(T));
+        try expect(result.imag == math.inf(T));
+    }
+}

lib/compiler_rt/divdc3.zig

@@ -0,0 +1,11 @@
+const common = @import("./common.zig");
+const divc3 = @import("./divc3.zig");
+const Complex = @import("./mulc3.zig").Complex;
+
+comptime {
+    @export(__divdc3, .{ .name = "__divdc3", .linkage = common.linkage });
+}
+
+pub fn __divdc3(a: f64, b: f64, c: f64, d: f64) callconv(.C) Complex(f64) {
+    return divc3.divc3(f64, a, b, c, d);
+}

lib/compiler_rt/divhc3.zig

@@ -0,0 +1,11 @@
+const common = @import("./common.zig");
+const divc3 = @import("./divc3.zig");
+const Complex = @import("./mulc3.zig").Complex;
+
+comptime {
+    @export(__divhc3, .{ .name = "__divhc3", .linkage = common.linkage });
+}
+
+pub fn __divhc3(a: f16, b: f16, c: f16, d: f16) callconv(.C) Complex(f16) {
+    return divc3.divc3(f16, a, b, c, d);
+}

lib/compiler_rt/divsc3.zig

@@ -0,0 +1,11 @@
+const common = @import("./common.zig");
+const divc3 = @import("./divc3.zig");
+const Complex = @import("./mulc3.zig").Complex;
+
+comptime {
+    @export(__divsc3, .{ .name = "__divsc3", .linkage = common.linkage });
+}
+
+pub fn __divsc3(a: f32, b: f32, c: f32, d: f32) callconv(.C) Complex(f32) {
+    return divc3.divc3(f32, a, b, c, d);
+}

lib/compiler_rt/divtc3.zig

@@ -0,0 +1,11 @@
+const common = @import("./common.zig");
+const divc3 = @import("./divc3.zig");
+const Complex = @import("./mulc3.zig").Complex;
+
+comptime {
+    @export(__divtc3, .{ .name = "__divtc3", .linkage = common.linkage });
+}
+
+pub fn __divtc3(a: f128, b: f128, c: f128, d: f128) callconv(.C) Complex(f128) {
+    return divc3.divc3(f128, a, b, c, d);
+}

lib/compiler_rt/divxc3.zig

@@ -0,0 +1,11 @@
+const common = @import("./common.zig");
+const divc3 = @import("./divc3.zig");
+const Complex = @import("./mulc3.zig").Complex;
+
+comptime {
+    @export(__divxc3, .{ .name = "__divxc3", .linkage = common.linkage });
+}
+
+pub fn __divxc3(a: f80, b: f80, c: f80, d: f80) callconv(.C) Complex(f80) {
+    return divc3.divc3(f80, a, b, c, d);
+}

lib/compiler_rt/mulc3.zig

@@ -0,0 +1,79 @@
+const std = @import("std");
+const isNan = std.math.isNan;
+const isInf = std.math.isInf;
+const copysign = std.math.copysign;
+
+pub fn Complex(comptime T: type) type {
+    return extern struct {
+        real: T,
+        imag: T,
+    };
+}
+
+/// Implementation based on Annex G of C17 Standard (N2176)
+pub inline fn mulc3(comptime T: type, a_in: T, b_in: T, c_in: T, d_in: T) Complex(T) {
+    var a = a_in;
+    var b = b_in;
+    var c = c_in;
+    var d = d_in;
+
+    const ac = a * c;
+    const bd = b * d;
+    const ad = a * d;
+    const bc = b * c;
+
+    const zero: T = 0.0;
+    const one: T = 1.0;
+
+    var z = Complex(T){
+        .real = ac - bd,
+        .imag = ad + bc,
+    };
+    if (isNan(z.real) and isNan(z.imag)) {
+        var recalc: bool = false;
+
+        if (isInf(a) or isInf(b)) { // (a + ib) is infinite
+
+            // "Box" the infinity (+/-inf goes to +/-1, all finite values go to 0)
+            a = copysign(if (isInf(a)) one else zero, a);
+            b = copysign(if (isInf(b)) one else zero, b);
+
+            // Replace NaNs in the other factor with (signed) 0
+            if (isNan(c)) c = copysign(zero, c);
+            if (isNan(d)) d = copysign(zero, d);
+
+            recalc = true;
+        }
+
+        if (isInf(c) or isInf(d)) { // (c + id) is infinite
+
+            // "Box" the infinity (+/-inf goes to +/-1, all finite values go to 0)
+            c = copysign(if (isInf(c)) one else zero, c);
+            d = copysign(if (isInf(d)) one else zero, d);
+
+            // Replace NaNs in the other factor with (signed) 0
+            if (isNan(a)) a = copysign(zero, a);
+            if (isNan(b)) b = copysign(zero, b);
+
+            recalc = true;
+        }
+
+        if (!recalc and (isInf(ac) or isInf(bd) or isInf(ad) or isInf(bc))) {
+
+            // Recover infinities from overflow by changing NaNs to 0
+            if (isNan(a)) a = copysign(zero, a);
+            if (isNan(b)) b = copysign(zero, b);
+            if (isNan(c)) c = copysign(zero, c);
+            if (isNan(d)) d = copysign(zero, d);
+
+            recalc = true;
+        }
+        if (recalc) {
+            return .{
+                .real = std.math.inf(T) * (a * c - b * d),
+                .imag = std.math.inf(T) * (a * d + b * c),
+            };
+        }
+    }
+    return z;
+}

lib/compiler_rt/mulc3_test.zig

@@ -0,0 +1,65 @@
+const std = @import("std");
+const math = std.math;
+const expect = std.testing.expect;
+
+const Complex = @import("./mulc3.zig").Complex;
+const __mulhc3 = @import("./mulhc3.zig").__mulhc3;
+const __mulsc3 = @import("./mulsc3.zig").__mulsc3;
+const __muldc3 = @import("./muldc3.zig").__muldc3;
+const __mulxc3 = @import("./mulxc3.zig").__mulxc3;
+const __multc3 = @import("./multc3.zig").__multc3;
+
+test {
+    try testMul(f16, __mulhc3);
+    try testMul(f32, __mulsc3);
+    try testMul(f64, __muldc3);
+    try testMul(f80, __mulxc3);
+    try testMul(f128, __multc3);
+}
+
+fn testMul(comptime T: type, comptime f: fn (T, T, T, T) callconv(.C) Complex(T)) !void {
+    {
+        var a: T = 1.0;
+        var b: T = 0.0;
+        var c: T = -1.0;
+        var d: T = 0.0;
+
+        const result = f(a, b, c, d);
+        try expect(result.real == -1.0);
+        try expect(result.imag == 0.0);
+    }
+    {
+        var a: T = 1.0;
+        var b: T = 0.0;
+        var c: T = -4.0;
+        var d: T = 0.0;
+
+        const result = f(a, b, c, d);
+        try expect(result.real == -4.0);
+        try expect(result.imag == 0.0);
+    }
+    {
+        // if one operand is an infinity and the other operand is a nonzero finite number or an infinity,
+        // then the result of the * operator is an infinity;
+        var a: T = math.inf(T);
+        var b: T = -math.inf(T);
+        var c: T = 1.0;
+        var d: T = 0.0;
+
+        const result = f(a, b, c, d);
+        try expect(result.real == math.inf(T));
+        try expect(result.imag == -math.inf(T));
+    }
+    {
+        // if one operand is an infinity and the other operand is a nonzero finite number or an infinity,
+        // then the result of the * operator is an infinity;
+        var a: T = math.inf(T);
+        var b: T = -1.0;
+        var c: T = 1.0;
+        var d: T = math.inf(T);
+
+        const result = f(a, b, c, d);
+        try expect(result.real == math.inf(T));
+        try expect(result.imag == math.inf(T));
+    }
+}

lib/compiler_rt/muldc3.zig

@@ -0,0 +1,12 @@
+const common = @import("./common.zig");
+const mulc3 = @import("./mulc3.zig");
+
+pub const panic = common.panic;
+
+comptime {
+    @export(__muldc3, .{ .name = "__muldc3", .linkage = common.linkage });
+}
+
+pub fn __muldc3(a: f64, b: f64, c: f64, d: f64) callconv(.C) mulc3.Complex(f64) {
+    return mulc3.mulc3(f64, a, b, c, d);
+}

lib/compiler_rt/mulhc3.zig

@@ -0,0 +1,12 @@
+const common = @import("./common.zig");
+const mulc3 = @import("./mulc3.zig");
+
+pub const panic = common.panic;
+
+comptime {
+    @export(__mulhc3, .{ .name = "__mulhc3", .linkage = common.linkage });
+}
+
+pub fn __mulhc3(a: f16, b: f16, c: f16, d: f16) callconv(.C) mulc3.Complex(f16) {
+    return mulc3.mulc3(f16, a, b, c, d);
+}

lib/compiler_rt/mulsc3.zig

@@ -0,0 +1,12 @@
+const common = @import("./common.zig");
+const mulc3 = @import("./mulc3.zig");
+
+pub const panic = common.panic;
+
+comptime {
+    @export(__mulsc3, .{ .name = "__mulsc3", .linkage = common.linkage });
+}
+
+pub fn __mulsc3(a: f32, b: f32, c: f32, d: f32) callconv(.C) mulc3.Complex(f32) {
+    return mulc3.mulc3(f32, a, b, c, d);
+}

lib/compiler_rt/multc3.zig

@@ -0,0 +1,12 @@
+const common = @import("./common.zig");
+const mulc3 = @import("./mulc3.zig");
+
+pub const panic = common.panic;
+
+comptime {
+    @export(__multc3, .{ .name = "__multc3", .linkage = common.linkage });
+}
+
+pub fn __multc3(a: f128, b: f128, c: f128, d: f128) callconv(.C) mulc3.Complex(f128) {
+    return mulc3.mulc3(f128, a, b, c, d);
+}

lib/compiler_rt/mulxc3.zig

@@ -0,0 +1,12 @@
+const common = @import("./common.zig");
+const mulc3 = @import("./mulc3.zig");
+
+pub const panic = common.panic;
+
+comptime {
+    @export(__mulxc3, .{ .name = "__mulxc3", .linkage = common.linkage });
+}
+
+pub fn __mulxc3(a: f80, b: f80, c: f80, d: f80) callconv(.C) mulc3.Complex(f80) {
+    return mulc3.mulc3(f80, a, b, c, d);
+}

test/standalone/c_compiler/test.c

@@ -1,4 +1,5 @@
 #include <assert.h>
+#include <complex.h>
 #include <stdio.h>
 #include <stdlib.h>
 
@@ -24,5 +25,30 @@ int main (int argc, char *argv[])
 
   if (!ok) abort();
 
+  // Test some basic arithmetic from compiler-rt
+  {
+    double complex z = 0.0 + I * 4.0;
+    double complex w = 0.0 + I * 16.0;
+    double complex product = z * w;
+    double complex quotient = z / w;
+
+    if (!(creal(product) == -64.0)) abort();
+    if (!(cimag(product) == 0.0)) abort();
+    if (!(creal(quotient) == 0.25)) abort();
+    if (!(cimag(quotient) == 0.0)) abort();
+  }
+
+  {
+    float complex z = 4.0 + I * 4.0;
+    float complex w = 2.0 - I * 2.0;
+    float complex product = z * w;
+    float complex quotient = z / w;
+
+    if (!(creal(product) == 16.0)) abort();
+    if (!(cimag(product) == 0.0)) abort();
+    if (!(creal(quotient) == 0.0)) abort();
+    if (!(cimag(quotient) == 2.0)) abort();
+  }
+
   return EXIT_SUCCESS;
 }