java - 在哪里可以找到 Java 的 Square Root 函数的源代码？

Question

<分区>

Answer 1

当您安装 JDK 时，可以在 src.zip 中找到标准库的源代码。不过，这对 StrictMath 没有帮助，因为 StrictMath.sqrt(double) 的实现方式如下:

public static native double sqrt(double a);

所以它实际上只是一个本地调用，Java 可能会在不同的平台上以不同的方式实现。

但是，正如 StrictMath 的文档所述:

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.

因此，通过找到适当版本的 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;