src/v8/src/third_party/fdlibm/fdlibm.cc

   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/v8.h"
  17
  18 #include "src/double.h"
  19 #include "src/third_party/fdlibm/fdlibm.h"
  20
  21
  22 namespace v8 {
  23 namespace fdlibm {
  24
  25 #ifdef _MSC_VER
  26 inline double scalbn(double x, int y) { return _scalb(x, y); }
  27 #endif  // _MSC_VER
  28
  29 const double MathConstants::constants[] = {
  30     6.36619772367581382433e-01,   // invpio2   0
  31     1.57079632673412561417e+00,   // pio2_1    1
  32     6.07710050650619224932e-11,   // pio2_1t   2
  33     6.07710050630396597660e-11,   // pio2_2    3
  34     2.02226624879595063154e-21,   // pio2_2t   4
  35     2.02226624871116645580e-21,   // pio2_3    5
  36     8.47842766036889956997e-32,   // pio2_3t   6
  37     -1.66666666666666324348e-01,  // S1        7  coefficients for sin
  38     8.33333333332248946124e-03,   //           8
  39     -1.98412698298579493134e-04,  //           9
  40     2.75573137070700676789e-06,   //          10
  41     -2.50507602534068634195e-08,  //          11
  42     1.58969099521155010221e-10,   // S6       12
  43     4.16666666666666019037e-02,   // C1       13  coefficients for cos
  44     -1.38888888888741095749e-03,  //          14
  45     2.48015872894767294178e-05,   //          15
  46     -2.75573143513906633035e-07,  //          16
  47     2.08757232129817482790e-09,   //          17
  48     -1.13596475577881948265e-11,  // C6       18
  49     3.33333333333334091986e-01,   // T0       19  coefficients for tan
  50     1.33333333333201242699e-01,   //          20
  51     5.39682539762260521377e-02,   //          21
  52     2.18694882948595424599e-02,   //          22
  53     8.86323982359930005737e-03,   //          23
  54     3.59207910759131235356e-03,   //          24
  55     1.45620945432529025516e-03,   //          25
  56     5.88041240820264096874e-04,   //          26
  57     2.46463134818469906812e-04,   //          27
  58     7.81794442939557092300e-05,   //          28
  59     7.14072491382608190305e-05,   //          29
  60     -1.85586374855275456654e-05,  //          30
  61     2.59073051863633712884e-05,   // T12      31
  62     7.85398163397448278999e-01,   // pio4     32
  63     3.06161699786838301793e-17,   // pio4lo   33
  64     6.93147180369123816490e-01,   // ln2_hi   34
  65     1.90821492927058770002e-10,   // ln2_lo   35
  66     1.80143985094819840000e+16,   // 2^54     36
  67     6.666666666666666666e-01,     // 2/3      37
  68     6.666666666666735130e-01,     // LP1      38  coefficients for log1p
  69     3.999999999940941908e-01,     //          39
  70     2.857142874366239149e-01,     //          40
  71     2.222219843214978396e-01,     //          41
  72     1.818357216161805012e-01,     //          42
  73     1.531383769920937332e-01,     //          43
  74     1.479819860511658591e-01,     // LP7      44
  75     7.09782712893383973096e+02,   //          45  overflow threshold for expm1
  76     1.44269504088896338700e+00,   // 1/ln2    46
  77     -3.33333333333331316428e-02,  // Q1       47  coefficients for expm1
  78     1.58730158725481460165e-03,   //          48
  79     -7.93650757867487942473e-05,  //          49
  80     4.00821782732936239552e-06,   //          50
  81     -2.01099218183624371326e-07,  // Q5       51
  82     710.4758600739439             //          52  overflow threshold sinh, cosh
  83 };
  84
  85
  86 // Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
  87 static const int two_over_pi[] = {
  88     0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, 0x95993C,
  89     0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, 0x424DD2, 0xE00649,
  90     0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, 0xA73EE8, 0x8235F5, 0x2EBB44,
  91     0x84E99C, 0x7026B4, 0x5F7E41, 0x3991D6, 0x398353, 0x39F49C, 0x845F8B,
  92     0xBDF928, 0x3B1FF8, 0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D,
  93     0x367ECF, 0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5,
  94     0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, 0x560330,
  95     0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, 0x91615E, 0xE61B08,
  96     0x659985, 0x5F14A0, 0x68408D, 0xFFD880, 0x4D7327, 0x310606, 0x1556CA,
  97     0x73A8C9, 0x60E27B, 0xC08C6B};
  98
  99 static const double zero = 0.0;
 100 static const double two24 = 1.6777216e+07;
 101 static const double one = 1.0;
 102 static const double twon24 = 5.9604644775390625e-08;
 103
 104 static const double PIo2[] = {
 105     1.57079625129699707031e+00,  // 0x3FF921FB, 0x40000000
 106     7.54978941586159635335e-08,  // 0x3E74442D, 0x00000000
 107     5.39030252995776476554e-15,  // 0x3CF84698, 0x80000000
 108     3.28200341580791294123e-22,  // 0x3B78CC51, 0x60000000
 109     1.27065575308067607349e-29,  // 0x39F01B83, 0x80000000
 110     1.22933308981111328932e-36,  // 0x387A2520, 0x40000000
 111     2.73370053816464559624e-44,  // 0x36E38222, 0x80000000
 112     2.16741683877804819444e-51   // 0x3569F31D, 0x00000000
 113 };
 114
 115
 116 int __kernel_rem_pio2(double* x, double* y, int e0, int nx) {
 117   static const int32_t jk = 3;
 118   double fw;
 119   int32_t jx = nx - 1;
 120   int32_t jv = (e0 - 3) / 24;
 121   if (jv < 0) jv = 0;
 122   int32_t q0 = e0 - 24 * (jv + 1);
 123   int32_t m = jx + jk;
 124
 125   double f[10];
 126   for (int i = 0, j = jv - jx; i <= m; i++, j++) {
 127     f[i] = (j < 0) ? zero : static_cast<double>(two_over_pi[j]);
 128   }
 129
 130   double q[10];
 131   for (int i = 0; i <= jk; i++) {
 132     fw = 0.0;
 133     for (int j = 0; j <= jx; j++) fw += x[j] * f[jx + i - j];
 134     q[i] = fw;
 135   }
 136
 137   int32_t jz = jk;
 138
 139 recompute:
 140
 141   int32_t iq[10];
 142   double z = q[jz];
 143   for (int i = 0, j = jz; j > 0; i++, j--) {
 144     fw = static_cast<double>(static_cast<int32_t>(twon24 * z));
 145     iq[i] = static_cast<int32_t>(z - two24 * fw);
 146     z = q[j - 1] + fw;
 147   }
 148
 149   z = scalbn(z, q0);
 150   z -= 8.0 * std::floor(z * 0.125);
 151   int32_t n = static_cast<int32_t>(z);
 152   z -= static_cast<double>(n);
 153   int32_t ih = 0;
 154   if (q0 > 0) {
 155     int32_t i = (iq[jz - 1] >> (24 - q0));
 156     n += i;
 157     iq[jz - 1] -= i << (24 - q0);
 158     ih = iq[jz - 1] >> (23 - q0);
 159   } else if (q0 == 0) {
 160     ih = iq[jz - 1] >> 23;
 161   } else if (z >= 0.5) {
 162     ih = 2;
 163   }
 164
 165   if (ih > 0) {
 166     n += 1;
 167     int32_t carry = 0;
 168     for (int i = 0; i < jz; i++) {
 169       int32_t j = iq[i];
 170       if (carry == 0) {
 171         if (j != 0) {
 172           carry = 1;
 173           iq[i] = 0x1000000 - j;
 174         }
 175       } else {
 176         iq[i] = 0xffffff - j;
 177       }
 178     }
 179     if (q0 == 1) {
 180       iq[jz - 1] &= 0x7fffff;
 181     } else if (q0 == 2) {
 182       iq[jz - 1] &= 0x3fffff;
 183     }
 184     if (ih == 2) {
 185       z = one - z;
 186       if (carry != 0) z -= scalbn(one, q0);
 187     }
 188   }
 189
 190   if (z == zero) {
 191     int32_t j = 0;
 192     for (int i = jz - 1; i >= jk; i--) j |= iq[i];
 193     if (j == 0) {
 194       int32_t k = 1;
 195       while (iq[jk - k] == 0) k++;
 196       for (int i = jz + 1; i <= jz + k; i++) {
 197         f[jx + i] = static_cast<double>(two_over_pi[jv + i]);
 198         for (j = 0, fw = 0.0; j <= jx; j++) fw += x[j] * f[jx + i - j];
 199         q[i] = fw;
 200       }
 201       jz += k;
 202       goto recompute;
 203     }
 204   }
 205
 206   if (z == 0.0) {
 207     jz -= 1;
 208     q0 -= 24;
 209     while (iq[jz] == 0) {
 210       jz--;
 211       q0 -= 24;
 212     }
 213   } else {
 214     z = scalbn(z, -q0);
 215     if (z >= two24) {
 216       fw = static_cast<double>(static_cast<int32_t>(twon24 * z));
 217       iq[jz] = static_cast<int32_t>(z - two24 * fw);
 218       jz += 1;
 219       q0 += 24;
 220       iq[jz] = static_cast<int32_t>(fw);
 221     } else {
 222       iq[jz] = static_cast<int32_t>(z);
 223     }
 224   }
 225
 226   fw = scalbn(one, q0);
 227   for (int i = jz; i >= 0; i--) {
 228     q[i] = fw * static_cast<double>(iq[i]);
 229     fw *= twon24;
 230   }
 231
 232   double fq[10];
 233   for (int i = jz; i >= 0; i--) {
 234     fw = 0.0;
 235     for (int k = 0; k <= jk && k <= jz - i; k++) fw += PIo2[k] * q[i + k];
 236     fq[jz - i] = fw;
 237   }
 238
 239   fw = 0.0;
 240   for (int i = jz; i >= 0; i--) fw += fq[i];
 241   y[0] = (ih == 0) ? fw : -fw;
 242   fw = fq[0] - fw;
 243   for (int i = 1; i <= jz; i++) fw += fq[i];
 244   y[1] = (ih == 0) ? fw : -fw;
 245   return n & 7;
 246 }
 247
 248
 249 int rempio2(double x, double* y) {
 250   int32_t hx = static_cast<int32_t>(internal::double_to_uint64(x) >> 32);
 251   int32_t ix = hx & 0x7fffffff;
 252
 253   if (ix >= 0x7ff00000) {
 254     *y = base::OS::nan_value();
 255     return 0;
 256   }
 257
 258   int32_t e0 = (ix >> 20) - 1046;
 259   uint64_t zi = internal::double_to_uint64(x) & 0xFFFFFFFFu;
 260   zi |= static_cast<uint64_t>(ix - (e0 << 20)) << 32;
 261   double z = internal::uint64_to_double(zi);
 262
 263   double tx[3];
 264   for (int i = 0; i < 2; i++) {
 265     tx[i] = static_cast<double>(static_cast<int32_t>(z));
 266     z = (z - tx[i]) * two24;
 267   }
 268   tx[2] = z;
 269
 270   int nx = 3;
 271   while (tx[nx - 1] == zero) nx--;
 272   int n = __kernel_rem_pio2(tx, y, e0, nx);
 273   if (hx < 0) {
 274     y[0] = -y[0];
 275     y[1] = -y[1];
 276     return -n;
 277   }
 278   return n;
 279 }
 280 }
 281 }  // namespace v8::internal