diff --git a/lib/compiler_rt.zig b/lib/compiler_rt.zig index 105c5ed7ad..d261c49ff1 100644 --- a/lib/compiler_rt.zig +++ b/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"); diff --git a/lib/compiler_rt/divc3.zig b/lib/compiler_rt/divc3.zig new file mode 100644 index 0000000000..4e4dba2856 --- /dev/null +++ b/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; +} diff --git a/lib/compiler_rt/divc3_test.zig b/lib/compiler_rt/divc3_test.zig new file mode 100644 index 0000000000..525e2f70fc --- /dev/null +++ b/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)); + } +} diff --git a/lib/compiler_rt/divdc3.zig b/lib/compiler_rt/divdc3.zig new file mode 100644 index 0000000000..f7592f639d --- /dev/null +++ b/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); +} diff --git a/lib/compiler_rt/divhc3.zig b/lib/compiler_rt/divhc3.zig new file mode 100644 index 0000000000..e2d682fe2d --- /dev/null +++ b/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); +} diff --git a/lib/compiler_rt/divsc3.zig b/lib/compiler_rt/divsc3.zig new file mode 100644 index 0000000000..a64f14629c --- /dev/null +++ b/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); +} diff --git a/lib/compiler_rt/divtc3.zig b/lib/compiler_rt/divtc3.zig new file mode 100644 index 0000000000..190df5c067 --- /dev/null +++ b/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); +} diff --git a/lib/compiler_rt/divxc3.zig b/lib/compiler_rt/divxc3.zig new file mode 100644 index 0000000000..32fb269ca5 --- /dev/null +++ b/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); +} diff --git a/lib/compiler_rt/mulc3.zig b/lib/compiler_rt/mulc3.zig new file mode 100644 index 0000000000..0033df70b6 --- /dev/null +++ b/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; +} diff --git a/lib/compiler_rt/mulc3_test.zig b/lib/compiler_rt/mulc3_test.zig new file mode 100644 index 0000000000..b07ecd12ae --- /dev/null +++ b/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)); + } +} diff --git a/lib/compiler_rt/muldc3.zig b/lib/compiler_rt/muldc3.zig new file mode 100644 index 0000000000..343ae5a064 --- /dev/null +++ b/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); +} diff --git a/lib/compiler_rt/mulhc3.zig b/lib/compiler_rt/mulhc3.zig new file mode 100644 index 0000000000..f1fad90aff --- /dev/null +++ b/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); +} diff --git a/lib/compiler_rt/mulsc3.zig b/lib/compiler_rt/mulsc3.zig new file mode 100644 index 0000000000..3ea055f9ff --- /dev/null +++ b/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); +} diff --git a/lib/compiler_rt/multc3.zig b/lib/compiler_rt/multc3.zig new file mode 100644 index 0000000000..c6afcdefb2 --- /dev/null +++ b/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); +} diff --git a/lib/compiler_rt/mulxc3.zig b/lib/compiler_rt/mulxc3.zig new file mode 100644 index 0000000000..decba868e8 --- /dev/null +++ b/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); +} diff --git a/test/standalone/c_compiler/test.c b/test/standalone/c_compiler/test.c index da0a137798..3beb1aa2ff 100644 --- a/test/standalone/c_compiler/test.c +++ b/test/standalone/c_compiler/test.c @@ -1,4 +1,5 @@ #include +#include #include #include @@ -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; }