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

      {
        if (y < 0)
          ax = 1 / ax;
        if (x < 0)
          {
            if (x == -1 && yisint == 0)
              ax = Double.NaN;
            else if (yisint == 1)
              ax = -ax;
          }
        return ax;
      }
    if (x < 0 && yisint == 0)
      return Double.NaN;

    // Now we can start!
    double t;
    double t1;
    double t2;
    double u;
    double v;
    double w;
    if (ay > TWO_31)
      {
        if (ay > TWO_64) // Automatic over/underflow.
          return ((ax < 1) ? y < 0 : y > 0) ? Double.POSITIVE_INFINITY : 0;
        // Over/underflow if x is not close to one.
        if (ax < 0.9999995231628418)
          return y < 0 ? Double.POSITIVE_INFINITY : 0;
        if (ax >= 1.0000009536743164)
          return y > 0 ? Double.POSITIVE_INFINITY : 0;
        // Now |1-x| is <= 2**-20, sufficient to compute
        // log(x) by x-x^2/2+x^3/3-x^4/4.
        t = x - 1;
        w = t * t * (0.5 - t * (1 / 3.0 - t * 0.25));
        u = INV_LN2_H * t;
        v = t * INV_LN2_L - w * INV_LN2;
        t1 = (float) (u + v);
        t2 = v - (t1 - u);
      }
    else
    {
      long bits = Double.doubleToLongBits(ax);
      int exp = (int) (bits >> 52);
      if (exp == 0) // Subnormal x.
        {
          ax *= TWO_54;
          bits = Double.doubleToLongBits(ax);
          exp = (int) (bits >> 52) - 54;
        }
      exp -= 1023; // Unbias exponent.
      ax = Double.longBitsToDouble((bits & 0x000fffffffffffffL)
                                   | 0x3ff0000000000000L);
      boolean k;
      if (ax < SQRT_1_5)  // |x|<sqrt(3/2).
        k = false;
      else if (ax < SQRT_3) // |x|<sqrt(3).
        k = true;
      else
        {
          k = false;
          ax *= 0.5;
          exp++;
        }

      // Compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5).
      u = ax - (k ? 1.5 : 1);
      v = 1 / (ax + (k ? 1.5 : 1));
      double s = u * v;
      double s_h = (float) s;
      double t_h = (float) (ax + (k ? 1.5 : 1));
      double t_l = ax - (t_h - (k ? 1.5 : 1));
      double s_l = v * ((u - s_h * t_h) - s_h * t_l);
      // Compute log(ax).
      double s2 = s * s;
      double r = s_l * (s_h + s) + s2 * s2
        * (L1 + s2 * (L2 + s2 * (L3 + s2 * (L4 + s2 * (L5 + s2 * L6)))));
      s2 = s_h * s_h;
      t_h = (float) (3.0 + s2 + r);
      t_l = r - (t_h - 3.0 - s2);
      // u+v = s*(1+...).
      u = s_h * t_h;
      v = s_l * t_h + t_l * s;
      // 2/(3log2)*(s+...).
      double p_h = (float) (u + v);
      double p_l = v - (p_h - u);
      double z_h = CP_H * p_h;
      double z_l = CP_L * p_h + p_l * CP + (k ? DP_L : 0);
      // log2(ax) = (s+..)*2/(3*log2) = exp + dp_h + z_h + z_l.
      t = exp;
      t1 = (float) (z_h + z_l + (k ? DP_H : 0) + t);
      t2 = z_l - (t1 - t - (k ? DP_H : 0) - z_h);
    }

    // Split up y into y1+y2 and compute (y1+y2)*(t1+t2).
    boolean negative = x < 0 && yisint == 1;
    double y1 = (float) y;
    double p_l = (y - y1) * t1 + y * t2;
    double p_h = y1 * t1;
    double z = p_l + p_h;
    if (z >= 1024) // Detect overflow.
      {
        if (z > 1024 || p_l + OVT > z - p_h)
          return negative ? Double.NEGATIVE_INFINITY
            : Double.POSITIVE_INFINITY;
      }
    else if (z <= -1075) // Detect underflow.
      {
        if (z < -1075 || p_l <= z - p_h)
          return negative ? -0.0 : 0;
      }

    // Compute 2**(p_h+p_l).
    int n = round((float) z);
    p_h -= n;
    t = (float) (p_l + p_h);
    u = t * LN2_H;
    v = (p_l - (t - p_h)) * LN2 + t * LN2_L;
    z = u + v;
    w = v - (z - u);
    t = z * z;
    t1 = z - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5))));
    double r = (z * t1) / (t1 - 2) - (w + z * w);
    z = scale(1 - (r - z), n);
    return negative ? -z : z;
  }

  /**
   * Get the IEEE 754 floating point remainder on two numbers. This is the
   * value of <code>x - y * <em>n</em></code>, where <em>n</em> is the closest
   * double to <code>x / y</code> (ties go to the even n); for a zero
   * remainder, the sign is that of <code>x</code>. If either argument is NaN,
   * the first argument is infinite, or the second argument is zero, the result
   * is NaN; if x is finite but y is infinte, the result is x.
   *
   * @param x the dividend (the top half)
   * @param y the divisor (the bottom half)
   * @return the IEEE 754-defined floating point remainder of x/y
   * @see #rint(double)
   */
  public static double IEEEremainder(double x, double y)
  {
    // Purge off exception values.
    if (x == Double.NEGATIVE_INFINITY || ! (x < Double.POSITIVE_INFINITY)
        || y == 0 || y != y)
      return Double.NaN;

    boolean negative = x < 0;
    x = abs(x);
    y = abs(y);
    if (x == y || x == 0)
      return 0 * x; // Get correct sign.

    // Achieve x < 2y, then take first shot at remainder.
    if (y < TWO_1023)
      x %= y + y;

    // Now adjust x to get correct precision.
    if (y < 4 / TWO_1023)
      {
        if (x + x > y)
          {
            x -= y;
            if (x + x >= y)
              x -= y;
          }
      }
    else
      {
        y *= 0.5;
        if (x > y)
          {
            x -= y;
            if (x >= y)
              x -= y;
          }
      }
    return negative ? -x : x;
  }

  /**
   * Take the nearest integer that is that is greater than or equal to the
   * argument. If the argument is NaN, infinite, or zero, the result is the
   * same; if the argument is between -1 and 0, the result is negative zero.
   * Note that <code>Math.ceil(x) == -Math.floor(-x)</code>.
   *
   * @param a the value to act upon
   * @return the nearest integer &gt;= <code>a</code>
   */
  public static double ceil(double a)
  {
    return -floor(-a);
  }

  /**
   * Take the nearest integer that is that is less than or equal to the
   * argument. If the argument is NaN, infinite, or zero, the result is the
   * same. Note that <code>Math.ceil(x) == -Math.floor(-x)</code>.
   *
   * @param a the value to act upon
   * @return the nearest integer &lt;= <code>a</code>
   */
  public static double floor(double a)
  {
    double x = abs(a);
    if (! (x < TWO_52) || (long) a == a)
      return a; // No fraction bits; includes NaN and infinity.
    if (x < 1)
      return a >= 0 ? 0 * a : -1; // Worry about signed zero.
    return a < 0 ? (long) a - 1.0 : (long) a; // Cast to long truncates.
  }

  /**
   * Take the nearest integer to the argument.  If it is exactly between
   * two integers, the even integer is taken. If the argument is NaN,
   * infinite, or zero, the result is the same.
   *
   * @param a the value to act upon
   * @return the nearest integer to <code>a</code>
   */
  public static double rint(double a)
  {
    double x = abs(a);
    if (! (x < TWO_52))
      return a; // No fraction bits; includes NaN and infinity.
    if (x <= 0.5)
      return 0 * a; // Worry about signed zero.
    if (x % 2 <= 0.5)
      return (long) a; // Catch round down to even.
    return (long) (a + (a < 0 ? -0.5 : 0.5)); // Cast to long truncates.
  }

  /**
   * Take the nearest integer to the argument.  This is equivalent to
   * <code>(int) Math.floor(f + 0.5f)</code>. If the argument is NaN, the
   * result is 0; otherwise if the argument is outside the range of int, the
   * result will be Integer.MIN_VALUE or Integer.MAX_VALUE, as appropriate.
   *
   * @param f the argument to round
   * @return the nearest integer to the argument
   * @see Integer#MIN_VALUE
   * @see Integer#MAX_VALUE
   */
  public static int round(float f)
  {
    return (int) floor(f + 0.5f);
  }

  /**
   * Take the nearest long to the argument.  This is equivalent to
   * <code>(long) Math.floor(d + 0.5)</code>. If the argument is NaN, the
   * result is 0; otherwise if the argument is outside the range of long, the
   * result will be Long.MIN_VALUE or Long.MAX_VALUE, as appropriate.
   *
   * @param d the argument to round
   * @return the nearest long to the argument
   * @see Long#MIN_VALUE
   * @see Long#MAX_VALUE
   */
  public static long round(double d)
  {
    return (long) floor(d + 0.5);
  }

  /**
   * Get a random number.  This behaves like Random.nextDouble(), seeded by
   * System.currentTimeMillis() when first called. In other words, the number
   * is from a pseudorandom sequence, and lies in the range [+0.0, 1.0).
   * This random sequence is only used by this method, and is threadsafe,
   * although you may want your own random number generator if it is shared
   * among threads.
   *
   * @return a random number
   * @see Random#nextDouble()
   * @see System#currentTimeMillis()
   */
  public static synchronized double random()
  {
    if (rand == null)
      rand = new Random();
    return rand.nextDouble();
  }

  /**
   * Convert from degrees to radians. The formula for this is
   * radians = degrees * (pi/180); however it is not always exact given the
   * limitations of floating point numbers.
   *
   * @param degrees an angle in degrees
   * @return the angle in radians
   */
  public static double toRadians(double degrees)
  {
    return degrees * (PI / 180);
  }

  /**
   * Convert from radians to degrees. The formula for this is
   * degrees = radians * (180/pi); however it is not always exact given the
   * limitations of floating point numbers.
   *
   * @param rads an angle in radians
   * @return the angle in degrees
   */
  public static double toDegrees(double rads)
  {
    return rads * (180 / PI);
  }

  /**
   * Constants for scaling and comparing doubles by powers of 2. The compiler
   * must automatically inline constructs like (1/TWO_54), so we don't list
   * negative powers of two here.
   */
  private static final double
    TWO_16 = 0x10000, // Long bits 0x40f0000000000000L.
    TWO_20 = 0x100000, // Long bits 0x4130000000000000L.
    TWO_24 = 0x1000000, // Long bits 0x4170000000000000L.
    TWO_27 = 0x8000000, // Long bits 0x41a0000000000000L.
    TWO_28 = 0x10000000, // Long bits 0x41b0000000000000L.
    TWO_29 = 0x20000000, // Long bits 0x41c0000000000000L.
    TWO_31 = 0x80000000L, // Long bits 0x41e0000000000000L.
    TWO_49 = 0x2000000000000L, // Long bits 0x4300000000000000L.
    TWO_52 = 0x10000000000000L, // Long bits 0x4330000000000000L.
    TWO_54 = 0x40000000000000L, // Long bits 0x4350000000000000L.
    TWO_57 = 0x200000000000000L, // Long bits 0x4380000000000000L.
    TWO_60 = 0x1000000000000000L, // Long bits 0x43b0000000000000L.
    TWO_64 = 1.8446744073709552e19, // Long bits 0x43f0000000000000L.
    TWO_66 = 7.378697629483821e19, // Long bits 0x4410000000000000L.
    TWO_1023 = 8.98846567431158e307; // Long bits 0x7fe0000000000000L.

  /**
   * Super precision for 2/pi in 24-bit chunks, for use in
   * {@link #remPiOver2()}.
   */
  private static final int TWO_OVER_PI[] = {
    0xa2f983, 0x6e4e44, 0x1529fc, 0x2757d1, 0xf534dd, 0xc0db62,
    0x95993c, 0x439041, 0xfe5163, 0xabdebb, 0xc561b7, 0x246e3a,
    0x424dd2, 0xe00649, 0x2eea09, 0xd1921c, 0xfe1deb, 0x1cb129,
    0xa73ee8, 0x8235f5, 0x2ebb44, 0x84e99c, 0x7026b4, 0x5f7e41,
    0x3991d6, 0x398353, 0x39f49c, 0x845f8b, 0xbdf928, 0x3b1ff8,
    0x97ffde, 0x05980f, 0xef2f11, 0x8b5a0a, 0x6d1f6d, 0x367ecf,
    0x27cb09, 0xb74f46, 0x3f669e, 0x5fea2d, 0x7527ba, 0xc7ebe5,
    0xf17b3d, 0x0739f7, 0x8a5292, 0xea6bfb, 0x5fb11f, 0x8d5d08,
    0x560330, 0x46fc7b, 0x6babf0, 0xcfbc20, 0x9af436, 0x1da9e3,
    0x91615e, 0xe61b08, 0x659985, 0x5f14a0, 0x68408d, 0xffd880,
    0x4d7327, 0x310606, 0x1556ca, 0x73a8c9, 0x60e27b, 0xc08c6b,
  };

  /**
   * Super precision for pi/2 in 24-bit chunks, for use in
   * {@link #remPiOver2()}.
   */
  private static final double PI_OVER_TWO[] = {
    1.570796251296997, // Long bits 0x3ff921fb40000000L.
    7.549789415861596e-8, // Long bits 0x3e74442d00000000L.
    5.390302529957765e-15, // Long bits 0x3cf8469880000000L.
    3.282003415807913e-22, // Long bits 0x3b78cc5160000000L.
    1.270655753080676e-29, // Long bits 0x39f01b8380000000L.
    1.2293330898111133e-36, // Long bits 0x387a252040000000L.
    2.7337005381646456e-44, // Long bits 0x36e3822280000000L.
    2.1674168387780482e-51, // Long bits 0x3569f31d00000000L.
  };

  /**
   * More constants related to pi, used in {@link #remPiOver2()} and
   * elsewhere.
   */
  private static final double
    PI_L = 1.2246467991473532e-16, // Long bits 0x3ca1a62633145c07L.
    PIO2_1 = 1.5707963267341256, // Long bits 0x3ff921fb54400000L.
    PIO2_1L = 6.077100506506192e-11, // Long bits 0x3dd0b4611a626331L.
    PIO2_2 = 6.077100506303966e-11, // Long bits 0x3dd0b4611a600000L.
    PIO2_2L = 2.0222662487959506e-21, // Long bits 0x3ba3198a2e037073L.
    PIO2_3 = 2.0222662487111665e-21, // Long bits 0x3ba3198a2e000000L.
    PIO2_3L = 8.4784276603689e-32; // Long bits 0x397b839a252049c1L.

  /**
   * Natural log and square root constants, for calculation of
   * {@link #exp(double)}, {@link #log(double)} and
   * {@link #power(double, double)}. CP is 2/(3*ln(2)).
   */
  private static final double
    SQRT_1_5 = 1.224744871391589, // Long bits 0x3ff3988e1409212eL.
    SQRT_2 = 1.4142135623730951, // Long bits 0x3ff6a09e667f3bcdL.
    SQRT_3 = 1.7320508075688772, // Long bits 0x3ffbb67ae8584caaL.
    EXP_LIMIT_H = 709.782712893384, // Long bits 0x40862e42fefa39efL.
    EXP_LIMIT_L = -745.1332191019411, // Long bits 0xc0874910d52d3051L.
    CP = 0.9617966939259756, // Long bits 0x3feec709dc3a03fdL.
    CP_H = 0.9617967009544373, // Long bits 0x3feec709e0000000L.
    CP_L = -7.028461650952758e-9, // Long bits 0xbe3e2fe0145b01f5L.
    LN2 = 0.6931471805599453, // Long bits 0x3fe62e42fefa39efL.
    LN2_H = 0.6931471803691238, // Long bits 0x3fe62e42fee00000L.
    LN2_L = 1.9082149292705877e-10, // Long bits 0x3dea39ef35793c76L.
    INV_LN2 = 1.4426950408889634, // Long bits 0x3ff71547652b82feL.
    INV_LN2_H = 1.4426950216293335, // Long bits 0x3ff7154760000000L.
    INV_LN2_L = 1.9259629911266175e-8; // Long bits 0x3e54ae0bf85ddf44L.

  /**
   * Constants for computing {@link #log(double)}.
   */
  private static final double
    LG1 = 0.6666666666666735, // Long bits 0x3fe5555555555593L.
    LG2 = 0.3999999999940942, // Long bits 0x3fd999999997fa04L.
    LG3 = 0.2857142874366239, // Long bits 0x3fd2492494229359L.
    LG4 = 0.22222198432149784, // Long bits 0x3fcc71c51d8e78afL.
    LG5 = 0.1818357216161805, // Long bits 0x3fc7466496cb03deL.
    LG6 = 0.15313837699209373, // Long bits 0x3fc39a09d078c69fL.
    LG7 = 0.14798198605116586; // Long bits 0x3fc2f112df3e5244L.

  /**
   * Constants for computing {@link #pow(double, double)}. L and P are
   * coefficients for series; OVT is -(1024-log2(ovfl+.5ulp)); and DP is ???.
   * The P coefficients also calculate {@link #exp(double)}.
   */
  private static final double
    L1 = 0.5999999999999946, // Long bits 0x3fe3333333333303L.
    L2 = 0.4285714285785502, // Long bits 0x3fdb6db6db6fabffL.
    L3 = 0.33333332981837743, // Long bits 0x3fd55555518f264dL.
    L4 = 0.272728123808534, // Long bits 0x3fd17460a91d4101L.
    L5 = 0.23066074577556175, // Long bits 0x3fcd864a93c9db65L.
    L6 = 0.20697501780033842, // Long bits 0x3fca7e284a454eefL.
    P1 = 0.16666666666666602, // Long bits 0x3fc555555555553eL.
    P2 = -2.7777777777015593e-3, // Long bits 0xbf66c16c16bebd93L.
    P3 = 6.613756321437934e-5, // Long bits 0x3f11566aaf25de2cL.
    P4 = -1.6533902205465252e-6, // Long bits 0xbebbbd41c5d26bf1L.
    P5 = 4.1381367970572385e-8, // Long bits 0x3e66376972bea4d0L.
    DP_H = 0.5849624872207642, // Long bits 0x3fe2b80340000000L.
    DP_L = 1.350039202129749e-8, // Long bits 0x3e4cfdeb43cfd006L.
    OVT = 8.008566259537294e-17; // Long bits 0x3c971547652b82feL.

  /**
   * Coefficients for computing {@link #sin(double)}.
   */
  private static final double
    S1 = -0.16666666666666632, // Long bits 0xbfc5555555555549L.
    S2 = 8.33333333332249e-3, // Long bits 0x3f8111111110f8a6L.
    S3 = -1.984126982985795e-4, // Long bits 0xbf2a01a019c161d5L.
    S4 = 2.7557313707070068e-6, // Long bits 0x3ec71de357b1fe7dL.
    S5 = -2.5050760253406863e-8, // Long bits 0xbe5ae5e68a2b9cebL.
    S6 = 1.58969099521155e-10; // Long bits 0x3de5d93a5acfd57cL.

  /**
   * Coefficients for computing {@link #cos(double)}.
   */
  private static final double
    C1 = 0.0416666666666666, // Long bits 0x3fa555555555554cL.
    C2 = -1.388888888887411e-3, // Long bits 0xbf56c16c16c15177L.
    C3 = 2.480158728947673e-5, // Long bits 0x3efa01a019cb1590L.
    C4 = -2.7557314351390663e-7, // Long bits 0xbe927e4f809c52adL.
    C5 = 2.087572321298175e-9, // Long bits 0x3e21ee9ebdb4b1c4L.
    C6 = -1.1359647557788195e-11; // Long bits 0xbda8fae9be8838d4L.

  /**
   * Coefficients for computing {@link #tan(double)}.
   */
  private static final double
    T0 = 0.3333333333333341, // Long bits 0x3fd5555555555563L.
    T1 = 0.13333333333320124, // Long bits 0x3fc111111110fe7aL.
    T2 = 0.05396825397622605, // Long bits 0x3faba1ba1bb341feL.
    T3 = 0.021869488294859542, // Long bits 0x3f9664f48406d637L.