strictmath.java

this gcc-g++-3.3.1.tar.gz is a source file of gcc, you can learn more about gcc through this codes f

    T4 = 8.8632398235993e-3, // Long bits 0x3f8226e3e96e8493L.
    T5 = 3.5920791075913124e-3, // Long bits 0x3f6d6d22c9560328L.
    T6 = 1.4562094543252903e-3, // Long bits 0x3f57dbc8fee08315L.
    T7 = 5.880412408202641e-4, // Long bits 0x3f4344d8f2f26501L.
    T8 = 2.464631348184699e-4, // Long bits 0x3f3026f71a8d1068L.
    T9 = 7.817944429395571e-5, // Long bits 0x3f147e88a03792a6L.
    T10 = 7.140724913826082e-5, // Long bits 0x3f12b80f32f0a7e9L.
    T11 = -1.8558637485527546e-5, // Long bits 0xbef375cbdb605373L.
    T12 = 2.590730518636337e-5; // Long bits 0x3efb2a7074bf7ad4L.

  /**
   * Coefficients for computing {@link #asin(double)} and
   * {@link #acos(double)}.
   */
  private static final double
    PS0 = 0.16666666666666666, // Long bits 0x3fc5555555555555L.
    PS1 = -0.3255658186224009, // Long bits 0xbfd4d61203eb6f7dL.
    PS2 = 0.20121253213486293, // Long bits 0x3fc9c1550e884455L.
    PS3 = -0.04005553450067941, // Long bits 0xbfa48228b5688f3bL.
    PS4 = 7.915349942898145e-4, // Long bits 0x3f49efe07501b288L.
    PS5 = 3.479331075960212e-5, // Long bits 0x3f023de10dfdf709L.
    QS1 = -2.403394911734414, // Long bits 0xc0033a271c8a2d4bL.
    QS2 = 2.0209457602335057, // Long bits 0x40002ae59c598ac8L.
    QS3 = -0.6882839716054533, // Long bits 0xbfe6066c1b8d0159L.
    QS4 = 0.07703815055590194; // Long bits 0x3fb3b8c5b12e9282L.

  /**
   * Coefficients for computing {@link #atan(double)}.
   */
  private static final double
    ATAN_0_5H = 0.4636476090008061, // Long bits 0x3fddac670561bb4fL.
    ATAN_0_5L = 2.2698777452961687e-17, // Long bits 0x3c7a2b7f222f65e2L.
    ATAN_1_5H = 0.982793723247329, // Long bits 0x3fef730bd281f69bL.
    ATAN_1_5L = 1.3903311031230998e-17, // Long bits 0x3c7007887af0cbbdL.
    AT0 = 0.3333333333333293, // Long bits 0x3fd555555555550dL.
    AT1 = -0.19999999999876483, // Long bits 0xbfc999999998ebc4L.
    AT2 = 0.14285714272503466, // Long bits 0x3fc24924920083ffL.
    AT3 = -0.11111110405462356, // Long bits 0xbfbc71c6fe231671L.
    AT4 = 0.09090887133436507, // Long bits 0x3fb745cdc54c206eL.
    AT5 = -0.0769187620504483, // Long bits 0xbfb3b0f2af749a6dL.
    AT6 = 0.06661073137387531, // Long bits 0x3fb10d66a0d03d51L.
    AT7 = -0.058335701337905735, // Long bits 0xbfadde2d52defd9aL.
    AT8 = 0.049768779946159324, // Long bits 0x3fa97b4b24760debL.
    AT9 = -0.036531572744216916, // Long bits 0xbfa2b4442c6a6c2fL.
    AT10 = 0.016285820115365782; // Long bits 0x3f90ad3ae322da11L.

  /**
   * Helper function for reducing an angle to a multiple of pi/2 within
   * [-pi/4, pi/4].
   *
   * @param x the angle; not infinity or NaN, and outside pi/4
   * @param y an array of 2 doubles modified to hold the remander x % pi/2
   * @return the quadrant of the result, mod 4: 0: [-pi/4, pi/4],
   *         1: [pi/4, 3*pi/4], 2: [3*pi/4, 5*pi/4], 3: [-3*pi/4, -pi/4]
   */
  private static int remPiOver2(double x, double[] y)
  {
    boolean negative = x < 0;
    x = abs(x);
    double z;
    int n;
    if (Configuration.DEBUG && (x <= PI / 4 || x != x
                                || x == Double.POSITIVE_INFINITY))
      throw new InternalError("Assertion failure");
    if (x < 3 * PI / 4) // If |x| is small.
      {
        z = x - PIO2_1;
        if ((float) x != (float) (PI / 2)) // 33+53 bit pi is good enough.
          {
            y[0] = z - PIO2_1L;
            y[1] = z - y[0] - PIO2_1L;
          }
        else // Near pi/2, use 33+33+53 bit pi.
          {
            z -= PIO2_2;
            y[0] = z - PIO2_2L;
            y[1] = z - y[0] - PIO2_2L;
          }
        n = 1;
      }
    else if (x <= TWO_20 * PI / 2) // Medium size.
      {
        n = (int) (2 / PI * x + 0.5);
        z = x - n * PIO2_1;
        double w = n * PIO2_1L; // First round good to 85 bits.
        y[0] = z - w;
        if (n >= 32 || (float) x == (float) (w))
          {
            if (x / y[0] >= TWO_16) // Second iteration, good to 118 bits.
              {
                double t = z;
                w = n * PIO2_2;
                z = t - w;
                w = n * PIO2_2L - (t - z - w);
                y[0] = z - w;
                if (x / y[0] >= TWO_49) // Third iteration, 151 bits accuracy.
                  {
                    t = z;
                    w = n * PIO2_3;
                    z = t - w;
                    w = n * PIO2_3L - (t - z - w);
                    y[0] = z - w;
                  }
              }
          }
        y[1] = z - y[0] - w;
      }
    else
      {
        // All other (large) arguments.
        int e0 = (int) (Double.doubleToLongBits(x) >> 52) - 1046;
        z = scale(x, -e0); // e0 = ilogb(z) - 23.
        double[] tx = new double[3];
        for (int i = 0; i < 2; i++)
          {
            tx[i] = (int) z;
            z = (z - tx[i]) * TWO_24;
          }
        tx[2] = z;
        int nx = 2;
        while (tx[nx] == 0)
          nx--;
        n = remPiOver2(tx, y, e0, nx);
      }
    if (negative)
      {
        y[0] = -y[0];
        y[1] = -y[1];
        return -n;
      }
    return n;
  }

  /**
   * Helper function for reducing an angle to a multiple of pi/2 within
   * [-pi/4, pi/4].
   *
   * @param x the positive angle, broken into 24-bit chunks
   * @param y an array of 2 doubles modified to hold the remander x % pi/2
   * @param e0 the exponent of x[0]
   * @param nx the last index used in x
   * @return the quadrant of the result, mod 4: 0: [-pi/4, pi/4],
   *         1: [pi/4, 3*pi/4], 2: [3*pi/4, 5*pi/4], 3: [-3*pi/4, -pi/4]
   */
  private static int remPiOver2(double[] x, double[] y, int e0, int nx)
  {
    int i;
    int ih;
    int n;
    double fw;
    double z;
    int[] iq = new int[20];
    double[] f = new double[20];
    double[] q = new double[20];
    boolean recompute = false;

    // Initialize jk, jz, jv, q0; note that 3>q0.
    int jk = 4;
    int jz = jk;
    int jv = max((e0 - 3) / 24, 0);
    int q0 = e0 - 24 * (jv + 1);

    // Set up f[0] to f[nx+jk] where f[nx+jk] = TWO_OVER_PI[jv+jk].
    int j = jv - nx;
    int m = nx + jk;
    for (i = 0; i <= m; i++, j++)
      f[i] = (j < 0) ? 0 : TWO_OVER_PI[j];

    // Compute q[0],q[1],...q[jk].
    for (i = 0; i <= jk; i++)
      {
        for (j = 0, fw = 0; j <= nx; j++)
          fw += x[j] * f[nx + i - j];
        q[i] = fw;
      }

    do
      {
        // Distill q[] into iq[] reversingly.
        for (i = 0, j = jz, z = q[jz]; j > 0; i++, j--)
          {
            fw = (int) (1 / TWO_24 * z);
            iq[i] = (int) (z - TWO_24 * fw);
            z = q[j - 1] + fw;
          }

        // Compute n.
        z = scale(z, q0);
        z -= 8 * floor(z * 0.125); // Trim off integer >= 8.
        n = (int) z;
        z -= n;
        ih = 0;
        if (q0 > 0) // Need iq[jz-1] to determine n.
          {
            i = iq[jz - 1] >> (24 - q0);
            n += i;
            iq[jz - 1] -= i << (24 - q0);
            ih = iq[jz - 1] >> (23 - q0);
          }
        else if (q0 == 0)
          ih = iq[jz - 1] >> 23;
        else if (z >= 0.5)
          ih = 2;

        if (ih > 0) // If q > 0.5.
          {
            n += 1;
            int carry = 0;
            for (i = 0; i < jz; i++) // Compute 1-q.
              {
                j = iq[i];
                if (carry == 0)
                  {
                    if (j != 0)
                      {
                        carry = 1;
                        iq[i] = 0x1000000 - j;
                      }
                  }
                else
                  iq[i] = 0xffffff - j;
              }
            switch (q0)
              {
              case 1: // Rare case: chance is 1 in 12 for non-default.
                iq[jz - 1] &= 0x7fffff;
                break;
              case 2:
                iq[jz - 1] &= 0x3fffff;
              }
            if (ih == 2)
              {
                z = 1 - z;
                if (carry != 0)
                  z -= scale(1, q0);
              }
          }

        // Check if recomputation is needed.
        if (z == 0)
          {
            j = 0;
            for (i = jz - 1; i >= jk; i--)
              j |= iq[i];
            if (j == 0) // Need recomputation.
              {
                int k;
                for (k = 1; iq[jk - k] == 0; k++); // k = no. of terms needed.

                for (i = jz + 1; i <= jz + k; i++) // Add q[jz+1] to q[jz+k].
                  {
                    f[nx + i] = TWO_OVER_PI[jv + i];
                    for (j = 0, fw = 0; j <= nx; j++)
                      fw += x[j] * f[nx + i - j];
                    q[i] = fw;
                  }
                jz += k;
                recompute = true;
              }
          }
      }
    while (recompute);

    // Chop off zero terms.
    if (z == 0)
      {
        jz--;
        q0 -= 24;
        while (iq[jz] == 0)
          {
            jz--;
            q0 -= 24;
          }
      }
    else // Break z into 24-bit if necessary.
      {
        z = scale(z, -q0);
        if (z >= TWO_24)
          {
            fw = (int) (1 / TWO_24 * z);
            iq[jz] = (int) (z - TWO_24 * fw);
            jz++;
            q0 += 24;
            iq[jz] = (int) fw;
          }
        else
          iq[jz] = (int) z;
      }

    // Convert integer "bit" chunk to floating-point value.
    fw = scale(1, q0);
    for (i = jz; i >= 0; i--)
      {
        q[i] = fw * iq[i];
        fw *= 1 / TWO_24;
      }

    // Compute PI_OVER_TWO[0,...,jk]*q[jz,...,0].
    double[] fq = new double[20];
    for (i = jz; i >= 0; i--)
      {
        fw = 0;
        for (int k = 0; k <= jk && k <= jz - i; k++)
          fw += PI_OVER_TWO[k] * q[i + k];
        fq[jz - i] = fw;
      }

    // Compress fq[] into y[].
    fw = 0;
    for (i = jz; i >= 0; i--)
      fw += fq[i];
    y[0] = (ih == 0) ? fw : -fw;
    fw = fq[0] - fw;
    for (i = 1; i <= jz; i++)
      fw += fq[i];
    y[1] = (ih == 0) ? fw : -fw;
    return n;
  }

  /**
   * Helper method for scaling a double by a power of 2.
   *
   * @param x the double
   * @param n the scale; |n| < 2048
   * @return x * 2**n
   */
  private static double scale(double x, int n)
  {
    if (Configuration.DEBUG && abs(n) >= 2048)
      throw new InternalError("Assertion failure");
    if (x == 0 || x == Double.NEGATIVE_INFINITY
        || ! (x < Double.POSITIVE_INFINITY) || n == 0)
      return x;
    long bits = Double.doubleToLongBits(x);
    int exp = (int) (bits >> 52) & 0x7ff;
    if (exp == 0) // Subnormal x.
      {
        x *= TWO_54;
        exp = ((int) (Double.doubleToLongBits(x) >> 52) & 0x7ff) - 54;
      }
    exp += n;
    if (exp > 0x7fe) // Overflow.
      return Double.POSITIVE_INFINITY * x;
    if (exp > 0) // Normal.
      return Double.longBitsToDouble((bits & 0x800fffffffffffffL)
                                     | ((long) exp << 52));
    if (exp <= -54)
      return 0 * x; // Underflow.
    exp += 54; // Subnormal result.
    x = Double.longBitsToDouble((bits & 0x800fffffffffffffL)
                                | ((long) exp << 52));
    return x * (1 / TWO_54);
  }

  /**
   * Helper trig function; computes sin in range [-pi/4, pi/4].
   *
   * @param x angle within about pi/4
   * @param y tail of x, created by remPiOver2
   * @return sin(x+y)
   */
  private static double sin(double x, double y)
  {
    if (Configuration.DEBUG && abs(x + y) > 0.7854)
      throw new InternalError("Assertion failure");
    if (abs(x) < 1 / TWO_27)
      return x;  // If |x| ~< 2**-27, already know answer.

    double z = x * x;
    double v = z * x;
    double r = S2 + z * (S3 + z * (S4 + z * (S5 + z * S6)));
    if (y == 0)
      return x + v * (S1 + z * r);
    return x - ((z * (0.5 * y - v * r) - y) - v * S1);
  }

  /**
   * Helper trig function; computes cos in range [-pi/4, pi/4].
   *
   * @param x angle within about pi/4
   * @param y tail of x, created by remPiOver2
   * @return cos(x+y)
   */
  private static double cos(double x, double y)
  {
    if (Configuration.DEBUG && abs(x + y) > 0.7854)
      throw new InternalError("Assertion failure");
    x = abs(x);
    if (x < 1 / TWO_27)
      return 1;  // If |x| ~< 2**-27, already know answer.

    double z = x * x;
    double r = z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * C6)))));

    if (x < 0.3)
      return 1 - (0.5 * z - (z * r - x * y));

    double qx = (x > 0.78125) ? 0.28125 : (x * 0.25);
    return 1 - qx - ((0.5 * z - qx) - (z * r - x * y));
  }

  /**
   * Helper trig function; computes tan in range [-pi/4, pi/4].
   *
   * @param x angle within about pi/4
   * @param y tail of x, created by remPiOver2
   * @param invert true iff -1/tan should be returned instead
   * @return tan(x+y)
   */
  private static double tan(double x, double y, boolean invert)
  {
    // PI/2 is irrational, so no double is a perfect multiple of it.
    if (Configuration.DEBUG && (abs(x + y) > 0.7854 || (x == 0 && invert)))
      throw new InternalError("Assertion failure");
    boolean negative = x < 0;
    if (negative)
      {
        x = -x;
        y = -y;
      }
    if (x < 1 / TWO_28) // If |x| ~< 2**-28, already know answer.
      return (negative ? -1 : 1) * (invert ? -1 / x : x);

    double z;
    double w;
    boolean large = x >= 0.6744;
    if (large)
      {
        z = PI / 4 - x;
        w = PI_L / 4 - y;
        x = z + w;
        y = 0;
      }
    z = x * x;
    w = z * z;
    // Break x**5*(T1+x**2*T2+...) into
    //   x**5(T1+x**4*T3+...+x**20*T11)
    // + x**5(x**2*(T2+x**4*T4+...+x**22*T12)).
    double r = T1 + w * (T3 + w * (T5 + w * (T7 + w * (T9 + w * T11))));
    double v = z * (T2 + w * (T4 + w * (T6 + w * (T8 + w * (T10 + w * T12)))));
    double s = z * x;
    r = y + z * (s * (r + v) + y);
    r += T0 * s;
    w = x + r;
    if (large)
      {
        v = invert ? -1 : 1;
        return (negative ? -1 : 1) * (v - 2 * (x - (w * w / (w + v) - r)));
      }
    if (! invert)
      return w;

    // Compute -1.0/(x+r) accurately.
    z = (float) w;
    v = r - (z - x);
    double a = -1 / w;
    double t = (float) a;
    return t + a * (1 + t * z + t * v);
  }
}