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/expf.c
5// https://git.musl-libc.org/cgit/musl/tree/src/math/exp.c
6
7const builtin = @import("builtin");
8const arch = builtin.cpu.arch;
9
10const std = @import("std");
11const math = std.math;
12const mem = std.mem;
13const expect = std.testing.expect;
14const expectEqual = std.testing.expectEqual;
15
16const compiler_rt = @import("../compiler_rt.zig");
17const symbol = compiler_rt.symbol;
18
19comptime {
20 symbol(&__exph, "__exph");
21 symbol(&expf, "expf");
22 symbol(&exp, "exp");
23 symbol(&__expx, "__expx");
24 symbol(&expq, "expf128");
25 symbol(&expl, "expl");
26}
27
28fn __exph(x: compiler_rt.f16.Abi) callconv(.c) compiler_rt.f16.Abi {
29 return compiler_rt.f16.toAbi(exp_f16(compiler_rt.f16.fromAbi(x)));
30}
31pub fn exp_f16(x: f16) f16 {
32 // TODO: more efficient implementation
33 return @floatCast(exp_f32(x));
34}
35
36fn expf(x: compiler_rt.f32.Abi) callconv(.c) compiler_rt.f32.Abi {
37 return compiler_rt.f32.toAbi(exp_f32(compiler_rt.f32.fromAbi(x)));
38}
39pub fn exp_f32(x_: f32) f32 {
40 const half = [_]f32{ 0.5, -0.5 };
41 const ln2hi = 6.9314575195e-1;
42 const ln2lo = 1.4286067653e-6;
43 const invln2 = 1.4426950216e+0;
44 const P1 = 1.6666625440e-1;
45 const P2 = -2.7667332906e-3;
46
47 var x = x_;
48 var hx: u32 = @bitCast(x);
49 const sign: i32 = @intCast(hx >> 31);
50 hx &= 0x7FFFFFFF;
51
52 if (math.isNan(x)) {
53 return x;
54 }
55
56 // |x| >= -87.33655 or nan
57 if (hx >= 0x42AEAC50) {
58 // nan
59 if (hx > 0x7F800000) {
60 return x;
61 }
62 // x >= 88.722839
63 if (hx >= 0x42b17218 and sign == 0) {
64 return x * 0x1.0p127;
65 }
66 if (sign != 0) {
67 if (compiler_rt.want_float_exceptions) mem.doNotOptimizeAway(-0x1.0p-149 / x); // overflow
68 // x <= -103.972084
69 if (hx >= 0x42CFF1B5) {
70 return 0;
71 }
72 }
73 }
74
75 var k: i32 = undefined;
76 var hi: f32 = undefined;
77 var lo: f32 = undefined;
78
79 // |x| > 0.5 * ln2
80 if (hx > 0x3EB17218) {
81 // |x| > 1.5 * ln2
82 if (hx > 0x3F851592) {
83 k = @intFromFloat(invln2 * x + half[@intCast(sign)]);
84 } else {
85 k = 1 - sign - sign;
86 }
87
88 const fk: f32 = @floatFromInt(k);
89 hi = x - fk * ln2hi;
90 lo = fk * ln2lo;
91 x = hi - lo;
92 }
93 // |x| > 2^(-14)
94 else if (hx > 0x39000000) {
95 k = 0;
96 hi = x;
97 lo = 0;
98 } else {
99 if (compiler_rt.want_float_exceptions) mem.doNotOptimizeAway(0x1.0p127 + x); // inexact
100 return 1 + x;
101 }
102
103 const xx = x * x;
104 const c = x - xx * (P1 + xx * P2);
105 const y = 1 + (x * c / (2 - c) - lo + hi);
106
107 if (k == 0) {
108 return y;
109 } else {
110 return math.scalbn(y, k);
111 }
112}
113
114fn exp(x: compiler_rt.f64.Abi) callconv(.c) compiler_rt.f64.Abi {
115 return compiler_rt.f64.toAbi(exp_f64(compiler_rt.f64.fromAbi(x)));
116}
117pub fn exp_f64(x_: f64) f64 {
118 const half = [_]f64{ 0.5, -0.5 };
119 const ln2hi: f64 = 6.93147180369123816490e-01;
120 const ln2lo: f64 = 1.90821492927058770002e-10;
121 const invln2: f64 = 1.44269504088896338700e+00;
122 const P1: f64 = 1.66666666666666019037e-01;
123 const P2: f64 = -2.77777777770155933842e-03;
124 const P3: f64 = 6.61375632143793436117e-05;
125 const P4: f64 = -1.65339022054652515390e-06;
126 const P5: f64 = 4.13813679705723846039e-08;
127
128 var x = x_;
129 const ux: u64 = @bitCast(x);
130 var hx = ux >> 32;
131 const sign: i32 = @intCast(hx >> 31);
132 hx &= 0x7FFFFFFF;
133
134 if (math.isNan(x)) {
135 return x;
136 }
137
138 // |x| >= 708.39 or nan
139 if (hx >= 0x4086232B) {
140 // nan
141 if (hx > 0x7FF00000) {
142 return x;
143 }
144 if (x > 709.782712893383973096) {
145 // overflow if x != inf
146 return if (compiler_rt.want_float_exceptions) x * 0x1p1023 else std.math.inf(f64);
147 }
148 if (x < -708.39641853226410622) {
149 // underflow if x != -inf
150 if (compiler_rt.want_float_exceptions) mem.doNotOptimizeAway(-0x0.0000000000001p-1022 / x);
151 if (x < -745.13321910194110842) {
152 return 0;
153 }
154 }
155 }
156
157 // argument reduction
158 var k: i32 = undefined;
159 var hi: f64 = undefined;
160 var lo: f64 = undefined;
161
162 // |x| > 0.5 * ln2
163 if (hx > 0x3FD62E42) {
164 // |x| >= 1.5 * ln2
165 if (hx > 0x3FF0A2B2) {
166 k = @intFromFloat(invln2 * x + half[@intCast(sign)]);
167 } else {
168 k = 1 - sign - sign;
169 }
170
171 const dk: f64 = @floatFromInt(k);
172 hi = x - dk * ln2hi;
173 lo = dk * ln2lo;
174 x = hi - lo;
175 }
176 // |x| > 2^(-28)
177 else if (hx > 0x3E300000) {
178 k = 0;
179 hi = x;
180 lo = 0;
181 } else {
182 // inexact if x != 0
183 if (compiler_rt.want_float_exceptions) mem.doNotOptimizeAway(0x1.0p1023 + x);
184 return 1 + x;
185 }
186
187 const xx = x * x;
188 const c = x - xx * (P1 + xx * (P2 + xx * (P3 + xx * (P4 + xx * P5))));
189 const y = 1 + (x * c / (2 - c) - lo + hi);
190
191 if (k == 0) {
192 return y;
193 } else {
194 return math.scalbn(y, k);
195 }
196}
197
198fn __expx(x: compiler_rt.f80.Abi) callconv(.c) compiler_rt.f80.Abi {
199 return compiler_rt.f80.toAbi(exp_f80(compiler_rt.f80.fromAbi(x)));
200}
201pub fn exp_f80(x: f80) f80 {
202 // TODO: more efficient implementation
203 return @floatCast(exp_f128(x));
204}
205
206fn expq(x: compiler_rt.f128.Abi) callconv(.c) compiler_rt.f128.Abi {
207 return compiler_rt.f128.toAbi(exp_f128(compiler_rt.f128.fromAbi(x)));
208}
209pub const exp_f128 = @import("exp_f128.zig").exp;
210
211pub fn expl(x: c_longdouble) callconv(.c) c_longdouble {
212 switch (@typeInfo(c_longdouble).float.bits) {
213 64 => return exp_f64(x),
214 80 => return exp_f80(x),
215 128 => return exp_f128(x),
216 else => comptime unreachable,
217 }
218}
219
220test "expf() special" {
221 try expectEqual(exp_f32(0.0), 1.0);
222 try expectEqual(exp_f32(-0.0), 1.0);
223 try expectEqual(exp_f32(1.0), math.e);
224 try expectEqual(exp_f32(math.ln2), 2.0);
225 try expectEqual(exp_f32(math.inf(f32)), math.inf(f32));
226 try expect(math.isPositiveZero(exp_f32(-math.inf(f32))));
227 try expect(math.isNan(exp_f32(math.nan(f32))));
228 try expect(math.isNan(exp_f32(math.snan(f32))));
229}
230
231test "expf() sanity" {
232 try expectEqual(exp_f32(-0x1.0223a0p+3), 0x1.490320p-12);
233 try expectEqual(exp_f32(0x1.161868p+2), 0x1.34712ap+6);
234 try expectEqual(exp_f32(-0x1.0c34b4p+3), 0x1.e06b1ap-13);
235 try expectEqual(exp_f32(-0x1.a206f0p+2), 0x1.7dd484p-10);
236 try expectEqual(exp_f32(0x1.288bbcp+3), 0x1.4abc80p+13);
237 try expectEqual(exp_f32(0x1.52efd0p-1), 0x1.f04a9cp+0);
238 try expectEqual(exp_f32(-0x1.a05cc8p-2), 0x1.54f1e0p-1);
239 try expectEqual(exp_f32(0x1.1f9efap-1), 0x1.c0f628p+0);
240 try expectEqual(exp_f32(0x1.8c5db0p-1), 0x1.1599b2p+1);
241 try expectEqual(exp_f32(-0x1.5b86eap-1), 0x1.03b572p-1);
242 try expectEqual(exp_f32(-0x1.57f25cp+2), 0x1.2fbea2p-8);
243 try expectEqual(exp_f32(0x1.c7d310p+3), 0x1.76eefp+20);
244 try expectEqual(exp_f32(0x1.19be70p+4), 0x1.52d3dep+25);
245 try expectEqual(exp_f32(-0x1.ab6d70p+3), 0x1.a88adep-20);
246 try expectEqual(exp_f32(-0x1.5ac18ep+2), 0x1.22b328p-8);
247 try expectEqual(exp_f32(-0x1.925982p-1), 0x1.d2acc0p-2);
248 try expectEqual(exp_f32(0x1.7221cep+3), 0x1.9c2ceap+16);
249 try expectEqual(exp_f32(0x1.11a0d4p+4), 0x1.980ee6p+24);
250 try expectEqual(exp_f32(-0x1.ae41a2p+1), 0x1.1c28d0p-5);
251 try expectEqual(exp_f32(-0x1.329154p+4), 0x1.47ef94p-28);
252}
253
254test "expf() boundary" {
255 try expectEqual(exp_f32(0x1.62e42ep+6), 0x1.ffff08p+127); // The last value before the result gets infinite
256 try expectEqual(exp_f32(0x1.62e430p+6), math.inf(f32)); // The first value that gives inf
257 try expectEqual(exp_f32(0x1.fffffep+127), math.inf(f32)); // Max input value
258 try expectEqual(exp_f32(0x1p-149), 1.0); // Min positive input value
259 try expectEqual(exp_f32(-0x1p-149), 1.0); // Min negative input value
260 try expectEqual(exp_f32(0x1p-126), 1.0); // First positive subnormal input
261 try expectEqual(exp_f32(-0x1p-126), 1.0); // First negative subnormal input
262 try expectEqual(exp_f32(-0x1.9fe368p+6), 0x1p-149); // The last value before the result flushes to zero
263 try expectEqual(exp_f32(-0x1.9fe36ap+6), 0.0); // The first value at which the result flushes to zero
264 try expectEqual(exp_f32(-0x1.5d589ep+6), 0x1.00004cp-126); // The last value before the result flushes to subnormal
265 try expectEqual(exp_f32(-0x1.5d58a0p+6), 0x1.ffff98p-127); // The first value for which the result flushes to subnormal
266}
267
268test "exp() special" {
269 try expectEqual(exp_f64(0.0), 1.0);
270 try expectEqual(exp_f64(-0.0), 1.0);
271 // TODO: Accuracy error - off in the last bit in 64-bit, disagreeing with GCC
272 // try expectEqual(exp(1.0), math.e);
273 try expectEqual(exp_f64(math.ln2), 2.0);
274 try expectEqual(exp_f64(math.inf(f64)), math.inf(f64));
275 try expect(math.isPositiveZero(exp_f64(-math.inf(f64))));
276 try expect(math.isNan(exp_f64(math.nan(f64))));
277 try expect(math.isNan(exp_f64(math.snan(f64))));
278}
279
280test "exp() sanity" {
281 try expectEqual(exp_f64(-0x1.02239f3c6a8f1p+3), 0x1.490327ea61235p-12);
282 try expectEqual(exp_f64(0x1.161868e18bc67p+2), 0x1.34712ed238c04p+6);
283 try expectEqual(exp_f64(-0x1.0c34b3e01e6e7p+3), 0x1.e06b1b6c18e64p-13);
284 try expectEqual(exp_f64(-0x1.a206f0a19dcc4p+2), 0x1.7dd47f810e68cp-10);
285 try expectEqual(exp_f64(0x1.288bbb0d6a1e6p+3), 0x1.4abc77496e07ep+13);
286 try expectEqual(exp_f64(0x1.52efd0cd80497p-1), 0x1.f04a9c1080500p+0);
287 try expectEqual(exp_f64(-0x1.a05cc754481d1p-2), 0x1.54f1e0fd3ea0dp-1);
288 try expectEqual(exp_f64(0x1.1f9ef934745cbp-1), 0x1.c0f6266a6a547p+0);
289 try expectEqual(exp_f64(0x1.8c5db097f7442p-1), 0x1.1599b1d4a25fbp+1);
290 try expectEqual(exp_f64(-0x1.5b86ea8118a0ep-1), 0x1.03b5728a00229p-1);
291 try expectEqual(exp_f64(-0x1.57f25b2b5006dp+2), 0x1.2fbea6a01cab9p-8);
292 try expectEqual(exp_f64(0x1.c7d30fb825911p+3), 0x1.76eeed45a0634p+20);
293 try expectEqual(exp_f64(0x1.19be709de7505p+4), 0x1.52d3eb7be6844p+25);
294 try expectEqual(exp_f64(-0x1.ab6d6fba96889p+3), 0x1.a88ae12f985d6p-20);
295 try expectEqual(exp_f64(-0x1.5ac18e27084ddp+2), 0x1.22b327da9cca6p-8);
296 try expectEqual(exp_f64(-0x1.925981b093c41p-1), 0x1.d2acc046b55f7p-2);
297 try expectEqual(exp_f64(0x1.7221cd18455f5p+3), 0x1.9c2cde8699cfbp+16);
298 try expectEqual(exp_f64(0x1.11a0d4a51b239p+4), 0x1.980ef612ff182p+24);
299 try expectEqual(exp_f64(-0x1.ae41a1079de4dp+1), 0x1.1c28d16bb3222p-5);
300 try expectEqual(exp_f64(-0x1.329153103b871p+4), 0x1.47efa6ddd0d22p-28);
301}
302
303test "exp() boundary" {
304 try expectEqual(exp_f64(0x1.62e42fefa39efp+9), 0x1.fffffffffff2ap+1023); // The last value before the result gets infinite
305 try expectEqual(exp_f64(0x1.62e42fefa39f0p+9), math.inf(f64)); // The first value that gives inf
306 try expectEqual(exp_f64(0x1.fffffffffffffp+1023), math.inf(f64)); // Max input value
307 try expectEqual(exp_f64(0x1p-1074), 1.0); // Min positive input value
308 try expectEqual(exp_f64(-0x1p-1074), 1.0); // Min negative input value
309 try expectEqual(exp_f64(0x1p-1022), 1.0); // First positive subnormal input
310 try expectEqual(exp_f64(-0x1p-1022), 1.0); // First negative subnormal input
311 try expectEqual(exp_f64(-0x1.74910d52d3051p+9), 0x1p-1074); // The last value before the result flushes to zero
312 try expectEqual(exp_f64(-0x1.74910d52d3052p+9), 0.0); // The first value at which the result flushes to zero
313 try expectEqual(exp_f64(-0x1.6232bdd7abcd2p+9), 0x1.000000000007cp-1022); // The last value before the result flushes to subnormal
314 try expectEqual(exp_f64(-0x1.6232bdd7abcd3p+9), 0x1.ffffffffffcf8p-1023); // The first value for which the result flushes to subnormal
315}