Best Practice Laufzeitvergleich Math.sqrt() aus Java und sqrt() aus C++

[OmegaPirat28]

Hallo,

Ich habe eine Simulation geschrieben und sie sowohl in Java als auch in C++ implementiert.
Als ich die Laufzeit verglich, erstaunte ich, dass das Java viermal so schnell war wie das C++ Programm.
Ich konnte jetzt genau einkesseln an welcher Stelle es liegt, dass das C++ Programm soviel langsamer ist und zwar müssen im Programm sehr viele Wurzeln berechnet werden.
Ich habe die Laufzeit der entsprechenden Bibliotheksfunktion von Java und C++ miteinander verglichen.
Auf meinem System werden in Java etwa 0,3 Sekunden benötigt um eine Milliarde Wurzeln zu berechnen, während die Funktion sqrt() aus C++ etwa 30 Sekunden braucht.
Das ist ein Faktor 100. Wieso ist die Javafunktion soviel schneller und was kann ich machen, um die Wurzelberechnungen in C++ zu optimieren?
Weiß eigentlich jemand wie die Math.sqrt() implementiert ist, weil die mir als verdammt schnell erscheint, auch mit eigenen Implementierungen zur Wurzelberechnung kam ich nicht annähernd an die Geschwindigkeit ran.

Vielen Dank im
Voraus

 

[tommysenf]

Math.sqrt leitet den Aufruf an StrictMath.sqrt weiter. Dies ist eine native Function und somit in C implementiert. Die Implementierung basiert auf fdlibm und kann hier gefunden werden (OpenJDK):
jdk/src/share/native/java/lang/fdlibm/src/e_sqrt.c

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;