strictmath.java

来自「纯java操作系统jnode,安装简单和操作简单的个人使用的Java操作系统」· Java 代码 · 共 1,754 行 · 第 1/4 页

		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.

	/**