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