mirror of
https://github.com/ziglang/zig.git
synced 2026-02-12 20:37:54 +00:00
Merge pull request #13117 from topolarity/compiler-rt-cmul
compiler-rt: Implement complex multiply/division
This commit is contained in:
Andrew Kelley
GitHub
commit
c3d67c5c4e
@ -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");
|
||||
|
||||
62
lib/compiler_rt/divc3.zig
Normal file
62
lib/compiler_rt/divc3.zig
Normal file
@ -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;
|
||||
}
|
||||
77
lib/compiler_rt/divc3_test.zig
Normal file
77
lib/compiler_rt/divc3_test.zig
Normal file
@ -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));
|
||||
}
|
||||
}
|
||||
11
lib/compiler_rt/divdc3.zig
Normal file
11
lib/compiler_rt/divdc3.zig
Normal file
@ -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);
|
||||
}
|
||||
11
lib/compiler_rt/divhc3.zig
Normal file
11
lib/compiler_rt/divhc3.zig
Normal file
@ -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);
|
||||
}
|
||||
11
lib/compiler_rt/divsc3.zig
Normal file
11
lib/compiler_rt/divsc3.zig
Normal file
@ -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);
|
||||
}
|
||||
11
lib/compiler_rt/divtc3.zig
Normal file
11
lib/compiler_rt/divtc3.zig
Normal file
@ -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);
|
||||
}
|
||||
11
lib/compiler_rt/divxc3.zig
Normal file
11
lib/compiler_rt/divxc3.zig
Normal file
@ -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);
|
||||
}
|
||||
@ -7,6 +7,6 @@ comptime {
|
||||
@export(__extenddfxf2, .{ .name = "__extenddfxf2", .linkage = common.linkage });
|
||||
}
|
||||
|
||||
fn __extenddfxf2(a: f64) callconv(.C) f80 {
|
||||
pub fn __extenddfxf2(a: f64) callconv(.C) f80 {
|
||||
return extend_f80(f64, @bitCast(u64, a));
|
||||
}
|
||||
|
||||
@ -92,6 +92,8 @@ pub inline fn extend_f80(comptime src_t: type, a: std.meta.Int(.unsigned, @typeI
|
||||
const src_qnan = 1 << (src_sig_bits - 1);
|
||||
const src_nan_code = src_qnan - 1;
|
||||
|
||||
const SrcShift = std.math.Log2Int(src_rep_t);
|
||||
|
||||
var dst: std.math.F80 = undefined;
|
||||
|
||||
// Break a into a sign and representation of the absolute value
|
||||
@ -124,7 +126,7 @@ pub inline fn extend_f80(comptime src_t: type, a: std.meta.Int(.unsigned, @typeI
|
||||
|
||||
dst.fraction = @as(u64, a_abs) << @intCast(u6, dst_sig_bits - src_sig_bits + scale);
|
||||
dst.fraction |= dst_int_bit; // bit 64 is always set for normal numbers
|
||||
dst.exp = @truncate(u16, a_abs >> @intCast(u4, src_sig_bits - scale));
|
||||
dst.exp = @truncate(u16, a_abs >> @intCast(SrcShift, src_sig_bits - scale));
|
||||
dst.exp ^= 1;
|
||||
dst.exp |= dst_exp_bias - src_exp_bias - scale + 1;
|
||||
} else {
|
||||
|
||||
@ -1,10 +1,27 @@
|
||||
const std = @import("std");
|
||||
const math = std.math;
|
||||
const builtin = @import("builtin");
|
||||
const __extendhfsf2 = @import("extendhfsf2.zig").__extendhfsf2;
|
||||
const __extendhftf2 = @import("extendhftf2.zig").__extendhftf2;
|
||||
const __extendsftf2 = @import("extendsftf2.zig").__extendsftf2;
|
||||
const __extenddftf2 = @import("extenddftf2.zig").__extenddftf2;
|
||||
const __extenddfxf2 = @import("extenddfxf2.zig").__extenddfxf2;
|
||||
const F16T = @import("./common.zig").F16T;
|
||||
|
||||
fn test__extenddfxf2(a: f64, expected: u80) !void {
|
||||
const x = __extenddfxf2(a);
|
||||
|
||||
const rep = @bitCast(u80, x);
|
||||
if (rep == expected)
|
||||
return;
|
||||
|
||||
// test other possible NaN representation(signal NaN)
|
||||
if (math.isNan(@bitCast(f80, expected)) and math.isNan(x))
|
||||
return;
|
||||
|
||||
@panic("__extenddfxf2 test failure");
|
||||
}
|
||||
|
||||
fn test__extenddftf2(a: f64, expected_hi: u64, expected_lo: u64) !void {
|
||||
const x = __extenddftf2(a);
|
||||
|
||||
@ -65,6 +82,33 @@ fn test__extendsftf2(a: f32, expected_hi: u64, expected_lo: u64) !void {
|
||||
return error.TestFailure;
|
||||
}
|
||||
|
||||
test "extenddfxf2" {
|
||||
// qNaN
|
||||
try test__extenddfxf2(makeQNaN64(), 0x7fffc000000000000000);
|
||||
|
||||
// NaN
|
||||
try test__extenddfxf2(makeNaN64(0x7100000000000), 0x7fffe080000000000000);
|
||||
// This is bad?
|
||||
|
||||
// inf
|
||||
try test__extenddfxf2(makeInf64(), 0x7fff8000000000000000);
|
||||
|
||||
// zero
|
||||
try test__extenddfxf2(0.0, 0x0);
|
||||
|
||||
try test__extenddfxf2(0x0.a3456789abcdefp+6, 0x4004a3456789abcdf000);
|
||||
|
||||
try test__extenddfxf2(0x0.edcba987654321fp-8, 0x3ff6edcba98765432000);
|
||||
|
||||
try test__extenddfxf2(0x0.a3456789abcdefp+46, 0x402ca3456789abcdf000);
|
||||
|
||||
try test__extenddfxf2(0x0.edcba987654321fp-44, 0x3fd2edcba98765432000);
|
||||
|
||||
// subnormal
|
||||
try test__extenddfxf2(0x1.8000000000001p-1022, 0x3c01c000000000000800);
|
||||
try test__extenddfxf2(0x1.8000000000002p-1023, 0x3c00c000000000001000);
|
||||
}
|
||||
|
||||
test "extenddftf2" {
|
||||
// qNaN
|
||||
try test__extenddftf2(makeQNaN64(), 0x7fff800000000000, 0x0);
|
||||
@ -85,6 +129,10 @@ test "extenddftf2" {
|
||||
try test__extenddftf2(0x1.23456789abcdefp+45, 0x402c23456789abcd, 0xf000000000000000);
|
||||
|
||||
try test__extenddftf2(0x1.edcba987654321fp-45, 0x3fd2edcba9876543, 0x2000000000000000);
|
||||
|
||||
// subnormal
|
||||
try test__extenddftf2(0x1.8p-1022, 0x3c01800000000000, 0x0);
|
||||
try test__extenddftf2(0x1.8p-1023, 0x3c00800000000000, 0x0);
|
||||
}
|
||||
|
||||
test "extendhfsf2" {
|
||||
|
||||
79
lib/compiler_rt/mulc3.zig
Normal file
79
lib/compiler_rt/mulc3.zig
Normal file
@ -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;
|
||||
}
|
||||
65
lib/compiler_rt/mulc3_test.zig
Normal file
65
lib/compiler_rt/mulc3_test.zig
Normal file
@ -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));
|
||||
}
|
||||
}
|
||||
12
lib/compiler_rt/muldc3.zig
Normal file
12
lib/compiler_rt/muldc3.zig
Normal file
@ -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);
|
||||
}
|
||||
12
lib/compiler_rt/mulhc3.zig
Normal file
12
lib/compiler_rt/mulhc3.zig
Normal file
@ -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);
|
||||
}
|
||||
12
lib/compiler_rt/mulsc3.zig
Normal file
12
lib/compiler_rt/mulsc3.zig
Normal file
@ -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);
|
||||
}
|
||||
12
lib/compiler_rt/multc3.zig
Normal file
12
lib/compiler_rt/multc3.zig
Normal file
@ -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);
|
||||
}
|
||||
12
lib/compiler_rt/mulxc3.zig
Normal file
12
lib/compiler_rt/mulxc3.zig
Normal file
@ -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);
|
||||
}
|
||||
@ -15,175 +15,157 @@ const minInt = std.math.minInt;
|
||||
///
|
||||
/// Special Cases:
|
||||
/// - ilogb(+-inf) = maxInt(i32)
|
||||
/// - ilogb(0) = maxInt(i32)
|
||||
/// - ilogb(nan) = maxInt(i32)
|
||||
/// - ilogb(+-0) = minInt(i32)
|
||||
/// - ilogb(nan) = minInt(i32)
|
||||
pub fn ilogb(x: anytype) i32 {
|
||||
const T = @TypeOf(x);
|
||||
return switch (T) {
|
||||
f32 => ilogb32(x),
|
||||
f64 => ilogb64(x),
|
||||
f128 => ilogb128(x),
|
||||
else => @compileError("ilogb not implemented for " ++ @typeName(T)),
|
||||
};
|
||||
return ilogbX(T, x);
|
||||
}
|
||||
|
||||
// TODO: unify these implementations with generics
|
||||
pub const fp_ilogbnan = minInt(i32);
|
||||
pub const fp_ilogb0 = minInt(i32);
|
||||
|
||||
// NOTE: Should these be exposed publicly?
|
||||
const fp_ilogbnan = -1 - @as(i32, maxInt(u32) >> 1);
|
||||
const fp_ilogb0 = fp_ilogbnan;
|
||||
fn ilogbX(comptime T: type, x: T) i32 {
|
||||
const typeWidth = @typeInfo(T).Float.bits;
|
||||
const significandBits = math.floatMantissaBits(T);
|
||||
const exponentBits = math.floatExponentBits(T);
|
||||
|
||||
fn ilogb32(x: f32) i32 {
|
||||
var u = @bitCast(u32, x);
|
||||
var e = @intCast(i32, (u >> 23) & 0xFF);
|
||||
const Z = std.meta.Int(.unsigned, typeWidth);
|
||||
|
||||
// TODO: We should be able to merge this with the lower check.
|
||||
if (math.isNan(x)) {
|
||||
return maxInt(i32);
|
||||
}
|
||||
const signBit = (@as(Z, 1) << (significandBits + exponentBits));
|
||||
const maxExponent = ((1 << exponentBits) - 1);
|
||||
const exponentBias = (maxExponent >> 1);
|
||||
|
||||
const absMask = signBit - 1;
|
||||
|
||||
var u = @bitCast(Z, x) & absMask;
|
||||
var e = @intCast(i32, u >> significandBits);
|
||||
|
||||
if (e == 0) {
|
||||
u <<= 9;
|
||||
if (u == 0) {
|
||||
math.raiseInvalid();
|
||||
return fp_ilogb0;
|
||||
}
|
||||
|
||||
// subnormal
|
||||
e = -0x7F;
|
||||
while (u >> 31 == 0) : (u <<= 1) {
|
||||
e -= 1;
|
||||
}
|
||||
return e;
|
||||
// offset sign bit, exponent bits, and integer bit (if present) + bias
|
||||
const offset = 1 + exponentBits + @boolToInt(T == f80) - exponentBias;
|
||||
return offset - @intCast(i32, @clz(u));
|
||||
}
|
||||
|
||||
if (e == 0xFF) {
|
||||
if (e == maxExponent) {
|
||||
math.raiseInvalid();
|
||||
if (u << 9 != 0) {
|
||||
return fp_ilogbnan;
|
||||
} else {
|
||||
return maxInt(i32);
|
||||
}
|
||||
if (u > @bitCast(Z, math.inf(T))) {
|
||||
return fp_ilogbnan; // u is a NaN
|
||||
} else return maxInt(i32);
|
||||
}
|
||||
|
||||
return e - 0x7F;
|
||||
}
|
||||
|
||||
fn ilogb64(x: f64) i32 {
|
||||
var u = @bitCast(u64, x);
|
||||
var e = @intCast(i32, (u >> 52) & 0x7FF);
|
||||
|
||||
if (math.isNan(x)) {
|
||||
return maxInt(i32);
|
||||
}
|
||||
|
||||
if (e == 0) {
|
||||
u <<= 12;
|
||||
if (u == 0) {
|
||||
math.raiseInvalid();
|
||||
return fp_ilogb0;
|
||||
}
|
||||
|
||||
// subnormal
|
||||
e = -0x3FF;
|
||||
while (u >> 63 == 0) : (u <<= 1) {
|
||||
e -= 1;
|
||||
}
|
||||
return e;
|
||||
}
|
||||
|
||||
if (e == 0x7FF) {
|
||||
math.raiseInvalid();
|
||||
if (u << 12 != 0) {
|
||||
return fp_ilogbnan;
|
||||
} else {
|
||||
return maxInt(i32);
|
||||
}
|
||||
}
|
||||
|
||||
return e - 0x3FF;
|
||||
}
|
||||
|
||||
fn ilogb128(x: f128) i32 {
|
||||
var u = @bitCast(u128, x);
|
||||
var e = @intCast(i32, (u >> 112) & 0x7FFF);
|
||||
|
||||
if (math.isNan(x)) {
|
||||
return maxInt(i32);
|
||||
}
|
||||
|
||||
if (e == 0) {
|
||||
u <<= 16;
|
||||
if (u == 0) {
|
||||
math.raiseInvalid();
|
||||
return fp_ilogb0;
|
||||
}
|
||||
|
||||
// subnormal x
|
||||
return ilogb128(x * 0x1p120) - 120;
|
||||
}
|
||||
|
||||
if (e == 0x7FFF) {
|
||||
math.raiseInvalid();
|
||||
if (u << 16 != 0) {
|
||||
return fp_ilogbnan;
|
||||
} else {
|
||||
return maxInt(i32);
|
||||
}
|
||||
}
|
||||
|
||||
return e - 0x3FFF;
|
||||
return e - exponentBias;
|
||||
}
|
||||
|
||||
test "type dispatch" {
|
||||
try expect(ilogb(@as(f32, 0.2)) == ilogb32(0.2));
|
||||
try expect(ilogb(@as(f64, 0.2)) == ilogb64(0.2));
|
||||
try expect(ilogb(@as(f32, 0.2)) == ilogbX(f32, 0.2));
|
||||
try expect(ilogb(@as(f64, 0.2)) == ilogbX(f64, 0.2));
|
||||
}
|
||||
|
||||
test "16" {
|
||||
try expect(ilogbX(f16, 0.0) == fp_ilogb0);
|
||||
try expect(ilogbX(f16, 0.5) == -1);
|
||||
try expect(ilogbX(f16, 0.8923) == -1);
|
||||
try expect(ilogbX(f16, 10.0) == 3);
|
||||
try expect(ilogbX(f16, -65504) == 15);
|
||||
try expect(ilogbX(f16, 2398.23) == 11);
|
||||
|
||||
try expect(ilogbX(f16, 0x1p-1) == -1);
|
||||
try expect(ilogbX(f16, 0x1p-17) == -17);
|
||||
try expect(ilogbX(f16, 0x1p-24) == -24);
|
||||
}
|
||||
|
||||
test "32" {
|
||||
try expect(ilogb32(0.0) == fp_ilogb0);
|
||||
try expect(ilogb32(0.5) == -1);
|
||||
try expect(ilogb32(0.8923) == -1);
|
||||
try expect(ilogb32(10.0) == 3);
|
||||
try expect(ilogb32(-123984) == 16);
|
||||
try expect(ilogb32(2398.23) == 11);
|
||||
try expect(ilogbX(f32, 0.0) == fp_ilogb0);
|
||||
try expect(ilogbX(f32, 0.5) == -1);
|
||||
try expect(ilogbX(f32, 0.8923) == -1);
|
||||
try expect(ilogbX(f32, 10.0) == 3);
|
||||
try expect(ilogbX(f32, -123984) == 16);
|
||||
try expect(ilogbX(f32, 2398.23) == 11);
|
||||
|
||||
try expect(ilogbX(f32, 0x1p-1) == -1);
|
||||
try expect(ilogbX(f32, 0x1p-122) == -122);
|
||||
try expect(ilogbX(f32, 0x1p-127) == -127);
|
||||
}
|
||||
|
||||
test "64" {
|
||||
try expect(ilogb64(0.0) == fp_ilogb0);
|
||||
try expect(ilogb64(0.5) == -1);
|
||||
try expect(ilogb64(0.8923) == -1);
|
||||
try expect(ilogb64(10.0) == 3);
|
||||
try expect(ilogb64(-123984) == 16);
|
||||
try expect(ilogb64(2398.23) == 11);
|
||||
try expect(ilogbX(f64, 0.0) == fp_ilogb0);
|
||||
try expect(ilogbX(f64, 0.5) == -1);
|
||||
try expect(ilogbX(f64, 0.8923) == -1);
|
||||
try expect(ilogbX(f64, 10.0) == 3);
|
||||
try expect(ilogbX(f64, -123984) == 16);
|
||||
try expect(ilogbX(f64, 2398.23) == 11);
|
||||
|
||||
try expect(ilogbX(f64, 0x1p-1) == -1);
|
||||
try expect(ilogbX(f64, 0x1p-127) == -127);
|
||||
try expect(ilogbX(f64, 0x1p-1012) == -1012);
|
||||
try expect(ilogbX(f64, 0x1p-1023) == -1023);
|
||||
}
|
||||
|
||||
test "80" {
|
||||
try expect(ilogbX(f80, 0.0) == fp_ilogb0);
|
||||
try expect(ilogbX(f80, 0.5) == -1);
|
||||
try expect(ilogbX(f80, 0.8923) == -1);
|
||||
try expect(ilogbX(f80, 10.0) == 3);
|
||||
try expect(ilogbX(f80, -123984) == 16);
|
||||
try expect(ilogbX(f80, 2398.23) == 11);
|
||||
|
||||
try expect(ilogbX(f80, 0x1p-1) == -1);
|
||||
try expect(ilogbX(f80, 0x1p-127) == -127);
|
||||
try expect(ilogbX(f80, 0x1p-1023) == -1023);
|
||||
try expect(ilogbX(f80, 0x1p-16383) == -16383);
|
||||
}
|
||||
|
||||
test "128" {
|
||||
try expect(ilogb128(0.0) == fp_ilogb0);
|
||||
try expect(ilogb128(0.5) == -1);
|
||||
try expect(ilogb128(0.8923) == -1);
|
||||
try expect(ilogb128(10.0) == 3);
|
||||
try expect(ilogb128(-123984) == 16);
|
||||
try expect(ilogb128(2398.23) == 11);
|
||||
try expect(ilogbX(f128, 0.0) == fp_ilogb0);
|
||||
try expect(ilogbX(f128, 0.5) == -1);
|
||||
try expect(ilogbX(f128, 0.8923) == -1);
|
||||
try expect(ilogbX(f128, 10.0) == 3);
|
||||
try expect(ilogbX(f128, -123984) == 16);
|
||||
try expect(ilogbX(f128, 2398.23) == 11);
|
||||
|
||||
try expect(ilogbX(f128, 0x1p-1) == -1);
|
||||
try expect(ilogbX(f128, 0x1p-127) == -127);
|
||||
try expect(ilogbX(f128, 0x1p-1023) == -1023);
|
||||
try expect(ilogbX(f128, 0x1p-16383) == -16383);
|
||||
}
|
||||
|
||||
test "16 special" {
|
||||
try expect(ilogbX(f16, math.inf(f16)) == maxInt(i32));
|
||||
try expect(ilogbX(f16, -math.inf(f16)) == maxInt(i32));
|
||||
try expect(ilogbX(f16, 0.0) == minInt(i32));
|
||||
try expect(ilogbX(f16, math.nan(f16)) == fp_ilogbnan);
|
||||
}
|
||||
|
||||
test "32 special" {
|
||||
try expect(ilogb32(math.inf(f32)) == maxInt(i32));
|
||||
try expect(ilogb32(-math.inf(f32)) == maxInt(i32));
|
||||
try expect(ilogb32(0.0) == minInt(i32));
|
||||
try expect(ilogb32(math.nan(f32)) == maxInt(i32));
|
||||
try expect(ilogbX(f32, math.inf(f32)) == maxInt(i32));
|
||||
try expect(ilogbX(f32, -math.inf(f32)) == maxInt(i32));
|
||||
try expect(ilogbX(f32, 0.0) == minInt(i32));
|
||||
try expect(ilogbX(f32, math.nan(f32)) == fp_ilogbnan);
|
||||
}
|
||||
|
||||
test "64 special" {
|
||||
try expect(ilogb64(math.inf(f64)) == maxInt(i32));
|
||||
try expect(ilogb64(-math.inf(f64)) == maxInt(i32));
|
||||
try expect(ilogb64(0.0) == minInt(i32));
|
||||
try expect(ilogb64(math.nan(f64)) == maxInt(i32));
|
||||
try expect(ilogbX(f64, math.inf(f64)) == maxInt(i32));
|
||||
try expect(ilogbX(f64, -math.inf(f64)) == maxInt(i32));
|
||||
try expect(ilogbX(f64, 0.0) == minInt(i32));
|
||||
try expect(ilogbX(f64, math.nan(f64)) == fp_ilogbnan);
|
||||
}
|
||||
|
||||
test "80 special" {
|
||||
try expect(ilogbX(f80, math.inf(f80)) == maxInt(i32));
|
||||
try expect(ilogbX(f80, -math.inf(f80)) == maxInt(i32));
|
||||
try expect(ilogbX(f80, 0.0) == minInt(i32));
|
||||
try expect(ilogbX(f80, math.nan(f80)) == fp_ilogbnan);
|
||||
}
|
||||
|
||||
test "128 special" {
|
||||
try expect(ilogb128(math.inf(f128)) == maxInt(i32));
|
||||
try expect(ilogb128(-math.inf(f128)) == maxInt(i32));
|
||||
try expect(ilogb128(0.0) == minInt(i32));
|
||||
try expect(ilogb128(math.nan(f128)) == maxInt(i32));
|
||||
try expect(ilogbX(f128, math.inf(f128)) == maxInt(i32));
|
||||
try expect(ilogbX(f128, -math.inf(f128)) == maxInt(i32));
|
||||
try expect(ilogbX(f128, 0.0) == minInt(i32));
|
||||
try expect(ilogbX(f128, math.nan(f128)) == fp_ilogbnan);
|
||||
}
|
||||
|
||||
@ -1,91 +1,148 @@
|
||||
// Ported from musl, which is licensed under the MIT license:
|
||||
// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
|
||||
//
|
||||
// https://git.musl-libc.org/cgit/musl/tree/src/math/ldexpf.c
|
||||
// https://git.musl-libc.org/cgit/musl/tree/src/math/ldexp.c
|
||||
|
||||
const std = @import("std");
|
||||
const math = std.math;
|
||||
const Log2Int = std.math.Log2Int;
|
||||
const assert = std.debug.assert;
|
||||
const expect = std.testing.expect;
|
||||
|
||||
/// Returns x * 2^n.
|
||||
pub fn ldexp(x: anytype, n: i32) @TypeOf(x) {
|
||||
var base = x;
|
||||
var shift = n;
|
||||
|
||||
const T = @TypeOf(base);
|
||||
const T = @TypeOf(x);
|
||||
const TBits = std.meta.Int(.unsigned, @typeInfo(T).Float.bits);
|
||||
|
||||
const exponent_bits = math.floatExponentBits(T);
|
||||
const mantissa_bits = math.floatMantissaBits(T);
|
||||
const exponent_min = math.floatExponentMin(T);
|
||||
const exponent_max = math.floatExponentMax(T);
|
||||
const fractional_bits = math.floatFractionalBits(T);
|
||||
|
||||
const exponent_bias = exponent_max;
|
||||
const max_biased_exponent = 2 * math.floatExponentMax(T);
|
||||
const mantissa_mask = @as(TBits, (1 << mantissa_bits) - 1);
|
||||
|
||||
// fix double rounding errors in subnormal ranges
|
||||
// https://git.musl-libc.org/cgit/musl/commit/src/math/ldexp.c?id=8c44a060243f04283ca68dad199aab90336141db
|
||||
const scale_min_expo = exponent_min + mantissa_bits + 1;
|
||||
const scale_min = @bitCast(T, @as(TBits, scale_min_expo + exponent_bias) << mantissa_bits);
|
||||
const scale_max = @bitCast(T, @intCast(TBits, exponent_max + exponent_bias) << mantissa_bits);
|
||||
const repr = @bitCast(TBits, x);
|
||||
const sign_bit = repr & (1 << (exponent_bits + mantissa_bits));
|
||||
|
||||
// scale `shift` within floating point limits, if possible
|
||||
// second pass is possible due to subnormal range
|
||||
// third pass always results in +/-0.0 or +/-inf
|
||||
if (shift > exponent_max) {
|
||||
base *= scale_max;
|
||||
shift -= exponent_max;
|
||||
if (shift > exponent_max) {
|
||||
base *= scale_max;
|
||||
shift -= exponent_max;
|
||||
if (shift > exponent_max) shift = exponent_max;
|
||||
if (math.isNan(x) or !math.isFinite(x))
|
||||
return x;
|
||||
|
||||
var exponent: i32 = @intCast(i32, (repr << 1) >> (mantissa_bits + 1));
|
||||
if (exponent == 0)
|
||||
exponent += (@as(i32, exponent_bits) + @boolToInt(T == f80)) - @clz(repr << 1);
|
||||
|
||||
if (n >= 0) {
|
||||
if (n > max_biased_exponent - exponent) {
|
||||
// Overflow. Return +/- inf
|
||||
return @bitCast(T, @bitCast(TBits, math.inf(T)) | sign_bit);
|
||||
} else if (exponent + n <= 0) {
|
||||
// Result is subnormal
|
||||
return @bitCast(T, (repr << @intCast(Log2Int(TBits), n)) | sign_bit);
|
||||
} else if (exponent <= 0) {
|
||||
// Result is normal, but needs shifting
|
||||
var result = @intCast(TBits, n + exponent) << mantissa_bits;
|
||||
result |= (repr << @intCast(Log2Int(TBits), 1 - exponent)) & mantissa_mask;
|
||||
return @bitCast(T, result | sign_bit);
|
||||
}
|
||||
} else if (shift < exponent_min) {
|
||||
base *= scale_min;
|
||||
shift -= scale_min_expo;
|
||||
if (shift < exponent_min) {
|
||||
base *= scale_min;
|
||||
shift -= scale_min_expo;
|
||||
if (shift < exponent_min) shift = exponent_min;
|
||||
|
||||
// Result needs no shifting
|
||||
return @bitCast(T, repr + (@intCast(TBits, n) << mantissa_bits));
|
||||
} else {
|
||||
if (n <= -exponent) {
|
||||
if (n < -(mantissa_bits + exponent))
|
||||
return @bitCast(T, sign_bit); // Severe underflow. Return +/- 0
|
||||
|
||||
// Result underflowed, we need to shift and round
|
||||
const shift = @intCast(Log2Int(TBits), math.min(-n, -(exponent + n) + 1));
|
||||
const exact_tie: bool = @ctz(repr) == shift - 1;
|
||||
var result = repr & mantissa_mask;
|
||||
|
||||
if (T != f80) // Include integer bit
|
||||
result |= @as(TBits, @boolToInt(exponent > 0)) << fractional_bits;
|
||||
result = @intCast(TBits, (result >> (shift - 1)));
|
||||
|
||||
// Round result, including round-to-even for exact ties
|
||||
result = ((result + 1) >> 1) & ~@as(TBits, @boolToInt(exact_tie));
|
||||
return @bitCast(T, result | sign_bit);
|
||||
}
|
||||
|
||||
// Result is exact, and needs no shifting
|
||||
return @bitCast(T, repr - (@intCast(TBits, -n) << mantissa_bits));
|
||||
}
|
||||
|
||||
return base * @bitCast(T, @intCast(TBits, shift + exponent_bias) << mantissa_bits);
|
||||
}
|
||||
|
||||
test "math.ldexp" {
|
||||
// TODO derive the various constants here with new maths API
|
||||
|
||||
// basic usage
|
||||
try expect(ldexp(@as(f16, 1.5), 4) == 24.0);
|
||||
try expect(ldexp(@as(f32, 1.5), 4) == 24.0);
|
||||
try expect(ldexp(@as(f64, 1.5), 4) == 24.0);
|
||||
try expect(ldexp(@as(f128, 1.5), 4) == 24.0);
|
||||
|
||||
// subnormals
|
||||
try expect(math.isNormal(ldexp(@as(f16, 1.0), -14)));
|
||||
try expect(!math.isNormal(ldexp(@as(f16, 1.0), -15)));
|
||||
try expect(math.isNormal(ldexp(@as(f32, 1.0), -126)));
|
||||
try expect(!math.isNormal(ldexp(@as(f32, 1.0), -127)));
|
||||
try expect(math.isNormal(ldexp(@as(f64, 1.0), -1022)));
|
||||
try expect(!math.isNormal(ldexp(@as(f64, 1.0), -1023)));
|
||||
try expect(math.isNormal(ldexp(@as(f128, 1.0), -16382)));
|
||||
try expect(!math.isNormal(ldexp(@as(f128, 1.0), -16383)));
|
||||
// unreliable due to lack of native f16 support, see talk on PR #8733
|
||||
// try expect(ldexp(@as(f16, 0x1.1FFp-1), -14 - 9) == math.floatTrueMin(f16));
|
||||
try expect(ldexp(@as(f16, 0x1.1FFp14), -14 - 9 - 15) == math.floatTrueMin(f16));
|
||||
try expect(ldexp(@as(f32, 0x1.3FFFFFp-1), -126 - 22) == math.floatTrueMin(f32));
|
||||
try expect(ldexp(@as(f64, 0x1.7FFFFFFFFFFFFp-1), -1022 - 51) == math.floatTrueMin(f64));
|
||||
try expect(ldexp(@as(f80, 0x1.7FFFFFFFFFFFFFFEp-1), -16382 - 62) == math.floatTrueMin(f80));
|
||||
try expect(ldexp(@as(f128, 0x1.7FFFFFFFFFFFFFFFFFFFFFFFFFFFp-1), -16382 - 111) == math.floatTrueMin(f128));
|
||||
|
||||
// float limits
|
||||
try expect(ldexp(math.floatMax(f32), -128 - 149) > 0.0);
|
||||
try expect(ldexp(math.floatMax(f32), -128 - 149 - 1) == 0.0);
|
||||
try expect(!math.isPositiveInf(ldexp(math.floatTrueMin(f16), 15 + 24)));
|
||||
try expect(math.isPositiveInf(ldexp(math.floatTrueMin(f16), 15 + 24 + 1)));
|
||||
try expect(!math.isPositiveInf(ldexp(math.floatTrueMin(f32), 127 + 149)));
|
||||
try expect(math.isPositiveInf(ldexp(math.floatTrueMin(f32), 127 + 149 + 1)));
|
||||
try expect(!math.isPositiveInf(ldexp(math.floatTrueMin(f64), 1023 + 1074)));
|
||||
try expect(math.isPositiveInf(ldexp(math.floatTrueMin(f64), 1023 + 1074 + 1)));
|
||||
try expect(!math.isPositiveInf(ldexp(math.floatTrueMin(f128), 16383 + 16494)));
|
||||
try expect(math.isPositiveInf(ldexp(math.floatTrueMin(f128), 16383 + 16494 + 1)));
|
||||
|
||||
@setEvalBranchQuota(10_000);
|
||||
|
||||
inline for ([_]type{ f16, f32, f64, f80, f128 }) |T| {
|
||||
const fractional_bits = math.floatFractionalBits(T);
|
||||
|
||||
const min_exponent = math.floatExponentMin(T);
|
||||
const max_exponent = math.floatExponentMax(T);
|
||||
const exponent_bias = max_exponent;
|
||||
|
||||
// basic usage
|
||||
try expect(ldexp(@as(T, 1.5), 4) == 24.0);
|
||||
|
||||
// normals -> subnormals
|
||||
try expect(math.isNormal(ldexp(@as(T, 1.0), min_exponent)));
|
||||
try expect(!math.isNormal(ldexp(@as(T, 1.0), min_exponent - 1)));
|
||||
|
||||
// normals -> zero
|
||||
try expect(ldexp(@as(T, 1.0), min_exponent - fractional_bits) > 0.0);
|
||||
try expect(ldexp(@as(T, 1.0), min_exponent - fractional_bits - 1) == 0.0);
|
||||
|
||||
// subnormals -> zero
|
||||
try expect(ldexp(math.floatTrueMin(T), 0) > 0.0);
|
||||
try expect(ldexp(math.floatTrueMin(T), -1) == 0.0);
|
||||
|
||||
// Multiplications might flush the denormals to zero, esp. at
|
||||
// runtime, so we manually construct the constants here instead.
|
||||
const Z = std.meta.Int(.unsigned, @bitSizeOf(T));
|
||||
const EightTimesTrueMin = @bitCast(T, @as(Z, 8));
|
||||
const TwoTimesTrueMin = @bitCast(T, @as(Z, 2));
|
||||
|
||||
// subnormals -> subnormals
|
||||
try expect(ldexp(math.floatTrueMin(T), 3) == EightTimesTrueMin);
|
||||
try expect(ldexp(EightTimesTrueMin, -2) == TwoTimesTrueMin);
|
||||
try expect(ldexp(EightTimesTrueMin, -3) == math.floatTrueMin(T));
|
||||
|
||||
// subnormals -> normals (+)
|
||||
try expect(ldexp(math.floatTrueMin(T), fractional_bits) == math.floatMin(T));
|
||||
try expect(ldexp(math.floatTrueMin(T), fractional_bits - 1) == math.floatMin(T) * 0.5);
|
||||
|
||||
// subnormals -> normals (-)
|
||||
try expect(ldexp(-math.floatTrueMin(T), fractional_bits) == -math.floatMin(T));
|
||||
try expect(ldexp(-math.floatTrueMin(T), fractional_bits - 1) == -math.floatMin(T) * 0.5);
|
||||
|
||||
// subnormals -> float limits (+inf)
|
||||
try expect(math.isFinite(ldexp(math.floatTrueMin(T), max_exponent + exponent_bias + fractional_bits - 1)));
|
||||
try expect(ldexp(math.floatTrueMin(T), max_exponent + exponent_bias + fractional_bits) == math.inf(T));
|
||||
|
||||
// subnormals -> float limits (-inf)
|
||||
try expect(math.isFinite(ldexp(-math.floatTrueMin(T), max_exponent + exponent_bias + fractional_bits - 1)));
|
||||
try expect(ldexp(-math.floatTrueMin(T), max_exponent + exponent_bias + fractional_bits) == -math.inf(T));
|
||||
|
||||
// infinity -> infinity
|
||||
try expect(ldexp(math.inf(T), math.maxInt(i32)) == math.inf(T));
|
||||
try expect(ldexp(math.inf(T), math.minInt(i32)) == math.inf(T));
|
||||
try expect(ldexp(math.inf(T), max_exponent) == math.inf(T));
|
||||
try expect(ldexp(math.inf(T), min_exponent) == math.inf(T));
|
||||
try expect(ldexp(-math.inf(T), math.maxInt(i32)) == -math.inf(T));
|
||||
try expect(ldexp(-math.inf(T), math.minInt(i32)) == -math.inf(T));
|
||||
|
||||
// extremely large n
|
||||
try expect(ldexp(math.floatMax(T), math.maxInt(i32)) == math.inf(T));
|
||||
try expect(ldexp(math.floatMax(T), -math.maxInt(i32)) == 0.0);
|
||||
try expect(ldexp(math.floatMax(T), math.minInt(i32)) == 0.0);
|
||||
try expect(ldexp(math.floatTrueMin(T), math.maxInt(i32)) == math.inf(T));
|
||||
try expect(ldexp(math.floatTrueMin(T), -math.maxInt(i32)) == 0.0);
|
||||
try expect(ldexp(math.floatTrueMin(T), math.minInt(i32)) == 0.0);
|
||||
}
|
||||
}
|
||||
|
||||
@ -2,6 +2,7 @@
|
||||
#include <stdlib.h>
|
||||
#include <stdbool.h>
|
||||
#include <string.h>
|
||||
#include <complex.h>
|
||||
|
||||
void zig_panic();
|
||||
|
||||
@ -50,6 +51,13 @@ void zig_ptr(void *);
|
||||
|
||||
void zig_bool(bool);
|
||||
|
||||
// Note: These two functions match the signature of __mulsc3 and __muldc3 in compiler-rt (and libgcc)
|
||||
float complex zig_cmultf_comp(float a_r, float a_i, float b_r, float b_i);
|
||||
double complex zig_cmultd_comp(double a_r, double a_i, double b_r, double b_i);
|
||||
|
||||
float complex zig_cmultf(float complex a, float complex b);
|
||||
double complex zig_cmultd(double complex a, double complex b);
|
||||
|
||||
struct BigStruct {
|
||||
uint64_t a;
|
||||
uint64_t b;
|
||||
@ -167,6 +175,43 @@ void run_c_tests(void) {
|
||||
|
||||
zig_bool(true);
|
||||
|
||||
// TODO: Resolve https://github.com/ziglang/zig/issues/8465
|
||||
//{
|
||||
// float complex a = 1.25f + I * 2.6f;
|
||||
// float complex b = 11.3f - I * 1.5f;
|
||||
// float complex z = zig_cmultf(a, b);
|
||||
// assert_or_panic(creal(z) == 1.5f);
|
||||
// assert_or_panic(cimag(z) == 13.5f);
|
||||
//}
|
||||
|
||||
{
|
||||
double complex a = 1.25 + I * 2.6;
|
||||
double complex b = 11.3 - I * 1.5;
|
||||
double complex z = zig_cmultd(a, b);
|
||||
assert_or_panic(creal(z) == 1.5);
|
||||
assert_or_panic(cimag(z) == 13.5);
|
||||
}
|
||||
|
||||
{
|
||||
float a_r = 1.25f;
|
||||
float a_i = 2.6f;
|
||||
float b_r = 11.3f;
|
||||
float b_i = -1.5f;
|
||||
float complex z = zig_cmultf_comp(a_r, a_i, b_r, b_i);
|
||||
assert_or_panic(creal(z) == 1.5f);
|
||||
assert_or_panic(cimag(z) == 13.5f);
|
||||
}
|
||||
|
||||
{
|
||||
double a_r = 1.25;
|
||||
double a_i = 2.6;
|
||||
double b_r = 11.3;
|
||||
double b_i = -1.5;
|
||||
double complex z = zig_cmultd_comp(a_r, a_i, b_r, b_i);
|
||||
assert_or_panic(creal(z) == 1.5);
|
||||
assert_or_panic(cimag(z) == 13.5);
|
||||
}
|
||||
|
||||
{
|
||||
struct BigStruct s = {1, 2, 3, 4, 5};
|
||||
zig_big_struct(s);
|
||||
@ -321,6 +366,42 @@ void c_five_floats(float a, float b, float c, float d, float e) {
|
||||
assert_or_panic(e == 5.0);
|
||||
}
|
||||
|
||||
float complex c_cmultf_comp(float a_r, float a_i, float b_r, float b_i) {
|
||||
assert_or_panic(a_r == 1.25f);
|
||||
assert_or_panic(a_i == 2.6f);
|
||||
assert_or_panic(b_r == 11.3f);
|
||||
assert_or_panic(b_i == -1.5f);
|
||||
|
||||
return 1.5f + I * 13.5f;
|
||||
}
|
||||
|
||||
double complex c_cmultd_comp(double a_r, double a_i, double b_r, double b_i) {
|
||||
assert_or_panic(a_r == 1.25);
|
||||
assert_or_panic(a_i == 2.6);
|
||||
assert_or_panic(b_r == 11.3);
|
||||
assert_or_panic(b_i == -1.5);
|
||||
|
||||
return 1.5 + I * 13.5;
|
||||
}
|
||||
|
||||
float complex c_cmultf(float complex a, float complex b) {
|
||||
assert_or_panic(creal(a) == 1.25f);
|
||||
assert_or_panic(cimag(a) == 2.6f);
|
||||
assert_or_panic(creal(b) == 11.3f);
|
||||
assert_or_panic(cimag(b) == -1.5f);
|
||||
|
||||
return 1.5f + I * 13.5f;
|
||||
}
|
||||
|
||||
double complex c_cmultd(double complex a, double complex b) {
|
||||
assert_or_panic(creal(a) == 1.25);
|
||||
assert_or_panic(cimag(a) == 2.6);
|
||||
assert_or_panic(creal(b) == 11.3);
|
||||
assert_or_panic(cimag(b) == -1.5);
|
||||
|
||||
return 1.5 + I * 13.5;
|
||||
}
|
||||
|
||||
void c_big_struct(struct BigStruct x) {
|
||||
assert_or_panic(x.a == 1);
|
||||
assert_or_panic(x.b == 2);
|
||||
|
||||
@ -145,6 +145,101 @@ export fn zig_bool(x: bool) void {
|
||||
expect(x) catch @panic("test failure: zig_bool");
|
||||
}
|
||||
|
||||
// TODO: Replace these with the correct types once we resolve
|
||||
// https://github.com/ziglang/zig/issues/8465
|
||||
//
|
||||
// For now, we have no way of referring to the _Complex C types from Zig,
|
||||
// so our ABI is unavoidably broken on some platforms (such as i386)
|
||||
const ComplexFloat = extern struct {
|
||||
real: f32,
|
||||
imag: f32,
|
||||
};
|
||||
const ComplexDouble = extern struct {
|
||||
real: f64,
|
||||
imag: f64,
|
||||
};
|
||||
|
||||
// Note: These two functions match the signature of __mulsc3 and __muldc3 in compiler-rt (and libgcc)
|
||||
extern fn c_cmultf_comp(a_r: f32, a_i: f32, b_r: f32, b_i: f32) ComplexFloat;
|
||||
extern fn c_cmultd_comp(a_r: f64, a_i: f64, b_r: f64, b_i: f64) ComplexDouble;
|
||||
|
||||
extern fn c_cmultf(a: ComplexFloat, b: ComplexFloat) ComplexFloat;
|
||||
extern fn c_cmultd(a: ComplexDouble, b: ComplexDouble) ComplexDouble;
|
||||
|
||||
test "C ABI complex float" {
|
||||
if (true) return error.SkipZigTest; // See https://github.com/ziglang/zig/issues/8465
|
||||
|
||||
const a = ComplexFloat{ .real = 1.25, .imag = 2.6 };
|
||||
const b = ComplexFloat{ .real = 11.3, .imag = -1.5 };
|
||||
|
||||
const z = c_cmultf(a, b);
|
||||
expect(z.real == 1.5) catch @panic("test failure: zig_complex_float 1");
|
||||
expect(z.imag == 13.5) catch @panic("test failure: zig_complex_float 2");
|
||||
}
|
||||
|
||||
test "C ABI complex float by component" {
|
||||
const a = ComplexFloat{ .real = 1.25, .imag = 2.6 };
|
||||
const b = ComplexFloat{ .real = 11.3, .imag = -1.5 };
|
||||
|
||||
const z2 = c_cmultf_comp(a.real, a.imag, b.real, b.imag);
|
||||
expect(z2.real == 1.5) catch @panic("test failure: zig_complex_float 3");
|
||||
expect(z2.imag == 13.5) catch @panic("test failure: zig_complex_float 4");
|
||||
}
|
||||
|
||||
test "C ABI complex double" {
|
||||
const a = ComplexDouble{ .real = 1.25, .imag = 2.6 };
|
||||
const b = ComplexDouble{ .real = 11.3, .imag = -1.5 };
|
||||
|
||||
const z = c_cmultd(a, b);
|
||||
expect(z.real == 1.5) catch @panic("test failure: zig_complex_double 1");
|
||||
expect(z.imag == 13.5) catch @panic("test failure: zig_complex_double 2");
|
||||
}
|
||||
|
||||
test "C ABI complex double by component" {
|
||||
const a = ComplexDouble{ .real = 1.25, .imag = 2.6 };
|
||||
const b = ComplexDouble{ .real = 11.3, .imag = -1.5 };
|
||||
|
||||
const z = c_cmultd_comp(a.real, a.imag, b.real, b.imag);
|
||||
expect(z.real == 1.5) catch @panic("test failure: zig_complex_double 3");
|
||||
expect(z.imag == 13.5) catch @panic("test failure: zig_complex_double 4");
|
||||
}
|
||||
|
||||
export fn zig_cmultf(a: ComplexFloat, b: ComplexFloat) ComplexFloat {
|
||||
expect(a.real == 1.25) catch @panic("test failure: zig_cmultf 1");
|
||||
expect(a.imag == 2.6) catch @panic("test failure: zig_cmultf 2");
|
||||
expect(b.real == 11.3) catch @panic("test failure: zig_cmultf 3");
|
||||
expect(b.imag == -1.5) catch @panic("test failure: zig_cmultf 4");
|
||||
|
||||
return .{ .real = 1.5, .imag = 13.5 };
|
||||
}
|
||||
|
||||
export fn zig_cmultd(a: ComplexDouble, b: ComplexDouble) ComplexDouble {
|
||||
expect(a.real == 1.25) catch @panic("test failure: zig_cmultd 1");
|
||||
expect(a.imag == 2.6) catch @panic("test failure: zig_cmultd 2");
|
||||
expect(b.real == 11.3) catch @panic("test failure: zig_cmultd 3");
|
||||
expect(b.imag == -1.5) catch @panic("test failure: zig_cmultd 4");
|
||||
|
||||
return .{ .real = 1.5, .imag = 13.5 };
|
||||
}
|
||||
|
||||
export fn zig_cmultf_comp(a_r: f32, a_i: f32, b_r: f32, b_i: f32) ComplexFloat {
|
||||
expect(a_r == 1.25) catch @panic("test failure: zig_cmultf_comp 1");
|
||||
expect(a_i == 2.6) catch @panic("test failure: zig_cmultf_comp 2");
|
||||
expect(b_r == 11.3) catch @panic("test failure: zig_cmultf_comp 3");
|
||||
expect(b_i == -1.5) catch @panic("test failure: zig_cmultf_comp 4");
|
||||
|
||||
return .{ .real = 1.5, .imag = 13.5 };
|
||||
}
|
||||
|
||||
export fn zig_cmultd_comp(a_r: f64, a_i: f64, b_r: f64, b_i: f64) ComplexDouble {
|
||||
expect(a_r == 1.25) catch @panic("test failure: zig_cmultd_comp 1");
|
||||
expect(a_i == 2.6) catch @panic("test failure: zig_cmultd_comp 2");
|
||||
expect(b_r == 11.3) catch @panic("test failure: zig_cmultd_comp 3");
|
||||
expect(b_i == -1.5) catch @panic("test failure: zig_cmultd_comp 4");
|
||||
|
||||
return .{ .real = 1.5, .imag = 13.5 };
|
||||
}
|
||||
|
||||
const BigStruct = extern struct {
|
||||
a: u64,
|
||||
b: u64,
|
||||
|
||||
@ -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;
|
||||
}
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user