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 16 #include "src/third_party/fdlibm/fdlibm.h" 17 18 #include <stdint.h> 19 #include <cmath> 20 #include <limits> 21 22 #include "src/base/macros.h" 23 #include "src/double.h" 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 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 63 INLINE(int __kernel_rem_pio2(double* x, double* y, int e0, int nx)) { 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 72 double f[20]; 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 77 double q[20]; 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 88 int32_t iq[20]; 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 179 double fq[20]; 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) { 201 *y = std::numeric_limits<double>::quiet_NaN(); 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 } 227 } // namespace internal 228 } // namespace v8 229
