这是一个有效/详细的实现:
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)
{
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;
}
从这里复制http://developer.classpath.org/doc/java/lang/StrictMath-source.html