Where can I find the source code for Java's Sq

Question

I know that Math.sqrt calls StrictMath.sqrt(double a) I was wanting to look at the actual code used to calculate it.

Answer 1

When you install a JDK the source code of the standard library can be found inside src.zip. This won't help you for StrictMath, though, as StrictMath.sqrt(double) is implemented as follows:

public static native double sqrt(double a);

So it's really just a native call and might be implemented differently on different platforms by Java.

However, as the documentation of StrictMath states:

To help ensure portability of Java programs, the definitions of some of the numeric functions in this package require that they produce the same results as certain published algorithms. These algorithms are available from the well-known network library netlib as the package "Freely Distributable Math Library," fdlibm. These algorithms, which are written in the C programming language, are then to be understood as executed with all floating-point operations following the rules of Java floating-point arithmetic.

The Java math library is defined with respect to fdlibm version 5.3. Where fdlibm provides more than one definition for a function (such as acos), use the "IEEE 754 core function" version (residing in a file whose name begins with the letter e). The methods which require fdlibm semantics are sin, cos, tan, asin, acos, atan, exp, log, log10, cbrt, atan2, pow, sinh, cosh, tanh, hypot, expm1, and log1p.

So by finding the appropriate version of the fdlibm source, you should also find the exact implementation used by Java (and mandated by the specification here).

The implementation used by fdlibm is

static const double one = 1.0, tiny=1.0e-300;

double z;
int sign = (int) 0x80000000; 
unsigned r, t1, s1, ix1, q1;
int ix0, s0, q, m, t, i;

ix0 = __HI(x); /* high word of x */
ix1 = __LO(x); /* low word of x */

/* take care of Inf and NaN */
if ((ix0 & 0x7ff00000) == 0x7ff00000) {            
    return x*x+x; /* sqrt(NaN) = NaN, 
                     sqrt(+inf) = +inf,
                     sqrt(-inf) = sNaN */
} 

/* take care of zero */
if (ix0 <= 0) {
    if (((ix0&(~sign)) | ix1) == 0) {
        return x; /* sqrt(+-0) = +-0 */
    } else if (ix0 < 0) {
        return (x-x) / (x-x); /* sqrt(-ve) = sNaN */
    }
}

/* normalize x */
m = (ix0 >> 20);
if (m == 0) { /* subnormal x */
    while (ix0==0) {
        m -= 21;
        ix0 |= (ix1 >> 11); ix1 <<= 21;
    }
    for (i=0; (ix0&0x00100000)==0; i++) {
        ix0 <<= 1;
    }
    m -= i-1;
    ix0 |= (ix1 >> (32-i));
    ix1 <<= i;
}

m -= 1023; /* unbias exponent */
ix0 = (ix0&0x000fffff)|0x00100000;
if (m&1) { /* odd m, double x to make it even */
    ix0 += ix0 + ((ix1&sign) >> 31);
    ix1 += ix1;
}

m >>= 1; /* m = [m/2] */

/* generate sqrt(x) bit by bit */
ix0 += ix0 + ((ix1 & sign)>>31);
ix1 += ix1;
q = q1 = s0 = s1 = 0; /* [q,q1] = sqrt(x) */
r = 0x00200000; /* r = moving bit from right to left */

while (r != 0) {
    t = s0 + r; 
    if (t <= ix0) { 
        s0 = t+r; 
        ix0 -= t; 
        q += r; 
    } 
    ix0 += ix0 + ((ix1&sign)>>31);
    ix1 += ix1;
    r>>=1;
}

r = sign;
while (r != 0) {
    t1 = s1+r; 
    t = s0;
    if ((t<ix0) || ((t == ix0) && (t1 <= ix1))) { 
        s1 = t1+r;
        if (((t1&sign) == sign) && (s1 & sign) == 0) {
            s0 += 1;
        }
        ix0 -= t;
        if (ix1 < t1) {
            ix0 -= 1;
        }
        ix1 -= t1;
        q1  += r;
    }
    ix0 += ix0 + ((ix1&sign) >> 31);
    ix1 += ix1;
    r >>= 1;
}

/* use floating add to find out rounding direction */
if((ix0 | ix1) != 0) {
    z = one - tiny; /* trigger inexact flag */
    if (z >= one) {
        z = one+tiny;
        if (q1 == (unsigned) 0xffffffff) { 
            q1=0; 
            q += 1;
        }
    } else if (z > one) {
        if (q1 == (unsigned) 0xfffffffe) {
            q+=1;
        }
        q1+=2; 
    } else
        q1 += (q1&1);
    }
}

ix0 = (q>>1) + 0x3fe00000;
ix1 =  q 1>> 1;
if ((q&1) == 1) ix1 |= sign;
ix0 += (m <<20);
__HI(z) = ix0;
__LO(z) = ix1;
return z;

Answer 2

I don't know exactly but I think you'll find Newton's algorithm at the end point.

UPD: as comments say concrete implementation depends on concrete java machine. For windows it's probably using assembler implementation, where the standard operator sqrt exists

