1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
3//
4// https://git.musl-libc.org/cgit/musl/tree/src/math/expmf.c
5// https://git.musl-libc.org/cgit/musl/tree/src/math/expm.c
6
7// TODO: Updated recently.
8
9const std = @import("../std.zig");
10const math = std.math;
11const mem = std.mem;
12const expect = std.testing.expect;
13const expectEqual = std.testing.expectEqual;
14
15/// Returns e raised to the power of x, minus 1 (e^x - 1). This is more accurate than exp(e, x) - 1
16/// when x is near 0.
17///
18/// Special Cases:
19/// - expm1(+inf) = +inf
20/// - expm1(-inf) = -1
21/// - expm1(nan) = nan
22pub fn expm1(x: anytype) @TypeOf(x) {
23 const T = @TypeOf(x);
24 return switch (T) {
25 f32 => expm1_32(x),
26 f64 => expm1_64(x),
27 else => @compileError("exp1m not implemented for " ++ @typeName(T)),
28 };
29}
30
31fn expm1_32(x_: f32) f32 {
32 if (math.isNan(x_))
33 return math.nan(f32);
34
35 const o_threshold: f32 = 8.8721679688e+01;
36 const ln2_hi: f32 = 6.9313812256e-01;
37 const ln2_lo: f32 = 9.0580006145e-06;
38 const invln2: f32 = 1.4426950216e+00;
39 const Q1: f32 = -3.3333212137e-2;
40 const Q2: f32 = 1.5807170421e-3;
41
42 var x = x_;
43 const ux: u32 = @bitCast(x);
44 const hx = ux & 0x7FFFFFFF;
45 const sign = ux >> 31;
46
47 // TODO: Shouldn't need this check explicitly.
48 if (math.isNegativeInf(x)) {
49 return -1.0;
50 }
51
52 // |x| >= 27 * ln2
53 if (hx >= 0x4195B844) {
54 // nan
55 if (hx > 0x7F800000) {
56 return x;
57 }
58 if (sign != 0) {
59 return -1;
60 }
61 if (x > o_threshold) {
62 x *= 0x1.0p127;
63 return x;
64 }
65 }
66
67 var hi: f32 = undefined;
68 var lo: f32 = undefined;
69 var c: f32 = undefined;
70 var k: i32 = undefined;
71
72 // |x| > 0.5 * ln2
73 if (hx > 0x3EB17218) {
74 // |x| < 1.5 * ln2
75 if (hx < 0x3F851592) {
76 if (sign == 0) {
77 hi = x - ln2_hi;
78 lo = ln2_lo;
79 k = 1;
80 } else {
81 hi = x + ln2_hi;
82 lo = -ln2_lo;
83 k = -1;
84 }
85 } else {
86 var kf = invln2 * x;
87 if (sign != 0) {
88 kf -= 0.5;
89 } else {
90 kf += 0.5;
91 }
92
93 k = @as(i32, @intFromFloat(kf));
94 const t = @as(f32, @floatFromInt(k));
95 hi = x - t * ln2_hi;
96 lo = t * ln2_lo;
97 }
98
99 x = hi - lo;
100 c = (hi - x) - lo;
101 }
102 // |x| < 2^(-25)
103 else if (hx < 0x33000000) {
104 if (hx < 0x00800000) {
105 mem.doNotOptimizeAway(x * x);
106 }
107 return x;
108 } else {
109 k = 0;
110 }
111
112 const hfx = 0.5 * x;
113 const hxs = x * hfx;
114 const r1 = 1.0 + hxs * (Q1 + hxs * Q2);
115 const t = 3.0 - r1 * hfx;
116 var e = hxs * ((r1 - t) / (6.0 - x * t));
117
118 // c is 0
119 if (k == 0) {
120 return x - (x * e - hxs);
121 }
122
123 e = x * (e - c) - c;
124 e -= hxs;
125
126 // exp(x) ~ 2^k (x_reduced - e + 1)
127 if (k == -1) {
128 return 0.5 * (x - e) - 0.5;
129 }
130 if (k == 1) {
131 if (x < -0.25) {
132 return -2.0 * (e - (x + 0.5));
133 } else {
134 return 1.0 + 2.0 * (x - e);
135 }
136 }
137
138 const twopk = @as(f32, @bitCast(@as(u32, @intCast((0x7F +% k) << 23))));
139
140 if (k < 0 or k > 56) {
141 var y = x - e + 1.0;
142 if (k == 128) {
143 y = y * 2.0 * 0x1.0p127;
144 } else {
145 y = y * twopk;
146 }
147
148 return y - 1.0;
149 }
150
151 const uf: f32 = @bitCast(@as(u32, @intCast(0x7F -% k)) << 23);
152 if (k < 23) {
153 return (x - e + (1 - uf)) * twopk;
154 } else {
155 return (x - (e + uf) + 1) * twopk;
156 }
157}
158
159fn expm1_64(x_: f64) f64 {
160 if (math.isNan(x_))
161 return math.nan(f64);
162
163 const o_threshold: f64 = 7.09782712893383973096e+02;
164 const ln2_hi: f64 = 6.93147180369123816490e-01;
165 const ln2_lo: f64 = 1.90821492927058770002e-10;
166 const invln2: f64 = 1.44269504088896338700e+00;
167 const Q1: f64 = -3.33333333333331316428e-02;
168 const Q2: f64 = 1.58730158725481460165e-03;
169 const Q3: f64 = -7.93650757867487942473e-05;
170 const Q4: f64 = 4.00821782732936239552e-06;
171 const Q5: f64 = -2.01099218183624371326e-07;
172
173 var x = x_;
174 const ux = @as(u64, @bitCast(x));
175 const hx = @as(u32, @intCast(ux >> 32)) & 0x7FFFFFFF;
176 const sign = ux >> 63;
177
178 if (math.isNegativeInf(x)) {
179 return -1.0;
180 }
181
182 // |x| >= 56 * ln2
183 if (hx >= 0x4043687A) {
184 // exp1md(nan) = nan
185 if (hx > 0x7FF00000) {
186 return x;
187 }
188 // exp1md(-ve) = -1
189 if (sign != 0) {
190 return -1;
191 }
192 if (x > o_threshold) {
193 math.raiseOverflow();
194 return math.inf(f64);
195 }
196 }
197
198 var hi: f64 = undefined;
199 var lo: f64 = undefined;
200 var c: f64 = undefined;
201 var k: i32 = undefined;
202
203 // |x| > 0.5 * ln2
204 if (hx > 0x3FD62E42) {
205 // |x| < 1.5 * ln2
206 if (hx < 0x3FF0A2B2) {
207 if (sign == 0) {
208 hi = x - ln2_hi;
209 lo = ln2_lo;
210 k = 1;
211 } else {
212 hi = x + ln2_hi;
213 lo = -ln2_lo;
214 k = -1;
215 }
216 } else {
217 var kf = invln2 * x;
218 if (sign != 0) {
219 kf -= 0.5;
220 } else {
221 kf += 0.5;
222 }
223
224 k = @as(i32, @intFromFloat(kf));
225 const t = @as(f64, @floatFromInt(k));
226 hi = x - t * ln2_hi;
227 lo = t * ln2_lo;
228 }
229
230 x = hi - lo;
231 c = (hi - x) - lo;
232 }
233 // |x| < 2^(-54)
234 else if (hx < 0x3C900000) {
235 if (hx < 0x00100000) {
236 mem.doNotOptimizeAway(@as(f32, @floatCast(x)));
237 }
238 return x;
239 } else {
240 k = 0;
241 }
242
243 const hfx = 0.5 * x;
244 const hxs = x * hfx;
245 const r1 = 1.0 + hxs * (Q1 + hxs * (Q2 + hxs * (Q3 + hxs * (Q4 + hxs * Q5))));
246 const t = 3.0 - r1 * hfx;
247 var e = hxs * ((r1 - t) / (6.0 - x * t));
248
249 // c is 0
250 if (k == 0) {
251 return x - (x * e - hxs);
252 }
253
254 e = x * (e - c) - c;
255 e -= hxs;
256
257 // exp(x) ~ 2^k (x_reduced - e + 1)
258 if (k == -1) {
259 return 0.5 * (x - e) - 0.5;
260 }
261 if (k == 1) {
262 if (x < -0.25) {
263 return -2.0 * (e - (x + 0.5));
264 } else {
265 return 1.0 + 2.0 * (x - e);
266 }
267 }
268
269 const twopk = @as(f64, @bitCast(@as(u64, @intCast(0x3FF +% k)) << 52));
270
271 if (k < 0 or k > 56) {
272 var y = x - e + 1.0;
273 if (k == 1024) {
274 y = y * 2.0 * 0x1.0p1023;
275 } else {
276 y = y * twopk;
277 }
278
279 return y - 1.0;
280 }
281
282 const uf = @as(f64, @bitCast(@as(u64, @intCast(0x3FF -% k)) << 52));
283 if (k < 20) {
284 return (x - e + (1 - uf)) * twopk;
285 } else {
286 return (x - (e + uf) + 1) * twopk;
287 }
288}
289
290test "expm1_32() special" {
291 try expect(math.isPositiveZero(expm1_32(0.0)));
292 try expect(math.isNegativeZero(expm1_32(-0.0)));
293 try expectEqual(expm1_32(math.ln2), 1.0);
294 try expectEqual(expm1_32(math.inf(f32)), math.inf(f32));
295 try expectEqual(expm1_32(-math.inf(f32)), -1.0);
296 try expect(math.isNan(expm1_32(math.nan(f32))));
297 try expect(math.isNan(expm1_32(math.snan(f32))));
298}
299
300test "expm1_32() sanity" {
301 try expectEqual(expm1_32(-0x1.0223a0p+3), -0x1.ffd6e0p-1);
302 try expectEqual(expm1_32(0x1.161868p+2), 0x1.30712ap+6);
303 try expectEqual(expm1_32(-0x1.0c34b4p+3), -0x1.ffe1fap-1);
304 try expectEqual(expm1_32(-0x1.a206f0p+2), -0x1.ff4116p-1);
305 try expectEqual(expm1_32(0x1.288bbcp+3), 0x1.4ab480p+13); // Disagrees with GCC in last bit
306 try expectEqual(expm1_32(0x1.52efd0p-1), 0x1.e09536p-1);
307 try expectEqual(expm1_32(-0x1.a05cc8p-2), -0x1.561c3ep-2);
308 try expectEqual(expm1_32(0x1.1f9efap-1), 0x1.81ec4ep-1);
309 try expectEqual(expm1_32(0x1.8c5db0p-1), 0x1.2b3364p+0);
310 try expectEqual(expm1_32(-0x1.5b86eap-1), -0x1.f8951ap-2);
311}
312
313test "expm1_32() boundary" {
314 // TODO: The last value before inf is actually 0x1.62e300p+6 -> 0x1.ff681ep+127
315 // try expectEqual(expm1_32(0x1.62e42ep+6), 0x1.ffff08p+127); // Last value before result is inf
316 try expectEqual(expm1_32(0x1.62e430p+6), math.inf(f32)); // First value that gives inf
317 try expectEqual(expm1_32(0x1.fffffep+127), math.inf(f32)); // Max input value
318 try expectEqual(expm1_32(0x1p-149), 0x1p-149); // Min positive input value
319 try expectEqual(expm1_32(-0x1p-149), -0x1p-149); // Min negative input value
320 try expectEqual(expm1_32(0x1p-126), 0x1p-126); // First positive subnormal input
321 try expectEqual(expm1_32(-0x1p-126), -0x1p-126); // First negative subnormal input
322 try expectEqual(expm1_32(0x1.fffffep-125), 0x1.fffffep-125); // Last positive value before subnormal
323 try expectEqual(expm1_32(-0x1.fffffep-125), -0x1.fffffep-125); // Last negative value before subnormal
324 try expectEqual(expm1_32(-0x1.154244p+4), -0x1.fffffep-1); // Last value before result is -1
325 try expectEqual(expm1_32(-0x1.154246p+4), -1); // First value where result is -1
326}
327
328test "expm1_64() special" {
329 try expect(math.isPositiveZero(expm1_64(0.0)));
330 try expect(math.isNegativeZero(expm1_64(-0.0)));
331 try expectEqual(expm1_64(math.ln2), 1.0);
332 try expectEqual(expm1_64(math.inf(f64)), math.inf(f64));
333 try expectEqual(expm1_64(-math.inf(f64)), -1.0);
334 try expect(math.isNan(expm1_64(math.nan(f64))));
335 try expect(math.isNan(expm1_64(math.snan(f64))));
336}
337
338test "expm1_64() sanity" {
339 try expectEqual(expm1_64(-0x1.02239f3c6a8f1p+3), -0x1.ffd6df9b02b3ep-1);
340 try expectEqual(expm1_64(0x1.161868e18bc67p+2), 0x1.30712ed238c04p+6);
341 try expectEqual(expm1_64(-0x1.0c34b3e01e6e7p+3), -0x1.ffe1f94e493e7p-1);
342 try expectEqual(expm1_64(-0x1.a206f0a19dcc4p+2), -0x1.ff4115c03f78dp-1);
343 try expectEqual(expm1_64(0x1.288bbb0d6a1e6p+3), 0x1.4ab477496e07ep+13);
344 try expectEqual(expm1_64(0x1.52efd0cd80497p-1), 0x1.e095382100a01p-1);
345 try expectEqual(expm1_64(-0x1.a05cc754481d1p-2), -0x1.561c3e0582be6p-2);
346 try expectEqual(expm1_64(0x1.1f9ef934745cbp-1), 0x1.81ec4cd4d4a8fp-1);
347 try expectEqual(expm1_64(0x1.8c5db097f7442p-1), 0x1.2b3363a944bf7p+0);
348 try expectEqual(expm1_64(-0x1.5b86ea8118a0ep-1), -0x1.f8951aebffbafp-2);
349}
350
351test "expm1_64() boundary" {
352 try expectEqual(expm1_64(0x1.62e42fefa39efp+9), 0x1.fffffffffff2ap+1023); // Last value before result is inf
353 try expectEqual(expm1_64(0x1.62e42fefa39f0p+9), math.inf(f64)); // First value that gives inf
354 try expectEqual(expm1_64(0x1.fffffffffffffp+1023), math.inf(f64)); // Max input value
355 try expectEqual(expm1_64(0x1p-1074), 0x1p-1074); // Min positive input value
356 try expectEqual(expm1_64(-0x1p-1074), -0x1p-1074); // Min negative input value
357 try expectEqual(expm1_64(0x1p-1022), 0x1p-1022); // First positive subnormal input
358 try expectEqual(expm1_64(-0x1p-1022), -0x1p-1022); // First negative subnormal input
359 try expectEqual(expm1_64(0x1.fffffffffffffp-1021), 0x1.fffffffffffffp-1021); // Last positive value before subnormal
360 try expectEqual(expm1_64(-0x1.fffffffffffffp-1021), -0x1.fffffffffffffp-1021); // Last negative value before subnormal
361 try expectEqual(expm1_64(-0x1.2b708872320e1p+5), -0x1.fffffffffffffp-1); // Last value before result is -1
362 try expectEqual(expm1_64(-0x1.2b708872320e2p+5), -1); // First value where result is -1
363}