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

                                                         * (PS4 + z * PS5)))));
        double q = 1 + z * (QS1 + z * (QS2 + z * (QS3 + z * QS4)));
        double s = sqrt(z);
        double w = p / q * s - PI_L / 2;
        return PI - 2 * (s + w);
      }
    double z = (1 - x) * 0.5; // x>0.5.
    double s = sqrt(z);
    double df = (float) s;
    double c = (z - df * df) / (s + df);
    double p = z * (PS0 + z * (PS1 + z * (PS2 + z * (PS3 + z
                                                     * (PS4 + z * PS5)))));
    double q = 1 + z * (QS1 + z * (QS2 + z * (QS3 + z * QS4)));
    double w = p / q * s + c;
    return 2 * (df + w);
  }

  /**
   * The trigonometric function <em>arcsin</em>. The range of angles returned
   * is -pi/2 to pi/2 radians (-90 to 90 degrees). If the argument is NaN, the
   * result is NaN; and the arctangent of 0 retains its sign.
   *
   * @param x the tan to turn back into an angle
   * @return arcsin(x)
   * @see #atan2(double, double)
   */
  public static double atan(double x)
  {
    double lo;
    double hi;
    boolean negative = x < 0;
    if (negative)
      x = -x;
    if (x >= TWO_66)
      return negative ? -PI / 2 : PI / 2;
    if (! (x >= 0.4375)) // |x|<7/16, or NaN.
      {
        if (! (x >= 1 / TWO_29)) // Small, or NaN.
          return negative ? -x : x;
        lo = hi = 0;
      }
    else if (x < 1.1875)
      {
        if (x < 0.6875) // 7/16<=|x|<11/16.
          {
            x = (2 * x - 1) / (2 + x);
            hi = ATAN_0_5H;
            lo = ATAN_0_5L;
          }
        else // 11/16<=|x|<19/16.
          {
            x = (x - 1) / (x + 1);
            hi = PI / 4;
            lo = PI_L / 4;
          }
      }
    else if (x < 2.4375) // 19/16<=|x|<39/16.
      {
        x = (x - 1.5) / (1 + 1.5 * x);
        hi = ATAN_1_5H;
        lo = ATAN_1_5L;
      }
    else // 39/16<=|x|<2**66.
      {
        x = -1 / x;
        hi = PI / 2;
        lo = PI_L / 2;
      }

    // Break sum from i=0 to 10 ATi*z**(i+1) into odd and even poly.
    double z = x * x;
    double w = z * z;
    double s1 = z * (AT0 + w * (AT2 + w * (AT4 + w * (AT6 + w
                                                      * (AT8 + w * AT10)))));
    double s2 = w * (AT1 + w * (AT3 + w * (AT5 + w * (AT7 + w * AT9))));
    if (hi == 0)
      return negative ? x * (s1 + s2) - x : x - x * (s1 + s2);
    z = hi - ((x * (s1 + s2) - lo) - x);
    return negative ? -z : z;
  }

  /**
   * A special version of the trigonometric function <em>arctan</em>, for
   * converting rectangular coordinates <em>(x, y)</em> to polar
   * <em>(r, theta)</em>. This computes the arctangent of x/y in the range
   * of -pi to pi radians (-180 to 180 degrees). Special cases:<ul>
   * <li>If either argument is NaN, the result is NaN.</li>
   * <li>If the first argument is positive zero and the second argument is
   * positive, or the first argument is positive and finite and the second
   * argument is positive infinity, then the result is positive zero.</li>
   * <li>If the first argument is negative zero and the second argument is
   * positive, or the first argument is negative and finite and the second
   * argument is positive infinity, then the result is negative zero.</li>
   * <li>If the first argument is positive zero and the second argument is
   * negative, or the first argument is positive and finite and the second
   * argument is negative infinity, then the result is the double value
   * closest to pi.</li>
   * <li>If the first argument is negative zero and the second argument is
   * negative, or the first argument is negative and finite and the second
   * argument is negative infinity, then the result is the double value
   * closest to -pi.</li>
   * <li>If the first argument is positive and the second argument is
   * positive zero or negative zero, or the first argument is positive
   * infinity and the second argument is finite, then the result is the
   * double value closest to pi/2.</li>
   * <li>If the first argument is negative and the second argument is
   * positive zero or negative zero, or the first argument is negative
   * infinity and the second argument is finite, then the result is the
   * double value closest to -pi/2.</li>
   * <li>If both arguments are positive infinity, then the result is the
   * double value closest to pi/4.</li>
   * <li>If the first argument is positive infinity and the second argument
   * is negative infinity, then the result is the double value closest to
   * 3*pi/4.</li>
   * <li>If the first argument is negative infinity and the second argument
   * is positive infinity, then the result is the double value closest to
   * -pi/4.</li>
   * <li>If both arguments are negative infinity, then the result is the
   * double value closest to -3*pi/4.</li>
   *
   * </ul><p>This returns theta, the angle of the point. To get r, albeit
   * slightly inaccurately, use sqrt(x*x+y*y).
   *
   * @param y the y position
   * @param x the x position
   * @return <em>theta</em> in the conversion of (x, y) to (r, theta)
   * @see #atan(double)
   */
  public static double atan2(double y, double x)
  {
    if (x != x || y != y)
      return Double.NaN;
    if (x == 1)
      return atan(y);
    if (x == Double.POSITIVE_INFINITY)
      {
        if (y == Double.POSITIVE_INFINITY)
          return PI / 4;
        if (y == Double.NEGATIVE_INFINITY)
          return -PI / 4;
        return 0 * y;
      }
    if (x == Double.NEGATIVE_INFINITY)
      {
        if (y == Double.POSITIVE_INFINITY)
          return 3 * PI / 4;
        if (y == Double.NEGATIVE_INFINITY)
          return -3 * PI / 4;
        return (1 / (0 * y) == Double.POSITIVE_INFINITY) ? PI : -PI;
      }
    if (y == 0)
      {
        if (1 / (0 * x) == Double.POSITIVE_INFINITY)
          return y;
        return (1 / y == Double.POSITIVE_INFINITY) ? PI : -PI;
      }
    if (y == Double.POSITIVE_INFINITY || y == Double.NEGATIVE_INFINITY
        || x == 0)
      return y < 0 ? -PI / 2 : PI / 2;

    double z = abs(y / x); // Safe to do y/x.
    if (z > TWO_60)
      z = PI / 2 + 0.5 * PI_L;
    else if (x < 0 && z < 1 / TWO_60)
      z = 0;
    else
      z = atan(z);
    if (x > 0)
      return y > 0 ? z : -z;
    return y > 0 ? PI - (z - PI_L) : z - PI_L - PI;
  }

  /**
   * Take <em>e</em><sup>a</sup>.  The opposite of <code>log()</code>. If the
   * argument is NaN, the result is NaN; if the argument is positive infinity,
   * the result is positive infinity; and if the argument is negative
   * infinity, the result is positive zero.
   *
   * @param x the number to raise to the power
   * @return the number raised to the power of <em>e</em>
   * @see #log(double)
   * @see #pow(double, double)
   */
  public static double exp(double x)
  {
    if (x != x)
      return x;
    if (x > EXP_LIMIT_H)
      return Double.POSITIVE_INFINITY;
    if (x < EXP_LIMIT_L)
      return 0;

    // Argument reduction.
    double hi;
    double lo;
    int k;
    double t = abs(x);
    if (t > 0.5 * LN2)
      {
        if (t < 1.5 * LN2)
          {
            hi = t - LN2_H;
            lo = LN2_L;
            k = 1;
          }
        else
          {
            k = (int) (INV_LN2 * t + 0.5);
            hi = t - k * LN2_H;
            lo = k * LN2_L;
          }
        if (x < 0)
          {
            hi = -hi;
            lo = -lo;
            k = -k;
          }
        x = hi - lo;
      }
    else if (t < 1 / TWO_28)
      return 1;
    else
      lo = hi = k = 0;

    // Now x is in primary range.
    t = x * x;
    double c = x - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5))));
    if (k == 0)
      return 1 - (x * c / (c - 2) - x);
    double y = 1 - (lo - x * c / (2 - c) - hi);
    return scale(y, k);
  }

  /**
   * Take ln(a) (the natural log).  The opposite of <code>exp()</code>. If the
   * argument is NaN or negative, the result is NaN; if the argument is
   * positive infinity, the result is positive infinity; and if the argument
   * is either zero, the result is negative infinity.
   *
   * <p>Note that the way to get log<sub>b</sub>(a) is to do this:
   * <code>ln(a) / ln(b)</code>.
   *
   * @param x the number to take the natural log of
   * @return the natural log of <code>a</code>
   * @see #exp(double)
   */
  public static double log(double x)
  {
    if (x == 0)
      return Double.NEGATIVE_INFINITY;
    if (x < 0)
      return Double.NaN;
    if (! (x < Double.POSITIVE_INFINITY))
      return x;

    // Normalize x.
    long bits = Double.doubleToLongBits(x);
    int exp = (int) (bits >> 52);
    if (exp == 0) // Subnormal x.
      {
        x *= TWO_54;
        bits = Double.doubleToLongBits(x);
        exp = (int) (bits >> 52) - 54;
      }
    exp -= 1023; // Unbias exponent.
    bits = (bits & 0x000fffffffffffffL) | 0x3ff0000000000000L;
    x = Double.longBitsToDouble(bits);
    if (x >= SQRT_2)
      {
        x *= 0.5;
        exp++;
      }
    x--;
    if (abs(x) < 1 / TWO_20)
      {
        if (x == 0)
          return exp * LN2_H + exp * LN2_L;
        double r = x * x * (0.5 - 1 / 3.0 * x);
        if (exp == 0)
          return x - r;
        return exp * LN2_H - ((r - exp * LN2_L) - x);
      }
    double s = x / (2 + x);
    double z = s * s;
    double w = z * z;
    double t1 = w * (LG2 + w * (LG4 + w * LG6));
    double t2 = z * (LG1 + w * (LG3 + w * (LG5 + w * LG7)));
    double r = t2 + t1;
    if (bits >= 0x3ff6174a00000000L && bits < 0x3ff6b85200000000L)
      {
        double h = 0.5 * x * x; // Need more accuracy for x near sqrt(2).
        if (exp == 0)
          return x - (h - s * (h + r));
        return exp * LN2_H - ((h - (s * (h + r) + exp * LN2_L)) - x);
      }
    if (exp == 0)
      return x - s * (x - r);
    return exp * LN2_H - ((s * (x - r) - exp * LN2_L) - x);
  }

  /**
   * Take a square root. If the argument is NaN or negative, the result is
   * NaN; if the argument is positive infinity, the result is positive
   * infinity; and if the result is either zero, the result is the same.
   *
   * <p>For other roots, use pow(x, 1/rootNumber).
   *
   * @param x the numeric argument
   * @return the square root of the argument
   * @see #pow(double, double)
   */
  public static double sqrt(double x)
  {
    if (x < 0)
      return Double.NaN;
    if (x == 0 || ! (x < Double.POSITIVE_INFINITY))
      return x;

    // Normalize x.
    long bits = Double.doubleToLongBits(x);
    int exp = (int) (bits >> 52);
    if (exp == 0) // Subnormal x.
      {
        x *= TWO_54;
        bits = Double.doubleToLongBits(x);
        exp = (int) (bits >> 52) - 54;
      }
    exp -= 1023; // Unbias exponent.
    bits = (bits & 0x000fffffffffffffL) | 0x0010000000000000L;
    if ((exp & 1) == 1) // Odd exp, double x to make it even.
      bits <<= 1;
    exp >>= 1;

    // Generate sqrt(x) bit by bit.
    bits <<= 1;
    long q = 0;
    long s = 0;
    long r = 0x0020000000000000L; // Move r right to left.
    while (r != 0)
      {
        long t = s + r;
        if (t <= bits)
          {
            s = t + r;
            bits -= t;
            q += r;
          }
        bits <<= 1;
        r >>= 1;
      }

    // Use floating add to round correctly.
    if (bits != 0)
      q += q & 1;
    return Double.longBitsToDouble((q >> 1) + ((exp + 1022L) << 52));
  }

  /**
   * Raise a number to a power. Special cases:<ul>
   * <li>If the second argument is positive or negative zero, then the result
   * is 1.0.</li>
   * <li>If the second argument is 1.0, then the result is the same as the
   * first argument.</li>
   * <li>If the second argument is NaN, then the result is NaN.</li>
   * <li>If the first argument is NaN and the second argument is nonzero,
   * then the result is NaN.</li>
   * <li>If the absolute value of the first argument is greater than 1 and
   * the second argument is positive infinity, or the absolute value of the
   * first argument is less than 1 and the second argument is negative
   * infinity, then the result is positive infinity.</li>
   * <li>If the absolute value of the first argument is greater than 1 and
   * the second argument is negative infinity, or the absolute value of the
   * first argument is less than 1 and the second argument is positive
   * infinity, then the result is positive zero.</li>
   * <li>If the absolute value of the first argument equals 1 and the second
   * argument is infinite, then the result is NaN.</li>
   * <li>If the first argument is positive zero and the second argument is
   * greater than zero, or the first argument is positive infinity and the
   * second argument is less than zero, then the result is positive zero.</li>
   * <li>If the first argument is positive zero and the second argument is
   * less than zero, or the first argument is positive infinity and the
   * second argument is greater than zero, then the result is positive
   * infinity.</li>
   * <li>If the first argument is negative zero and the second argument is
   * greater than zero but not a finite odd integer, or the first argument is
   * negative infinity and the second argument is less than zero but not a
   * finite odd integer, then the result is positive zero.</li>
   * <li>If the first argument is negative zero and the second argument is a
   * positive finite odd integer, or the first argument is negative infinity
   * and the second argument is a negative finite odd integer, then the result
   * is negative zero.</li>
   * <li>If the first argument is negative zero and the second argument is
   * less than zero but not a finite odd integer, or the first argument is
   * negative infinity and the second argument is greater than zero but not a
   * finite odd integer, then the result is positive infinity.</li>
   * <li>If the first argument is negative zero and the second argument is a
   * negative finite odd integer, or the first argument is negative infinity
   * and the second argument is a positive finite odd integer, then the result
   * is negative infinity.</li>
   * <li>If the first argument is less than zero and the second argument is a
   * finite even integer, then the result is equal to the result of raising
   * the absolute value of the first argument to the power of the second
   * argument.</li>
   * <li>If the first argument is less than zero and the second argument is a
   * finite odd integer, then the result is equal to the negative of the
   * result of raising the absolute value of the first argument to the power
   * of the second argument.</li>
   * <li>If the first argument is finite and less than zero and the second
   * argument is finite and not an integer, then the result is NaN.</li>
   * <li>If both arguments are integers, then the result is exactly equal to
   * the mathematical result of raising the first argument to the power of
   * the second argument if that result can in fact be represented exactly as
   * a double value.</li>
   *
   * </ul><p>(In the foregoing descriptions, a floating-point value is
   * considered to be an integer if and only if it is a fixed point of the
   * method {@link #ceil(double)} or, equivalently, a fixed point of the
   * method {@link #floor(double)}. A value is a fixed point of a one-argument
   * method if and only if the result of applying the method to the value is
   * equal to the value.)
   *
   * @param x the number to raise
   * @param y the power to raise it to
   * @return x<sup>y</sup>
   */
  public static double pow(double x, double y)
  {
    // Special cases first.
    if (y == 0)
      return 1;
    if (y == 1)
      return x;
    if (y == -1)
      return 1 / x;
    if (x != x || y != y)
      return Double.NaN;

    // When x < 0, yisint tells if y is not an integer (0), even(1),
    // or odd (2).
    int yisint = 0;
    if (x < 0 && floor(y) == y)
      yisint = (y % 2 == 0) ? 2 : 1;
    double ax = abs(x);
    double ay = abs(y);

    // More special cases, of y.
    if (ay == Double.POSITIVE_INFINITY)
      {
        if (ax == 1)
          return Double.NaN;
        if (ax > 1)
          return y > 0 ? y : 0;
        return y < 0 ? -y : 0;
      }
    if (y == 2)
      return x * x;
    if (y == 0.5)
      return sqrt(x);

    // More special cases, of x.
    if (x == 0 || ax == Double.POSITIVE_INFINITY || ax == 1)