strictmathtest.java
来自「纯java操作系统jnode,安装简单和操作简单的个人使用的Java操作系统」· Java 代码 · 共 1,689 行 · 第 1/4 页
JAVA
1,689 行
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);
}
}
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?