Где я могу найти исходный код функции Square Root Java?

Question

Я знаю, что Math.sqrt называет StrictMath.sqrt (double a) Я хотел посмотреть на фактический код, используемый для его расчета.

Answer 1

Загрузите исходный код из OpenJDK.

Answer 2

Я точно не знаю, но, думаю, вы найдете алгоритм Ньютона в конечной точке.

UPD: как указано, конкретная реализация зависит от конкретной машины java. Для окон, вероятно, используется реализация ассемблера, где существует стандартный оператор sqrt

Answer 3

При установке JDK исходный код стандартной библиотеки можно найти внутри src.zip. Это не поможет для StrictMath, хотя, как StrictMath.sqrt(double) реализуется следующим образом:

public static native double sqrt(double a); 

Так что это на самом деле просто родной вызов и может быть реализован по-разному на разных платформах с помощью Java.

Однако, как документации StrictMath состояний:

Чтобы помочь обеспечить переносимость программ Java, определение некоторых числовых функций в этом пакете требует, чтобы они производят одни и те же результаты, что и некоторые опубликованные алгоритмов. Эти алгоритмы доступны из известной сетевой библиотеки netlib в виде пакета «Свободно распределяемая математическая библиотека», fdlibm. Эти алгоритмы, написанные на языке программирования C, следует понимать как выполненные со всеми операциями с плавающей запятой, которые следуют правилам арифметики с плавающей запятой Java.

Java-математическая библиотека определена в отношении версии fdlibm версии 5.3. Если fdlibm предоставляет более одного определения для функции (например, acos), используйте версию «Основная функция IEEE 754» (находящаяся в файле, имя которого начинается с буквы e). Методами, которые требуют семантики fdlibm, являются sin, cos, tan, asin, acos, atan, exp, log, log10, cbrt, atan2, pow, sinh, cosh, tanh, hypot, expm1 и log1p.

Поэтому, найдя подходящую версию источника fdlibm, вы также должны найти точную реализацию, используемую Java (и указанную здесь спецификацией).

Реализация используется fdlibm является

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 4

Поскольку я посчастливилось иметь OpenJDK валяется, я покажу здесь ее реализацию.

В JDK/SRC/доля/родной/Java/языки/StrictMath.c:

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

jsqrt определяется как sqrt в JDK/SRC/доля/родной/Java/языки/fdlibm/SRC/w_sqrt.c (имя изменено через препроцессор):

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

И __ieee754_sqrt определяется в JDK/SRC/доля/родной/Java/языки/fdlibm/SRC/e_sqrt.с как:

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

Есть обильные комментарии в файле, объясняющие используемые методы, которые я опущенные для (пола) краткости. Here's the file in Mercurial (Я надеюсь, что это правильный способ связаться с ним).