Merge pull request #11491 from poinwn/copysign-signbit

std.math: simplify signbit and copysign
This commit is contained in:
Andrew Kelley 2022-05-19 21:11:45 -04:00 committed by GitHub
commit 85a9bd932f
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
14 changed files with 67 additions and 191 deletions

View File

@ -57,7 +57,7 @@ pub fn fma(x: f64, y: f64, z: f64) callconv(.C) f64 {
if (spread <= 53 * 2) {
zs = math.scalbn(zs, -spread);
} else {
zs = math.copysign(f64, math.floatMin(f64), zs);
zs = math.copysign(math.floatMin(f64), zs);
}
const xy = dd_mul(xs, ys);
@ -116,7 +116,7 @@ pub fn fmaq(x: f128, y: f128, z: f128) callconv(.C) f128 {
if (spread <= 113 * 2) {
zs = math.scalbn(zs, -spread);
} else {
zs = math.copysign(f128, math.floatMin(f128), zs);
zs = math.copysign(math.floatMin(f128), zs);
}
const xy = dd_mul128(xs, ys);

View File

@ -33,7 +33,7 @@ fn atanh_32(x: f32) f32 {
var y = @bitCast(f32, i); // |x|
if (y == 1.0) {
return math.copysign(f32, math.inf(f32), x);
return math.copysign(math.inf(f32), x);
}
if (u < 0x3F800000 - (1 << 23)) {
@ -62,7 +62,7 @@ fn atanh_64(x: f64) f64 {
var y = @bitCast(f64, u & (maxInt(u64) >> 1)); // |x|
if (y == 1.0) {
return math.copysign(f64, math.inf(f64), x);
return math.copysign(math.inf(f64), x);
}
if (e < 0x3FF - 1) {

View File

@ -45,13 +45,13 @@ fn cosh32(z: Complex(f32)) Complex(f32) {
if (ix < 0x42b17218) {
// x < 88.7: exp(|x|) won't overflow
const h = @exp(@fabs(x)) * 0.5;
return Complex(f32).init(math.copysign(f32, h, x) * @cos(y), h * @sin(y));
return Complex(f32).init(math.copysign(h, x) * @cos(y), h * @sin(y));
}
// x < 192.7: scale to avoid overflow
else if (ix < 0x4340b1e7) {
const v = Complex(f32).init(@fabs(x), y);
const r = ldexp_cexp(v, -1);
return Complex(f32).init(r.re, r.im * math.copysign(f32, 1, x));
return Complex(f32).init(r.re, r.im * math.copysign(@as(f32, 1.0), x));
}
// x >= 192.7: result always overflows
else {
@ -61,14 +61,14 @@ fn cosh32(z: Complex(f32)) Complex(f32) {
}
if (ix == 0 and iy >= 0x7f800000) {
return Complex(f32).init(y - y, math.copysign(f32, 0, x * (y - y)));
return Complex(f32).init(y - y, math.copysign(@as(f32, 0.0), x * (y - y)));
}
if (iy == 0 and ix >= 0x7f800000) {
if (hx & 0x7fffff == 0) {
return Complex(f32).init(x * x, math.copysign(f32, 0, x) * y);
return Complex(f32).init(x * x, math.copysign(@as(f32, 0.0), x) * y);
}
return Complex(f32).init(x, math.copysign(f32, 0, (x + x) * y));
return Complex(f32).init(x, math.copysign(@as(f32, 0.0), (x + x) * y));
}
if (ix < 0x7f800000 and iy >= 0x7f800000) {
@ -113,13 +113,13 @@ fn cosh64(z: Complex(f64)) Complex(f64) {
if (ix < 0x40862e42) {
// x < 710: exp(|x|) won't overflow
const h = @exp(@fabs(x)) * 0.5;
return Complex(f64).init(h * @cos(y), math.copysign(f64, h, x) * @sin(y));
return Complex(f64).init(h * @cos(y), math.copysign(h, x) * @sin(y));
}
// x < 1455: scale to avoid overflow
else if (ix < 0x4096bbaa) {
const v = Complex(f64).init(@fabs(x), y);
const r = ldexp_cexp(v, -1);
return Complex(f64).init(r.re, r.im * math.copysign(f64, 1, x));
return Complex(f64).init(r.re, r.im * math.copysign(@as(f64, 1.0), x));
}
// x >= 1455: result always overflows
else {
@ -129,14 +129,14 @@ fn cosh64(z: Complex(f64)) Complex(f64) {
}
if (ix | lx == 0 and iy >= 0x7ff00000) {
return Complex(f64).init(y - y, math.copysign(f64, 0, x * (y - y)));
return Complex(f64).init(y - y, math.copysign(@as(f64, 0.0), x * (y - y)));
}
if (iy | ly == 0 and ix >= 0x7ff00000) {
if ((hx & 0xfffff) | lx == 0) {
return Complex(f64).init(x * x, math.copysign(f64, 0, x) * y);
return Complex(f64).init(x * x, math.copysign(@as(f64, 0.0), x) * y);
}
return Complex(f64).init(x * x, math.copysign(f64, 0, (x + x) * y));
return Complex(f64).init(x * x, math.copysign(@as(f64, 0.0), (x + x) * y));
}
if (ix < 0x7ff00000 and iy >= 0x7ff00000) {

View File

@ -9,7 +9,7 @@ pub fn proj(z: anytype) Complex(@TypeOf(z.re)) {
const T = @TypeOf(z.re);
if (math.isInf(z.re) or math.isInf(z.im)) {
return Complex(T).init(math.inf(T), math.copysign(T, 0, z.re));
return Complex(T).init(math.inf(T), math.copysign(@as(T, 0.0), z.re));
}
return Complex(T).init(z.re, z.im);

View File

@ -45,13 +45,13 @@ fn sinh32(z: Complex(f32)) Complex(f32) {
if (ix < 0x42b17218) {
// x < 88.7: exp(|x|) won't overflow
const h = @exp(@fabs(x)) * 0.5;
return Complex(f32).init(math.copysign(f32, h, x) * @cos(y), h * @sin(y));
return Complex(f32).init(math.copysign(h, x) * @cos(y), h * @sin(y));
}
// x < 192.7: scale to avoid overflow
else if (ix < 0x4340b1e7) {
const v = Complex(f32).init(@fabs(x), y);
const r = ldexp_cexp(v, -1);
return Complex(f32).init(r.re * math.copysign(f32, 1, x), r.im);
return Complex(f32).init(r.re * math.copysign(@as(f32, 1.0), x), r.im);
}
// x >= 192.7: result always overflows
else {
@ -61,14 +61,14 @@ fn sinh32(z: Complex(f32)) Complex(f32) {
}
if (ix == 0 and iy >= 0x7f800000) {
return Complex(f32).init(math.copysign(f32, 0, x * (y - y)), y - y);
return Complex(f32).init(math.copysign(@as(f32, 0.0), x * (y - y)), y - y);
}
if (iy == 0 and ix >= 0x7f800000) {
if (hx & 0x7fffff == 0) {
return Complex(f32).init(x, y);
}
return Complex(f32).init(x, math.copysign(f32, 0, y));
return Complex(f32).init(x, math.copysign(@as(f32, 0.0), y));
}
if (ix < 0x7f800000 and iy >= 0x7f800000) {
@ -112,13 +112,13 @@ fn sinh64(z: Complex(f64)) Complex(f64) {
if (ix < 0x40862e42) {
// x < 710: exp(|x|) won't overflow
const h = @exp(@fabs(x)) * 0.5;
return Complex(f64).init(math.copysign(f64, h, x) * @cos(y), h * @sin(y));
return Complex(f64).init(math.copysign(h, x) * @cos(y), h * @sin(y));
}
// x < 1455: scale to avoid overflow
else if (ix < 0x4096bbaa) {
const v = Complex(f64).init(@fabs(x), y);
const r = ldexp_cexp(v, -1);
return Complex(f64).init(r.re * math.copysign(f64, 1, x), r.im);
return Complex(f64).init(r.re * math.copysign(@as(f64, 1.0), x), r.im);
}
// x >= 1455: result always overflows
else {
@ -128,14 +128,14 @@ fn sinh64(z: Complex(f64)) Complex(f64) {
}
if (ix | lx == 0 and iy >= 0x7ff00000) {
return Complex(f64).init(math.copysign(f64, 0, x * (y - y)), y - y);
return Complex(f64).init(math.copysign(@as(f64, 0.0), x * (y - y)), y - y);
}
if (iy | ly == 0 and ix >= 0x7ff00000) {
if ((hx & 0xfffff) | lx == 0) {
return Complex(f64).init(x, y);
}
return Complex(f64).init(x, math.copysign(f64, 0, y));
return Complex(f64).init(x, math.copysign(@as(f64, 0.0), y));
}
if (ix < 0x7ff00000 and iy >= 0x7ff00000) {

View File

@ -43,9 +43,9 @@ fn sqrt32(z: Complex(f32)) Complex(f32) {
// sqrt(-inf + i nan) = nan +- inf i
// sqrt(-inf + iy) = 0 + inf i
if (math.signbit(x)) {
return Complex(f32).init(@fabs(x - y), math.copysign(f32, x, y));
return Complex(f32).init(@fabs(x - y), math.copysign(x, y));
} else {
return Complex(f32).init(x, math.copysign(f32, y - y, y));
return Complex(f32).init(x, math.copysign(y - y, y));
}
}
@ -65,7 +65,7 @@ fn sqrt32(z: Complex(f32)) Complex(f32) {
const t = @sqrt((-dx + math.hypot(f64, dx, dy)) * 0.5);
return Complex(f32).init(
@floatCast(f32, @fabs(y) / (2.0 * t)),
@floatCast(f32, math.copysign(f64, t, y)),
@floatCast(f32, math.copysign(t, y)),
);
}
}
@ -94,9 +94,9 @@ fn sqrt64(z: Complex(f64)) Complex(f64) {
// sqrt(-inf + i nan) = nan +- inf i
// sqrt(-inf + iy) = 0 + inf i
if (math.signbit(x)) {
return Complex(f64).init(@fabs(x - y), math.copysign(f64, x, y));
return Complex(f64).init(@fabs(x - y), math.copysign(x, y));
} else {
return Complex(f64).init(x, math.copysign(f64, y - y, y));
return Complex(f64).init(x, math.copysign(y - y, y));
}
}
@ -116,7 +116,7 @@ fn sqrt64(z: Complex(f64)) Complex(f64) {
result = Complex(f64).init(t, y / (2.0 * t));
} else {
const t = @sqrt((-x + math.hypot(f64, x, y)) * 0.5);
result = Complex(f64).init(@fabs(y) / (2.0 * t), math.copysign(f64, t, y));
result = Complex(f64).init(@fabs(y) / (2.0 * t), math.copysign(t, y));
}
if (scale) {

View File

@ -34,7 +34,7 @@ fn tanh32(z: Complex(f32)) Complex(f32) {
}
const xx = @bitCast(f32, hx - 0x40000000);
const r = if (math.isInf(y)) y else @sin(y) * @cos(y);
return Complex(f32).init(xx, math.copysign(f32, 0, r));
return Complex(f32).init(xx, math.copysign(@as(f32, 0.0), r));
}
if (!math.isFinite(y)) {
@ -45,7 +45,7 @@ fn tanh32(z: Complex(f32)) Complex(f32) {
// x >= 11
if (ix >= 0x41300000) {
const exp_mx = @exp(-@fabs(x));
return Complex(f32).init(math.copysign(f32, 1, x), 4 * @sin(y) * @cos(y) * exp_mx * exp_mx);
return Complex(f32).init(math.copysign(@as(f32, 1.0), x), 4 * @sin(y) * @cos(y) * exp_mx * exp_mx);
}
// Kahan's algorithm
@ -77,7 +77,7 @@ fn tanh64(z: Complex(f64)) Complex(f64) {
const xx = @bitCast(f64, (@as(u64, hx - 0x40000000) << 32) | lx);
const r = if (math.isInf(y)) y else @sin(y) * @cos(y);
return Complex(f64).init(xx, math.copysign(f64, 0, r));
return Complex(f64).init(xx, math.copysign(@as(f64, 0.0), r));
}
if (!math.isFinite(y)) {
@ -88,7 +88,7 @@ fn tanh64(z: Complex(f64)) Complex(f64) {
// x >= 22
if (ix >= 0x40360000) {
const exp_mx = @exp(-@fabs(x));
return Complex(f64).init(math.copysign(f64, 1, x), 4 * @sin(y) * @cos(y) * exp_mx * exp_mx);
return Complex(f64).init(math.copysign(@as(f64, 1.0), x), 4 * @sin(y) * @cos(y) * exp_mx * exp_mx);
}
// Kahan's algorithm

View File

@ -1,92 +1,25 @@
// 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/copysignf.c
// https://git.musl-libc.org/cgit/musl/tree/src/math/copysign.c
const std = @import("../std.zig");
const math = std.math;
const expect = std.testing.expect;
const maxInt = std.math.maxInt;
/// Returns a value with the magnitude of x and the sign of y.
pub fn copysign(comptime T: type, x: T, y: T) T {
return switch (T) {
f16 => copysign16(x, y),
f32 => copysign32(x, y),
f64 => copysign64(x, y),
f128 => copysign128(x, y),
else => @compileError("copysign not implemented for " ++ @typeName(T)),
};
}
fn copysign16(x: f16, y: f16) f16 {
const ux = @bitCast(u16, x);
const uy = @bitCast(u16, y);
const h1 = ux & (maxInt(u16) / 2);
const h2 = uy & (@as(u16, 1) << 15);
return @bitCast(f16, h1 | h2);
}
fn copysign32(x: f32, y: f32) f32 {
const ux = @bitCast(u32, x);
const uy = @bitCast(u32, y);
const h1 = ux & (maxInt(u32) / 2);
const h2 = uy & (@as(u32, 1) << 31);
return @bitCast(f32, h1 | h2);
}
fn copysign64(x: f64, y: f64) f64 {
const ux = @bitCast(u64, x);
const uy = @bitCast(u64, y);
const h1 = ux & (maxInt(u64) / 2);
const h2 = uy & (@as(u64, 1) << 63);
return @bitCast(f64, h1 | h2);
}
fn copysign128(x: f128, y: f128) f128 {
const ux = @bitCast(u128, x);
const uy = @bitCast(u128, y);
const h1 = ux & (maxInt(u128) / 2);
const h2 = uy & (@as(u128, 1) << 127);
return @bitCast(f128, h1 | h2);
/// Returns a value with the magnitude of `magnitude` and the sign of `sign`.
pub fn copysign(magnitude: anytype, sign: @TypeOf(magnitude)) @TypeOf(magnitude) {
const T = @TypeOf(magnitude);
const TBits = std.meta.Int(.unsigned, @typeInfo(T).Float.bits);
const sign_bit_mask = @as(TBits, 1) << (@bitSizeOf(T) - 1);
const mag = @bitCast(TBits, magnitude) & ~sign_bit_mask;
const sgn = @bitCast(TBits, sign) & sign_bit_mask;
return @bitCast(T, mag | sgn);
}
test "math.copysign" {
try expect(copysign(f16, 1.0, 1.0) == copysign16(1.0, 1.0));
try expect(copysign(f32, 1.0, 1.0) == copysign32(1.0, 1.0));
try expect(copysign(f64, 1.0, 1.0) == copysign64(1.0, 1.0));
try expect(copysign(f128, 1.0, 1.0) == copysign128(1.0, 1.0));
}
test "math.copysign16" {
try expect(copysign16(5.0, 1.0) == 5.0);
try expect(copysign16(5.0, -1.0) == -5.0);
try expect(copysign16(-5.0, -1.0) == -5.0);
try expect(copysign16(-5.0, 1.0) == 5.0);
}
test "math.copysign32" {
try expect(copysign32(5.0, 1.0) == 5.0);
try expect(copysign32(5.0, -1.0) == -5.0);
try expect(copysign32(-5.0, -1.0) == -5.0);
try expect(copysign32(-5.0, 1.0) == 5.0);
}
test "math.copysign64" {
try expect(copysign64(5.0, 1.0) == 5.0);
try expect(copysign64(5.0, -1.0) == -5.0);
try expect(copysign64(-5.0, -1.0) == -5.0);
try expect(copysign64(-5.0, 1.0) == 5.0);
}
test "math.copysign128" {
try expect(copysign128(5.0, 1.0) == 5.0);
try expect(copysign128(5.0, -1.0) == -5.0);
try expect(copysign128(-5.0, -1.0) == -5.0);
try expect(copysign128(-5.0, 1.0) == 5.0);
inline for ([_]type{ f16, f32, f64, f80, f128 }) |T| {
try expect(copysign(@as(T, 1.0), @as(T, 1.0)) == 1.0);
try expect(copysign(@as(T, 2.0), @as(T, -2.0)) == -2.0);
try expect(copysign(@as(T, -3.0), @as(T, 3.0)) == 3.0);
try expect(copysign(@as(T, -4.0), @as(T, -4.0)) == -4.0);
try expect(copysign(@as(T, 5.0), @as(T, -500.0)) == -5.0);
try expect(copysign(math.inf(T), @as(T, -0.0)) == -math.inf(T));
try expect(copysign(@as(T, 6.0), -math.nan(T)) == -6.0);
}
}

View File

@ -5,10 +5,7 @@ const expect = std.testing.expect;
/// Returns whether x is a finite value.
pub fn isFinite(x: anytype) bool {
const T = @TypeOf(x);
const TBits = std.meta.Int(.unsigned, @bitSizeOf(T));
if (@typeInfo(T) != .Float) {
@compileError("isFinite not implemented for " ++ @typeName(T));
}
const TBits = std.meta.Int(.unsigned, @typeInfo(T).Float.bits);
const remove_sign = ~@as(TBits, 0) >> 1;
return @bitCast(TBits, x) & remove_sign < @bitCast(TBits, math.inf(T));
}

View File

@ -5,10 +5,7 @@ const expect = std.testing.expect;
/// Returns whether x is an infinity, ignoring sign.
pub fn isInf(x: anytype) bool {
const T = @TypeOf(x);
const TBits = std.meta.Int(.unsigned, @bitSizeOf(T));
if (@typeInfo(T) != .Float) {
@compileError("isInf not implemented for " ++ @typeName(T));
}
const TBits = std.meta.Int(.unsigned, @typeInfo(T).Float.bits);
const remove_sign = ~@as(TBits, 0) >> 1;
return @bitCast(TBits, x) & remove_sign == @bitCast(TBits, math.inf(T));
}

View File

@ -5,10 +5,7 @@ const expect = std.testing.expect;
/// Returns whether x is neither zero, subnormal, infinity, or NaN.
pub fn isNormal(x: anytype) bool {
const T = @TypeOf(x);
const TBits = std.meta.Int(.unsigned, @bitSizeOf(T));
if (@typeInfo(T) != .Float) {
@compileError("isNormal not implemented for " ++ @typeName(T));
}
const TBits = std.meta.Int(.unsigned, @typeInfo(T).Float.bits);
const increment_exp = 1 << math.floatMantissaBits(T);
const remove_sign = ~@as(TBits, 0) >> 1;

View File

@ -15,10 +15,7 @@ pub fn ldexp(x: anytype, n: i32) @TypeOf(x) {
var shift = n;
const T = @TypeOf(base);
const TBits = std.meta.Int(.unsigned, @bitSizeOf(T));
if (@typeInfo(T) != .Float) {
@compileError("ldexp not implemented for " ++ @typeName(T));
}
const TBits = std.meta.Int(.unsigned, @typeInfo(T).Float.bits);
const mantissa_bits = math.floatMantissaBits(T);
const exponent_min = math.floatExponentMin(T);

View File

@ -60,7 +60,7 @@ pub fn pow(comptime T: type, x: T, y: T) T {
if (y < 0) {
// pow(+-0, y) = +- 0 for y an odd integer
if (isOddInteger(y)) {
return math.copysign(T, math.inf(T), x);
return math.copysign(math.inf(T), x);
}
// pow(+-0, y) = +inf for y an even integer
else {

View File

@ -5,64 +5,19 @@ const expect = std.testing.expect;
/// Returns whether x is negative or negative 0.
pub fn signbit(x: anytype) bool {
const T = @TypeOf(x);
return switch (T) {
f16 => signbit16(x),
f32 => signbit32(x),
f64 => signbit64(x),
f80 => signbit80(x),
f128 => signbit128(x),
else => @compileError("signbit not implemented for " ++ @typeName(T)),
};
}
fn signbit16(x: f16) bool {
const bits = @bitCast(u16, x);
return bits >> 15 != 0;
}
fn signbit32(x: f32) bool {
const bits = @bitCast(u32, x);
return bits >> 31 != 0;
}
fn signbit64(x: f64) bool {
const bits = @bitCast(u64, x);
return bits >> 63 != 0;
}
fn signbit80(x: f80) bool {
const bits = @bitCast(u80, x);
return bits >> 79 != 0;
}
fn signbit128(x: f128) bool {
const bits = @bitCast(u128, x);
return bits >> 127 != 0;
const TBits = std.meta.Int(.unsigned, @typeInfo(T).Float.bits);
return @bitCast(TBits, x) >> (@bitSizeOf(T) - 1) != 0;
}
test "math.signbit" {
try expect(signbit(@as(f16, 4.0)) == signbit16(4.0));
try expect(signbit(@as(f32, 4.0)) == signbit32(4.0));
try expect(signbit(@as(f64, 4.0)) == signbit64(4.0));
try expect(signbit(@as(f128, 4.0)) == signbit128(4.0));
}
test "math.signbit16" {
try expect(!signbit16(4.0));
try expect(signbit16(-3.0));
}
test "math.signbit32" {
try expect(!signbit32(4.0));
try expect(signbit32(-3.0));
}
test "math.signbit64" {
try expect(!signbit64(4.0));
try expect(signbit64(-3.0));
}
test "math.signbit128" {
try expect(!signbit128(4.0));
try expect(signbit128(-3.0));
inline for ([_]type{ f16, f32, f64, f80, f128 }) |T| {
try expect(!signbit(@as(T, 0.0)));
try expect(!signbit(@as(T, 1.0)));
try expect(signbit(@as(T, -2.0)));
try expect(signbit(@as(T, -0.0)));
try expect(!signbit(math.inf(T)));
try expect(signbit(-math.inf(T)));
try expect(!signbit(math.nan(T)));
try expect(signbit(-math.nan(T)));
}
}