| author | |
| committer | |
| log | 50422d5c37e85f08da7579cbffb4c0281a5cb2a8 |
| tree | f9f64279b97584ae6ca3abacc05e0621a30e9682 |
| parent | 962903e0fa222ee002355832e17144f66ec8d6ec |
| parent | 81097669270c9f525eb0556915bdbcc3492fcec3 |
Reviewed-on: https://codeberg.org/ziglang/zig/pulls/30648
Reviewed-by: Andrew Kelley <andrewrk@noreply.codeberg.org>
Reviewed-by: Alex Rønne Petersen <alex@alexrp.com>25 files changed, 357 insertions(+), 1335 deletions(-)
lib/compiler_rt/exp.zig+3-6| ... | @@ -134,14 +134,11 @@ pub fn exp(x_: f64) callconv(.c) f64 { | ... | @@ -134,14 +134,11 @@ pub fn exp(x_: f64) callconv(.c) f64 { |
| 134 | } | 134 | } |
| 135 | if (x > 709.782712893383973096) { | 135 | if (x > 709.782712893383973096) { |
| 136 | // overflow if x != inf | 136 | // overflow if x != inf |
| 137 | if (!math.isInf(x)) { | 137 | return if (common.want_float_exceptions) x * 0x1p1023 else std.math.inf(f64); |
| 138 | math.raiseOverflow(); | ||
| 139 | } | ||
| 140 | return math.inf(f64); | ||
| 141 | } | 138 | } |
| 142 | if (x < -708.39641853226410622) { | 139 | if (x < -708.39641853226410622) { |
| 143 | // underflow if x != -inf | 140 | // underflow if x != -inf |
| 144 | // if (common.want_float_exceptions) mem.doNotOptimizeAway(@as(f32, -0x1.0p-149 / x)); | 141 | if (common.want_float_exceptions) mem.doNotOptimizeAway(-0x0.0000000000001p-1022 / x); |
| 145 | if (x < -745.13321910194110842) { | 142 | if (x < -745.13321910194110842) { |
| 146 | return 0; | 143 | return 0; |
| 147 | } | 144 | } |
| ... | @@ -174,7 +171,7 @@ pub fn exp(x_: f64) callconv(.c) f64 { | ... | @@ -174,7 +171,7 @@ pub fn exp(x_: f64) callconv(.c) f64 { |
| 174 | lo = 0; | 171 | lo = 0; |
| 175 | } else { | 172 | } else { |
| 176 | // inexact if x != 0 | 173 | // inexact if x != 0 |
| 177 | // if (common.want_float_exceptions) mem.doNotOptimizeAway(0x1.0p1023 + x); | 174 | if (common.want_float_exceptions) mem.doNotOptimizeAway(0x1.0p1023 + x); |
| 178 | return 1 + x; | 175 | return 1 + x; |
| 179 | } | 176 | } |
| 180 | 177 |
lib/compiler_rt/exp2.zig+1-2| ... | @@ -108,8 +108,7 @@ pub fn exp2(x: f64) callconv(.c) f64 { | ... | @@ -108,8 +108,7 @@ pub fn exp2(x: f64) callconv(.c) f64 { |
| 108 | if (ix >= 0x408FF000) { | 108 | if (ix >= 0x408FF000) { |
| 109 | // x >= 1024 or nan | 109 | // x >= 1024 or nan |
| 110 | if (ix >= 0x40900000 and ux >> 63 == 0) { | 110 | if (ix >= 0x40900000 and ux >> 63 == 0) { |
| 111 | math.raiseOverflow(); | 111 | return if (common.want_float_exceptions) x * 0x1p1023 else std.math.inf(f64); |
| 112 | return math.inf(f64); | ||
| 113 | } | 112 | } |
| 114 | // -inf or -nan | 113 | // -inf or -nan |
| 115 | if (ix >= 0x7FF00000) { | 114 | if (ix >= 0x7FF00000) { |
lib/compiler_rt/log.zig+345-50| ... | @@ -1,8 +1,8 @@ | ... | @@ -1,8 +1,8 @@ |
| 1 | //! Ported from musl, which is licensed under the MIT license: | 1 | //! Ported from musl, which is licensed under the MIT license: |
| 2 | //! https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT | 2 | //! https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT |
| 3 | //! | 3 | //! |
| 4 | //! https://git.musl-libc.org/cgit/musl/tree/src/math/lnf.c | 4 | //! https://git.musl-libc.org/cgit/musl/tree/src/math/logf.c?h=2d7d05f031e014068a61d3076c6178513395d2ae |
| 5 | //! https://git.musl-libc.org/cgit/musl/tree/src/math/ln.c | 5 | //! https://git.musl-libc.org/cgit/musl/tree/src/math/log.c?h=1b76ff0767d01df72f692806ee5adee13c67ef88 |
| 6 | 6 | ||
| 7 | const std = @import("std"); | 7 | const std = @import("std"); |
| 8 | const builtin = @import("builtin"); | 8 | const builtin = @import("builtin"); |
| ... | @@ -45,11 +45,11 @@ pub fn logf(x_: f32) callconv(.c) f32 { | ... | @@ -45,11 +45,11 @@ pub fn logf(x_: f32) callconv(.c) f32 { |
| 45 | if (ix < 0x00800000 or ix >> 31 != 0) { | 45 | if (ix < 0x00800000 or ix >> 31 != 0) { |
| 46 | // log(+-0) = -inf | 46 | // log(+-0) = -inf |
| 47 | if (ix << 1 == 0) { | 47 | if (ix << 1 == 0) { |
| 48 | return -math.inf(f32); | 48 | return if (common.want_float_exceptions) -1 / (x * x) else -std.math.inf(f64); |
| 49 | } | 49 | } |
| 50 | // log(-#) = nan | 50 | // log(-#) = nan |
| 51 | if (ix >> 31 != 0) { | 51 | if (ix >> 31 != 0) { |
| 52 | return math.nan(f32); | 52 | return if (common.want_float_exceptions) (x - x) / 0.0 else math.nan(f64); |
| 53 | } | 53 | } |
| 54 | 54 | ||
| 55 | // subnormal, scale x | 55 | // subnormal, scale x |
| ... | @@ -81,60 +81,355 @@ pub fn logf(x_: f32) callconv(.c) f32 { | ... | @@ -81,60 +81,355 @@ pub fn logf(x_: f32) callconv(.c) f32 { |
| 81 | return s * (hfsq + R) + dk * ln2_lo - hfsq + f + dk * ln2_hi; | 81 | return s * (hfsq + R) + dk * ln2_lo - hfsq + f + dk * ln2_hi; |
| 82 | } | 82 | } |
| 83 | 83 | ||
| 84 | pub fn log(x_: f64) callconv(.c) f64 { | 84 | pub fn log(x: f64) callconv(.c) f64 { |
| 85 | const ln2_hi: f64 = 6.93147180369123816490e-01; | 85 | const poly1 = [_]f64{ |
| 86 | const ln2_lo: f64 = 1.90821492927058770002e-10; | 86 | -0x1p-1, |
| 87 | const Lg1: f64 = 6.666666666666735130e-01; | 87 | 0x1.5555555555577p-2, |
| 88 | const Lg2: f64 = 3.999999999940941908e-01; | 88 | -0x1.ffffffffffdcbp-3, |
| 89 | const Lg3: f64 = 2.857142874366239149e-01; | 89 | 0x1.999999995dd0cp-3, |
| 90 | const Lg4: f64 = 2.222219843214978396e-01; | 90 | -0x1.55555556745a7p-3, |
| 91 | const Lg5: f64 = 1.818357216161805012e-01; | 91 | 0x1.24924a344de3p-3, |
| 92 | const Lg6: f64 = 1.531383769920937332e-01; | 92 | -0x1.fffffa4423d65p-4, |
| 93 | const Lg7: f64 = 1.479819860511658591e-01; | 93 | 0x1.c7184282ad6cap-4, |
| 94 | -0x1.999eb43b068ffp-4, | ||
| 95 | 0x1.78182f7afd085p-4, | ||
| 96 | -0x1.5521375d145cdp-4, | ||
| 97 | }; | ||
| 94 | 98 | ||
| 95 | var x = x_; | 99 | const poly = [_]f64{ |
| 96 | var ix: u64 = @bitCast(x); | 100 | -0x1.0000000000001p-1, |
| 97 | var hx: u32 = @intCast(ix >> 32); | 101 | 0x1.555555551305bp-2, |
| 98 | var k: i32 = 0; | 102 | -0x1.fffffffeb459p-3, |
| 103 | 0x1.999b324f10111p-3, | ||
| 104 | -0x1.55575e506c89fp-3, | ||
| 105 | }; | ||
| 99 | 106 | ||
| 100 | if (hx < 0x00100000 or hx >> 31 != 0) { | 107 | const tab = [128]struct { invc: f64, logc: f64 }{ |
| 101 | // log(+-0) = -inf | 108 | .{ .invc = 0x1.734f0c3e0de9fp+0, .logc = -0x1.7cc7f79e69000p-2 }, |
| 102 | if (ix << 1 == 0) { | 109 | .{ .invc = 0x1.713786a2ce91fp+0, .logc = -0x1.76feec20d0000p-2 }, |
| 103 | return -math.inf(f64); | 110 | .{ .invc = 0x1.6f26008fab5a0p+0, .logc = -0x1.713e31351e000p-2 }, |
| 104 | } | 111 | .{ .invc = 0x1.6d1a61f138c7dp+0, .logc = -0x1.6b85b38287800p-2 }, |
| 105 | // log(-#) = nan | 112 | .{ .invc = 0x1.6b1490bc5b4d1p+0, .logc = -0x1.65d5590807800p-2 }, |
| 106 | if (hx >> 31 != 0) { | 113 | .{ .invc = 0x1.69147332f0cbap+0, .logc = -0x1.602d076180000p-2 }, |
| 107 | return math.nan(f64); | 114 | .{ .invc = 0x1.6719f18224223p+0, .logc = -0x1.5a8ca86909000p-2 }, |
| 115 | .{ .invc = 0x1.6524f99a51ed9p+0, .logc = -0x1.54f4356035000p-2 }, | ||
| 116 | .{ .invc = 0x1.63356aa8f24c4p+0, .logc = -0x1.4f637c36b4000p-2 }, | ||
| 117 | .{ .invc = 0x1.614b36b9ddc14p+0, .logc = -0x1.49da7fda85000p-2 }, | ||
| 118 | .{ .invc = 0x1.5f66452c65c4cp+0, .logc = -0x1.445923989a800p-2 }, | ||
| 119 | .{ .invc = 0x1.5d867b5912c4fp+0, .logc = -0x1.3edf439b0b800p-2 }, | ||
| 120 | .{ .invc = 0x1.5babccb5b90dep+0, .logc = -0x1.396ce448f7000p-2 }, | ||
| 121 | .{ .invc = 0x1.59d61f2d91a78p+0, .logc = -0x1.3401e17bda000p-2 }, | ||
| 122 | .{ .invc = 0x1.5805612465687p+0, .logc = -0x1.2e9e2ef468000p-2 }, | ||
| 123 | .{ .invc = 0x1.56397cee76bd3p+0, .logc = -0x1.2941b3830e000p-2 }, | ||
| 124 | .{ .invc = 0x1.54725e2a77f93p+0, .logc = -0x1.23ec58cda8800p-2 }, | ||
| 125 | .{ .invc = 0x1.52aff42064583p+0, .logc = -0x1.1e9e129279000p-2 }, | ||
| 126 | .{ .invc = 0x1.50f22dbb2bddfp+0, .logc = -0x1.1956d2b48f800p-2 }, | ||
| 127 | .{ .invc = 0x1.4f38f4734ded7p+0, .logc = -0x1.141679ab9f800p-2 }, | ||
| 128 | .{ .invc = 0x1.4d843cfde2840p+0, .logc = -0x1.0edd094ef9800p-2 }, | ||
| 129 | .{ .invc = 0x1.4bd3ec078a3c8p+0, .logc = -0x1.09aa518db1000p-2 }, | ||
| 130 | .{ .invc = 0x1.4a27fc3e0258ap+0, .logc = -0x1.047e65263b800p-2 }, | ||
| 131 | .{ .invc = 0x1.4880524d48434p+0, .logc = -0x1.feb224586f000p-3 }, | ||
| 132 | .{ .invc = 0x1.46dce1b192d0bp+0, .logc = -0x1.f474a7517b000p-3 }, | ||
| 133 | .{ .invc = 0x1.453d9d3391854p+0, .logc = -0x1.ea4443d103000p-3 }, | ||
| 134 | .{ .invc = 0x1.43a2744b4845ap+0, .logc = -0x1.e020d44e9b000p-3 }, | ||
| 135 | .{ .invc = 0x1.420b54115f8fbp+0, .logc = -0x1.d60a22977f000p-3 }, | ||
| 136 | .{ .invc = 0x1.40782da3ef4b1p+0, .logc = -0x1.cc00104959000p-3 }, | ||
| 137 | .{ .invc = 0x1.3ee8f5d57fe8fp+0, .logc = -0x1.c202956891000p-3 }, | ||
| 138 | .{ .invc = 0x1.3d5d9a00b4ce9p+0, .logc = -0x1.b81178d811000p-3 }, | ||
| 139 | .{ .invc = 0x1.3bd60c010c12bp+0, .logc = -0x1.ae2c9ccd3d000p-3 }, | ||
| 140 | .{ .invc = 0x1.3a5242b75dab8p+0, .logc = -0x1.a45402e129000p-3 }, | ||
| 141 | .{ .invc = 0x1.38d22cd9fd002p+0, .logc = -0x1.9a877681df000p-3 }, | ||
| 142 | .{ .invc = 0x1.3755bc5847a1cp+0, .logc = -0x1.90c6d69483000p-3 }, | ||
| 143 | .{ .invc = 0x1.35dce49ad36e2p+0, .logc = -0x1.87120a645c000p-3 }, | ||
| 144 | .{ .invc = 0x1.34679984dd440p+0, .logc = -0x1.7d68fb4143000p-3 }, | ||
| 145 | .{ .invc = 0x1.32f5cceffcb24p+0, .logc = -0x1.73cb83c627000p-3 }, | ||
| 146 | .{ .invc = 0x1.3187775a10d49p+0, .logc = -0x1.6a39a9b376000p-3 }, | ||
| 147 | .{ .invc = 0x1.301c8373e3990p+0, .logc = -0x1.60b3154b7a000p-3 }, | ||
| 148 | .{ .invc = 0x1.2eb4ebb95f841p+0, .logc = -0x1.5737d76243000p-3 }, | ||
| 149 | .{ .invc = 0x1.2d50a0219a9d1p+0, .logc = -0x1.4dc7b8fc23000p-3 }, | ||
| 150 | .{ .invc = 0x1.2bef9a8b7fd2ap+0, .logc = -0x1.4462c51d20000p-3 }, | ||
| 151 | .{ .invc = 0x1.2a91c7a0c1babp+0, .logc = -0x1.3b08abc830000p-3 }, | ||
| 152 | .{ .invc = 0x1.293726014b530p+0, .logc = -0x1.31b996b490000p-3 }, | ||
| 153 | .{ .invc = 0x1.27dfa5757a1f5p+0, .logc = -0x1.2875490a44000p-3 }, | ||
| 154 | .{ .invc = 0x1.268b39b1d3bbfp+0, .logc = -0x1.1f3b9f879a000p-3 }, | ||
| 155 | .{ .invc = 0x1.2539d838ff5bdp+0, .logc = -0x1.160c8252ca000p-3 }, | ||
| 156 | .{ .invc = 0x1.23eb7aac9083bp+0, .logc = -0x1.0ce7f57f72000p-3 }, | ||
| 157 | .{ .invc = 0x1.22a012ba940b6p+0, .logc = -0x1.03cdc49fea000p-3 }, | ||
| 158 | .{ .invc = 0x1.2157996cc4132p+0, .logc = -0x1.f57bdbc4b8000p-4 }, | ||
| 159 | .{ .invc = 0x1.201201dd2fc9bp+0, .logc = -0x1.e370896404000p-4 }, | ||
| 160 | .{ .invc = 0x1.1ecf4494d480bp+0, .logc = -0x1.d17983ef94000p-4 }, | ||
| 161 | .{ .invc = 0x1.1d8f5528f6569p+0, .logc = -0x1.bf9674ed8a000p-4 }, | ||
| 162 | .{ .invc = 0x1.1c52311577e7cp+0, .logc = -0x1.adc79202f6000p-4 }, | ||
| 163 | .{ .invc = 0x1.1b17c74cb26e9p+0, .logc = -0x1.9c0c3e7288000p-4 }, | ||
| 164 | .{ .invc = 0x1.19e010c2c1ab6p+0, .logc = -0x1.8a646b372c000p-4 }, | ||
| 165 | .{ .invc = 0x1.18ab07bb670bdp+0, .logc = -0x1.78d01b3ac0000p-4 }, | ||
| 166 | .{ .invc = 0x1.1778a25efbcb6p+0, .logc = -0x1.674f145380000p-4 }, | ||
| 167 | .{ .invc = 0x1.1648d354c31dap+0, .logc = -0x1.55e0e6d878000p-4 }, | ||
| 168 | .{ .invc = 0x1.151b990275fddp+0, .logc = -0x1.4485cdea1e000p-4 }, | ||
| 169 | .{ .invc = 0x1.13f0ea432d24cp+0, .logc = -0x1.333d94d6aa000p-4 }, | ||
| 170 | .{ .invc = 0x1.12c8b7210f9dap+0, .logc = -0x1.22079f8c56000p-4 }, | ||
| 171 | .{ .invc = 0x1.11a3028ecb531p+0, .logc = -0x1.10e4698622000p-4 }, | ||
| 172 | .{ .invc = 0x1.107fbda8434afp+0, .logc = -0x1.ffa6c6ad20000p-5 }, | ||
| 173 | .{ .invc = 0x1.0f5ee0f4e6bb3p+0, .logc = -0x1.dda8d4a774000p-5 }, | ||
| 174 | .{ .invc = 0x1.0e4065d2a9fcep+0, .logc = -0x1.bbcece4850000p-5 }, | ||
| 175 | .{ .invc = 0x1.0d244632ca521p+0, .logc = -0x1.9a1894012c000p-5 }, | ||
| 176 | .{ .invc = 0x1.0c0a77ce2981ap+0, .logc = -0x1.788583302c000p-5 }, | ||
| 177 | .{ .invc = 0x1.0af2f83c636d1p+0, .logc = -0x1.5715e67d68000p-5 }, | ||
| 178 | .{ .invc = 0x1.09ddb98a01339p+0, .logc = -0x1.35c8a49658000p-5 }, | ||
| 179 | .{ .invc = 0x1.08cabaf52e7dfp+0, .logc = -0x1.149e364154000p-5 }, | ||
| 180 | .{ .invc = 0x1.07b9f2f4e28fbp+0, .logc = -0x1.e72c082eb8000p-6 }, | ||
| 181 | .{ .invc = 0x1.06ab58c358f19p+0, .logc = -0x1.a55f152528000p-6 }, | ||
| 182 | .{ .invc = 0x1.059eea5ecf92cp+0, .logc = -0x1.63d62cf818000p-6 }, | ||
| 183 | .{ .invc = 0x1.04949cdd12c90p+0, .logc = -0x1.228fb8caa0000p-6 }, | ||
| 184 | .{ .invc = 0x1.038c6c6f0ada9p+0, .logc = -0x1.c317b20f90000p-7 }, | ||
| 185 | .{ .invc = 0x1.02865137932a9p+0, .logc = -0x1.419355daa0000p-7 }, | ||
| 186 | .{ .invc = 0x1.0182427ea7348p+0, .logc = -0x1.81203c2ec0000p-8 }, | ||
| 187 | .{ .invc = 0x1.008040614b195p+0, .logc = -0x1.0040979240000p-9 }, | ||
| 188 | .{ .invc = 0x1.fe01ff726fa1ap-1, .logc = 0x1.feff384900000p-9 }, | ||
| 189 | .{ .invc = 0x1.fa11cc261ea74p-1, .logc = 0x1.7dc41353d0000p-7 }, | ||
| 190 | .{ .invc = 0x1.f6310b081992ep-1, .logc = 0x1.3cea3c4c28000p-6 }, | ||
| 191 | .{ .invc = 0x1.f25f63ceeadcdp-1, .logc = 0x1.b9fc114890000p-6 }, | ||
| 192 | .{ .invc = 0x1.ee9c8039113e7p-1, .logc = 0x1.1b0d8ce110000p-5 }, | ||
| 193 | .{ .invc = 0x1.eae8078cbb1abp-1, .logc = 0x1.58a5bd001c000p-5 }, | ||
| 194 | .{ .invc = 0x1.e741aa29d0c9bp-1, .logc = 0x1.95c8340d88000p-5 }, | ||
| 195 | .{ .invc = 0x1.e3a91830a99b5p-1, .logc = 0x1.d276aef578000p-5 }, | ||
| 196 | .{ .invc = 0x1.e01e009609a56p-1, .logc = 0x1.07598e598c000p-4 }, | ||
| 197 | .{ .invc = 0x1.dca01e577bb98p-1, .logc = 0x1.253f5e30d2000p-4 }, | ||
| 198 | .{ .invc = 0x1.d92f20b7c9103p-1, .logc = 0x1.42edd8b380000p-4 }, | ||
| 199 | .{ .invc = 0x1.d5cac66fb5ccep-1, .logc = 0x1.606598757c000p-4 }, | ||
| 200 | .{ .invc = 0x1.d272caa5ede9dp-1, .logc = 0x1.7da76356a0000p-4 }, | ||
| 201 | .{ .invc = 0x1.cf26e3e6b2ccdp-1, .logc = 0x1.9ab434e1c6000p-4 }, | ||
| 202 | .{ .invc = 0x1.cbe6da2a77902p-1, .logc = 0x1.b78c7bb0d6000p-4 }, | ||
| 203 | .{ .invc = 0x1.c8b266d37086dp-1, .logc = 0x1.d431332e72000p-4 }, | ||
| 204 | .{ .invc = 0x1.c5894bd5d5804p-1, .logc = 0x1.f0a3171de6000p-4 }, | ||
| 205 | .{ .invc = 0x1.c26b533bb9f8cp-1, .logc = 0x1.067152b914000p-3 }, | ||
| 206 | .{ .invc = 0x1.bf583eeece73fp-1, .logc = 0x1.147858292b000p-3 }, | ||
| 207 | .{ .invc = 0x1.bc4fd75db96c1p-1, .logc = 0x1.2266ecdca3000p-3 }, | ||
| 208 | .{ .invc = 0x1.b951e0c864a28p-1, .logc = 0x1.303d7a6c55000p-3 }, | ||
| 209 | .{ .invc = 0x1.b65e2c5ef3e2cp-1, .logc = 0x1.3dfc33c331000p-3 }, | ||
| 210 | .{ .invc = 0x1.b374867c9888bp-1, .logc = 0x1.4ba366b7a8000p-3 }, | ||
| 211 | .{ .invc = 0x1.b094b211d304ap-1, .logc = 0x1.5933928d1f000p-3 }, | ||
| 212 | .{ .invc = 0x1.adbe885f2ef7ep-1, .logc = 0x1.66acd2418f000p-3 }, | ||
| 213 | .{ .invc = 0x1.aaf1d31603da2p-1, .logc = 0x1.740f8ec669000p-3 }, | ||
| 214 | .{ .invc = 0x1.a82e63fd358a7p-1, .logc = 0x1.815c0f51af000p-3 }, | ||
| 215 | .{ .invc = 0x1.a5740ef09738bp-1, .logc = 0x1.8e92954f68000p-3 }, | ||
| 216 | .{ .invc = 0x1.a2c2a90ab4b27p-1, .logc = 0x1.9bb3602f84000p-3 }, | ||
| 217 | .{ .invc = 0x1.a01a01393f2d1p-1, .logc = 0x1.a8bed1c2c0000p-3 }, | ||
| 218 | .{ .invc = 0x1.9d79f24db3c1bp-1, .logc = 0x1.b5b515c01d000p-3 }, | ||
| 219 | .{ .invc = 0x1.9ae2505c7b190p-1, .logc = 0x1.c2967ccbcc000p-3 }, | ||
| 220 | .{ .invc = 0x1.9852ef297ce2fp-1, .logc = 0x1.cf635d5486000p-3 }, | ||
| 221 | .{ .invc = 0x1.95cbaeea44b75p-1, .logc = 0x1.dc1bd3446c000p-3 }, | ||
| 222 | .{ .invc = 0x1.934c69de74838p-1, .logc = 0x1.e8c01b8cfe000p-3 }, | ||
| 223 | .{ .invc = 0x1.90d4f2f6752e6p-1, .logc = 0x1.f5509c0179000p-3 }, | ||
| 224 | .{ .invc = 0x1.8e6528effd79dp-1, .logc = 0x1.00e6c121fb800p-2 }, | ||
| 225 | .{ .invc = 0x1.8bfce9fcc007cp-1, .logc = 0x1.071b80e93d000p-2 }, | ||
| 226 | .{ .invc = 0x1.899c0dabec30ep-1, .logc = 0x1.0d46b9e867000p-2 }, | ||
| 227 | .{ .invc = 0x1.87427aa2317fbp-1, .logc = 0x1.13687334bd000p-2 }, | ||
| 228 | .{ .invc = 0x1.84f00acb39a08p-1, .logc = 0x1.1980d67234800p-2 }, | ||
| 229 | .{ .invc = 0x1.82a49e8653e55p-1, .logc = 0x1.1f8ffe0cc8000p-2 }, | ||
| 230 | .{ .invc = 0x1.8060195f40260p-1, .logc = 0x1.2595fd7636800p-2 }, | ||
| 231 | .{ .invc = 0x1.7e22563e0a329p-1, .logc = 0x1.2b9300914a800p-2 }, | ||
| 232 | .{ .invc = 0x1.7beb377dcb5adp-1, .logc = 0x1.3187210436000p-2 }, | ||
| 233 | .{ .invc = 0x1.79baa679725c2p-1, .logc = 0x1.377266dec1800p-2 }, | ||
| 234 | .{ .invc = 0x1.77907f2170657p-1, .logc = 0x1.3d54ffbaf3000p-2 }, | ||
| 235 | .{ .invc = 0x1.756cadbd6130cp-1, .logc = 0x1.432eee32fe000p-2 }, | ||
| 236 | }; | ||
| 237 | |||
| 238 | const tab2 = [128]struct { chi: f64, clo: f64 }{ | ||
| 239 | .{ .chi = 0x1.61000014fb66bp-1, .clo = 0x1.e026c91425b3cp-56 }, | ||
| 240 | .{ .chi = 0x1.63000034db495p-1, .clo = 0x1.dbfea48005d41p-55 }, | ||
| 241 | .{ .chi = 0x1.650000d94d478p-1, .clo = 0x1.e7fa786d6a5b7p-55 }, | ||
| 242 | .{ .chi = 0x1.67000074e6fadp-1, .clo = 0x1.1fcea6b54254cp-57 }, | ||
| 243 | .{ .chi = 0x1.68ffffedf0faep-1, .clo = -0x1.c7e274c590efdp-56 }, | ||
| 244 | .{ .chi = 0x1.6b0000763c5bcp-1, .clo = -0x1.ac16848dcda01p-55 }, | ||
| 245 | .{ .chi = 0x1.6d0001e5cc1f6p-1, .clo = 0x1.33f1c9d499311p-55 }, | ||
| 246 | .{ .chi = 0x1.6efffeb05f63ep-1, .clo = -0x1.e80041ae22d53p-56 }, | ||
| 247 | .{ .chi = 0x1.710000e86978p-1, .clo = 0x1.bff6671097952p-56 }, | ||
| 248 | .{ .chi = 0x1.72ffffc67e912p-1, .clo = 0x1.c00e226bd8724p-55 }, | ||
| 249 | .{ .chi = 0x1.74fffdf81116ap-1, .clo = -0x1.e02916ef101d2p-57 }, | ||
| 250 | .{ .chi = 0x1.770000f679c9p-1, .clo = -0x1.7fc71cd549c74p-57 }, | ||
| 251 | .{ .chi = 0x1.78ffffa7ec835p-1, .clo = 0x1.1bec19ef50483p-55 }, | ||
| 252 | .{ .chi = 0x1.7affffe20c2e6p-1, .clo = -0x1.07e1729cc6465p-56 }, | ||
| 253 | .{ .chi = 0x1.7cfffed3fc9p-1, .clo = -0x1.08072087b8b1cp-55 }, | ||
| 254 | .{ .chi = 0x1.7efffe9261a76p-1, .clo = 0x1.dc0286d9df9aep-55 }, | ||
| 255 | .{ .chi = 0x1.81000049ca3e8p-1, .clo = 0x1.97fd251e54c33p-55 }, | ||
| 256 | .{ .chi = 0x1.8300017932c8fp-1, .clo = -0x1.afee9b630f381p-55 }, | ||
| 257 | .{ .chi = 0x1.850000633739cp-1, .clo = 0x1.9bfbf6b6535bcp-55 }, | ||
| 258 | .{ .chi = 0x1.87000204289c6p-1, .clo = -0x1.bbf65f3117b75p-55 }, | ||
| 259 | .{ .chi = 0x1.88fffebf57904p-1, .clo = -0x1.9006ea23dcb57p-55 }, | ||
| 260 | .{ .chi = 0x1.8b00022bc04dfp-1, .clo = -0x1.d00df38e04b0ap-56 }, | ||
| 261 | .{ .chi = 0x1.8cfffe50c1b8ap-1, .clo = -0x1.8007146ff9f05p-55 }, | ||
| 262 | .{ .chi = 0x1.8effffc918e43p-1, .clo = 0x1.3817bd07a7038p-55 }, | ||
| 263 | .{ .chi = 0x1.910001efa5fc7p-1, .clo = 0x1.93e9176dfb403p-55 }, | ||
| 264 | .{ .chi = 0x1.9300013467bb9p-1, .clo = 0x1.f804e4b980276p-56 }, | ||
| 265 | .{ .chi = 0x1.94fffe6ee076fp-1, .clo = -0x1.f7ef0d9ff622ep-55 }, | ||
| 266 | .{ .chi = 0x1.96fffde3c12d1p-1, .clo = -0x1.082aa962638bap-56 }, | ||
| 267 | .{ .chi = 0x1.98ffff4458a0dp-1, .clo = -0x1.7801b9164a8efp-55 }, | ||
| 268 | .{ .chi = 0x1.9afffdd982e3ep-1, .clo = -0x1.740e08a5a9337p-55 }, | ||
| 269 | .{ .chi = 0x1.9cfffed49fb66p-1, .clo = 0x1.fce08c19bep-60 }, | ||
| 270 | .{ .chi = 0x1.9f00020f19c51p-1, .clo = -0x1.a3faa27885b0ap-55 }, | ||
| 271 | .{ .chi = 0x1.a10001145b006p-1, .clo = 0x1.4ff489958da56p-56 }, | ||
| 272 | .{ .chi = 0x1.a300007bbf6fap-1, .clo = 0x1.cbeab8a2b6d18p-55 }, | ||
| 273 | .{ .chi = 0x1.a500010971d79p-1, .clo = 0x1.8fecadd78793p-55 }, | ||
| 274 | .{ .chi = 0x1.a70001df52e48p-1, .clo = -0x1.f41763dd8abdbp-55 }, | ||
| 275 | .{ .chi = 0x1.a90001c593352p-1, .clo = -0x1.ebf0284c27612p-55 }, | ||
| 276 | .{ .chi = 0x1.ab0002a4f3e4bp-1, .clo = -0x1.9fd043cff3f5fp-57 }, | ||
| 277 | .{ .chi = 0x1.acfffd7ae1ed1p-1, .clo = -0x1.23ee7129070b4p-55 }, | ||
| 278 | .{ .chi = 0x1.aefffee510478p-1, .clo = 0x1.a063ee00edea3p-57 }, | ||
| 279 | .{ .chi = 0x1.b0fffdb650d5bp-1, .clo = 0x1.a06c8381f0ab9p-58 }, | ||
| 280 | .{ .chi = 0x1.b2ffffeaaca57p-1, .clo = -0x1.9011e74233c1dp-56 }, | ||
| 281 | .{ .chi = 0x1.b4fffd995badcp-1, .clo = -0x1.9ff1068862a9fp-56 }, | ||
| 282 | .{ .chi = 0x1.b7000249e659cp-1, .clo = 0x1.aff45d0864f3ep-55 }, | ||
| 283 | .{ .chi = 0x1.b8ffff987164p-1, .clo = 0x1.cfe7796c2c3f9p-56 }, | ||
| 284 | .{ .chi = 0x1.bafffd204cb4fp-1, .clo = -0x1.3ff27eef22bc4p-57 }, | ||
| 285 | .{ .chi = 0x1.bcfffd2415c45p-1, .clo = -0x1.cffb7ee3bea21p-57 }, | ||
| 286 | .{ .chi = 0x1.beffff86309dfp-1, .clo = -0x1.14103972e0b5cp-55 }, | ||
| 287 | .{ .chi = 0x1.c0fffe1b57653p-1, .clo = 0x1.bc16494b76a19p-55 }, | ||
| 288 | .{ .chi = 0x1.c2ffff1fa57e3p-1, .clo = -0x1.4feef8d30c6edp-57 }, | ||
| 289 | .{ .chi = 0x1.c4fffdcbfe424p-1, .clo = -0x1.43f68bcec4775p-55 }, | ||
| 290 | .{ .chi = 0x1.c6fffed54b9f7p-1, .clo = 0x1.47ea3f053e0ecp-55 }, | ||
| 291 | .{ .chi = 0x1.c8fffeb998fd5p-1, .clo = 0x1.383068df992f1p-56 }, | ||
| 292 | .{ .chi = 0x1.cb0002125219ap-1, .clo = -0x1.8fd8e64180e04p-57 }, | ||
| 293 | .{ .chi = 0x1.ccfffdd94469cp-1, .clo = 0x1.e7ebe1cc7ea72p-55 }, | ||
| 294 | .{ .chi = 0x1.cefffeafdc476p-1, .clo = 0x1.ebe39ad9f88fep-55 }, | ||
| 295 | .{ .chi = 0x1.d1000169af82bp-1, .clo = 0x1.57d91a8b95a71p-56 }, | ||
| 296 | .{ .chi = 0x1.d30000d0ff71dp-1, .clo = 0x1.9c1906970c7dap-55 }, | ||
| 297 | .{ .chi = 0x1.d4fffea790fc4p-1, .clo = -0x1.80e37c558fe0cp-58 }, | ||
| 298 | .{ .chi = 0x1.d70002edc87e5p-1, .clo = -0x1.f80d64dc10f44p-56 }, | ||
| 299 | .{ .chi = 0x1.d900021dc82aap-1, .clo = -0x1.47c8f94fd5c5cp-56 }, | ||
| 300 | .{ .chi = 0x1.dafffd86b0283p-1, .clo = 0x1.c7f1dc521617ep-55 }, | ||
| 301 | .{ .chi = 0x1.dd000296c4739p-1, .clo = 0x1.8019eb2ffb153p-55 }, | ||
| 302 | .{ .chi = 0x1.defffe54490f5p-1, .clo = 0x1.e00d2c652cc89p-57 }, | ||
| 303 | .{ .chi = 0x1.e0fffcdabf694p-1, .clo = -0x1.f8340202d69d2p-56 }, | ||
| 304 | .{ .chi = 0x1.e2fffdb52c8ddp-1, .clo = 0x1.b00c1ca1b0864p-56 }, | ||
| 305 | .{ .chi = 0x1.e4ffff24216efp-1, .clo = 0x1.2ffa8b094ab51p-56 }, | ||
| 306 | .{ .chi = 0x1.e6fffe88a5e11p-1, .clo = -0x1.7f673b1efbe59p-58 }, | ||
| 307 | .{ .chi = 0x1.e9000119eff0dp-1, .clo = -0x1.4808d5e0bc801p-55 }, | ||
| 308 | .{ .chi = 0x1.eafffdfa51744p-1, .clo = 0x1.80006d54320b5p-56 }, | ||
| 309 | .{ .chi = 0x1.ed0001a127fa1p-1, .clo = -0x1.002f860565c92p-58 }, | ||
| 310 | .{ .chi = 0x1.ef00007babcc4p-1, .clo = -0x1.540445d35e611p-55 }, | ||
| 311 | .{ .chi = 0x1.f0ffff57a8d02p-1, .clo = -0x1.ffb3139ef9105p-59 }, | ||
| 312 | .{ .chi = 0x1.f30001ee58ac7p-1, .clo = 0x1.a81acf2731155p-55 }, | ||
| 313 | .{ .chi = 0x1.f4ffff5823494p-1, .clo = 0x1.a3f41d4d7c743p-55 }, | ||
| 314 | .{ .chi = 0x1.f6ffffca94c6bp-1, .clo = -0x1.202f41c987875p-57 }, | ||
| 315 | .{ .chi = 0x1.f8fffe1f9c441p-1, .clo = 0x1.77dd1f477e74bp-56 }, | ||
| 316 | .{ .chi = 0x1.fafffd2e0e37ep-1, .clo = -0x1.f01199a7ca331p-57 }, | ||
| 317 | .{ .chi = 0x1.fd0001c77e49ep-1, .clo = 0x1.181ee4bceacb1p-56 }, | ||
| 318 | .{ .chi = 0x1.feffff7e0c331p-1, .clo = -0x1.e05370170875ap-57 }, | ||
| 319 | .{ .chi = 0x1.00ffff465606ep+0, .clo = -0x1.a7ead491c0adap-55 }, | ||
| 320 | .{ .chi = 0x1.02ffff3867a58p+0, .clo = -0x1.77f69c3fcb2ep-54 }, | ||
| 321 | .{ .chi = 0x1.04ffffdfc0d17p+0, .clo = 0x1.7bffe34cb945bp-54 }, | ||
| 322 | .{ .chi = 0x1.0700003cd4d82p+0, .clo = 0x1.20083c0e456cbp-55 }, | ||
| 323 | .{ .chi = 0x1.08ffff9f2cbe8p+0, .clo = -0x1.dffdfbe37751ap-57 }, | ||
| 324 | .{ .chi = 0x1.0b000010cda65p+0, .clo = -0x1.13f7faee626ebp-54 }, | ||
| 325 | .{ .chi = 0x1.0d00001a4d338p+0, .clo = 0x1.07dfa79489ff7p-55 }, | ||
| 326 | .{ .chi = 0x1.0effffadafdfdp+0, .clo = -0x1.7040570d66bcp-56 }, | ||
| 327 | .{ .chi = 0x1.110000bbafd96p+0, .clo = 0x1.e80d4846d0b62p-55 }, | ||
| 328 | .{ .chi = 0x1.12ffffae5f45dp+0, .clo = 0x1.dbffa64fd36efp-54 }, | ||
| 329 | .{ .chi = 0x1.150000dd59ad9p+0, .clo = 0x1.a0077701250aep-54 }, | ||
| 330 | .{ .chi = 0x1.170000f21559ap+0, .clo = 0x1.dfdf9e2e3deeep-55 }, | ||
| 331 | .{ .chi = 0x1.18ffffc275426p+0, .clo = 0x1.10030dc3b7273p-54 }, | ||
| 332 | .{ .chi = 0x1.1b000123d3c59p+0, .clo = 0x1.97f7980030188p-54 }, | ||
| 333 | .{ .chi = 0x1.1cffff8299eb7p+0, .clo = -0x1.5f932ab9f8c67p-57 }, | ||
| 334 | .{ .chi = 0x1.1effff48ad4p+0, .clo = 0x1.37fbf9da75bebp-54 }, | ||
| 335 | .{ .chi = 0x1.210000c8b86a4p+0, .clo = 0x1.f806b91fd5b22p-54 }, | ||
| 336 | .{ .chi = 0x1.2300003854303p+0, .clo = 0x1.3ffc2eb9fbf33p-54 }, | ||
| 337 | .{ .chi = 0x1.24fffffbcf684p+0, .clo = 0x1.601e77e2e2e72p-56 }, | ||
| 338 | .{ .chi = 0x1.26ffff52921d9p+0, .clo = 0x1.ffcbb767f0c61p-56 }, | ||
| 339 | .{ .chi = 0x1.2900014933a3cp+0, .clo = -0x1.202ca3c02412bp-56 }, | ||
| 340 | .{ .chi = 0x1.2b00014556313p+0, .clo = -0x1.2808233f21f02p-54 }, | ||
| 341 | .{ .chi = 0x1.2cfffebfe523bp+0, .clo = -0x1.8ff7e384fdcf2p-55 }, | ||
| 342 | .{ .chi = 0x1.2f0000bb8ad96p+0, .clo = -0x1.5ff51503041c5p-55 }, | ||
| 343 | .{ .chi = 0x1.30ffffb7ae2afp+0, .clo = -0x1.10071885e289dp-55 }, | ||
| 344 | .{ .chi = 0x1.32ffffeac5f7fp+0, .clo = -0x1.1ff5d3fb7b715p-54 }, | ||
| 345 | .{ .chi = 0x1.350000ca66756p+0, .clo = 0x1.57f82228b82bdp-54 }, | ||
| 346 | .{ .chi = 0x1.3700011fbf721p+0, .clo = 0x1.000bac40dd5ccp-55 }, | ||
| 347 | .{ .chi = 0x1.38ffff9592fb9p+0, .clo = -0x1.43f9d2db2a751p-54 }, | ||
| 348 | .{ .chi = 0x1.3b00004ddd242p+0, .clo = 0x1.57f6b707638e1p-55 }, | ||
| 349 | .{ .chi = 0x1.3cffff5b2c957p+0, .clo = 0x1.a023a10bf1231p-56 }, | ||
| 350 | .{ .chi = 0x1.3efffeab0b418p+0, .clo = 0x1.87f6d66b152bp-54 }, | ||
| 351 | .{ .chi = 0x1.410001532aff4p+0, .clo = 0x1.7f8375f198524p-57 }, | ||
| 352 | .{ .chi = 0x1.4300017478b29p+0, .clo = 0x1.301e672dc5143p-55 }, | ||
| 353 | .{ .chi = 0x1.44fffe795b463p+0, .clo = 0x1.9ff69b8b2895ap-55 }, | ||
| 354 | .{ .chi = 0x1.46fffe80475ep+0, .clo = -0x1.5c0b19bc2f254p-54 }, | ||
| 355 | .{ .chi = 0x1.48fffef6fc1e7p+0, .clo = 0x1.b4009f23a2a72p-54 }, | ||
| 356 | .{ .chi = 0x1.4afffe5bea704p+0, .clo = -0x1.4ffb7bf0d7d45p-54 }, | ||
| 357 | .{ .chi = 0x1.4d000171027dep+0, .clo = -0x1.9c06471dc6a3dp-54 }, | ||
| 358 | .{ .chi = 0x1.4f0000ff03ee2p+0, .clo = 0x1.77f890b85531cp-54 }, | ||
| 359 | .{ .chi = 0x1.5100012dc4bd1p+0, .clo = 0x1.004657166a436p-57 }, | ||
| 360 | .{ .chi = 0x1.530001605277ap+0, .clo = -0x1.6bfcece233209p-54 }, | ||
| 361 | .{ .chi = 0x1.54fffecdb704cp+0, .clo = -0x1.902720505a1d7p-55 }, | ||
| 362 | .{ .chi = 0x1.56fffef5f54a9p+0, .clo = 0x1.bbfe60ec96412p-54 }, | ||
| 363 | .{ .chi = 0x1.5900017e61012p+0, .clo = 0x1.87ec581afef9p-55 }, | ||
| 364 | .{ .chi = 0x1.5b00003c93e92p+0, .clo = -0x1.f41080abf0ccp-54 }, | ||
| 365 | .{ .chi = 0x1.5d0001d4919bcp+0, .clo = -0x1.8812afb254729p-54 }, | ||
| 366 | .{ .chi = 0x1.5efffe7b87a89p+0, .clo = -0x1.47eb780ed6904p-54 }, | ||
| 367 | }; | ||
| 368 | |||
| 369 | var ix: i64 = @bitCast(x); | ||
| 370 | |||
| 371 | const LO: i64 = @bitCast(@as(f64, 1.0 - 0x1p-4)); | ||
| 372 | const HI: i64 = @bitCast(@as(f64, 1.0 + 0x1.09p-4)); | ||
| 373 | if (LO <= ix and ix < HI) { | ||
| 374 | @branchHint(.unlikely); | ||
| 375 | if (ix == @as(i64, @bitCast(@as(f64, 1)))) { | ||
| 376 | @branchHint(.unlikely); | ||
| 377 | return 0; | ||
| 108 | } | 378 | } |
| 109 | 379 | ||
| 110 | // subnormal, scale x | 380 | const r = x - 1; |
| 111 | k -= 54; | 381 | const r2 = r * r; |
| 112 | x *= 0x1p54; | 382 | const r3 = r * r2; |
| 113 | hx = @intCast(@as(u64, @bitCast(x)) >> 32); | 383 | |
| 114 | } else if (hx >= 0x7FF00000) { | 384 | const y = r3 * (poly1[1] + r * poly1[2] + r2 * poly1[3] + |
| 115 | return x; | 385 | r3 * (poly1[4] + r * poly1[5] + r2 * poly1[6] + |
| 116 | } else if (hx == 0x3FF00000 and ix << 32 == 0) { | 386 | r3 * (poly1[7] + r * poly1[8] + r2 * poly1[9] + r3 * poly1[10]))); |
| 117 | return 0; | 387 | |
| 388 | var w = r * 0x1p27; | ||
| 389 | const rhi = r + w - w; | ||
| 390 | const rlo = r - rhi; | ||
| 391 | w = rhi * rhi * poly1[0]; | ||
| 392 | const hi = r + w; | ||
| 393 | const lo = r - hi + w + poly1[0] * rlo * (rhi + r); | ||
| 394 | return y + lo + hi; | ||
| 118 | } | 395 | } |
| 396 | const top = @as(u64, @bitCast(ix)) >> 48; | ||
| 397 | if (top < 0x0010 or 0x7ff0 <= top) { | ||
| 398 | @branchHint(.unlikely); | ||
| 119 | 399 | ||
| 120 | // x into [sqrt(2) / 2, sqrt(2)] | 400 | if (ix << 1 == 0) |
| 121 | hx += 0x3FF00000 - 0x3FE6A09E; | 401 | return if (common.want_float_exceptions) -1 / (x * x) else -std.math.inf(f64); |
| 122 | k += @as(i32, @intCast(hx >> 20)) - 0x3FF; | ||
| 123 | hx = (hx & 0x000FFFFF) + 0x3FE6A09E; | ||
| 124 | ix = (@as(u64, hx) << 32) | (ix & 0xFFFFFFFF); | ||
| 125 | x = @bitCast(ix); | ||
| 126 | 402 | ||
| 127 | const f = x - 1.0; | 403 | if (ix == @as(i64, @bitCast(std.math.inf(f64)))) |
| 128 | const hfsq = 0.5 * f * f; | 404 | return x; |
| 129 | const s = f / (2.0 + f); | ||
| 130 | const z = s * s; | ||
| 131 | const w = z * z; | ||
| 132 | const t1 = w * (Lg2 + w * (Lg4 + w * Lg6)); | ||
| 133 | const t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7))); | ||
| 134 | const R = t2 + t1; | ||
| 135 | const dk: f64 = @floatFromInt(k); | ||
| 136 | 405 | ||
| 137 | return s * (hfsq + R) + dk * ln2_lo - hfsq + f + dk * ln2_hi; | 406 | if (top & 0x8000 != 0 or top & 0x7ff0 == 0x7ff0) |
| 407 | return if (common.want_float_exceptions) (x - x) / 0.0 else math.nan(f64); | ||
| 408 | |||
| 409 | ix = @as(i64, @bitCast(x * 0x1p52)) - (52 << 52); | ||
| 410 | } | ||
| 411 | |||
| 412 | const tmp: packed struct(i64) { unused: u45, i: u7, k: i12 } = @bitCast(ix - 0x3fe6000000000000); | ||
| 413 | const i = tmp.i; | ||
| 414 | const k = tmp.k; | ||
| 415 | const iz = ix - (@as(i64, tmp.k) << 52); | ||
| 416 | const invc = tab[i].invc; | ||
| 417 | const logc = tab[i].logc; | ||
| 418 | const z: f64 = @bitCast(iz); | ||
| 419 | |||
| 420 | // maybe use fma as musl does | ||
| 421 | const r = (z - tab2[i].chi - tab2[i].clo) * invc; | ||
| 422 | |||
| 423 | // https://github.com/ziglang/zig/issues/18614 | ||
| 424 | const kd: f64 = @floatFromInt(k); | ||
| 425 | const w = kd * 0x1.62e42fefa3800p-1 + logc; | ||
| 426 | const hi = w + r; | ||
| 427 | const lo = w - hi + r + kd * 0x1.ef35793c76730p-45; | ||
| 428 | |||
| 429 | const r2 = r * r; | ||
| 430 | |||
| 431 | const y = lo + r2 * poly[0] + r * r2 * (poly[1] + r * poly[2] + r2 * (poly[3] + r * poly[4])) + hi; | ||
| 432 | return @bitCast(y); | ||
| 138 | } | 433 | } |
| 139 | 434 | ||
| 140 | pub fn __logx(a: f80) callconv(.c) f80 { | 435 | pub fn __logx(a: f80) callconv(.c) f80 { |
lib/compiler_rt/log10.zig+4-4| ... | @@ -49,11 +49,11 @@ pub fn log10f(x_: f32) callconv(.c) f32 { | ... | @@ -49,11 +49,11 @@ pub fn log10f(x_: f32) callconv(.c) f32 { |
| 49 | if (ix < 0x00800000 or ix >> 31 != 0) { | 49 | if (ix < 0x00800000 or ix >> 31 != 0) { |
| 50 | // log(+-0) = -inf | 50 | // log(+-0) = -inf |
| 51 | if (ix << 1 == 0) { | 51 | if (ix << 1 == 0) { |
| 52 | return -math.inf(f32); | 52 | return if (common.want_float_exceptions) -1 / (x * x) else -std.math.inf(f64); |
| 53 | } | 53 | } |
| 54 | // log(-#) = nan | 54 | // log(-#) = nan |
| 55 | if (ix >> 31 != 0) { | 55 | if (ix >> 31 != 0) { |
| 56 | return math.nan(f32); | 56 | return if (common.want_float_exceptions) (x - x) / 0.0 else math.nan(f64); |
| 57 | } | 57 | } |
| 58 | 58 | ||
| 59 | k -= 25; | 59 | k -= 25; |
| ... | @@ -111,11 +111,11 @@ pub fn log10(x_: f64) callconv(.c) f64 { | ... | @@ -111,11 +111,11 @@ pub fn log10(x_: f64) callconv(.c) f64 { |
| 111 | if (hx < 0x00100000 or hx >> 31 != 0) { | 111 | if (hx < 0x00100000 or hx >> 31 != 0) { |
| 112 | // log(+-0) = -inf | 112 | // log(+-0) = -inf |
| 113 | if (ix << 1 == 0) { | 113 | if (ix << 1 == 0) { |
| 114 | return -math.inf(f64); | 114 | return if (common.want_float_exceptions) -1 / (x * x) else -std.math.inf(f64); |
| 115 | } | 115 | } |
| 116 | // log(-#) = nan | 116 | // log(-#) = nan |
| 117 | if (hx >> 31 != 0) { | 117 | if (hx >> 31 != 0) { |
| 118 | return math.nan(f64); | 118 | return if (common.want_float_exceptions) (x - x) / 0.0 else math.nan(f64); |
| 119 | } | 119 | } |
| 120 | 120 | ||
| 121 | // subnormal, scale x | 121 | // subnormal, scale x |
lib/compiler_rt/log2.zig+4-4| ... | @@ -47,11 +47,11 @@ pub fn log2f(x_: f32) callconv(.c) f32 { | ... | @@ -47,11 +47,11 @@ pub fn log2f(x_: f32) callconv(.c) f32 { |
| 47 | if (ix < 0x00800000 or ix >> 31 != 0) { | 47 | if (ix < 0x00800000 or ix >> 31 != 0) { |
| 48 | // log(+-0) = -inf | 48 | // log(+-0) = -inf |
| 49 | if (ix << 1 == 0) { | 49 | if (ix << 1 == 0) { |
| 50 | return -math.inf(f32); | 50 | return if (common.want_float_exceptions) -1 / (x * x) else -std.math.inf(f64); |
| 51 | } | 51 | } |
| 52 | // log(-#) = nan | 52 | // log(-#) = nan |
| 53 | if (ix >> 31 != 0) { | 53 | if (ix >> 31 != 0) { |
| 54 | return math.nan(f32); | 54 | return if (common.want_float_exceptions) (x - x) / 0.0 else math.nan(f64); |
| 55 | } | 55 | } |
| 56 | 56 | ||
| 57 | k -= 25; | 57 | k -= 25; |
| ... | @@ -105,11 +105,11 @@ pub fn log2(x_: f64) callconv(.c) f64 { | ... | @@ -105,11 +105,11 @@ pub fn log2(x_: f64) callconv(.c) f64 { |
| 105 | if (hx < 0x00100000 or hx >> 31 != 0) { | 105 | if (hx < 0x00100000 or hx >> 31 != 0) { |
| 106 | // log(+-0) = -inf | 106 | // log(+-0) = -inf |
| 107 | if (ix << 1 == 0) { | 107 | if (ix << 1 == 0) { |
| 108 | return -math.inf(f64); | 108 | return if (common.want_float_exceptions) -1 / (x * x) else -std.math.inf(f64); |
| 109 | } | 109 | } |
| 110 | // log(-#) = nan | 110 | // log(-#) = nan |
| 111 | if (hx >> 31 != 0) { | 111 | if (hx >> 31 != 0) { |
| 112 | return math.nan(f64); | 112 | return if (common.want_float_exceptions) (x - x) / 0.0 else math.nan(f64); |
| 113 | } | 113 | } |
| 114 | 114 | ||
| 115 | // subnormal, scale x | 115 | // subnormal, scale x |
lib/libc/mingw/math/expf.c deleted-11| ... | @@ -1,11 +0,0 @@ | ||
| 1 | /** | ||
| 2 | * This file has no copyright assigned and is placed in the Public Domain. | ||
| 3 | * This file is part of the mingw-w64 runtime package. | ||
| 4 | * No warranty is given; refer to the file DISCLAIMER.PD within this package. | ||
| 5 | */ | ||
| 6 | #include <math.h> | ||
| 7 | |||
| 8 | float expf (float x) | ||
| 9 | { | ||
| 10 | return (float) exp (x); | ||
| 11 | } | ||
lib/libc/mingw/math/log10f.c deleted-11| ... | @@ -1,11 +0,0 @@ | ||
| 1 | /** | ||
| 2 | * This file has no copyright assigned and is placed in the Public Domain. | ||
| 3 | * This file is part of the mingw-w64 runtime package. | ||
| 4 | * No warranty is given; refer to the file DISCLAIMER.PD within this package. | ||
| 5 | */ | ||
| 6 | #include <math.h> | ||
| 7 | |||
| 8 | float log10f(float _X) | ||
| 9 | { | ||
| 10 | return ((float)log10((double)_X)); | ||
| 11 | } | ||
lib/libc/mingw/math/logf.c deleted-11| ... | @@ -1,11 +0,0 @@ | ||
| 1 | /** | ||
| 2 | * This file has no copyright assigned and is placed in the Public Domain. | ||
| 3 | * This file is part of the mingw-w64 runtime package. | ||
| 4 | * No warranty is given; refer to the file DISCLAIMER.PD within this package. | ||
| 5 | */ | ||
| 6 | #include <math.h> | ||
| 7 | |||
| 8 | float logf(float _X) | ||
| 9 | { | ||
| 10 | return ((float)log((double)_X)); | ||
| 11 | } | ||
lib/libc/musl/src/math/exp.c deleted-134| ... | @@ -1,134 +0,0 @@ | ||
| 1 | /* | ||
| 2 | * Double-precision e^x function. | ||
| 3 | * | ||
| 4 | * Copyright (c) 2018, Arm Limited. | ||
| 5 | * SPDX-License-Identifier: MIT | ||
| 6 | */ | ||
| 7 | |||
| 8 | #include <math.h> | ||
| 9 | #include <stdint.h> | ||
| 10 | #include "libm.h" | ||
| 11 | #include "exp_data.h" | ||
| 12 | |||
| 13 | #define N (1 << EXP_TABLE_BITS) | ||
| 14 | #define InvLn2N __exp_data.invln2N | ||
| 15 | #define NegLn2hiN __exp_data.negln2hiN | ||
| 16 | #define NegLn2loN __exp_data.negln2loN | ||
| 17 | #define Shift __exp_data.shift | ||
| 18 | #define T __exp_data.tab | ||
| 19 | #define C2 __exp_data.poly[5 - EXP_POLY_ORDER] | ||
| 20 | #define C3 __exp_data.poly[6 - EXP_POLY_ORDER] | ||
| 21 | #define C4 __exp_data.poly[7 - EXP_POLY_ORDER] | ||
| 22 | #define C5 __exp_data.poly[8 - EXP_POLY_ORDER] | ||
| 23 | |||
| 24 | /* Handle cases that may overflow or underflow when computing the result that | ||
| 25 | is scale*(1+TMP) without intermediate rounding. The bit representation of | ||
| 26 | scale is in SBITS, however it has a computed exponent that may have | ||
| 27 | overflown into the sign bit so that needs to be adjusted before using it as | ||
| 28 | a double. (int32_t)KI is the k used in the argument reduction and exponent | ||
| 29 | adjustment of scale, positive k here means the result may overflow and | ||
| 30 | negative k means the result may underflow. */ | ||
| 31 | static inline double specialcase(double_t tmp, uint64_t sbits, uint64_t ki) | ||
| 32 | { | ||
| 33 | 	double_t scale, y; | ||
| 34 | |||
| 35 | 	if ((ki & 0x80000000) == 0) { | ||
| 36 | 		/* k > 0, the exponent of scale might have overflowed by <= 460. */ | ||
| 37 | 		sbits -= 1009ull << 52; | ||
| 38 | 		scale = asdouble(sbits); | ||
| 39 | 		y = 0x1p1009 * (scale + scale * tmp); | ||
| 40 | 		return eval_as_double(y); | ||
| 41 | 	} | ||
| 42 | 	/* k < 0, need special care in the subnormal range. */ | ||
| 43 | 	sbits += 1022ull << 52; | ||
| 44 | 	scale = asdouble(sbits); | ||
| 45 | 	y = scale + scale * tmp; | ||
| 46 | 	if (y < 1.0) { | ||
| 47 | 		/* Round y to the right precision before scaling it into the subnormal | ||
| 48 | 		 range to avoid double rounding that can cause 0.5+E/2 ulp error where | ||
| 49 | 		 E is the worst-case ulp error outside the subnormal range. So this | ||
| 50 | 		 is only useful if the goal is better than 1 ulp worst-case error. */ | ||
| 51 | 		double_t hi, lo; | ||
| 52 | 		lo = scale - y + scale * tmp; | ||
| 53 | 		hi = 1.0 + y; | ||
| 54 | 		lo = 1.0 - hi + y + lo; | ||
| 55 | 		y = eval_as_double(hi + lo) - 1.0; | ||
| 56 | 		/* Avoid -0.0 with downward rounding. */ | ||
| 57 | 		if (WANT_ROUNDING && y == 0.0) | ||
| 58 | 			y = 0.0; | ||
| 59 | 		/* The underflow exception needs to be signaled explicitly. */ | ||
| 60 | 		fp_force_eval(fp_barrier(0x1p-1022) * 0x1p-1022); | ||
| 61 | 	} | ||
| 62 | 	y = 0x1p-1022 * y; | ||
| 63 | 	return eval_as_double(y); | ||
| 64 | } | ||
| 65 | |||
| 66 | /* Top 12 bits of a double (sign and exponent bits). */ | ||
| 67 | static inline uint32_t top12(double x) | ||
| 68 | { | ||
| 69 | 	return asuint64(x) >> 52; | ||
| 70 | } | ||
| 71 | |||
| 72 | double exp(double x) | ||
| 73 | { | ||
| 74 | 	uint32_t abstop; | ||
| 75 | 	uint64_t ki, idx, top, sbits; | ||
| 76 | 	double_t kd, z, r, r2, scale, tail, tmp; | ||
| 77 | |||
| 78 | 	abstop = top12(x) & 0x7ff; | ||
| 79 | 	if (predict_false(abstop - top12(0x1p-54) >= top12(512.0) - top12(0x1p-54))) { | ||
| 80 | 		if (abstop - top12(0x1p-54) >= 0x80000000) | ||
| 81 | 			/* Avoid spurious underflow for tiny x. */ | ||
| 82 | 			/* Note: 0 is common input. */ | ||
| 83 | 			return WANT_ROUNDING ? 1.0 + x : 1.0; | ||
| 84 | 		if (abstop >= top12(1024.0)) { | ||
| 85 | 			if (asuint64(x) == asuint64(-INFINITY)) | ||
| 86 | 				return 0.0; | ||
| 87 | 			if (abstop >= top12(INFINITY)) | ||
| 88 | 				return 1.0 + x; | ||
| 89 | 			if (asuint64(x) >> 63) | ||
| 90 | 				return __math_uflow(0); | ||
| 91 | 			else | ||
| 92 | 				return __math_oflow(0); | ||
| 93 | 		} | ||
| 94 | 		/* Large x is special cased below. */ | ||
| 95 | 		abstop = 0; | ||
| 96 | 	} | ||
| 97 | |||
| 98 | 	/* exp(x) = 2^(k/N) * exp(r), with exp(r) in [2^(-1/2N),2^(1/2N)]. */ | ||
| 99 | 	/* x = ln2/N*k + r, with int k and r in [-ln2/2N, ln2/2N]. */ | ||
| 100 | 	z = InvLn2N * x; | ||
| 101 | #if TOINT_INTRINSICS | ||
| 102 | 	kd = roundtoint(z); | ||
| 103 | 	ki = converttoint(z); | ||
| 104 | #elif EXP_USE_TOINT_NARROW | ||
| 105 | 	/* z - kd is in [-0.5-2^-16, 0.5] in all rounding modes. */ | ||
| 106 | 	kd = eval_as_double(z + Shift); | ||
| 107 | 	ki = asuint64(kd) >> 16; | ||
| 108 | 	kd = (double_t)(int32_t)ki; | ||
| 109 | #else | ||
| 110 | 	/* z - kd is in [-1, 1] in non-nearest rounding modes. */ | ||
| 111 | 	kd = eval_as_double(z + Shift); | ||
| 112 | 	ki = asuint64(kd); | ||
| 113 | 	kd -= Shift; | ||
| 114 | #endif | ||
| 115 | 	r = x + kd * NegLn2hiN + kd * NegLn2loN; | ||
| 116 | 	/* 2^(k/N) ~= scale * (1 + tail). */ | ||
| 117 | 	idx = 2 * (ki % N); | ||
| 118 | 	top = ki << (52 - EXP_TABLE_BITS); | ||
| 119 | 	tail = asdouble(T[idx]); | ||
| 120 | 	/* This is only a valid scale when -1023*N < k < 1024*N. */ | ||
| 121 | 	sbits = T[idx + 1] + top; | ||
| 122 | 	/* exp(x) = 2^(k/N) * exp(r) ~= scale + scale * (tail + exp(r) - 1). */ | ||
| 123 | 	/* Evaluation is optimized assuming superscalar pipelined execution. */ | ||
| 124 | 	r2 = r * r; | ||
| 125 | 	/* Without fma the worst case error is 0.25/N ulp larger. */ | ||
| 126 | 	/* Worst case error is less than 0.5+1.11/N+(abs poly error * 2^53) ulp. */ | ||
| 127 | 	tmp = tail + r + r2 * (C2 + r * C3) + r2 * r2 * (C4 + r * C5); | ||
| 128 | 	if (predict_false(abstop == 0)) | ||
| 129 | 		return specialcase(tmp, sbits, ki); | ||
| 130 | 	scale = asdouble(sbits); | ||
| 131 | 	/* Note: tmp == 0 or |tmp| > 2^-200 and scale > 2^-739, so there | ||
| 132 | 	 is no spurious underflow here even without fma. */ | ||
| 133 | 	return eval_as_double(scale + scale * tmp); | ||
| 134 | } | ||
lib/libc/musl/src/math/exp2.c deleted-121| ... | @@ -1,121 +0,0 @@ | ||
| 1 | /* | ||
| 2 | * Double-precision 2^x function. | ||
| 3 | * | ||
| 4 | * Copyright (c) 2018, Arm Limited. | ||
| 5 | * SPDX-License-Identifier: MIT | ||
| 6 | */ | ||
| 7 | |||
| 8 | #include <math.h> | ||
| 9 | #include <stdint.h> | ||
| 10 | #include "libm.h" | ||
| 11 | #include "exp_data.h" | ||
| 12 | |||
| 13 | #define N (1 << EXP_TABLE_BITS) | ||
| 14 | #define Shift __exp_data.exp2_shift | ||
| 15 | #define T __exp_data.tab | ||
| 16 | #define C1 __exp_data.exp2_poly[0] | ||
| 17 | #define C2 __exp_data.exp2_poly[1] | ||
| 18 | #define C3 __exp_data.exp2_poly[2] | ||
| 19 | #define C4 __exp_data.exp2_poly[3] | ||
| 20 | #define C5 __exp_data.exp2_poly[4] | ||
| 21 | |||
| 22 | /* Handle cases that may overflow or underflow when computing the result that | ||
| 23 | is scale*(1+TMP) without intermediate rounding. The bit representation of | ||
| 24 | scale is in SBITS, however it has a computed exponent that may have | ||
| 25 | overflown into the sign bit so that needs to be adjusted before using it as | ||
| 26 | a double. (int32_t)KI is the k used in the argument reduction and exponent | ||
| 27 | adjustment of scale, positive k here means the result may overflow and | ||
| 28 | negative k means the result may underflow. */ | ||
| 29 | static inline double specialcase(double_t tmp, uint64_t sbits, uint64_t ki) | ||
| 30 | { | ||
| 31 | 	double_t scale, y; | ||
| 32 | |||
| 33 | 	if ((ki & 0x80000000) == 0) { | ||
| 34 | 		/* k > 0, the exponent of scale might have overflowed by 1. */ | ||
| 35 | 		sbits -= 1ull << 52; | ||
| 36 | 		scale = asdouble(sbits); | ||
| 37 | 		y = 2 * (scale + scale * tmp); | ||
| 38 | 		return eval_as_double(y); | ||
| 39 | 	} | ||
| 40 | 	/* k < 0, need special care in the subnormal range. */ | ||
| 41 | 	sbits += 1022ull << 52; | ||
| 42 | 	scale = asdouble(sbits); | ||
| 43 | 	y = scale + scale * tmp; | ||
| 44 | 	if (y < 1.0) { | ||
| 45 | 		/* Round y to the right precision before scaling it into the subnormal | ||
| 46 | 		 range to avoid double rounding that can cause 0.5+E/2 ulp error where | ||
| 47 | 		 E is the worst-case ulp error outside the subnormal range. So this | ||
| 48 | 		 is only useful if the goal is better than 1 ulp worst-case error. */ | ||
| 49 | 		double_t hi, lo; | ||
| 50 | 		lo = scale - y + scale * tmp; | ||
| 51 | 		hi = 1.0 + y; | ||
| 52 | 		lo = 1.0 - hi + y + lo; | ||
| 53 | 		y = eval_as_double(hi + lo) - 1.0; | ||
| 54 | 		/* Avoid -0.0 with downward rounding. */ | ||
| 55 | 		if (WANT_ROUNDING && y == 0.0) | ||
| 56 | 			y = 0.0; | ||
| 57 | 		/* The underflow exception needs to be signaled explicitly. */ | ||
| 58 | 		fp_force_eval(fp_barrier(0x1p-1022) * 0x1p-1022); | ||
| 59 | 	} | ||
| 60 | 	y = 0x1p-1022 * y; | ||
| 61 | 	return eval_as_double(y); | ||
| 62 | } | ||
| 63 | |||
| 64 | /* Top 12 bits of a double (sign and exponent bits). */ | ||
| 65 | static inline uint32_t top12(double x) | ||
| 66 | { | ||
| 67 | 	return asuint64(x) >> 52; | ||
| 68 | } | ||
| 69 | |||
| 70 | double exp2(double x) | ||
| 71 | { | ||
| 72 | 	uint32_t abstop; | ||
| 73 | 	uint64_t ki, idx, top, sbits; | ||
| 74 | 	double_t kd, r, r2, scale, tail, tmp; | ||
| 75 | |||
| 76 | 	abstop = top12(x) & 0x7ff; | ||
| 77 | 	if (predict_false(abstop - top12(0x1p-54) >= top12(512.0) - top12(0x1p-54))) { | ||
| 78 | 		if (abstop - top12(0x1p-54) >= 0x80000000) | ||
| 79 | 			/* Avoid spurious underflow for tiny x. */ | ||
| 80 | 			/* Note: 0 is common input. */ | ||
| 81 | 			return WANT_ROUNDING ? 1.0 + x : 1.0; | ||
| 82 | 		if (abstop >= top12(1024.0)) { | ||
| 83 | 			if (asuint64(x) == asuint64(-INFINITY)) | ||
| 84 | 				return 0.0; | ||
| 85 | 			if (abstop >= top12(INFINITY)) | ||
| 86 | 				return 1.0 + x; | ||
| 87 | 			if (!(asuint64(x) >> 63)) | ||
| 88 | 				return __math_oflow(0); | ||
| 89 | 			else if (asuint64(x) >= asuint64(-1075.0)) | ||
| 90 | 				return __math_uflow(0); | ||
| 91 | 		} | ||
| 92 | 		if (2 * asuint64(x) > 2 * asuint64(928.0)) | ||
| 93 | 			/* Large x is special cased below. */ | ||
| 94 | 			abstop = 0; | ||
| 95 | 	} | ||
| 96 | |||
| 97 | 	/* exp2(x) = 2^(k/N) * 2^r, with 2^r in [2^(-1/2N),2^(1/2N)]. */ | ||
| 98 | 	/* x = k/N + r, with int k and r in [-1/2N, 1/2N]. */ | ||
| 99 | 	kd = eval_as_double(x + Shift); | ||
| 100 | 	ki = asuint64(kd); /* k. */ | ||
| 101 | 	kd -= Shift; /* k/N for int k. */ | ||
| 102 | 	r = x - kd; | ||
| 103 | 	/* 2^(k/N) ~= scale * (1 + tail). */ | ||
| 104 | 	idx = 2 * (ki % N); | ||
| 105 | 	top = ki << (52 - EXP_TABLE_BITS); | ||
| 106 | 	tail = asdouble(T[idx]); | ||
| 107 | 	/* This is only a valid scale when -1023*N < k < 1024*N. */ | ||
| 108 | 	sbits = T[idx + 1] + top; | ||
| 109 | 	/* exp2(x) = 2^(k/N) * 2^r ~= scale + scale * (tail + 2^r - 1). */ | ||
| 110 | 	/* Evaluation is optimized assuming superscalar pipelined execution. */ | ||
| 111 | 	r2 = r * r; | ||
| 112 | 	/* Without fma the worst case error is 0.5/N ulp larger. */ | ||
| 113 | 	/* Worst case error is less than 0.5+0.86/N+(abs poly error * 2^53) ulp. */ | ||
| 114 | 	tmp = tail + r * C1 + r2 * (C2 + r * C3) + r2 * r2 * (C4 + r * C5); | ||
| 115 | 	if (predict_false(abstop == 0)) | ||
| 116 | 		return specialcase(tmp, sbits, ki); | ||
| 117 | 	scale = asdouble(sbits); | ||
| 118 | 	/* Note: tmp == 0 or |tmp| > 2^-65 and scale > 2^-928, so there | ||
| 119 | 	 is no spurious underflow here even without fma. */ | ||
| 120 | 	return eval_as_double(scale + scale * tmp); | ||
| 121 | } | ||
lib/libc/musl/src/math/exp2f.c deleted-69| ... | @@ -1,69 +0,0 @@ | ||
| 1 | /* | ||
| 2 | * Single-precision 2^x function. | ||
| 3 | * | ||
| 4 | * Copyright (c) 2017-2018, Arm Limited. | ||
| 5 | * SPDX-License-Identifier: MIT | ||
| 6 | */ | ||
| 7 | |||
| 8 | #include <math.h> | ||
| 9 | #include <stdint.h> | ||
| 10 | #include "libm.h" | ||
| 11 | #include "exp2f_data.h" | ||
| 12 | |||
| 13 | /* | ||
| 14 | EXP2F_TABLE_BITS = 5 | ||
| 15 | EXP2F_POLY_ORDER = 3 | ||
| 16 | |||
| 17 | ULP error: 0.502 (nearest rounding.) | ||
| 18 | Relative error: 1.69 * 2^-34 in [-1/64, 1/64] (before rounding.) | ||
| 19 | Wrong count: 168353 (all nearest rounding wrong results with fma.) | ||
| 20 | Non-nearest ULP error: 1 (rounded ULP error) | ||
| 21 | */ | ||
| 22 | |||
| 23 | #define N (1 << EXP2F_TABLE_BITS) | ||
| 24 | #define T __exp2f_data.tab | ||
| 25 | #define C __exp2f_data.poly | ||
| 26 | #define SHIFT __exp2f_data.shift_scaled | ||
| 27 | |||
| 28 | static inline uint32_t top12(float x) | ||
| 29 | { | ||
| 30 | 	return asuint(x) >> 20; | ||
| 31 | } | ||
| 32 | |||
| 33 | float exp2f(float x) | ||
| 34 | { | ||
| 35 | 	uint32_t abstop; | ||
| 36 | 	uint64_t ki, t; | ||
| 37 | 	double_t kd, xd, z, r, r2, y, s; | ||
| 38 | |||
| 39 | 	xd = (double_t)x; | ||
| 40 | 	abstop = top12(x) & 0x7ff; | ||
| 41 | 	if (predict_false(abstop >= top12(128.0f))) { | ||
| 42 | 		/* |x| >= 128 or x is nan. */ | ||
| 43 | 		if (asuint(x) == asuint(-INFINITY)) | ||
| 44 | 			return 0.0f; | ||
| 45 | 		if (abstop >= top12(INFINITY)) | ||
| 46 | 			return x + x; | ||
| 47 | 		if (x > 0.0f) | ||
| 48 | 			return __math_oflowf(0); | ||
| 49 | 		if (x <= -150.0f) | ||
| 50 | 			return __math_uflowf(0); | ||
| 51 | 	} | ||
| 52 | |||
| 53 | 	/* x = k/N + r with r in [-1/(2N), 1/(2N)] and int k. */ | ||
| 54 | 	kd = eval_as_double(xd + SHIFT); | ||
| 55 | 	ki = asuint64(kd); | ||
| 56 | 	kd -= SHIFT; /* k/N for int k. */ | ||
| 57 | 	r = xd - kd; | ||
| 58 | |||
| 59 | 	/* exp2(x) = 2^(k/N) * 2^r ~= s * (C0*r^3 + C1*r^2 + C2*r + 1) */ | ||
| 60 | 	t = T[ki % N]; | ||
| 61 | 	t += ki << (52 - EXP2F_TABLE_BITS); | ||
| 62 | 	s = asdouble(t); | ||
| 63 | 	z = C[0] * r + C[1]; | ||
| 64 | 	r2 = r * r; | ||
| 65 | 	y = C[2] * r + 1; | ||
| 66 | 	y = z * r2 + y; | ||
| 67 | 	y = y * s; | ||
| 68 | 	return eval_as_float(y); | ||
| 69 | } | ||
lib/libc/musl/src/math/expf.c deleted-80| ... | @@ -1,80 +0,0 @@ | ||
| 1 | /* | ||
| 2 | * Single-precision e^x function. | ||
| 3 | * | ||
| 4 | * Copyright (c) 2017-2018, Arm Limited. | ||
| 5 | * SPDX-License-Identifier: MIT | ||
| 6 | */ | ||
| 7 | |||
| 8 | #include <math.h> | ||
| 9 | #include <stdint.h> | ||
| 10 | #include "libm.h" | ||
| 11 | #include "exp2f_data.h" | ||
| 12 | |||
| 13 | /* | ||
| 14 | EXP2F_TABLE_BITS = 5 | ||
| 15 | EXP2F_POLY_ORDER = 3 | ||
| 16 | |||
| 17 | ULP error: 0.502 (nearest rounding.) | ||
| 18 | Relative error: 1.69 * 2^-34 in [-ln2/64, ln2/64] (before rounding.) | ||
| 19 | Wrong count: 170635 (all nearest rounding wrong results with fma.) | ||
| 20 | Non-nearest ULP error: 1 (rounded ULP error) | ||
| 21 | */ | ||
| 22 | |||
| 23 | #define N (1 << EXP2F_TABLE_BITS) | ||
| 24 | #define InvLn2N __exp2f_data.invln2_scaled | ||
| 25 | #define T __exp2f_data.tab | ||
| 26 | #define C __exp2f_data.poly_scaled | ||
| 27 | |||
| 28 | static inline uint32_t top12(float x) | ||
| 29 | { | ||
| 30 | 	return asuint(x) >> 20; | ||
| 31 | } | ||
| 32 | |||
| 33 | float expf(float x) | ||
| 34 | { | ||
| 35 | 	uint32_t abstop; | ||
| 36 | 	uint64_t ki, t; | ||
| 37 | 	double_t kd, xd, z, r, r2, y, s; | ||
| 38 | |||
| 39 | 	xd = (double_t)x; | ||
| 40 | 	abstop = top12(x) & 0x7ff; | ||
| 41 | 	if (predict_false(abstop >= top12(88.0f))) { | ||
| 42 | 		/* |x| >= 88 or x is nan. */ | ||
| 43 | 		if (asuint(x) == asuint(-INFINITY)) | ||
| 44 | 			return 0.0f; | ||
| 45 | 		if (abstop >= top12(INFINITY)) | ||
| 46 | 			return x + x; | ||
| 47 | 		if (x > 0x1.62e42ep6f) /* x > log(0x1p128) ~= 88.72 */ | ||
| 48 | 			return __math_oflowf(0); | ||
| 49 | 		if (x < -0x1.9fe368p6f) /* x < log(0x1p-150) ~= -103.97 */ | ||
| 50 | 			return __math_uflowf(0); | ||
| 51 | 	} | ||
| 52 | |||
| 53 | 	/* x*N/Ln2 = k + r with r in [-1/2, 1/2] and int k. */ | ||
| 54 | 	z = InvLn2N * xd; | ||
| 55 | |||
| 56 | 	/* Round and convert z to int, the result is in [-150*N, 128*N] and | ||
| 57 | 	 ideally ties-to-even rule is used, otherwise the magnitude of r | ||
| 58 | 	 can be bigger which gives larger approximation error. */ | ||
| 59 | #if TOINT_INTRINSICS | ||
| 60 | 	kd = roundtoint(z); | ||
| 61 | 	ki = converttoint(z); | ||
| 62 | #else | ||
| 63 | # define SHIFT __exp2f_data.shift | ||
| 64 | 	kd = eval_as_double(z + SHIFT); | ||
| 65 | 	ki = asuint64(kd); | ||
| 66 | 	kd -= SHIFT; | ||
| 67 | #endif | ||
| 68 | 	r = z - kd; | ||
| 69 | |||
| 70 | 	/* exp(x) = 2^(k/N) * 2^(r/N) ~= s * (C0*r^3 + C1*r^2 + C2*r + 1) */ | ||
| 71 | 	t = T[ki % N]; | ||
| 72 | 	t += ki << (52 - EXP2F_TABLE_BITS); | ||
| 73 | 	s = asdouble(t); | ||
| 74 | 	z = C[0] * r + C[1]; | ||
| 75 | 	r2 = r * r; | ||
| 76 | 	y = C[2] * r + 1; | ||
| 77 | 	y = z * r2 + y; | ||
| 78 | 	y = y * s; | ||
| 79 | 	return eval_as_float(y); | ||
| 80 | } | ||
lib/libc/musl/src/math/fmod.c deleted-68| ... | @@ -1,68 +0,0 @@ | ||
| 1 | #include <math.h> | ||
| 2 | #include <stdint.h> | ||
| 3 | |||
| 4 | double fmod(double x, double y) | ||
| 5 | { | ||
| 6 | 	union {double f; uint64_t i;} ux = {x}, uy = {y}; | ||
| 7 | 	int ex = ux.i>>52 & 0x7ff; | ||
| 8 | 	int ey = uy.i>>52 & 0x7ff; | ||
| 9 | 	int sx = ux.i>>63; | ||
| 10 | 	uint64_t i; | ||
| 11 | |||
| 12 | 	/* in the followings uxi should be ux.i, but then gcc wrongly adds */ | ||
| 13 | 	/* float load/store to inner loops ruining performance and code size */ | ||
| 14 | 	uint64_t uxi = ux.i; | ||
| 15 | |||
| 16 | 	if (uy.i<<1 == 0 || isnan(y) || ex == 0x7ff) | ||
| 17 | 		return (x*y)/(x*y); | ||
| 18 | 	if (uxi<<1 <= uy.i<<1) { | ||
| 19 | 		if (uxi<<1 == uy.i<<1) | ||
| 20 | 			return 0*x; | ||
| 21 | 		return x; | ||
| 22 | 	} | ||
| 23 | |||
| 24 | 	/* normalize x and y */ | ||
| 25 | 	if (!ex) { | ||
| 26 | 		for (i = uxi<<12; i>>63 == 0; ex--, i <<= 1); | ||
| 27 | 		uxi <<= -ex + 1; | ||
| 28 | 	} else { | ||
| 29 | 		uxi &= -1ULL >> 12; | ||
| 30 | 		uxi |= 1ULL << 52; | ||
| 31 | 	} | ||
| 32 | 	if (!ey) { | ||
| 33 | 		for (i = uy.i<<12; i>>63 == 0; ey--, i <<= 1); | ||
| 34 | 		uy.i <<= -ey + 1; | ||
| 35 | 	} else { | ||
| 36 | 		uy.i &= -1ULL >> 12; | ||
| 37 | 		uy.i |= 1ULL << 52; | ||
| 38 | 	} | ||
| 39 | |||
| 40 | 	/* x mod y */ | ||
| 41 | 	for (; ex > ey; ex--) { | ||
| 42 | 		i = uxi - uy.i; | ||
| 43 | 		if (i >> 63 == 0) { | ||
| 44 | 			if (i == 0) | ||
| 45 | 				return 0*x; | ||
| 46 | 			uxi = i; | ||
| 47 | 		} | ||
| 48 | 		uxi <<= 1; | ||
| 49 | 	} | ||
| 50 | 	i = uxi - uy.i; | ||
| 51 | 	if (i >> 63 == 0) { | ||
| 52 | 		if (i == 0) | ||
| 53 | 			return 0*x; | ||
| 54 | 		uxi = i; | ||
| 55 | 	} | ||
| 56 | 	for (; uxi>>52 == 0; uxi <<= 1, ex--); | ||
| 57 | |||
| 58 | 	/* scale result */ | ||
| 59 | 	if (ex > 0) { | ||
| 60 | 		uxi -= 1ULL << 52; | ||
| 61 | 		uxi |= (uint64_t)ex << 52; | ||
| 62 | 	} else { | ||
| 63 | 		uxi >>= -ex + 1; | ||
| 64 | 	} | ||
| 65 | 	uxi |= (uint64_t)sx << 63; | ||
| 66 | 	ux.i = uxi; | ||
| 67 | 	return ux.f; | ||
| 68 | } | ||
lib/libc/musl/src/math/fmodf.c deleted-65| ... | @@ -1,65 +0,0 @@ | ||
| 1 | #include <math.h> | ||
| 2 | #include <stdint.h> | ||
| 3 | |||
| 4 | float fmodf(float x, float y) | ||
| 5 | { | ||
| 6 | 	union {float f; uint32_t i;} ux = {x}, uy = {y}; | ||
| 7 | 	int ex = ux.i>>23 & 0xff; | ||
| 8 | 	int ey = uy.i>>23 & 0xff; | ||
| 9 | 	uint32_t sx = ux.i & 0x80000000; | ||
| 10 | 	uint32_t i; | ||
| 11 | 	uint32_t uxi = ux.i; | ||
| 12 | |||
| 13 | 	if (uy.i<<1 == 0 || isnan(y) || ex == 0xff) | ||
| 14 | 		return (x*y)/(x*y); | ||
| 15 | 	if (uxi<<1 <= uy.i<<1) { | ||
| 16 | 		if (uxi<<1 == uy.i<<1) | ||
| 17 | 			return 0*x; | ||
| 18 | 		return x; | ||
| 19 | 	} | ||
| 20 | |||
| 21 | 	/* normalize x and y */ | ||
| 22 | 	if (!ex) { | ||
| 23 | 		for (i = uxi<<9; i>>31 == 0; ex--, i <<= 1); | ||
| 24 | 		uxi <<= -ex + 1; | ||
| 25 | 	} else { | ||
| 26 | 		uxi &= -1U >> 9; | ||
| 27 | 		uxi |= 1U << 23; | ||
| 28 | 	} | ||
| 29 | 	if (!ey) { | ||
| 30 | 		for (i = uy.i<<9; i>>31 == 0; ey--, i <<= 1); | ||
| 31 | 		uy.i <<= -ey + 1; | ||
| 32 | 	} else { | ||
| 33 | 		uy.i &= -1U >> 9; | ||
| 34 | 		uy.i |= 1U << 23; | ||
| 35 | 	} | ||
| 36 | |||
| 37 | 	/* x mod y */ | ||
| 38 | 	for (; ex > ey; ex--) { | ||
| 39 | 		i = uxi - uy.i; | ||
| 40 | 		if (i >> 31 == 0) { | ||
| 41 | 			if (i == 0) | ||
| 42 | 				return 0*x; | ||
| 43 | 			uxi = i; | ||
| 44 | 		} | ||
| 45 | 		uxi <<= 1; | ||
| 46 | 	} | ||
| 47 | 	i = uxi - uy.i; | ||
| 48 | 	if (i >> 31 == 0) { | ||
| 49 | 		if (i == 0) | ||
| 50 | 			return 0*x; | ||
| 51 | 		uxi = i; | ||
| 52 | 	} | ||
| 53 | 	for (; uxi>>23 == 0; uxi <<= 1, ex--); | ||
| 54 | |||
| 55 | 	/* scale result up */ | ||
| 56 | 	if (ex > 0) { | ||
| 57 | 		uxi -= 1U << 23; | ||
| 58 | 		uxi |= (uint32_t)ex << 23; | ||
| 59 | 	} else { | ||
| 60 | 		uxi >>= -ex + 1; | ||
| 61 | 	} | ||
| 62 | 	uxi |= sx; | ||
| 63 | 	ux.i = uxi; | ||
| 64 | 	return ux.f; | ||
| 65 | } | ||
lib/libc/musl/src/math/fmodl.c deleted-105| ... | @@ -1,105 +0,0 @@ | ||
| 1 | #include "libm.h" | ||
| 2 | |||
| 3 | #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 | ||
| 4 | long double fmodl(long double x, long double y) | ||
| 5 | { | ||
| 6 | 	return fmod(x, y); | ||
| 7 | } | ||
| 8 | #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 | ||
| 9 | long double fmodl(long double x, long double y) | ||
| 10 | { | ||
| 11 | 	union ldshape ux = {x}, uy = {y}; | ||
| 12 | 	int ex = ux.i.se & 0x7fff; | ||
| 13 | 	int ey = uy.i.se & 0x7fff; | ||
| 14 | 	int sx = ux.i.se & 0x8000; | ||
| 15 | |||
| 16 | 	if (y == 0 || isnan(y) || ex == 0x7fff) | ||
| 17 | 		return (x*y)/(x*y); | ||
| 18 | 	ux.i.se = ex; | ||
| 19 | 	uy.i.se = ey; | ||
| 20 | 	if (ux.f <= uy.f) { | ||
| 21 | 		if (ux.f == uy.f) | ||
| 22 | 			return 0*x; | ||
| 23 | 		return x; | ||
| 24 | 	} | ||
| 25 | |||
| 26 | 	/* normalize x and y */ | ||
| 27 | 	if (!ex) { | ||
| 28 | 		ux.f *= 0x1p120f; | ||
| 29 | 		ex = ux.i.se - 120; | ||
| 30 | 	} | ||
| 31 | 	if (!ey) { | ||
| 32 | 		uy.f *= 0x1p120f; | ||
| 33 | 		ey = uy.i.se - 120; | ||
| 34 | 	} | ||
| 35 | |||
| 36 | 	/* x mod y */ | ||
| 37 | #if LDBL_MANT_DIG == 64 | ||
| 38 | 	uint64_t i, mx, my; | ||
| 39 | 	mx = ux.i.m; | ||
| 40 | 	my = uy.i.m; | ||
| 41 | 	for (; ex > ey; ex--) { | ||
| 42 | 		i = mx - my; | ||
| 43 | 		if (mx >= my) { | ||
| 44 | 			if (i == 0) | ||
| 45 | 				return 0*x; | ||
| 46 | 			mx = 2*i; | ||
| 47 | 		} else if (2*mx < mx) { | ||
| 48 | 			mx = 2*mx - my; | ||
| 49 | 		} else { | ||
| 50 | 			mx = 2*mx; | ||
| 51 | 		} | ||
| 52 | 	} | ||
| 53 | 	i = mx - my; | ||
| 54 | 	if (mx >= my) { | ||
| 55 | 		if (i == 0) | ||
| 56 | 			return 0*x; | ||
| 57 | 		mx = i; | ||
| 58 | 	} | ||
| 59 | 	for (; mx >> 63 == 0; mx *= 2, ex--); | ||
| 60 | 	ux.i.m = mx; | ||
| 61 | #elif LDBL_MANT_DIG == 113 | ||
| 62 | 	uint64_t hi, lo, xhi, xlo, yhi, ylo; | ||
| 63 | 	xhi = (ux.i2.hi & -1ULL>>16) | 1ULL<<48; | ||
| 64 | 	yhi = (uy.i2.hi & -1ULL>>16) | 1ULL<<48; | ||
| 65 | 	xlo = ux.i2.lo; | ||
| 66 | 	ylo = uy.i2.lo; | ||
| 67 | 	for (; ex > ey; ex--) { | ||
| 68 | 		hi = xhi - yhi; | ||
| 69 | 		lo = xlo - ylo; | ||
| 70 | 		if (xlo < ylo) | ||
| 71 | 			hi -= 1; | ||
| 72 | 		if (hi >> 63 == 0) { | ||
| 73 | 			if ((hi|lo) == 0) | ||
| 74 | 				return 0*x; | ||
| 75 | 			xhi = 2*hi + (lo>>63); | ||
| 76 | 			xlo = 2*lo; | ||
| 77 | 		} else { | ||
| 78 | 			xhi = 2*xhi + (xlo>>63); | ||
| 79 | 			xlo = 2*xlo; | ||
| 80 | 		} | ||
| 81 | 	} | ||
| 82 | 	hi = xhi - yhi; | ||
| 83 | 	lo = xlo - ylo; | ||
| 84 | 	if (xlo < ylo) | ||
| 85 | 		hi -= 1; | ||
| 86 | 	if (hi >> 63 == 0) { | ||
| 87 | 		if ((hi|lo) == 0) | ||
| 88 | 			return 0*x; | ||
| 89 | 		xhi = hi; | ||
| 90 | 		xlo = lo; | ||
| 91 | 	} | ||
| 92 | 	for (; xhi >> 48 == 0; xhi = 2*xhi + (xlo>>63), xlo = 2*xlo, ex--); | ||
| 93 | 	ux.i2.hi = xhi; | ||
| 94 | 	ux.i2.lo = xlo; | ||
| 95 | #endif | ||
| 96 | |||
| 97 | 	/* scale result */ | ||
| 98 | 	if (ex <= 0) { | ||
| 99 | 		ux.i.se = (ex+120)|sx; | ||
| 100 | 		ux.f *= 0x1p-120f; | ||
| 101 | 	} else | ||
| 102 | 		ux.i.se = ex|sx; | ||
| 103 | 	return ux.f; | ||
| 104 | } | ||
| 105 | #endif | ||
lib/libc/musl/src/math/log.c deleted-112| ... | @@ -1,112 +0,0 @@ | ||
| 1 | /* | ||
| 2 | * Double-precision log(x) function. | ||
| 3 | * | ||
| 4 | * Copyright (c) 2018, Arm Limited. | ||
| 5 | * SPDX-License-Identifier: MIT | ||
| 6 | */ | ||
| 7 | |||
| 8 | #include <math.h> | ||
| 9 | #include <stdint.h> | ||
| 10 | #include "libm.h" | ||
| 11 | #include "log_data.h" | ||
| 12 | |||
| 13 | #define T __log_data.tab | ||
| 14 | #define T2 __log_data.tab2 | ||
| 15 | #define B __log_data.poly1 | ||
| 16 | #define A __log_data.poly | ||
| 17 | #define Ln2hi __log_data.ln2hi | ||
| 18 | #define Ln2lo __log_data.ln2lo | ||
| 19 | #define N (1 << LOG_TABLE_BITS) | ||
| 20 | #define OFF 0x3fe6000000000000 | ||
| 21 | |||
| 22 | /* Top 16 bits of a double. */ | ||
| 23 | static inline uint32_t top16(double x) | ||
| 24 | { | ||
| 25 | 	return asuint64(x) >> 48; | ||
| 26 | } | ||
| 27 | |||
| 28 | double log(double x) | ||
| 29 | { | ||
| 30 | 	double_t w, z, r, r2, r3, y, invc, logc, kd, hi, lo; | ||
| 31 | 	uint64_t ix, iz, tmp; | ||
| 32 | 	uint32_t top; | ||
| 33 | 	int k, i; | ||
| 34 | |||
| 35 | 	ix = asuint64(x); | ||
| 36 | 	top = top16(x); | ||
| 37 | #define LO asuint64(1.0 - 0x1p-4) | ||
| 38 | #define HI asuint64(1.0 + 0x1.09p-4) | ||
| 39 | 	if (predict_false(ix - LO < HI - LO)) { | ||
| 40 | 		/* Handle close to 1.0 inputs separately. */ | ||
| 41 | 		/* Fix sign of zero with downward rounding when x==1. */ | ||
| 42 | 		if (WANT_ROUNDING && predict_false(ix == asuint64(1.0))) | ||
| 43 | 			return 0; | ||
| 44 | 		r = x - 1.0; | ||
| 45 | 		r2 = r * r; | ||
| 46 | 		r3 = r * r2; | ||
| 47 | 		y = r3 * | ||
| 48 | 		 (B[1] + r * B[2] + r2 * B[3] + | ||
| 49 | 		 r3 * (B[4] + r * B[5] + r2 * B[6] + | ||
| 50 | 			 r3 * (B[7] + r * B[8] + r2 * B[9] + r3 * B[10]))); | ||
| 51 | 		/* Worst-case error is around 0.507 ULP. */ | ||
| 52 | 		w = r * 0x1p27; | ||
| 53 | 		double_t rhi = r + w - w; | ||
| 54 | 		double_t rlo = r - rhi; | ||
| 55 | 		w = rhi * rhi * B[0]; /* B[0] == -0.5. */ | ||
| 56 | 		hi = r + w; | ||
| 57 | 		lo = r - hi + w; | ||
| 58 | 		lo += B[0] * rlo * (rhi + r); | ||
| 59 | 		y += lo; | ||
| 60 | 		y += hi; | ||
| 61 | 		return eval_as_double(y); | ||
| 62 | 	} | ||
| 63 | 	if (predict_false(top - 0x0010 >= 0x7ff0 - 0x0010)) { | ||
| 64 | 		/* x < 0x1p-1022 or inf or nan. */ | ||
| 65 | 		if (ix * 2 == 0) | ||
| 66 | 			return __math_divzero(1); | ||
| 67 | 		if (ix == asuint64(INFINITY)) /* log(inf) == inf. */ | ||
| 68 | 			return x; | ||
| 69 | 		if ((top & 0x8000) || (top & 0x7ff0) == 0x7ff0) | ||
| 70 | 			return __math_invalid(x); | ||
| 71 | 		/* x is subnormal, normalize it. */ | ||
| 72 | 		ix = asuint64(x * 0x1p52); | ||
| 73 | 		ix -= 52ULL << 52; | ||
| 74 | 	} | ||
| 75 | |||
| 76 | 	/* x = 2^k z; where z is in range [OFF,2*OFF) and exact. | ||
| 77 | 	 The range is split into N subintervals. | ||
| 78 | 	 The ith subinterval contains z and c is near its center. */ | ||
| 79 | 	tmp = ix - OFF; | ||
| 80 | 	i = (tmp >> (52 - LOG_TABLE_BITS)) % N; | ||
| 81 | 	k = (int64_t)tmp >> 52; /* arithmetic shift */ | ||
| 82 | 	iz = ix - (tmp & 0xfffULL << 52); | ||
| 83 | 	invc = T[i].invc; | ||
| 84 | 	logc = T[i].logc; | ||
| 85 | 	z = asdouble(iz); | ||
| 86 | |||
| 87 | 	/* log(x) = log1p(z/c-1) + log(c) + k*Ln2. */ | ||
| 88 | 	/* r ~= z/c - 1, |r| < 1/(2*N). */ | ||
| 89 | #if __FP_FAST_FMA | ||
| 90 | 	/* rounding error: 0x1p-55/N. */ | ||
| 91 | 	r = __builtin_fma(z, invc, -1.0); | ||
| 92 | #else | ||
| 93 | 	/* rounding error: 0x1p-55/N + 0x1p-66. */ | ||
| 94 | 	r = (z - T2[i].chi - T2[i].clo) * invc; | ||
| 95 | #endif | ||
| 96 | 	kd = (double_t)k; | ||
| 97 | |||
| 98 | 	/* hi + lo = r + log(c) + k*Ln2. */ | ||
| 99 | 	w = kd * Ln2hi + logc; | ||
| 100 | 	hi = w + r; | ||
| 101 | 	lo = w - hi + r + kd * Ln2lo; | ||
| 102 | |||
| 103 | 	/* log(x) = lo + (log1p(r) - r) + hi. */ | ||
| 104 | 	r2 = r * r; /* rounding error: 0x1p-54/N^2. */ | ||
| 105 | 	/* Worst case error if |y| > 0x1p-5: | ||
| 106 | 	 0.5 + 4.13/N + abs-poly-error*2^57 ULP (+ 0.002 ULP without fma) | ||
| 107 | 	 Worst case error if |y| > 0x1p-4: | ||
| 108 | 	 0.5 + 2.06/N + abs-poly-error*2^56 ULP (+ 0.001 ULP without fma). */ | ||
| 109 | 	y = lo + r2 * A[0] + | ||
| 110 | 	 r * r2 * (A[1] + r * A[2] + r2 * (A[3] + r * A[4])) + hi; | ||
| 111 | 	return eval_as_double(y); | ||
| 112 | } | ||
lib/libc/musl/src/math/log10.c deleted-101| ... | @@ -1,101 +0,0 @@ | ||
| 1 | /* origin: FreeBSD /usr/src/lib/msun/src/e_log10.c */ | ||
| 2 | /* | ||
| 3 | * ==================================================== | ||
| 4 | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. | ||
| 5 | * | ||
| 6 | * Developed at SunSoft, a Sun Microsystems, Inc. business. | ||
| 7 | * Permission to use, copy, modify, and distribute this | ||
| 8 | * software is freely granted, provided that this notice | ||
| 9 | * is preserved. | ||
| 10 | * ==================================================== | ||
| 11 | */ | ||
| 12 | /* | ||
| 13 | * Return the base 10 logarithm of x. See log.c for most comments. | ||
| 14 | * | ||
| 15 | * Reduce x to 2^k (1+f) and calculate r = log(1+f) - f + f*f/2 | ||
| 16 | * as in log.c, then combine and scale in extra precision: | ||
| 17 | * log10(x) = (f - f*f/2 + r)/log(10) + k*log10(2) | ||
| 18 | */ | ||
| 19 | |||
| 20 | #include <math.h> | ||
| 21 | #include <stdint.h> | ||
| 22 | |||
| 23 | static const double | ||
| 24 | ivln10hi = 4.34294481878168880939e-01, /* 0x3fdbcb7b, 0x15200000 */ | ||
| 25 | ivln10lo = 2.50829467116452752298e-11, /* 0x3dbb9438, 0xca9aadd5 */ | ||
| 26 | log10_2hi = 3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */ | ||
| 27 | log10_2lo = 3.69423907715893078616e-13, /* 0x3D59FEF3, 0x11F12B36 */ | ||
| 28 | Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */ | ||
| 29 | Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */ | ||
| 30 | Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */ | ||
| 31 | Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */ | ||
| 32 | Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */ | ||
| 33 | Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ | ||
| 34 | Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ | ||
| 35 | |||
| 36 | double log10(double x) | ||
| 37 | { | ||
| 38 | 	union {double f; uint64_t i;} u = {x}; | ||
| 39 | 	double_t hfsq,f,s,z,R,w,t1,t2,dk,y,hi,lo,val_hi,val_lo; | ||
| 40 | 	uint32_t hx; | ||
| 41 | 	int k; | ||
| 42 | |||
| 43 | 	hx = u.i>>32; | ||
| 44 | 	k = 0; | ||
| 45 | 	if (hx < 0x00100000 || hx>>31) { | ||
| 46 | 		if (u.i<<1 == 0) | ||
| 47 | 			return -1/(x*x); /* log(+-0)=-inf */ | ||
| 48 | 		if (hx>>31) | ||
| 49 | 			return (x-x)/0.0; /* log(-#) = NaN */ | ||
| 50 | 		/* subnormal number, scale x up */ | ||
| 51 | 		k -= 54; | ||
| 52 | 		x *= 0x1p54; | ||
| 53 | 		u.f = x; | ||
| 54 | 		hx = u.i>>32; | ||
| 55 | 	} else if (hx >= 0x7ff00000) { | ||
| 56 | 		return x; | ||
| 57 | 	} else if (hx == 0x3ff00000 && u.i<<32 == 0) | ||
| 58 | 		return 0; | ||
| 59 | |||
| 60 | 	/* reduce x into [sqrt(2)/2, sqrt(2)] */ | ||
| 61 | 	hx += 0x3ff00000 - 0x3fe6a09e; | ||
| 62 | 	k += (int)(hx>>20) - 0x3ff; | ||
| 63 | 	hx = (hx&0x000fffff) + 0x3fe6a09e; | ||
| 64 | 	u.i = (uint64_t)hx<<32 | (u.i&0xffffffff); | ||
| 65 | 	x = u.f; | ||
| 66 | |||
| 67 | 	f = x - 1.0; | ||
| 68 | 	hfsq = 0.5*f*f; | ||
| 69 | 	s = f/(2.0+f); | ||
| 70 | 	z = s*s; | ||
| 71 | 	w = z*z; | ||
| 72 | 	t1 = w*(Lg2+w*(Lg4+w*Lg6)); | ||
| 73 | 	t2 = z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); | ||
| 74 | 	R = t2 + t1; | ||
| 75 | |||
| 76 | 	/* See log2.c for details. */ | ||
| 77 | 	/* hi+lo = f - hfsq + s*(hfsq+R) ~ log(1+f) */ | ||
| 78 | 	hi = f - hfsq; | ||
| 79 | 	u.f = hi; | ||
| 80 | 	u.i &= (uint64_t)-1<<32; | ||
| 81 | 	hi = u.f; | ||
| 82 | 	lo = f - hi - hfsq + s*(hfsq+R); | ||
| 83 | |||
| 84 | 	/* val_hi+val_lo ~ log10(1+f) + k*log10(2) */ | ||
| 85 | 	val_hi = hi*ivln10hi; | ||
| 86 | 	dk = k; | ||
| 87 | 	y = dk*log10_2hi; | ||
| 88 | 	val_lo = dk*log10_2lo + (lo+hi)*ivln10lo + lo*ivln10hi; | ||
| 89 | |||
| 90 | 	/* | ||
| 91 | 	 * Extra precision in for adding y is not strictly needed | ||
| 92 | 	 * since there is no very large cancellation near x = sqrt(2) or | ||
| 93 | 	 * x = 1/sqrt(2), but we do it anyway since it costs little on CPUs | ||
| 94 | 	 * with some parallelism and it reduces the error for many args. | ||
| 95 | 	 */ | ||
| 96 | 	w = y + val_hi; | ||
| 97 | 	val_lo += (y - w) + val_hi; | ||
| 98 | 	val_hi = w; | ||
| 99 | |||
| 100 | 	return val_lo + val_hi; | ||
| 101 | } | ||
lib/libc/musl/src/math/log10f.c deleted-77| ... | @@ -1,77 +0,0 @@ | ||
| 1 | /* origin: FreeBSD /usr/src/lib/msun/src/e_log10f.c */ | ||
| 2 | /* | ||
| 3 | * ==================================================== | ||
| 4 | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. | ||
| 5 | * | ||
| 6 | * Developed at SunPro, a Sun Microsystems, Inc. business. | ||
| 7 | * Permission to use, copy, modify, and distribute this | ||
| 8 | * software is freely granted, provided that this notice | ||
| 9 | * is preserved. | ||
| 10 | * ==================================================== | ||
| 11 | */ | ||
| 12 | /* | ||
| 13 | * See comments in log10.c. | ||
| 14 | */ | ||
| 15 | |||
| 16 | #include <math.h> | ||
| 17 | #include <stdint.h> | ||
| 18 | |||
| 19 | static const float | ||
| 20 | ivln10hi = 4.3432617188e-01, /* 0x3ede6000 */ | ||
| 21 | ivln10lo = -3.1689971365e-05, /* 0xb804ead9 */ | ||
| 22 | log10_2hi = 3.0102920532e-01, /* 0x3e9a2080 */ | ||
| 23 | log10_2lo = 7.9034151668e-07, /* 0x355427db */ | ||
| 24 | /* |(log(1+s)-log(1-s))/s - Lg(s)| < 2**-34.24 (~[-4.95e-11, 4.97e-11]). */ | ||
| 25 | Lg1 = 0xaaaaaa.0p-24, /* 0.66666662693 */ | ||
| 26 | Lg2 = 0xccce13.0p-25, /* 0.40000972152 */ | ||
| 27 | Lg3 = 0x91e9ee.0p-25, /* 0.28498786688 */ | ||
| 28 | Lg4 = 0xf89e26.0p-26; /* 0.24279078841 */ | ||
| 29 | |||
| 30 | float log10f(float x) | ||
| 31 | { | ||
| 32 | 	union {float f; uint32_t i;} u = {x}; | ||
| 33 | 	float_t hfsq,f,s,z,R,w,t1,t2,dk,hi,lo; | ||
| 34 | 	uint32_t ix; | ||
| 35 | 	int k; | ||
| 36 | |||
| 37 | 	ix = u.i; | ||
| 38 | 	k = 0; | ||
| 39 | 	if (ix < 0x00800000 || ix>>31) { /* x < 2**-126 */ | ||
| 40 | 		if (ix<<1 == 0) | ||
| 41 | 			return -1/(x*x); /* log(+-0)=-inf */ | ||
| 42 | 		if (ix>>31) | ||
| 43 | 			return (x-x)/0.0f; /* log(-#) = NaN */ | ||
| 44 | 		/* subnormal number, scale up x */ | ||
| 45 | 		k -= 25; | ||
| 46 | 		x *= 0x1p25f; | ||
| 47 | 		u.f = x; | ||
| 48 | 		ix = u.i; | ||
| 49 | 	} else if (ix >= 0x7f800000) { | ||
| 50 | 		return x; | ||
| 51 | 	} else if (ix == 0x3f800000) | ||
| 52 | 		return 0; | ||
| 53 | |||
| 54 | 	/* reduce x into [sqrt(2)/2, sqrt(2)] */ | ||
| 55 | 	ix += 0x3f800000 - 0x3f3504f3; | ||
| 56 | 	k += (int)(ix>>23) - 0x7f; | ||
| 57 | 	ix = (ix&0x007fffff) + 0x3f3504f3; | ||
| 58 | 	u.i = ix; | ||
| 59 | 	x = u.f; | ||
| 60 | |||
| 61 | 	f = x - 1.0f; | ||
| 62 | 	s = f/(2.0f + f); | ||
| 63 | 	z = s*s; | ||
| 64 | 	w = z*z; | ||
| 65 | 	t1= w*(Lg2+w*Lg4); | ||
| 66 | 	t2= z*(Lg1+w*Lg3); | ||
| 67 | 	R = t2 + t1; | ||
| 68 | 	hfsq = 0.5f*f*f; | ||
| 69 | |||
| 70 | 	hi = f - hfsq; | ||
| 71 | 	u.f = hi; | ||
| 72 | 	u.i &= 0xfffff000; | ||
| 73 | 	hi = u.f; | ||
| 74 | 	lo = f - hi - hfsq + s*(hfsq+R); | ||
| 75 | 	dk = k; | ||
| 76 | 	return dk*log10_2lo + (lo+hi)*ivln10lo + lo*ivln10hi + hi*ivln10hi + dk*log10_2hi; | ||
| 77 | } | ||
lib/libc/musl/src/math/log2.c deleted-122| ... | @@ -1,122 +0,0 @@ | ||
| 1 | /* | ||
| 2 | * Double-precision log2(x) function. | ||
| 3 | * | ||
| 4 | * Copyright (c) 2018, Arm Limited. | ||
| 5 | * SPDX-License-Identifier: MIT | ||
| 6 | */ | ||
| 7 | |||
| 8 | #include <math.h> | ||
| 9 | #include <stdint.h> | ||
| 10 | #include "libm.h" | ||
| 11 | #include "log2_data.h" | ||
| 12 | |||
| 13 | #define T __log2_data.tab | ||
| 14 | #define T2 __log2_data.tab2 | ||
| 15 | #define B __log2_data.poly1 | ||
| 16 | #define A __log2_data.poly | ||
| 17 | #define InvLn2hi __log2_data.invln2hi | ||
| 18 | #define InvLn2lo __log2_data.invln2lo | ||
| 19 | #define N (1 << LOG2_TABLE_BITS) | ||
| 20 | #define OFF 0x3fe6000000000000 | ||
| 21 | |||
| 22 | /* Top 16 bits of a double. */ | ||
| 23 | static inline uint32_t top16(double x) | ||
| 24 | { | ||
| 25 | 	return asuint64(x) >> 48; | ||
| 26 | } | ||
| 27 | |||
| 28 | double log2(double x) | ||
| 29 | { | ||
| 30 | 	double_t z, r, r2, r4, y, invc, logc, kd, hi, lo, t1, t2, t3, p; | ||
| 31 | 	uint64_t ix, iz, tmp; | ||
| 32 | 	uint32_t top; | ||
| 33 | 	int k, i; | ||
| 34 | |||
| 35 | 	ix = asuint64(x); | ||
| 36 | 	top = top16(x); | ||
| 37 | #define LO asuint64(1.0 - 0x1.5b51p-5) | ||
| 38 | #define HI asuint64(1.0 + 0x1.6ab2p-5) | ||
| 39 | 	if (predict_false(ix - LO < HI - LO)) { | ||
| 40 | 		/* Handle close to 1.0 inputs separately. */ | ||
| 41 | 		/* Fix sign of zero with downward rounding when x==1. */ | ||
| 42 | 		if (WANT_ROUNDING && predict_false(ix == asuint64(1.0))) | ||
| 43 | 			return 0; | ||
| 44 | 		r = x - 1.0; | ||
| 45 | #if __FP_FAST_FMA | ||
| 46 | 		hi = r * InvLn2hi; | ||
| 47 | 		lo = r * InvLn2lo + __builtin_fma(r, InvLn2hi, -hi); | ||
| 48 | #else | ||
| 49 | 		double_t rhi, rlo; | ||
| 50 | 		rhi = asdouble(asuint64(r) & -1ULL << 32); | ||
| 51 | 		rlo = r - rhi; | ||
| 52 | 		hi = rhi * InvLn2hi; | ||
| 53 | 		lo = rlo * InvLn2hi + r * InvLn2lo; | ||
| 54 | #endif | ||
| 55 | 		r2 = r * r; /* rounding error: 0x1p-62. */ | ||
| 56 | 		r4 = r2 * r2; | ||
| 57 | 		/* Worst-case error is less than 0.54 ULP (0.55 ULP without fma). */ | ||
| 58 | 		p = r2 * (B[0] + r * B[1]); | ||
| 59 | 		y = hi + p; | ||
| 60 | 		lo += hi - y + p; | ||
| 61 | 		lo += r4 * (B[2] + r * B[3] + r2 * (B[4] + r * B[5]) + | ||
| 62 | 			 r4 * (B[6] + r * B[7] + r2 * (B[8] + r * B[9]))); | ||
| 63 | 		y += lo; | ||
| 64 | 		return eval_as_double(y); | ||
| 65 | 	} | ||
| 66 | 	if (predict_false(top - 0x0010 >= 0x7ff0 - 0x0010)) { | ||
| 67 | 		/* x < 0x1p-1022 or inf or nan. */ | ||
| 68 | 		if (ix * 2 == 0) | ||
| 69 | 			return __math_divzero(1); | ||
| 70 | 		if (ix == asuint64(INFINITY)) /* log(inf) == inf. */ | ||
| 71 | 			return x; | ||
| 72 | 		if ((top & 0x8000) || (top & 0x7ff0) == 0x7ff0) | ||
| 73 | 			return __math_invalid(x); | ||
| 74 | 		/* x is subnormal, normalize it. */ | ||
| 75 | 		ix = asuint64(x * 0x1p52); | ||
| 76 | 		ix -= 52ULL << 52; | ||
| 77 | 	} | ||
| 78 | |||
| 79 | 	/* x = 2^k z; where z is in range [OFF,2*OFF) and exact. | ||
| 80 | 	 The range is split into N subintervals. | ||
| 81 | 	 The ith subinterval contains z and c is near its center. */ | ||
| 82 | 	tmp = ix - OFF; | ||
| 83 | 	i = (tmp >> (52 - LOG2_TABLE_BITS)) % N; | ||
| 84 | 	k = (int64_t)tmp >> 52; /* arithmetic shift */ | ||
| 85 | 	iz = ix - (tmp & 0xfffULL << 52); | ||
| 86 | 	invc = T[i].invc; | ||
| 87 | 	logc = T[i].logc; | ||
| 88 | 	z = asdouble(iz); | ||
| 89 | 	kd = (double_t)k; | ||
| 90 | |||
| 91 | 	/* log2(x) = log2(z/c) + log2(c) + k. */ | ||
| 92 | 	/* r ~= z/c - 1, |r| < 1/(2*N). */ | ||
| 93 | #if __FP_FAST_FMA | ||
| 94 | 	/* rounding error: 0x1p-55/N. */ | ||
| 95 | 	r = __builtin_fma(z, invc, -1.0); | ||
| 96 | 	t1 = r * InvLn2hi; | ||
| 97 | 	t2 = r * InvLn2lo + __builtin_fma(r, InvLn2hi, -t1); | ||
| 98 | #else | ||
| 99 | 	double_t rhi, rlo; | ||
| 100 | 	/* rounding error: 0x1p-55/N + 0x1p-65. */ | ||
| 101 | 	r = (z - T2[i].chi - T2[i].clo) * invc; | ||
| 102 | 	rhi = asdouble(asuint64(r) & -1ULL << 32); | ||
| 103 | 	rlo = r - rhi; | ||
| 104 | 	t1 = rhi * InvLn2hi; | ||
| 105 | 	t2 = rlo * InvLn2hi + r * InvLn2lo; | ||
| 106 | #endif | ||
| 107 | |||
| 108 | 	/* hi + lo = r/ln2 + log2(c) + k. */ | ||
| 109 | 	t3 = kd + logc; | ||
| 110 | 	hi = t3 + t1; | ||
| 111 | 	lo = t3 - hi + t1 + t2; | ||
| 112 | |||
| 113 | 	/* log2(r+1) = r/ln2 + r^2*poly(r). */ | ||
| 114 | 	/* Evaluation is optimized assuming superscalar pipelined execution. */ | ||
| 115 | 	r2 = r * r; /* rounding error: 0x1p-54/N^2. */ | ||
| 116 | 	r4 = r2 * r2; | ||
| 117 | 	/* Worst-case error if |y| > 0x1p-4: 0.547 ULP (0.550 ULP without fma). | ||
| 118 | 	 ~ 0.5 + 2/N/ln2 + abs-poly-error*0x1p56 ULP (+ 0.003 ULP without fma). */ | ||
| 119 | 	p = A[0] + r * A[1] + r2 * (A[2] + r * A[3]) + r4 * (A[4] + r * A[5]); | ||
| 120 | 	y = lo + r2 * p + hi; | ||
| 121 | 	return eval_as_double(y); | ||
| 122 | } | ||
lib/libc/musl/src/math/log2f.c deleted-72| ... | @@ -1,72 +0,0 @@ | ||
| 1 | /* | ||
| 2 | * Single-precision log2 function. | ||
| 3 | * | ||
| 4 | * Copyright (c) 2017-2018, Arm Limited. | ||
| 5 | * SPDX-License-Identifier: MIT | ||
| 6 | */ | ||
| 7 | |||
| 8 | #include <math.h> | ||
| 9 | #include <stdint.h> | ||
| 10 | #include "libm.h" | ||
| 11 | #include "log2f_data.h" | ||
| 12 | |||
| 13 | /* | ||
| 14 | LOG2F_TABLE_BITS = 4 | ||
| 15 | LOG2F_POLY_ORDER = 4 | ||
| 16 | |||
| 17 | ULP error: 0.752 (nearest rounding.) | ||
| 18 | Relative error: 1.9 * 2^-26 (before rounding.) | ||
| 19 | */ | ||
| 20 | |||
| 21 | #define N (1 << LOG2F_TABLE_BITS) | ||
| 22 | #define T __log2f_data.tab | ||
| 23 | #define A __log2f_data.poly | ||
| 24 | #define OFF 0x3f330000 | ||
| 25 | |||
| 26 | float log2f(float x) | ||
| 27 | { | ||
| 28 | 	double_t z, r, r2, p, y, y0, invc, logc; | ||
| 29 | 	uint32_t ix, iz, top, tmp; | ||
| 30 | 	int k, i; | ||
| 31 | |||
| 32 | 	ix = asuint(x); | ||
| 33 | 	/* Fix sign of zero with downward rounding when x==1. */ | ||
| 34 | 	if (WANT_ROUNDING && predict_false(ix == 0x3f800000)) | ||
| 35 | 		return 0; | ||
| 36 | 	if (predict_false(ix - 0x00800000 >= 0x7f800000 - 0x00800000)) { | ||
| 37 | 		/* x < 0x1p-126 or inf or nan. */ | ||
| 38 | 		if (ix * 2 == 0) | ||
| 39 | 			return __math_divzerof(1); | ||
| 40 | 		if (ix == 0x7f800000) /* log2(inf) == inf. */ | ||
| 41 | 			return x; | ||
| 42 | 		if ((ix & 0x80000000) || ix * 2 >= 0xff000000) | ||
| 43 | 			return __math_invalidf(x); | ||
| 44 | 		/* x is subnormal, normalize it. */ | ||
| 45 | 		ix = asuint(x * 0x1p23f); | ||
| 46 | 		ix -= 23 << 23; | ||
| 47 | 	} | ||
| 48 | |||
| 49 | 	/* x = 2^k z; where z is in range [OFF,2*OFF] and exact. | ||
| 50 | 	 The range is split into N subintervals. | ||
| 51 | 	 The ith subinterval contains z and c is near its center. */ | ||
| 52 | 	tmp = ix - OFF; | ||
| 53 | 	i = (tmp >> (23 - LOG2F_TABLE_BITS)) % N; | ||
| 54 | 	top = tmp & 0xff800000; | ||
| 55 | 	iz = ix - top; | ||
| 56 | 	k = (int32_t)tmp >> 23; /* arithmetic shift */ | ||
| 57 | 	invc = T[i].invc; | ||
| 58 | 	logc = T[i].logc; | ||
| 59 | 	z = (double_t)asfloat(iz); | ||
| 60 | |||
| 61 | 	/* log2(x) = log1p(z/c-1)/ln2 + log2(c) + k */ | ||
| 62 | 	r = z * invc - 1; | ||
| 63 | 	y0 = logc + (double_t)k; | ||
| 64 | |||
| 65 | 	/* Pipelined polynomial evaluation to approximate log1p(r)/ln2. */ | ||
| 66 | 	r2 = r * r; | ||
| 67 | 	y = A[1] * r + A[2]; | ||
| 68 | 	y = A[0] * r2 + y; | ||
| 69 | 	p = A[3] * r + y0; | ||
| 70 | 	y = y * r2 + p; | ||
| 71 | 	return eval_as_float(y); | ||
| 72 | } | ||
lib/libc/musl/src/math/logf.c deleted-71| ... | @@ -1,71 +0,0 @@ | ||
| 1 | /* | ||
| 2 | * Single-precision log function. | ||
| 3 | * | ||
| 4 | * Copyright (c) 2017-2018, Arm Limited. | ||
| 5 | * SPDX-License-Identifier: MIT | ||
| 6 | */ | ||
| 7 | |||
| 8 | #include <math.h> | ||
| 9 | #include <stdint.h> | ||
| 10 | #include "libm.h" | ||
| 11 | #include "logf_data.h" | ||
| 12 | |||
| 13 | /* | ||
| 14 | LOGF_TABLE_BITS = 4 | ||
| 15 | LOGF_POLY_ORDER = 4 | ||
| 16 | |||
| 17 | ULP error: 0.818 (nearest rounding.) | ||
| 18 | Relative error: 1.957 * 2^-26 (before rounding.) | ||
| 19 | */ | ||
| 20 | |||
| 21 | #define T __logf_data.tab | ||
| 22 | #define A __logf_data.poly | ||
| 23 | #define Ln2 __logf_data.ln2 | ||
| 24 | #define N (1 << LOGF_TABLE_BITS) | ||
| 25 | #define OFF 0x3f330000 | ||
| 26 | |||
| 27 | float logf(float x) | ||
| 28 | { | ||
| 29 | 	double_t z, r, r2, y, y0, invc, logc; | ||
| 30 | 	uint32_t ix, iz, tmp; | ||
| 31 | 	int k, i; | ||
| 32 | |||
| 33 | 	ix = asuint(x); | ||
| 34 | 	/* Fix sign of zero with downward rounding when x==1. */ | ||
| 35 | 	if (WANT_ROUNDING && predict_false(ix == 0x3f800000)) | ||
| 36 | 		return 0; | ||
| 37 | 	if (predict_false(ix - 0x00800000 >= 0x7f800000 - 0x00800000)) { | ||
| 38 | 		/* x < 0x1p-126 or inf or nan. */ | ||
| 39 | 		if (ix * 2 == 0) | ||
| 40 | 			return __math_divzerof(1); | ||
| 41 | 		if (ix == 0x7f800000) /* log(inf) == inf. */ | ||
| 42 | 			return x; | ||
| 43 | 		if ((ix & 0x80000000) || ix * 2 >= 0xff000000) | ||
| 44 | 			return __math_invalidf(x); | ||
| 45 | 		/* x is subnormal, normalize it. */ | ||
| 46 | 		ix = asuint(x * 0x1p23f); | ||
| 47 | 		ix -= 23 << 23; | ||
| 48 | 	} | ||
| 49 | |||
| 50 | 	/* x = 2^k z; where z is in range [OFF,2*OFF] and exact. | ||
| 51 | 	 The range is split into N subintervals. | ||
| 52 | 	 The ith subinterval contains z and c is near its center. */ | ||
| 53 | 	tmp = ix - OFF; | ||
| 54 | 	i = (tmp >> (23 - LOGF_TABLE_BITS)) % N; | ||
| 55 | 	k = (int32_t)tmp >> 23; /* arithmetic shift */ | ||
| 56 | 	iz = ix - (tmp & 0xff800000); | ||
| 57 | 	invc = T[i].invc; | ||
| 58 | 	logc = T[i].logc; | ||
| 59 | 	z = (double_t)asfloat(iz); | ||
| 60 | |||
| 61 | 	/* log(x) = log1p(z/c-1) + log(c) + k*Ln2 */ | ||
| 62 | 	r = z * invc - 1; | ||
| 63 | 	y0 = logc + (double_t)k * Ln2; | ||
| 64 | |||
| 65 | 	/* Pipelined polynomial evaluation to approximate log1p(r). */ | ||
| 66 | 	r2 = r * r; | ||
| 67 | 	y = A[1] * r + A[2]; | ||
| 68 | 	y = A[0] * r2 + y; | ||
| 69 | 	y = y * r2 + (y0 + r); | ||
| 70 | 	return eval_as_float(y); | ||
| 71 | } | ||
lib/libc/musl/src/string/bcmp.c deleted-8| ... | @@ -1,8 +0,0 @@ | ||
| 1 | #define _BSD_SOURCE | ||
| 2 | #include <string.h> | ||
| 3 | #include <strings.h> | ||
| 4 | |||
| 5 | int bcmp(const void *s1, const void *s2, size_t n) | ||
| 6 | { | ||
| 7 | 	return memcmp(s1, s2, n); | ||
| 8 | } | ||
src/libs/mingw.zig-3| ... | @@ -998,9 +998,6 @@ const mingw32_x86_src = [_][]const u8{ | ... | @@ -998,9 +998,6 @@ const mingw32_x86_src = [_][]const u8{ |
| 998 | const mingw32_x86_32_src = [_][]const u8{ | 998 | const mingw32_x86_32_src = [_][]const u8{ |
| 999 | // ucrtbase | 999 | // ucrtbase |
| 1000 | "math" ++ path.sep_str ++ "coshf.c", | 1000 | "math" ++ path.sep_str ++ "coshf.c", |
| 1001 | "math" ++ path.sep_str ++ "expf.c", | ||
| 1002 | "math" ++ path.sep_str ++ "log10f.c", | ||
| 1003 | "math" ++ path.sep_str ++ "logf.c", | ||
| 1004 | "math" ++ path.sep_str ++ "modff.c", | 1001 | "math" ++ path.sep_str ++ "modff.c", |
| 1005 | "math" ++ path.sep_str ++ "powf.c", | 1002 | "math" ++ path.sep_str ++ "powf.c", |
| 1006 | "math" ++ path.sep_str ++ "sinhf.c", | 1003 | "math" ++ path.sep_str ++ "sinhf.c", |
src/libs/musl.zig-14| ... | @@ -900,13 +900,9 @@ const src_files = [_][]const u8{ | ... | @@ -900,13 +900,9 @@ const src_files = [_][]const u8{ |
| 900 | "musl/src/math/exp10.c", | 900 | "musl/src/math/exp10.c", |
| 901 | "musl/src/math/exp10f.c", | 901 | "musl/src/math/exp10f.c", |
| 902 | "musl/src/math/exp10l.c", | 902 | "musl/src/math/exp10l.c", |
| 903 | "musl/src/math/exp2.c", | ||
| 904 | "musl/src/math/exp2f.c", | ||
| 905 | "musl/src/math/exp2f_data.c", | 903 | "musl/src/math/exp2f_data.c", |
| 906 | "musl/src/math/exp2l.c", | 904 | "musl/src/math/exp2l.c", |
| 907 | "musl/src/math/exp.c", | ||
| 908 | "musl/src/math/exp_data.c", | 905 | "musl/src/math/exp_data.c", |
| 909 | "musl/src/math/expf.c", | ||
| 910 | "musl/src/math/expl.c", | 906 | "musl/src/math/expl.c", |
| 911 | "musl/src/math/expm1.c", | 907 | "musl/src/math/expm1.c", |
| 912 | "musl/src/math/expm1f.c", | 908 | "musl/src/math/expm1f.c", |
| ... | @@ -927,9 +923,6 @@ const src_files = [_][]const u8{ | ... | @@ -927,9 +923,6 @@ const src_files = [_][]const u8{ |
| 927 | "musl/src/math/fmin.c", | 923 | "musl/src/math/fmin.c", |
| 928 | "musl/src/math/fminf.c", | 924 | "musl/src/math/fminf.c", |
| 929 | "musl/src/math/fminl.c", | 925 | "musl/src/math/fminl.c", |
| 930 | "musl/src/math/fmod.c", | ||
| 931 | "musl/src/math/fmodf.c", | ||
| 932 | "musl/src/math/fmodl.c", | ||
| 933 | "musl/src/math/__fpclassify.c", | 926 | "musl/src/math/__fpclassify.c", |
| 934 | "musl/src/math/__fpclassifyf.c", | 927 | "musl/src/math/__fpclassifyf.c", |
| 935 | "musl/src/math/__fpclassifyl.c", | 928 | "musl/src/math/__fpclassifyl.c", |
| ... | @@ -1024,23 +1017,17 @@ const src_files = [_][]const u8{ | ... | @@ -1024,23 +1017,17 @@ const src_files = [_][]const u8{ |
| 1024 | "musl/src/math/llround.c", | 1017 | "musl/src/math/llround.c", |
| 1025 | "musl/src/math/llroundf.c", | 1018 | "musl/src/math/llroundf.c", |
| 1026 | "musl/src/math/llroundl.c", | 1019 | "musl/src/math/llroundl.c", |
| 1027 | "musl/src/math/log10.c", | ||
| 1028 | "musl/src/math/log10f.c", | ||
| 1029 | "musl/src/math/log10l.c", | 1020 | "musl/src/math/log10l.c", |
| 1030 | "musl/src/math/log1p.c", | 1021 | "musl/src/math/log1p.c", |
| 1031 | "musl/src/math/log1pf.c", | 1022 | "musl/src/math/log1pf.c", |
| 1032 | "musl/src/math/log1pl.c", | 1023 | "musl/src/math/log1pl.c", |
| 1033 | "musl/src/math/log2.c", | ||
| 1034 | "musl/src/math/log2_data.c", | 1024 | "musl/src/math/log2_data.c", |
| 1035 | "musl/src/math/log2f.c", | ||
| 1036 | "musl/src/math/log2f_data.c", | 1025 | "musl/src/math/log2f_data.c", |
| 1037 | "musl/src/math/log2l.c", | 1026 | "musl/src/math/log2l.c", |
| 1038 | "musl/src/math/logb.c", | 1027 | "musl/src/math/logb.c", |
| 1039 | "musl/src/math/logbf.c", | 1028 | "musl/src/math/logbf.c", |
| 1040 | "musl/src/math/logbl.c", | 1029 | "musl/src/math/logbl.c", |
| 1041 | "musl/src/math/log.c", | ||
| 1042 | "musl/src/math/log_data.c", | 1030 | "musl/src/math/log_data.c", |
| 1043 | "musl/src/math/logf.c", | ||
| 1044 | "musl/src/math/logf_data.c", | 1031 | "musl/src/math/logf_data.c", |
| 1045 | "musl/src/math/logl.c", | 1032 | "musl/src/math/logl.c", |
| 1046 | "musl/src/math/lrint.c", | 1033 | "musl/src/math/lrint.c", |
| ... | @@ -1752,7 +1739,6 @@ const src_files = [_][]const u8{ | ... | @@ -1752,7 +1739,6 @@ const src_files = [_][]const u8{ |
| 1752 | "musl/src/stdlib/strtol.c", | 1739 | "musl/src/stdlib/strtol.c", |
| 1753 | "musl/src/stdlib/wcstod.c", | 1740 | "musl/src/stdlib/wcstod.c", |
| 1754 | "musl/src/stdlib/wcstol.c", | 1741 | "musl/src/stdlib/wcstol.c", |
| 1755 | "musl/src/string/bcmp.c", | ||
| 1756 | "musl/src/string/bcopy.c", | 1742 | "musl/src/string/bcopy.c", |
| 1757 | "musl/src/string/explicit_bzero.c", | 1743 | "musl/src/string/explicit_bzero.c", |
| 1758 | "musl/src/string/index.c", | 1744 | "musl/src/string/index.c", |
src/libs/wasi_libc.zig-14| ... | @@ -739,13 +739,9 @@ const libc_top_half_src_files = [_][]const u8{ | ... | @@ -739,13 +739,9 @@ const libc_top_half_src_files = [_][]const u8{ |
| 739 | "musl/src/math/exp10.c", | 739 | "musl/src/math/exp10.c", |
| 740 | "musl/src/math/exp10f.c", | 740 | "musl/src/math/exp10f.c", |
| 741 | "musl/src/math/exp10l.c", | 741 | "musl/src/math/exp10l.c", |
| 742 | "musl/src/math/exp2.c", | ||
| 743 | "musl/src/math/exp2f.c", | ||
| 744 | "musl/src/math/exp2f_data.c", | 742 | "musl/src/math/exp2f_data.c", |
| 745 | "musl/src/math/exp2l.c", | 743 | "musl/src/math/exp2l.c", |
| 746 | "musl/src/math/exp.c", | ||
| 747 | "musl/src/math/exp_data.c", | 744 | "musl/src/math/exp_data.c", |
| 748 | "musl/src/math/expf.c", | ||
| 749 | "musl/src/math/expl.c", | 745 | "musl/src/math/expl.c", |
| 750 | "musl/src/math/expm1.c", | 746 | "musl/src/math/expm1.c", |
| 751 | "musl/src/math/expm1f.c", | 747 | "musl/src/math/expm1f.c", |
| ... | @@ -759,9 +755,6 @@ const libc_top_half_src_files = [_][]const u8{ | ... | @@ -759,9 +755,6 @@ const libc_top_half_src_files = [_][]const u8{ |
| 759 | "musl/src/math/fmaf.c", | 755 | "musl/src/math/fmaf.c", |
| 760 | "musl/src/math/fmaxl.c", | 756 | "musl/src/math/fmaxl.c", |
| 761 | "musl/src/math/fminl.c", | 757 | "musl/src/math/fminl.c", |
| 762 | "musl/src/math/fmod.c", | ||
| 763 | "musl/src/math/fmodf.c", | ||
| 764 | "musl/src/math/fmodl.c", | ||
| 765 | "musl/src/math/frexp.c", | 758 | "musl/src/math/frexp.c", |
| 766 | "musl/src/math/frexpf.c", | 759 | "musl/src/math/frexpf.c", |
| 767 | "musl/src/math/frexpl.c", | 760 | "musl/src/math/frexpl.c", |
| ... | @@ -792,23 +785,17 @@ const libc_top_half_src_files = [_][]const u8{ | ... | @@ -792,23 +785,17 @@ const libc_top_half_src_files = [_][]const u8{ |
| 792 | "musl/src/math/llround.c", | 785 | "musl/src/math/llround.c", |
| 793 | "musl/src/math/llroundf.c", | 786 | "musl/src/math/llroundf.c", |
| 794 | "musl/src/math/llroundl.c", | 787 | "musl/src/math/llroundl.c", |
| 795 | "musl/src/math/log10.c", | ||
| 796 | "musl/src/math/log10f.c", | ||
| 797 | "musl/src/math/log10l.c", | 788 | "musl/src/math/log10l.c", |
| 798 | "musl/src/math/log1p.c", | 789 | "musl/src/math/log1p.c", |
| 799 | "musl/src/math/log1pf.c", | 790 | "musl/src/math/log1pf.c", |
| 800 | "musl/src/math/log1pl.c", | 791 | "musl/src/math/log1pl.c", |
| 801 | "musl/src/math/log2.c", | ||
| 802 | "musl/src/math/log2_data.c", | 792 | "musl/src/math/log2_data.c", |
| 803 | "musl/src/math/log2f.c", | ||
| 804 | "musl/src/math/log2f_data.c", | 793 | "musl/src/math/log2f_data.c", |
| 805 | "musl/src/math/log2l.c", | 794 | "musl/src/math/log2l.c", |
| 806 | "musl/src/math/logb.c", | 795 | "musl/src/math/logb.c", |
| 807 | "musl/src/math/logbf.c", | 796 | "musl/src/math/logbf.c", |
| 808 | "musl/src/math/logbl.c", | 797 | "musl/src/math/logbl.c", |
| 809 | "musl/src/math/log.c", | ||
| 810 | "musl/src/math/log_data.c", | 798 | "musl/src/math/log_data.c", |
| 811 | "musl/src/math/logf.c", | ||
| 812 | "musl/src/math/logf_data.c", | 799 | "musl/src/math/logf_data.c", |
| 813 | "musl/src/math/logl.c", | 800 | "musl/src/math/logl.c", |
| 814 | "musl/src/math/lrint.c", | 801 | "musl/src/math/lrint.c", |
| ... | @@ -1031,7 +1018,6 @@ const libc_top_half_src_files = [_][]const u8{ | ... | @@ -1031,7 +1018,6 @@ const libc_top_half_src_files = [_][]const u8{ |
| 1031 | "musl/src/stdlib/qsort.c", | 1018 | "musl/src/stdlib/qsort.c", |
| 1032 | "musl/src/stdlib/qsort_nr.c", | 1019 | "musl/src/stdlib/qsort_nr.c", |
| 1033 | "musl/src/stdlib/strtol.c", | 1020 | "musl/src/stdlib/strtol.c", |
| 1034 | "musl/src/string/bcmp.c", | ||
| 1035 | "musl/src/string/bcopy.c", | 1021 | "musl/src/string/bcopy.c", |
| 1036 | "musl/src/string/explicit_bzero.c", | 1022 | "musl/src/string/explicit_bzero.c", |
| 1037 | "musl/src/string/index.c", | 1023 | "musl/src/string/index.c", |