authorgravatar for thatlemon@gmail.comLemonBoy <thatlemon@gmail.com> 2021-06-11 09:23:48+02:00
committergravatar for git@vexu.euVeikka Tuominen <git@vexu.eu> 2021-06-11 14:04:46+03:00
log56c0a66ce30094d199420a3f2b2d2e759c1b0613
tree6886f18e18b63384d5345ccd226dcb28cd2a0820
parent2ce033f4159409db66866e28a3591256ea0dfadc

std: Fix tanh for negative inputs

It turns out the code was not ported correctly from C and produced wrong results for negative input values. As a bonus fix the NaN codepath by adding yet another missing piece of code. Spotted in #9047

1 files changed, 25 insertions(+), 33 deletions(-)

lib/std/math/tanh.zig+25-33
...@@ -37,89 +37,75 @@ fn tanh32(x: f32) f32 {...@@ -37,89 +37,75 @@ fn tanh32(x: f32) f32 {
37 const u = @bitCast(u32, x);37 const u = @bitCast(u32, x);
38 const ux = u & 0x7FFFFFFF;38 const ux = u & 0x7FFFFFFF;
39 const ax = @bitCast(f32, ux);39 const ax = @bitCast(f32, ux);
40 const sign = (u >> 31) != 0;
4041
41 var t: f32 = undefined;42 var t: f32 = undefined;
4243
43 if (x == 0.0 or math.isNan(x)) {
44 return x;
45 }
46
47 // |x| < log(3) / 2 ~= 0.5493 or nan44 // |x| < log(3) / 2 ~= 0.5493 or nan
48 if (ux > 0x3F0C9F54) {45 if (ux > 0x3F0C9F54) {
49 // |x| > 1046 // |x| > 10
50 if (ux > 0x41200000) {47 if (ux > 0x41200000) {
51 t = 1.0;48 t = 1.0 + 0 / x;
52 } else {49 } else {
53 t = math.expm1(2 * x);50 t = math.expm1(2 * ax);
54 t = 1 - 2 / (t + 2);51 t = 1 - 2 / (t + 2);
55 }52 }
56 }53 }
57 // |x| > log(5 / 3) / 2 ~= 0.255454 // |x| > log(5 / 3) / 2 ~= 0.2554
58 else if (ux > 0x3E82C578) {55 else if (ux > 0x3E82C578) {
59 t = math.expm1(2 * x);56 t = math.expm1(2 * ax);
60 t = t / (t + 2);57 t = t / (t + 2);
61 }58 }
62 // |x| >= 0x1.0p-12659 // |x| >= 0x1.0p-126
63 else if (ux >= 0x00800000) {60 else if (ux >= 0x00800000) {
64 t = math.expm1(-2 * x);61 t = math.expm1(-2 * ax);
65 t = -t / (t + 2);62 t = -t / (t + 2);
66 }63 }
67 // |x| is subnormal64 // |x| is subnormal
68 else {65 else {
69 math.doNotOptimizeAway(x * x);66 math.doNotOptimizeAway(ax * ax);
70 t = x;67 t = ax;
71 }68 }
7269
73 if (u >> 31 != 0) {70 return if (sign) -t else t;
74 return -t;
75 } else {
76 return t;
77 }
78}71}
7972
80fn tanh64(x: f64) f64 {73fn tanh64(x: f64) f64 {
81 const u = @bitCast(u64, x);74 const u = @bitCast(u64, x);
82 const w = @intCast(u32, u >> 32);75 const ux = u & 0x7FFFFFFFFFFFFFFF;
83 const ax = @bitCast(f64, u & (maxInt(u64) >> 1));76 const w = @intCast(u32, ux >> 32);
77 const ax = @bitCast(f64, ux);
78 const sign = (u >> 63) != 0;
8479
85 var t: f64 = undefined;80 var t: f64 = undefined;
8681
87 // TODO: Shouldn't need these checks.
88 if (x == 0.0 or math.isNan(x)) {
89 return x;
90 }
91
92 // |x| < log(3) / 2 ~= 0.5493 or nan82 // |x| < log(3) / 2 ~= 0.5493 or nan
93 if (w > 0x3FE193EA) {83 if (w > 0x3FE193EA) {
94 // |x| > 20 or nan84 // |x| > 20 or nan
95 if (w > 0x40340000) {85 if (w > 0x40340000) {
96 t = 1.0;86 t = 1.0 - 0 / ax;
97 } else {87 } else {
98 t = math.expm1(2 * x);88 t = math.expm1(2 * ax);
99 t = 1 - 2 / (t + 2);89 t = 1 - 2 / (t + 2);
100 }90 }
101 }91 }
102 // |x| > log(5 / 3) / 2 ~= 0.255492 // |x| > log(5 / 3) / 2 ~= 0.2554
103 else if (w > 0x3FD058AE) {93 else if (w > 0x3FD058AE) {
104 t = math.expm1(2 * x);94 t = math.expm1(2 * ax);
105 t = t / (t + 2);95 t = t / (t + 2);
106 }96 }
107 // |x| >= 0x1.0p-102297 // |x| >= 0x1.0p-1022
108 else if (w >= 0x00100000) {98 else if (w >= 0x00100000) {
109 t = math.expm1(-2 * x);99 t = math.expm1(-2 * ax);
110 t = -t / (t + 2);100 t = -t / (t + 2);
111 }101 }
112 // |x| is subnormal102 // |x| is subnormal
113 else {103 else {
114 math.doNotOptimizeAway(@floatCast(f32, x));104 math.doNotOptimizeAway(@floatCast(f32, ax));
115 t = x;105 t = ax;
116 }106 }
117107
118 if (u >> 63 != 0) {108 return if (sign) -t else t;
119 return -t;
120 } else {
121 return t;
122 }
123}109}
124110
125test "math.tanh" {111test "math.tanh" {
...@@ -135,6 +121,9 @@ test "math.tanh32" {...@@ -135,6 +121,9 @@ test "math.tanh32" {
135 try expect(math.approxEqAbs(f32, tanh32(0.8923), 0.712528, epsilon));121 try expect(math.approxEqAbs(f32, tanh32(0.8923), 0.712528, epsilon));
136 try expect(math.approxEqAbs(f32, tanh32(1.5), 0.905148, epsilon));122 try expect(math.approxEqAbs(f32, tanh32(1.5), 0.905148, epsilon));
137 try expect(math.approxEqAbs(f32, tanh32(37.45), 1.0, epsilon));123 try expect(math.approxEqAbs(f32, tanh32(37.45), 1.0, epsilon));
124 try expect(math.approxEqAbs(f32, tanh32(-0.8923), -0.712528, epsilon));
125 try expect(math.approxEqAbs(f32, tanh32(-1.5), -0.905148, epsilon));
126 try expect(math.approxEqAbs(f32, tanh32(-37.45), -1.0, epsilon));
138}127}
139128
140test "math.tanh64" {129test "math.tanh64" {
...@@ -145,6 +134,9 @@ test "math.tanh64" {...@@ -145,6 +134,9 @@ test "math.tanh64" {
145 try expect(math.approxEqAbs(f64, tanh64(0.8923), 0.712528, epsilon));134 try expect(math.approxEqAbs(f64, tanh64(0.8923), 0.712528, epsilon));
146 try expect(math.approxEqAbs(f64, tanh64(1.5), 0.905148, epsilon));135 try expect(math.approxEqAbs(f64, tanh64(1.5), 0.905148, epsilon));
147 try expect(math.approxEqAbs(f64, tanh64(37.45), 1.0, epsilon));136 try expect(math.approxEqAbs(f64, tanh64(37.45), 1.0, epsilon));
137 try expect(math.approxEqAbs(f64, tanh64(-0.8923), -0.712528, epsilon));
138 try expect(math.approxEqAbs(f64, tanh64(-1.5), -0.905148, epsilon));
139 try expect(math.approxEqAbs(f64, tanh64(-37.45), -1.0, epsilon));
148}140}
149141
150test "math.tanh32.special" {142test "math.tanh32.special" {