Math in Java What is the algorithm behind pow()?

I've just started using Java. As my first project, I'm writing a program to find the root of a given number (in this case, the cube root) At present, I am trying to achieve this goal This 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));


    }
}

This can only reach 11 digits because it is too large for the IDE But if I just use math Pow (number, 1.0 / 3.0), which can be used for larger numbers and can be calculated at any time

So, math What algorithm does pow () use to give an immediate answer? I know my method depends on guessing, I guess math Pow () may actually calculate the answer, but how?

Solution

This is an interesting question If you look at the source code of the Java math class, you will find that it calls strictmath Pow (double1, double2), the signature of strictmath is public static native pow (double a, double B);

So, in the end, this is a real local phone, which may vary depending on the platform However, having an implementation somewhere is not easy to see Here is the description of the function and the code of the function itself:

be careful

Looking at mathematics and trying to understand it may inevitably lead to more problems However, you can better understand the native function by searching the GitHub in Java math function source code and scanning the mathematical summary Happy exploration:)

Method description

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)

exceptional case

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,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

accuracy

pow(x,y) returns x**y nearly rounded. In particular
                      pow(integer,integer)
       always returns the correct integer provided it is
       representable.

source code

#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;
}