| ... | @@ -28,30 +28,25 @@ pub fn asinh(x: anytype) @TypeOf(x) { | ... | @@ -28,30 +28,25 @@ pub fn asinh(x: anytype) @TypeOf(x) { |
| 28 | fn asinh32(x: f32) f32 { | 28 | fn asinh32(x: f32) f32 { |
| 29 | const u = @as(u32, @bitCast(x)); | 29 | const u = @as(u32, @bitCast(x)); |
| 30 | const i = u & 0x7FFFFFFF; | 30 | const i = u & 0x7FFFFFFF; |
| 31 | const s = i >> 31; | 31 | const s = u >> 31; |
| 32 | | 32 | |
| 33 | var rx = @as(f32, @bitCast(i)); // |x| | 33 | var rx = @as(f32, @bitCast(i)); // |x| |
| 34 | | 34 | |
| 35 | // TODO: Shouldn't need this explicit check. | | |
| 36 | if (math.isNegativeInf(x)) { | | |
| 37 | return x; | | |
| 38 | } | | |
| 39 | | | |
| 40 | // |x| >= 0x1p12 or inf or nan | 35 | // |x| >= 0x1p12 or inf or nan |
| 41 | if (i >= 0x3F800000 + (12 << 23)) { | 36 | if (i >= 0x3F800000 + (12 << 23)) { |
| 42 | rx = @log(rx) + 0.69314718055994530941723212145817656; | 37 | rx = @log(rx) + 0.69314718055994530941723212145817656; |
| 43 | } | 38 | } |
| 44 | // |x| >= 2 | 39 | // |x| >= 2 |
| 45 | else if (i >= 0x3F800000 + (1 << 23)) { | 40 | else if (i >= 0x3F800000 + (1 << 23)) { |
| 46 | rx = @log(2 * x + 1 / (@sqrt(x * x + 1) + x)); | 41 | rx = @log(2 * rx + 1 / (@sqrt(rx * rx + 1) + rx)); |
| 47 | } | 42 | } |
| 48 | // |x| >= 0x1p-12, up to 1.6ulp error | 43 | // |x| >= 0x1p-12, up to 1.6ulp error |
| 49 | else if (i >= 0x3F800000 - (12 << 23)) { | 44 | else if (i >= 0x3F800000 - (12 << 23)) { |
| 50 | rx = math.log1p(x + x * x / (@sqrt(x * x + 1) + 1)); | 45 | rx = math.log1p(rx + rx * rx / (@sqrt(rx * rx + 1) + 1)); |
| 51 | } | 46 | } |
| 52 | // |x| < 0x1p-12, inexact if x != 0 | 47 | // |x| < 0x1p-12, inexact if x != 0 |
| 53 | else { | 48 | else { |
| 54 | math.doNotOptimizeAway(x + 0x1.0p120); | 49 | math.doNotOptimizeAway(rx + 0x1.0p120); |
| 55 | } | 50 | } |
| 56 | | 51 | |
| 57 | return if (s != 0) -rx else rx; | 52 | return if (s != 0) -rx else rx; |
| ... | @@ -60,29 +55,25 @@ fn asinh32(x: f32) f32 { | ... | @@ -60,29 +55,25 @@ fn asinh32(x: f32) f32 { |
| 60 | fn asinh64(x: f64) f64 { | 55 | fn asinh64(x: f64) f64 { |
| 61 | const u = @as(u64, @bitCast(x)); | 56 | const u = @as(u64, @bitCast(x)); |
| 62 | const e = (u >> 52) & 0x7FF; | 57 | const e = (u >> 52) & 0x7FF; |
| 63 | const s = e >> 63; | 58 | const s = u >> 63; |
| 64 | | 59 | |
| 65 | var rx = @as(f64, @bitCast(u & (maxInt(u64) >> 1))); // |x| | 60 | var rx = @as(f64, @bitCast(u & (maxInt(u64) >> 1))); // |x| |
| 66 | | 61 | |
| 67 | if (math.isNegativeInf(x)) { | | |
| 68 | return x; | | |
| 69 | } | | |
| 70 | | | |
| 71 | // |x| >= 0x1p26 or inf or nan | 62 | // |x| >= 0x1p26 or inf or nan |
| 72 | if (e >= 0x3FF + 26) { | 63 | if (e >= 0x3FF + 26) { |
| 73 | rx = @log(rx) + 0.693147180559945309417232121458176568; | 64 | rx = @log(rx) + 0.693147180559945309417232121458176568; |
| 74 | } | 65 | } |
| 75 | // |x| >= 2 | 66 | // |x| >= 2 |
| 76 | else if (e >= 0x3FF + 1) { | 67 | else if (e >= 0x3FF + 1) { |
| 77 | rx = @log(2 * x + 1 / (@sqrt(x * x + 1) + x)); | 68 | rx = @log(2 * rx + 1 / (@sqrt(rx * rx + 1) + rx)); |
| 78 | } | 69 | } |
| 79 | // |x| >= 0x1p-12, up to 1.6ulp error | 70 | // |x| >= 0x1p-12, up to 1.6ulp error |
| 80 | else if (e >= 0x3FF - 26) { | 71 | else if (e >= 0x3FF - 26) { |
| 81 | rx = math.log1p(x + x * x / (@sqrt(x * x + 1) + 1)); | 72 | rx = math.log1p(rx + rx * rx / (@sqrt(rx * rx + 1) + 1)); |
| 82 | } | 73 | } |
| 83 | // |x| < 0x1p-12, inexact if x != 0 | 74 | // |x| < 0x1p-12, inexact if x != 0 |
| 84 | else { | 75 | else { |
| 85 | math.doNotOptimizeAway(x + 0x1.0p120); | 76 | math.doNotOptimizeAway(rx + 0x1.0p120); |
| 86 | } | 77 | } |
| 87 | | 78 | |
| 88 | return if (s != 0) -rx else rx; | 79 | return if (s != 0) -rx else rx; |
| ... | @@ -121,7 +112,7 @@ test "math.asinh64" { | ... | @@ -121,7 +112,7 @@ test "math.asinh64" { |
| 121 | | 112 | |
| 122 | test "math.asinh32.special" { | 113 | test "math.asinh32.special" { |
| 123 | try expect(asinh32(0.0) == 0.0); | 114 | try expect(asinh32(0.0) == 0.0); |
| 124 | try expect(asinh32(-0.0) == -0.0); | 115 | try expect(@as(u32, @bitCast(asinh32(-0.0))) == @as(u32, 2147483648)); |
| 125 | try expect(math.isPositiveInf(asinh32(math.inf(f32)))); | 116 | try expect(math.isPositiveInf(asinh32(math.inf(f32)))); |
| 126 | try expect(math.isNegativeInf(asinh32(-math.inf(f32)))); | 117 | try expect(math.isNegativeInf(asinh32(-math.inf(f32)))); |
| 127 | try expect(math.isNan(asinh32(math.nan(f32)))); | 118 | try expect(math.isNan(asinh32(math.nan(f32)))); |
| ... | @@ -129,7 +120,7 @@ test "math.asinh32.special" { | ... | @@ -129,7 +120,7 @@ test "math.asinh32.special" { |
| 129 | | 120 | |
| 130 | test "math.asinh64.special" { | 121 | test "math.asinh64.special" { |
| 131 | try expect(asinh64(0.0) == 0.0); | 122 | try expect(asinh64(0.0) == 0.0); |
| 132 | try expect(asinh64(-0.0) == -0.0); | 123 | try expect(@as(u64, @bitCast(asinh64(-0.0))) == @as(u64, 9223372036854775808)); |
| 133 | try expect(math.isPositiveInf(asinh64(math.inf(f64)))); | 124 | try expect(math.isPositiveInf(asinh64(math.inf(f64)))); |
| 134 | try expect(math.isNegativeInf(asinh64(-math.inf(f64)))); | 125 | try expect(math.isNegativeInf(asinh64(-math.inf(f64)))); |
| 135 | try expect(math.isNan(asinh64(math.nan(f64)))); | 126 | try expect(math.isNan(asinh64(math.nan(f64)))); |