log.zig (8304B) - 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/lnf.c 5 //! https://git.musl-libc.org/cgit/musl/tree/src/math/ln.c 6 7 const std = @import("std"); 8 const builtin = @import("builtin"); 9 const math = std.math; 10 const expect = std.testing.expect; 11 const expectEqual = std.testing.expectEqual; 12 const arch = builtin.cpu.arch; 13 const common = @import("common.zig"); 14 15 pub const panic = common.panic; 16 17 comptime { 18 @export(&__logh, .{ .name = "__logh", .linkage = common.linkage, .visibility = common.visibility }); 19 @export(&logf, .{ .name = "logf", .linkage = common.linkage, .visibility = common.visibility }); 20 @export(&log, .{ .name = "log", .linkage = common.linkage, .visibility = common.visibility }); 21 @export(&__logx, .{ .name = "__logx", .linkage = common.linkage, .visibility = common.visibility }); 22 if (common.want_ppc_abi) { 23 @export(&logq, .{ .name = "logf128", .linkage = common.linkage, .visibility = common.visibility }); 24 } 25 @export(&logq, .{ .name = "logq", .linkage = common.linkage, .visibility = common.visibility }); 26 @export(&logl, .{ .name = "logl", .linkage = common.linkage, .visibility = common.visibility }); 27 } 28 29 pub fn __logh(a: f16) callconv(.c) f16 { 30 // TODO: more efficient implementation 31 return @floatCast(logf(a)); 32 } 33 34 pub fn logf(x_: f32) callconv(.c) f32 { 35 const ln2_hi: f32 = 6.9313812256e-01; 36 const ln2_lo: f32 = 9.0580006145e-06; 37 const Lg1: f32 = 0xaaaaaa.0p-24; 38 const Lg2: f32 = 0xccce13.0p-25; 39 const Lg3: f32 = 0x91e9ee.0p-25; 40 const Lg4: f32 = 0xf89e26.0p-26; 41 42 var x = x_; 43 var ix: u32 = @bitCast(x); 44 var k: i32 = 0; 45 46 // x < 2^(-126) 47 if (ix < 0x00800000 or ix >> 31 != 0) { 48 // log(+-0) = -inf 49 if (ix << 1 == 0) { 50 return -math.inf(f32); 51 } 52 // log(-#) = nan 53 if (ix >> 31 != 0) { 54 return math.nan(f32); 55 } 56 57 // subnormal, scale x 58 k -= 25; 59 x *= 0x1.0p25; 60 ix = @bitCast(x); 61 } else if (ix >= 0x7F800000) { 62 return x; 63 } else if (ix == 0x3F800000) { 64 return 0; 65 } 66 67 // x into [sqrt(2) / 2, sqrt(2)] 68 ix += 0x3F800000 - 0x3F3504F3; 69 k += @as(i32, @intCast(ix >> 23)) - 0x7F; 70 ix = (ix & 0x007FFFFF) + 0x3F3504F3; 71 x = @bitCast(ix); 72 73 const f = x - 1.0; 74 const s = f / (2.0 + f); 75 const z = s * s; 76 const w = z * z; 77 const t1 = w * (Lg2 + w * Lg4); 78 const t2 = z * (Lg1 + w * Lg3); 79 const R = t2 + t1; 80 const hfsq = 0.5 * f * f; 81 const dk: f32 = @floatFromInt(k); 82 83 return s * (hfsq + R) + dk * ln2_lo - hfsq + f + dk * ln2_hi; 84 } 85 86 pub fn log(x_: f64) callconv(.c) f64 { 87 const ln2_hi: f64 = 6.93147180369123816490e-01; 88 const ln2_lo: f64 = 1.90821492927058770002e-10; 89 const Lg1: f64 = 6.666666666666735130e-01; 90 const Lg2: f64 = 3.999999999940941908e-01; 91 const Lg3: f64 = 2.857142874366239149e-01; 92 const Lg4: f64 = 2.222219843214978396e-01; 93 const Lg5: f64 = 1.818357216161805012e-01; 94 const Lg6: f64 = 1.531383769920937332e-01; 95 const Lg7: f64 = 1.479819860511658591e-01; 96 97 var x = x_; 98 var ix: u64 = @bitCast(x); 99 var hx: u32 = @intCast(ix >> 32); 100 var k: i32 = 0; 101 102 if (hx < 0x00100000 or hx >> 31 != 0) { 103 // log(+-0) = -inf 104 if (ix << 1 == 0) { 105 return -math.inf(f64); 106 } 107 // log(-#) = nan 108 if (hx >> 31 != 0) { 109 return math.nan(f64); 110 } 111 112 // subnormal, scale x 113 k -= 54; 114 x *= 0x1p54; 115 hx = @intCast(@as(u64, @bitCast(x)) >> 32); 116 } else if (hx >= 0x7FF00000) { 117 return x; 118 } else if (hx == 0x3FF00000 and ix << 32 == 0) { 119 return 0; 120 } 121 122 // x into [sqrt(2) / 2, sqrt(2)] 123 hx += 0x3FF00000 - 0x3FE6A09E; 124 k += @as(i32, @intCast(hx >> 20)) - 0x3FF; 125 hx = (hx & 0x000FFFFF) + 0x3FE6A09E; 126 ix = (@as(u64, hx) << 32) | (ix & 0xFFFFFFFF); 127 x = @bitCast(ix); 128 129 const f = x - 1.0; 130 const hfsq = 0.5 * f * f; 131 const s = f / (2.0 + f); 132 const z = s * s; 133 const w = z * z; 134 const t1 = w * (Lg2 + w * (Lg4 + w * Lg6)); 135 const t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7))); 136 const R = t2 + t1; 137 const dk: f64 = @floatFromInt(k); 138 139 return s * (hfsq + R) + dk * ln2_lo - hfsq + f + dk * ln2_hi; 140 } 141 142 pub fn __logx(a: f80) callconv(.c) f80 { 143 // TODO: more efficient implementation 144 return @floatCast(logq(a)); 145 } 146 147 pub fn logq(a: f128) callconv(.c) f128 { 148 // TODO: more correct implementation 149 return log(@floatCast(a)); 150 } 151 152 pub fn logl(x: c_longdouble) callconv(.c) c_longdouble { 153 switch (@typeInfo(c_longdouble).float.bits) { 154 16 => return __logh(x), 155 32 => return logf(x), 156 64 => return log(x), 157 80 => return __logx(x), 158 128 => return logq(x), 159 else => @compileError("unreachable"), 160 } 161 } 162 163 test "logf() special" { 164 try expectEqual(logf(0.0), -math.inf(f32)); 165 try expectEqual(logf(-0.0), -math.inf(f32)); 166 try expect(math.isPositiveZero(logf(1.0))); 167 try expectEqual(logf(math.e), 1.0); 168 try expectEqual(logf(math.inf(f32)), math.inf(f32)); 169 try expect(math.isNan(logf(-1.0))); 170 try expect(math.isNan(logf(-math.inf(f32)))); 171 try expect(math.isNan(logf(math.nan(f32)))); 172 try expect(math.isNan(logf(math.snan(f32)))); 173 } 174 175 test "logf() sanity" { 176 try expect(math.isNan(logf(-0x1.0223a0p+3))); 177 try expectEqual(logf(0x1.161868p+2), 0x1.7815b0p+0); 178 try expect(math.isNan(logf(-0x1.0c34b4p+3))); 179 try expect(math.isNan(logf(-0x1.a206f0p+2))); 180 try expectEqual(logf(0x1.288bbcp+3), 0x1.1cfcd6p+1); 181 try expectEqual(logf(0x1.52efd0p-1), -0x1.a6694cp-2); 182 try expect(math.isNan(logf(-0x1.a05cc8p-2))); 183 try expectEqual(logf(0x1.1f9efap-1), -0x1.2742bap-1); 184 try expectEqual(logf(0x1.8c5db0p-1), -0x1.062160p-2); 185 try expect(math.isNan(logf(-0x1.5b86eap-1))); 186 } 187 188 test "logf() boundary" { 189 try expectEqual(logf(0x1.fffffep+127), 0x1.62e430p+6); // Max input value 190 try expectEqual(logf(0x1p-149), -0x1.9d1da0p+6); // Min positive input value 191 try expect(math.isNan(logf(-0x1p-149))); // Min negative input value 192 try expectEqual(logf(0x1.000002p+0), 0x1.fffffep-24); // Last value before result reaches +0 193 try expectEqual(logf(0x1.fffffep-1), -0x1p-24); // Last value before result reaches -0 194 try expectEqual(logf(0x1p-126), -0x1.5d58a0p+6); // First subnormal 195 try expect(math.isNan(logf(-0x1p-126))); // First negative subnormal 196 } 197 198 test "log() special" { 199 try expectEqual(log(0.0), -math.inf(f64)); 200 try expectEqual(log(-0.0), -math.inf(f64)); 201 try expect(math.isPositiveZero(log(1.0))); 202 try expectEqual(log(math.e), 1.0); 203 try expectEqual(log(math.inf(f64)), math.inf(f64)); 204 try expect(math.isNan(log(-1.0))); 205 try expect(math.isNan(log(-math.inf(f64)))); 206 try expect(math.isNan(log(math.nan(f64)))); 207 try expect(math.isNan(log(math.snan(f64)))); 208 } 209 210 test "log() sanity" { 211 try expect(math.isNan(log(-0x1.02239f3c6a8f1p+3))); 212 try expectEqual(log(0x1.161868e18bc67p+2), 0x1.7815b08f99c65p+0); 213 try expect(math.isNan(log(-0x1.0c34b3e01e6e7p+3))); 214 try expect(math.isNan(log(-0x1.a206f0a19dcc4p+2))); 215 try expectEqual(log(0x1.288bbb0d6a1e6p+3), 0x1.1cfcd53d72604p+1); 216 try expectEqual(log(0x1.52efd0cd80497p-1), -0x1.a6694a4a85621p-2); 217 try expect(math.isNan(log(-0x1.a05cc754481d1p-2))); 218 try expectEqual(log(0x1.1f9ef934745cbp-1), -0x1.2742bc03d02ddp-1); 219 try expectEqual(log(0x1.8c5db097f7442p-1), -0x1.06215de4a3f92p-2); 220 try expect(math.isNan(log(-0x1.5b86ea8118a0ep-1))); 221 } 222 223 test "log() boundary" { 224 try expectEqual(log(0x1.fffffffffffffp+1023), 0x1.62e42fefa39efp+9); // Max input value 225 try expectEqual(log(0x1p-1074), -0x1.74385446d71c3p+9); // Min positive input value 226 try expect(math.isNan(log(-0x1p-1074))); // Min negative input value 227 try expectEqual(log(0x1.0000000000001p+0), 0x1.fffffffffffffp-53); // Last value before result reaches +0 228 try expectEqual(log(0x1.fffffffffffffp-1), -0x1p-53); // Last value before result reaches -0 229 try expectEqual(log(0x1p-1022), -0x1.6232bdd7abcd2p+9); // First subnormal 230 try expect(math.isNan(log(-0x1p-1022))); // First negative subnormal 231 }