Ben Murdoch | b8a8cc1 | 2014-11-26 15:28:44 +0000 | [diff] [blame] | 1 | // The following is adapted from fdlibm (http://www.netlib.org/fdlibm). |
| 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 | // The original source code covered by the above license above has been |
| 13 | // modified significantly by Google Inc. |
| 14 | // Copyright 2014 the V8 project authors. All rights reserved. |
| 15 | |
Emily Bernier | d0a1eb7 | 2015-03-24 16:35:39 -0400 | [diff] [blame] | 16 | #include "src/third_party/fdlibm/fdlibm.h" |
Ben Murdoch | b8a8cc1 | 2014-11-26 15:28:44 +0000 | [diff] [blame] | 17 | |
Ben Murdoch | 4a90d5f | 2016-03-22 12:00:34 +0000 | [diff] [blame^] | 18 | #include <stdint.h> |
| 19 | #include <cmath> |
| 20 | #include <limits> |
| 21 | |
| 22 | #include "src/base/macros.h" |
| 23 | #include "src/double.h" |
Ben Murdoch | b8a8cc1 | 2014-11-26 15:28:44 +0000 | [diff] [blame] | 24 | |
| 25 | namespace v8 { |
| 26 | namespace fdlibm { |
| 27 | |
| 28 | #ifdef _MSC_VER |
| 29 | inline double scalbn(double x, int y) { return _scalb(x, y); } |
| 30 | #endif // _MSC_VER |
| 31 | |
Ben Murdoch | b8a8cc1 | 2014-11-26 15:28:44 +0000 | [diff] [blame] | 32 | |
| 33 | // Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi |
| 34 | static const int two_over_pi[] = { |
| 35 | 0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, 0x95993C, |
| 36 | 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, 0x424DD2, 0xE00649, |
| 37 | 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, 0xA73EE8, 0x8235F5, 0x2EBB44, |
| 38 | 0x84E99C, 0x7026B4, 0x5F7E41, 0x3991D6, 0x398353, 0x39F49C, 0x845F8B, |
| 39 | 0xBDF928, 0x3B1FF8, 0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, |
| 40 | 0x367ECF, 0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, |
| 41 | 0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, 0x560330, |
| 42 | 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, 0x91615E, 0xE61B08, |
| 43 | 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, 0x4D7327, 0x310606, 0x1556CA, |
| 44 | 0x73A8C9, 0x60E27B, 0xC08C6B}; |
| 45 | |
| 46 | static const double zero = 0.0; |
| 47 | static const double two24 = 1.6777216e+07; |
| 48 | static const double one = 1.0; |
| 49 | static const double twon24 = 5.9604644775390625e-08; |
| 50 | |
| 51 | static const double PIo2[] = { |
| 52 | 1.57079625129699707031e+00, // 0x3FF921FB, 0x40000000 |
| 53 | 7.54978941586159635335e-08, // 0x3E74442D, 0x00000000 |
| 54 | 5.39030252995776476554e-15, // 0x3CF84698, 0x80000000 |
| 55 | 3.28200341580791294123e-22, // 0x3B78CC51, 0x60000000 |
| 56 | 1.27065575308067607349e-29, // 0x39F01B83, 0x80000000 |
| 57 | 1.22933308981111328932e-36, // 0x387A2520, 0x40000000 |
| 58 | 2.73370053816464559624e-44, // 0x36E38222, 0x80000000 |
| 59 | 2.16741683877804819444e-51 // 0x3569F31D, 0x00000000 |
| 60 | }; |
| 61 | |
| 62 | |
Ben Murdoch | 4a90d5f | 2016-03-22 12:00:34 +0000 | [diff] [blame^] | 63 | INLINE(int __kernel_rem_pio2(double* x, double* y, int e0, int nx)) { |
Ben Murdoch | b8a8cc1 | 2014-11-26 15:28:44 +0000 | [diff] [blame] | 64 | static const int32_t jk = 3; |
| 65 | double fw; |
| 66 | int32_t jx = nx - 1; |
| 67 | int32_t jv = (e0 - 3) / 24; |
| 68 | if (jv < 0) jv = 0; |
| 69 | int32_t q0 = e0 - 24 * (jv + 1); |
| 70 | int32_t m = jx + jk; |
| 71 | |
Ben Murdoch | 4a90d5f | 2016-03-22 12:00:34 +0000 | [diff] [blame^] | 72 | double f[20]; |
Ben Murdoch | b8a8cc1 | 2014-11-26 15:28:44 +0000 | [diff] [blame] | 73 | for (int i = 0, j = jv - jx; i <= m; i++, j++) { |
| 74 | f[i] = (j < 0) ? zero : static_cast<double>(two_over_pi[j]); |
| 75 | } |
| 76 | |
Ben Murdoch | 4a90d5f | 2016-03-22 12:00:34 +0000 | [diff] [blame^] | 77 | double q[20]; |
Ben Murdoch | b8a8cc1 | 2014-11-26 15:28:44 +0000 | [diff] [blame] | 78 | for (int i = 0; i <= jk; i++) { |
| 79 | fw = 0.0; |
| 80 | for (int j = 0; j <= jx; j++) fw += x[j] * f[jx + i - j]; |
| 81 | q[i] = fw; |
| 82 | } |
| 83 | |
| 84 | int32_t jz = jk; |
| 85 | |
| 86 | recompute: |
| 87 | |
Ben Murdoch | 4a90d5f | 2016-03-22 12:00:34 +0000 | [diff] [blame^] | 88 | int32_t iq[20]; |
Ben Murdoch | b8a8cc1 | 2014-11-26 15:28:44 +0000 | [diff] [blame] | 89 | double z = q[jz]; |
| 90 | for (int i = 0, j = jz; j > 0; i++, j--) { |
| 91 | fw = static_cast<double>(static_cast<int32_t>(twon24 * z)); |
| 92 | iq[i] = static_cast<int32_t>(z - two24 * fw); |
| 93 | z = q[j - 1] + fw; |
| 94 | } |
| 95 | |
| 96 | z = scalbn(z, q0); |
| 97 | z -= 8.0 * std::floor(z * 0.125); |
| 98 | int32_t n = static_cast<int32_t>(z); |
| 99 | z -= static_cast<double>(n); |
| 100 | int32_t ih = 0; |
| 101 | if (q0 > 0) { |
| 102 | int32_t i = (iq[jz - 1] >> (24 - q0)); |
| 103 | n += i; |
| 104 | iq[jz - 1] -= i << (24 - q0); |
| 105 | ih = iq[jz - 1] >> (23 - q0); |
| 106 | } else if (q0 == 0) { |
| 107 | ih = iq[jz - 1] >> 23; |
| 108 | } else if (z >= 0.5) { |
| 109 | ih = 2; |
| 110 | } |
| 111 | |
| 112 | if (ih > 0) { |
| 113 | n += 1; |
| 114 | int32_t carry = 0; |
| 115 | for (int i = 0; i < jz; i++) { |
| 116 | int32_t j = iq[i]; |
| 117 | if (carry == 0) { |
| 118 | if (j != 0) { |
| 119 | carry = 1; |
| 120 | iq[i] = 0x1000000 - j; |
| 121 | } |
| 122 | } else { |
| 123 | iq[i] = 0xffffff - j; |
| 124 | } |
| 125 | } |
| 126 | if (q0 == 1) { |
| 127 | iq[jz - 1] &= 0x7fffff; |
| 128 | } else if (q0 == 2) { |
| 129 | iq[jz - 1] &= 0x3fffff; |
| 130 | } |
| 131 | if (ih == 2) { |
| 132 | z = one - z; |
| 133 | if (carry != 0) z -= scalbn(one, q0); |
| 134 | } |
| 135 | } |
| 136 | |
| 137 | if (z == zero) { |
| 138 | int32_t j = 0; |
| 139 | for (int i = jz - 1; i >= jk; i--) j |= iq[i]; |
| 140 | if (j == 0) { |
| 141 | int32_t k = 1; |
| 142 | while (iq[jk - k] == 0) k++; |
| 143 | for (int i = jz + 1; i <= jz + k; i++) { |
| 144 | f[jx + i] = static_cast<double>(two_over_pi[jv + i]); |
| 145 | for (j = 0, fw = 0.0; j <= jx; j++) fw += x[j] * f[jx + i - j]; |
| 146 | q[i] = fw; |
| 147 | } |
| 148 | jz += k; |
| 149 | goto recompute; |
| 150 | } |
| 151 | } |
| 152 | |
| 153 | if (z == 0.0) { |
| 154 | jz -= 1; |
| 155 | q0 -= 24; |
| 156 | while (iq[jz] == 0) { |
| 157 | jz--; |
| 158 | q0 -= 24; |
| 159 | } |
| 160 | } else { |
| 161 | z = scalbn(z, -q0); |
| 162 | if (z >= two24) { |
| 163 | fw = static_cast<double>(static_cast<int32_t>(twon24 * z)); |
| 164 | iq[jz] = static_cast<int32_t>(z - two24 * fw); |
| 165 | jz += 1; |
| 166 | q0 += 24; |
| 167 | iq[jz] = static_cast<int32_t>(fw); |
| 168 | } else { |
| 169 | iq[jz] = static_cast<int32_t>(z); |
| 170 | } |
| 171 | } |
| 172 | |
| 173 | fw = scalbn(one, q0); |
| 174 | for (int i = jz; i >= 0; i--) { |
| 175 | q[i] = fw * static_cast<double>(iq[i]); |
| 176 | fw *= twon24; |
| 177 | } |
| 178 | |
Ben Murdoch | 4a90d5f | 2016-03-22 12:00:34 +0000 | [diff] [blame^] | 179 | double fq[20]; |
Ben Murdoch | b8a8cc1 | 2014-11-26 15:28:44 +0000 | [diff] [blame] | 180 | for (int i = jz; i >= 0; i--) { |
| 181 | fw = 0.0; |
| 182 | for (int k = 0; k <= jk && k <= jz - i; k++) fw += PIo2[k] * q[i + k]; |
| 183 | fq[jz - i] = fw; |
| 184 | } |
| 185 | |
| 186 | fw = 0.0; |
| 187 | for (int i = jz; i >= 0; i--) fw += fq[i]; |
| 188 | y[0] = (ih == 0) ? fw : -fw; |
| 189 | fw = fq[0] - fw; |
| 190 | for (int i = 1; i <= jz; i++) fw += fq[i]; |
| 191 | y[1] = (ih == 0) ? fw : -fw; |
| 192 | return n & 7; |
| 193 | } |
| 194 | |
| 195 | |
| 196 | int rempio2(double x, double* y) { |
| 197 | int32_t hx = static_cast<int32_t>(internal::double_to_uint64(x) >> 32); |
| 198 | int32_t ix = hx & 0x7fffffff; |
| 199 | |
| 200 | if (ix >= 0x7ff00000) { |
Ben Murdoch | 4a90d5f | 2016-03-22 12:00:34 +0000 | [diff] [blame^] | 201 | *y = std::numeric_limits<double>::quiet_NaN(); |
Ben Murdoch | b8a8cc1 | 2014-11-26 15:28:44 +0000 | [diff] [blame] | 202 | return 0; |
| 203 | } |
| 204 | |
| 205 | int32_t e0 = (ix >> 20) - 1046; |
| 206 | uint64_t zi = internal::double_to_uint64(x) & 0xFFFFFFFFu; |
| 207 | zi |= static_cast<uint64_t>(ix - (e0 << 20)) << 32; |
| 208 | double z = internal::uint64_to_double(zi); |
| 209 | |
| 210 | double tx[3]; |
| 211 | for (int i = 0; i < 2; i++) { |
| 212 | tx[i] = static_cast<double>(static_cast<int32_t>(z)); |
| 213 | z = (z - tx[i]) * two24; |
| 214 | } |
| 215 | tx[2] = z; |
| 216 | |
| 217 | int nx = 3; |
| 218 | while (tx[nx - 1] == zero) nx--; |
| 219 | int n = __kernel_rem_pio2(tx, y, e0, nx); |
| 220 | if (hx < 0) { |
| 221 | y[0] = -y[0]; |
| 222 | y[1] = -y[1]; |
| 223 | return -n; |
| 224 | } |
| 225 | return n; |
| 226 | } |
Ben Murdoch | 4a90d5f | 2016-03-22 12:00:34 +0000 | [diff] [blame^] | 227 | } // namespace internal |
| 228 | } // namespace v8 |