strictmath.java

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

	 * 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.
		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 (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 (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 (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) {
		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.
		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);
	}
}