Answer 3

Since I happen to have OpenJDK lying around, I'll show its implementation here.

In jdk/src/share/native/java/lang/StrictMath.c:

JNIEXPORT jdouble JNICALL
Java_java_lang_StrictMath_sqrt(JNIEnv *env, jclass unused, jdouble d)
{
    return (jdouble) jsqrt((double)d);
}

jsqrt is defined as sqrt in jdk/src/share/native/java/lang/fdlibm/src/w_sqrt.c (the name is changed through the preprocessor):

#ifdef __STDC__
        double sqrt(double x)           /* wrapper sqrt */
#else
        double sqrt(x)                  /* wrapper sqrt */
        double x;
#endif
{
#ifdef _IEEE_LIBM
        return __ieee754_sqrt(x);
#else
        double z;
        z = __ieee754_sqrt(x);
        if(_LIB_VERSION == _IEEE_ || isnan(x)) return z;
        if(x<0.0) {
            return __kernel_standard(x,x,26); /* sqrt(negative) */
        } else
            return z;
#endif
}

And __ieee754_sqrt is defined in jdk/src/share/native/java/lang/fdlibm/src/e_sqrt.c as:

#ifdef __STDC__
static  const double    one     = 1.0, tiny=1.0e-300;
#else
static  double  one     = 1.0, tiny=1.0e-300;
#endif

#ifdef __STDC__
        double __ieee754_sqrt(double x)
#else
        double __ieee754_sqrt(x)
        double x;
#endif
{
        double z;
        int     sign = (int)0x80000000;
        unsigned r,t1,s1,ix1,q1;
        int ix0,s0,q,m,t,i;

        ix0 = __HI(x);                  /* high word of x */
        ix1 = __LO(x);          /* low word of x */

    /* take care of Inf and NaN */
        if((ix0&0x7ff00000)==0x7ff00000) {
            return x*x+x;               /* sqrt(NaN)=NaN, sqrt(+inf)=+inf
                                           sqrt(-inf)=sNaN */
        }
    /* take care of zero */
        if(ix0<=0) {
            if(((ix0&(~sign))|ix1)==0) return x;/* sqrt(+-0) = +-0 */
            else if(ix0<0)
                return (x-x)/(x-x);             /* sqrt(-ve) = sNaN */
        }
    /* normalize x */
        m = (ix0>>20);
        if(m==0) {                              /* subnormal x */
            while(ix0==0) {
                m -= 21;
                ix0 |= (ix1>>11); ix1 <<= 21;
            }
            for(i=0;(ix0&0x00100000)==0;i++) ix0<<=1;
            m -= i-1;
            ix0 |= (ix1>>(32-i));
            ix1 <<= i;
        }
        m -= 1023;      /* unbias exponent */
        ix0 = (ix0&0x000fffff)|0x00100000;
        if(m&1){        /* odd m, double x to make it even */
            ix0 += ix0 + ((ix1&sign)>>31);
            ix1 += ix1;
        }
        m >>= 1;        /* m = [m/2] */

    /* generate sqrt(x) bit by bit */
        ix0 += ix0 + ((ix1&sign)>>31);
        ix1 += ix1;
        q = q1 = s0 = s1 = 0;   /* [q,q1] = sqrt(x) */
        r = 0x00200000;         /* r = moving bit from right to left */

        while(r!=0) {
            t = s0+r;
            if(t<=ix0) {
                s0   = t+r;
                ix0 -= t;
                q   += r;
            }
            ix0 += ix0 + ((ix1&sign)>>31);
            ix1 += ix1;
            r>>=1;
        }

        r = sign;
        while(r!=0) {
            t1 = s1+r;
            t  = s0;
            if((t<ix0)||((t==ix0)&&(t1<=ix1))) {
                s1  = t1+r;
                if(((t1&sign)==sign)&&(s1&sign)==0) s0 += 1;
                ix0 -= t;
                if (ix1 < t1) ix0 -= 1;
                ix1 -= t1;
                q1  += r;
            }
            ix0 += ix0 + ((ix1&sign)>>31);
            ix1 += ix1;
            r>>=1;
        }

    /* use floating add to find out rounding direction */
        if((ix0|ix1)!=0) {
            z = one-tiny; /* trigger inexact flag */
            if (z>=one) {
                z = one+tiny;
                if (q1==(unsigned)0xffffffff) { q1=0; q += 1;}
                else if (z>one) {
                    if (q1==(unsigned)0xfffffffe) q+=1;
                    q1+=2;
                } else
                    q1 += (q1&1);
            }
        }
        ix0 = (q>>1)+0x3fe00000;
        ix1 =  q1>>1;
        if ((q&1)==1) ix1 |= sign;
        ix0 += (m <<20);
        __HI(z) = ix0;
        __LO(z) = ix1;
        return z;
}

There are copious comments in the file explaining the methods used, which I have omitted for (semi-) brevity. Here's the file in Mercurial (I hope this is the right way to link to it).

Answer 4

Download the sourcecode from the OpenJDK.