Java中Math.pow()背后的算法是什么 [英] What's the algorithm behind Math.pow() in Java
问题描述
我刚开始使用Java，并且作为第一个项目，我正在编写一个程序，该程序查找给定数字的根(在这种情况下为立方根).目前，我正在尝试Newton-Ralphson来实现这一目标.这是代码-
I just started with Java and as a first project I am writing a program that finds roots(cube root in this case) of a given number. Presently I am trying out Newton-Ralphson to achieve this. Here is the code-
import java.util.Scanner;
import static java.lang.Math.abs;
public class newClass {
public static void main(String[] args) {
Scanner input = new Scanner(System.in);
System.out.println("Number whose cube root u wanna find:");
Double number = input.nextDouble();
Double epsilon = 0.0001;
Double ans = number/2.00;
while (abs((abs(number) - abs(Math.pow(ans,3))))>epsilon){
System.out.println("in loop");
ans = ans - ((Math.pow(ans,3) - number)/(3*Math.pow(ans,2)));
System.out.println(ans);
if ((number - ans)<=epsilon){
System.out.println(ans);
}
}
//System.out.println(Math.pow(number,1.0/3.0));
}
}
此功能最多只能使用11位数字，因为它太大了以至于IDE无法处理.但是，如果我只使用Math.pow(number,1.0/3.0)，它就可以处理更大的数字，并且可以立即进行计算.
this works only upto 11 digit numbers after that it gets too big for IDE to handle. But if I simply use Math.pow(number,1.0/3.0) it works for much bigger numbers and computes it in no time.
那么，Math.pow()使用什么算法可以立即给出答案?
我了解我的方法依赖于猜测，我猜math.pow()可能实际上是在计算答案，但是如何?
So, what is the algorithm that Math.pow() uses that gives an instant answer?
I understand that my method relies on guessing and I guess math.pow() may actually be calculating the answer, but how?
推荐答案
这是一个有趣的问题.如果查看Java Math类的源代码，您会发现它调用StrictMath.pow(double1，double2)，并且StrictMath的签名为public static native double pow(double a, double b);
This is a fun question. If you look into the the source code for Java's Math class, you will find that it calls StrictMath.pow(double1, double2), and StrictMath's signature is public static native double pow(double a, double b);
因此，最后，这是一个真正的本地调用，可能会因平台而异.但是，某处有一个实现，而且看起来并不容易.这是该函数的说明以及该函数本身的代码:
So, in the end, it is a truly native call that might differ depending on the platform. However, there is an implementation somewhere, and it isn't very easy to look at. Here is the description of the function and the code for the function itself:
注意
仔细研究数学，试图理解它可能不可避免地导致更多问题.但是，通过在 Java Math上搜索此Github函数源代码，并浏览一下数学摘要，您肯定可以更好地理解本机函数.快乐探索:)
Looking through the math, trying to understand it might inevitably lead to even more questions. But, by searching through this Github on Java Math Function Source Code and glancing out the mathematical summaries, you can definitely understand the native functions better. Happy Exploring :)
方法说明
Method: Let x = 2 * (1+f)
1. Compute and return log2(x) in two pieces:
log2(x) = w1 + w2,
where w1 has 53-24 = 29 bit trailing zeros.
2. Perform y*log2(x) = n+y' by simulating muti-precision
arithmetic, where |y'|<=0.5.
3. Return x**y = 2**n*exp(y'*log2)
特殊情况
1. (anything) ** 0 is 1
2. (anything) ** 1 is itself
3. (anything) ** NAN is NAN
4. NAN ** (anything except 0) is NAN
5. +-(|x| > 1) ** +INF is +INF
6. +-(|x| > 1) ** -INF is +0
7. +-(|x| < 1) ** +INF is +0
8. +-(|x| < 1) ** -INF is +INF
9. +-1 ** +-INF is NAN
10. +0 ** (+anything except 0, NAN) is +0
11. -0 ** (+anything except 0, NAN, odd integer) is +0
12. +0 ** (-anything except 0, NAN) is +INF
13. -0 ** (-anything except 0, NAN, odd integer) is +INF
14. -0 ** (odd integer) = -( +0 ** (odd integer) )
15. +INF ** (+anything except 0,NAN) is +INF
16. +INF ** (-anything except 0,NAN) is +0
17. -INF ** (anything) = -0 ** (-anything)
18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
19. (-anything except 0 and inf) ** (non-integer) is NAN
准确性
pow(x,y) returns x**y nearly rounded. In particular
pow(integer,integer)
always returns the correct integer provided it is
representable.
源代码
#ifdef __STDC__
double __ieee754_pow(double x, double y)
#else
double __ieee754_pow(x,y)
double x, y;
#endif
{
double z,ax,z_h,z_l,p_h,p_l;
double y1,t1,t2,r,s,t,u,v,w;
int i0,i1,i,j,k,yisint,n;
int hx,hy,ix,iy;
unsigned lx,ly;
i0 = ((*(int*)&one)>>29)^1; i1=1-i0;
hx = __HI(x); lx = __LO(x);
hy = __HI(y); ly = __LO(y);
ix = hx&0x7fffffff; iy = hy&0x7fffffff;
/* y==zero: x**0 = 1 */
if((iy|ly)==0) return one;
/* +-NaN return x+y */
if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) ||
iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0)))
return x+y;
/* determine if y is an odd int when x < 0
* yisint = 0 ... y is not an integer
* yisint = 1 ... y is an odd int
* yisint = 2 ... y is an even int
*/
yisint = 0;
if(hx<0) {
if(iy>=0x43400000) yisint = 2; /* even integer y */
else if(iy>=0x3ff00000) {
k = (iy>>20)-0x3ff; /* exponent */
if(k>20) {
j = ly>>(52-k);
if((j<<(52-k))==ly) yisint = 2-(j&1);
} else if(ly==0) {
j = iy>>(20-k);
if((j<<(20-k))==iy) yisint = 2-(j&1);
}
}
}
/* special value of y */
if(ly==0) {
if (iy==0x7ff00000) { /* y is +-inf */
if(((ix-0x3ff00000)|lx)==0)
return y - y; /* inf**+-1 is NaN */
else if (ix >= 0x3ff00000)/* (|x|>1)**+-inf = inf,0 */
return (hy>=0)? y: zero;
else /* (|x|<1)**-,+inf = inf,0 */
return (hy<0)?-y: zero;
}
if(iy==0x3ff00000) { /* y is +-1 */
if(hy<0) return one/x; else return x;
}
if(hy==0x40000000) return x*x; /* y is 2 */
if(hy==0x3fe00000) { /* y is 0.5 */
if(hx>=0) /* x >= +0 */
return sqrt(x);
}
}
ax = fabs(x);
/* special value of x */
if(lx==0) {
if(ix==0x7ff00000||ix==0||ix==0x3ff00000){
z = ax; /*x is +-0,+-inf,+-1*/
if(hy<0) z = one/z; /* z = (1/|x|) */
if(hx<0) {
if(((ix-0x3ff00000)|yisint)==0) {
z = (z-z)/(z-z); /* (-1)**non-int is NaN */
} else if(yisint==1)
z = -1.0*z; /* (x<0)**odd = -(|x|**odd) */
}
return z;
}
}
n = (hx>>31)+1;
/* (x<0)**(non-int) is NaN */
if((n|yisint)==0) return (x-x)/(x-x);
s = one; /* s (sign of result -ve**odd) = -1 else = 1 */
if((n|(yisint-1))==0) s = -one;/* (-ve)**(odd int) */
/* |y| is huge */
if(iy>0x41e00000) { /* if |y| > 2**31 */
if(iy>0x43f00000){ /* if |y| > 2**64, must o/uflow */
if(ix<=0x3fefffff) return (hy<0)? huge*huge:tiny*tiny;
if(ix>=0x3ff00000) return (hy>0)? huge*huge:tiny*tiny;
}
/* over/underflow if x is not close to one */
if(ix<0x3fefffff) return (hy<0)? s*huge*huge:s*tiny*tiny;
if(ix>0x3ff00000) return (hy>0)? s*huge*huge:s*tiny*tiny;
/* now |1-x| is tiny <= 2**-20, suffice to compute
log(x) by x-x^2/2+x^3/3-x^4/4 */
t = ax-one; /* t has 20 trailing zeros */
w = (t*t)*(0.5-t*(0.3333333333333333333333-t*0.25));
u = ivln2_h*t; /* ivln2_h has 21 sig. bits */
v = t*ivln2_l-w*ivln2;
t1 = u+v;
__LO(t1) = 0;
t2 = v-(t1-u);
} else {
double ss,s2,s_h,s_l,t_h,t_l;
n = 0;
/* take care subnormal number */
if(ix<0x00100000)
{ax *= two53; n -= 53; ix = __HI(ax); }
n += ((ix)>>20)-0x3ff;
j = ix&0x000fffff;
/* determine interval */
ix = j|0x3ff00000; /* normalize ix */
if(j<=0x3988E) k=0; /* |x|<sqrt(3/2) */
else if(j<0xBB67A) k=1; /* |x|<sqrt(3) */
else {k=0;n+=1;ix -= 0x00100000;}
__HI(ax) = ix;
/* compute ss = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
u = ax-bp[k]; /* bp[0]=1.0, bp[1]=1.5 */
v = one/(ax+bp[k]);
ss = u*v;
s_h = ss;
__LO(s_h) = 0;
/* t_h=ax+bp[k] High */
t_h = zero;
__HI(t_h)=((ix>>1)|0x20000000)+0x00080000+(k<<18);
t_l = ax - (t_h-bp[k]);
s_l = v*((u-s_h*t_h)-s_h*t_l);
/* compute log(ax) */
s2 = ss*ss;
r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6)))));
r += s_l*(s_h+ss);
s2 = s_h*s_h;
t_h = 3.0+s2+r;
__LO(t_h) = 0;
t_l = r-((t_h-3.0)-s2);
/* u+v = ss*(1+...) */
u = s_h*t_h;
v = s_l*t_h+t_l*ss;
/* 2/(3log2)*(ss+...) */
p_h = u+v;
__LO(p_h) = 0;
p_l = v-(p_h-u);
z_h = cp_h*p_h; /* cp_h+cp_l = 2/(3*log2) */
z_l = cp_l*p_h+p_l*cp+dp_l[k];
/* log2(ax) = (ss+..)*2/(3*log2) = n + dp_h + z_h + z_l */
t = (double)n;
t1 = (((z_h+z_l)+dp_h[k])+t);
__LO(t1) = 0;
t2 = z_l-(((t1-t)-dp_h[k])-z_h);
}
/* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
y1 = y;
__LO(y1) = 0;
p_l = (y-y1)*t1+y*t2;
p_h = y1*t1;
z = p_l+p_h;
j = __HI(z);
i = __LO(z);
if (j>=0x40900000) { /* z >= 1024 */
if(((j-0x40900000)|i)!=0) /* if z > 1024 */
return s*huge*huge; /* overflow */
else {
if(p_l+ovt>z-p_h) return s*huge*huge; /* overflow */
}
} else if((j&0x7fffffff)>=0x4090cc00 ) { /* z <= -1075 */
if(((j-0xc090cc00)|i)!=0) /* z < -1075 */
return s*tiny*tiny; /* underflow */
else {
if(p_l<=z-p_h) return s*tiny*tiny; /* underflow */
}
}
/*
* compute 2**(p_h+p_l)
*/
i = j&0x7fffffff;
k = (i>>20)-0x3ff;
n = 0;
if(i>0x3fe00000) { /* if |z| > 0.5, set n = [z+0.5] */
n = j+(0x00100000>>(k+1));
k = ((n&0x7fffffff)>>20)-0x3ff; /* new k for n */
t = zero;
__HI(t) = (n&~(0x000fffff>>k));
n = ((n&0x000fffff)|0x00100000)>>(20-k);
if(j<0) n = -n;
p_h -= t;
}
t = p_l+p_h;
__LO(t) = 0;
u = t*lg2_h;
v = (p_l-(t-p_h))*lg2+t*lg2_l;
z = u+v;
w = v-(z-u);
t = z*z;
t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
r = (z*t1)/(t1-two)-(w+z*w);
z = one-(r-z);
j = __HI(z);
j += (n<<20);
if((j>>20)<=0) z = scalbn(z,n); /* subnormal output */
else __HI(z) += (n<<20);
return s*z;
}
这篇关于Java中Math.pow()背后的算法是什么的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
