Where can I find the source code for Java's Square Root function?

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

From stackoverflow

• Download the sourcecode from the OpenJDK.

• 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

  Joey : Assembler opcodes are not OS-dependent, so that has nothing to do with Windows. But yes, the JVM will favor a native instruction as detailed in the various comments of the C source.

• 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 =  q1>>1;
  if ((q&1)==1) ix1 |= sign;
  ix0 += (m <<20);
  __HI(z) = ix0;
  __LO(z) = ix1;
  return z;
  

  Brian Agnew : Mmmm. It's almost too simple :-)

• 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).