转载

卡马克是最快的开根号方法吗

intro1 ：卡马克算法时最快的开根号方式吗？C/Java语言本身是怎么实现开根号的？

intro2 : java的内置sqrt和c的内置sqrt哪个更快？Java的编译/JIT优化和GCC的编译优化是否有不足之处？

intro3 : java的C2优化效果一定比C1效果好吗(对性能而言)？

问题或现象

这是最近看一位博主解Leet Code题想到的，原题简化一下是：给一个正整数(32位int)开根号后得到x，再对x取整返回。

博文的解法是使用二分查找，Java代码实现，不过这里想对该题再讨论几点

1）二分查找也可以优化下，建立一个简单的范围表，再二分查找，某几个国产IP库查询也是该做法(因为比BTree省太多内存)。

2）其实还可以用牛顿切线法，每个ACMer入门练手时都会碰到的算法。

3）因为本题目只是要求返回正整数，那么如果我在牛顿切分法阈值判断的时候，再加一个条件，判断本轮的整数部分和上一轮的整数部分是否相等，会不会更快？

4）开根号怎能少了卡马克算法，要知道在关于开根号方法中祭出卡马克算法，也就基本意味着本次交谈该结束了，但是卡马克是最快的方法吗？

所以下文中我写了几行代码来简单验证下(C/Java版)

先了解下牛顿切分法是什么：

许多网站链接有介绍，这里简单描述下，中学课本里的二次方程求根公式正式提出距今不过1000年多，高次方程求根公式知道高斯/阿贝尔/伽罗华才告一段落，在此之前物理学家们怎么做的呢？牛顿提出了牛顿迭代法，但这是在实/复数域上求解方程的近似根，思想就是通过直线逼近曲线(世纪更早的我国数学家刘徽也提出该思想并求出圆周率近似值)，只不过那时候没有系统的考虑连续和收敛的问题。

牛顿切分转化为牛顿迭代就是比如对方程 f(n) => x(n+1)=F(x(n))转换，也就是从当前态计算出下一个状态，这是适合人们手工去计算的，更是非常适合计算机代码和执行。

我们把上述 f(n) 用于开根号情况，也就是求方程 X^2-c=0 的解，借助于大学里的泰勒展开，我们可以得到：

X(n+1) = ( X(n) + C/X(n) )/2

当X(n+1)与X(n) 差值足够小的时候，我们就认为 X(n+1) 接近于真实解了。

至于卡马克算法讨论需要更多篇幅，网上有更多介绍。简单来说是基于牛顿迭代使用一个“魔数”免于多次迭代近似对单精度数字求根。

下面代码是C实现。

其中：

1) sqrt_newton：牛顿迭代，sqrt_newton_int int比较的牛顿迭代，InvSqrt 卡马克算法开根号

2) 在main函数里，我使用简单的for循环来统计耗时，暂不考虑精度以及其他问题，且先相信只要循环次数大，耗时差距大，就可以得出足够可信的结论。

#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#include<time.h>

double FLT_MIN = 1e-7;

double sqrt_newton(int x){
    if(x <= 0)  return 0;
    double res, lastres;
    res = x;    //初始值，可以为任意非0的值
    do{
        lastres = res;
        res = (res + x/res)/2;
    }while(fabs(lastres-res) > FLT_MIN);
    return res;
}

double sqrt_newton_int(int x){
    if(x <= 0)  return 0;
    double res, lastres;
    res = x;
    int last = 0;
    do{
        lastres = res;
        res = (res + x/res)/2;
        int cur = (int)res;
        if (last == cur){
            return res;
        }
        last = cur;
    }while(fabs(lastres-res) > FLT_MIN);
    return res;
}

float InvSqrt(float x){
    float xhalf = 0.5f * x;
    int i = *(int*)&x;
    i = 0x5f375a86 - (i>>1);
    x = *(float*)&i;
    x = x*(1.5f-xhalf*x*x);
    x = x*(1.5f-xhalf*x*x);
    x = x*(1.5f-xhalf*x*x);
    return 1/x;
}

int main(int argc, char **argv){
    clock_t begin, end;
    int num = atoi(argv[1]);
    double res1=0,res2=0,res4=0;
    float res3=0;
    int i;
    float num_f;
    int loopcnts = 1000000;
    if(sizeof(argv)>2){
        loopcnts = atoi(argv[2]);
    }
 
    begin = clock();
    for(i = 0; i < loopcnts; i++){
        res1+=sqrt_newton_int(num+i);
        res2+=sqrt_newton(num+i);
        res3+=InvSqrt(num+i);
        res4+=sqrt(num+i);
    }
    end = clock();
    printf("hot %f/t:%f,%f,%f,%f/n", (double)(end-begin),res1,res2,res3,res4);

    begin = clock();
    for(i = 0; i < loopcnts; i++)
        res1 = sqrt_newton(num);
    end = clock();
    printf(" newton_cos(%d) = %f, /t/tcost: %f/n", num, res1, (double)(end-begin));
 
    begin = clock();
    for(i = 0; i < loopcnts; i++)
        res2 = sqrt_newton_int(num);
    end = clock();
    printf("newton_int(%d) = %f, /t/tcost: %f/n", num, res2, (double)(end-begin));
 
    num_f=1.0f*num;
    begin = clock();
    for(i = 0; i < loopcnts; i++)
        res3 = InvSqrt(num_f);
    end = clock();
    printf(" InvSqrt(%d) = %f, /t/tcost: %f/n", num, res3, (double)(end-begin));
 
    begin = clock();
    for(i = 0; i < loopcnts; i++)
        res4 = sqrt(num);
    end = clock();
    printf("sys sqrt(%d) = %f, /t/tcost: %f/n", num, res4, (double)(end-begin));
    // printf("res:%f, %f, %f, %f", res1,res2,res3,res4);
}

因int范围限制，故每轮采取循环1亿次，多运行几次，发现输出下面结果

➜  Test gcc quic_newton.c -o quic_newton     
➜  Test ./quic_newton 2855 100000000   
hot 40405059.000000 :666695110416.069336,666695110197.906128,274877906944.000000,666695110197.906128
 newton_cos(2855) = 53.432200,    cost: 11803758.000000
newton_int(2855) = 53.432201,     cost: 10048574.000000
 InvSqrt(2855) = 53.432198,     cost: 967581.000000
sys sqrt(2855) = 53.432200,     cost: 614861.000000
➜  Test
➜  Test gcc -O3  quic_newton.c -o quic_newton
➜  Test ./quic_newton 2855 100000000         
hot 14997438.000000 :666695110416.069336,666695110197.906128,274877906944.000000,666695110197.906128
 newton_cos(2855) = 53.432200,    cost: 3292126.000000
newton_int(2855) = 53.432201,     cost: 2845886.000000
 InvSqrt(2855) = 53.432198,     cost: 1.000000
sys sqrt(2855) = 53.432200,     cost: 0.000000
...
➜  Test gcc -O1  quic_newton.c -o quic_newton
➜  Test ./quic_newton 2855 100000000         
hot 19265085.000000 :666695110416.069336,666695110197.906128,274877906944.000000,666695110197.906128
 newton_cos(2855) = 53.432200,    cost: 1.000000
newton_int(2855) = 53.432201,     cost: 0.000000
 InvSqrt(2855) = 53.432198,     cost: 0.000000
sys sqrt(2855) = 53.432200,     cost: 0.000000
➜  Test

1) 忽略上述InvSqrt float数值不正确的情况，这里是为了避免类型转换，减少对结果的影响。

2) 虽有for循环，以及类型/加操作，但如果测试下来耗时是一半，那有理由相信性能其实是大于一倍的。

3) 可发现第一次编译后，运行几次，那段卡马克效果都不如使用C系统库的sqrt函数，相差0.5倍左右，不过已经比牛顿迭代少一个数量级了。

4) 当开启 -O3 优化的时候，出现耗时为1，这是GCC编译优化的缘故，也可看到优化后调用牛顿迭代耗时下降。

5) 但如果看 hot 的总耗时，还是可见-O3 效果优于 -O1 优化的。

6) 开启 -O1 优化时候，看到耗时都为1。因为优化实际上相当只运行最后一次循环结果，这点看汇编代码可以看到。不过这里有个问题是单是对于main函数而言，-O3优化效果反而不如 -O1。

所以为了避免上述问题，这里需要改进下代码:见附件quic_newton.c

1) 上文中 res4 = sqrt(num); 都类似对应的改为 res4 += sqrt(num+i)格式

2) hot之后将 res1-4都置零一下。

3) int num 改为double num，并且相应printf语句里 %d改为%f。

（部分可能受制于发生隐式类型转换，不过因为已经是double/float，对于再加int性能消耗是可忽略，或可能通过 -O1优化掉。查计算机指令周期可知，相对于前三个计算方法这个时间影响也很小）

➜  Test gcc quic_newton.c -o quic_newton
➜  Test ./quic_newton 2500 100000000
hot 35850465.000000 :666691578733.306519,666691578514.627563,274877906944.000000,666691578514.627686
 newton_inl(2500.000000) = 666691578514.627563,     cost: 17829585.000000
 newton_cos(2500.000000) = 666691578514.627563,     cost: 17944181.000000
newton_int(2500.000000) = 666691578733.306519,    cost: 18475511.000000
 InvSqrt(2500.000000) = 274877906944.000000,    cost: 3108139.000000
sys sqrt(2500.000000) = 666691578514.627686,    cost: 252962.000000
➜  Test gcc -O3 quic_newton.c -o quic_newton
➜  Test ./quic_newton 2500 100000000       
hot 15424798.000000 :666691578733.306519,666691578514.627563,274877906944.000000,666691578514.627686
 newton_inl(2500.000000) = 666691578514.627563,     cost: 7050614.000000
 newton_cos(2500.000000) = 666691578514.627563,     cost: 7113696.000000
newton_int(2500.000000) = 666691578733.306519,    cost: 7668225.000000
 InvSqrt(2500.000000) = 274877906944.000000,    cost: 417802.000000
sys sqrt(2500.000000) = 666691578514.627686,    cost: 168760.000000
➜  Test ./quic_newton 2500 100000000
hot 15394188.000000 :666691578733.306519,666691578514.627563,274877906944.000000,666691578514.627686
 newton_inl(2500.000000) = 666691578514.627563,     cost: 7063698.000000
 newton_cos(2500.000000) = 666691578514.627563,     cost: 7111073.000000
newton_int(2500.000000) = 666691578733.306519,    cost: 7663540.000000
 InvSqrt(2500.000000) = 274877906944.000000,    cost: 416777.000000
sys sqrt(2500.000000) = 666691578514.627686,    cost: 167816.000000
➜  Test

请注意上述命令执行时的参数和输出，可以看到InvSqrt接近一个数量级的提升，sys sqrt 40%提升。这里的优化显然是减少了循环体内的指令。

比如：

begin = clock();
for(i = 0; i < loopcnts; i++)
    res4 += sqrt(num+i);
end = clock();

对比一下 -O0和-O1关键部分

# -O0
  callq _clock
  movq  %rax, -24(%rbp)
  movl  $0, -72(%rbp)
LBB3_24:                                ## =>This Inner Loop Header: Depth=1
  movl  -72(%rbp), %eax
  cmpl  -80(%rbp), %eax
  jge LBB3_27
## %bb.25:                              ##   in Loop: Header=BB3_24 Depth=1
  movsd -40(%rbp), %xmm0        ## xmm0 = mem[0],zero
  movl  -72(%rbp), %eax
  cvtsi2sdl %eax, %xmm1
  addsd %xmm1, %xmm0
  sqrtsd  %xmm0, %xmm0
  addsd -64(%rbp), %xmm0
  movsd %xmm0, -64(%rbp)
## %bb.26:                              ##   in Loop: Header=BB3_24 Depth=1
  movl  -72(%rbp), %eax
  addl  $1, %eax
  movl  %eax, -72(%rbp)
  jmp LBB3_24
LBB3_27:
  callq _clock


#-O1
  callq _clock
  movsd -48(%rbp), %xmm3        ## 8-byte Reload
                                        ## xmm3 = mem[0],zero
  movq  %rax, %r14
  testl %r15d, %r15d
  xorpd %xmm4, %xmm4
  jle LBB3_22
## %bb.20:
  xorpd %xmm0, %xmm0
  movsd LCPI3_0(%rip), %xmm1    ## xmm1 = mem[0],zero
  xorpd %xmm4, %xmm4
.p2align  4, 0x90
LBB3_21:                                ## =>This Inner Loop Header: Depth=1
  movapd  %xmm3, %xmm2
  addsd %xmm0, %xmm2
  sqrtsd  %xmm2, %xmm2
  addsd %xmm2, %xmm4
  addsd %xmm1, %xmm0
  decl  %r15d
  jne LBB3_21
LBB3_22:
  movsd %xmm4, -40(%rbp)        ## 8-byte Spill
  callq _clock

-O1 优化掉了一个int转浮点型，而且还巧妙的通过自增的方式减少了指令。即优化后相当于

double initial=2500f;
for(int i=10...0;i>0;i--){
  sqrt(initial);
  initial=initial+1.0f
}

上面这个优化有趣，我在后面的Java版代码里也会讨论下这一点。

需要说明一点，改进后的quic_newton.c耗时统计，跟你添加的gcc参数相关(O1/3支持参数多样来控制你的优化)，当然，也跟你的机器，系统，硬件性能相关。比如我在一台旧的intel xeon+centos5上编译后运行，InvSqrt也在4.1秒左右，sys sqrt在1.8秒了，是mac机结果十倍了，这个时候是看不到 -O3和-O1效果的差异的。

至此，要说的其实就是，现代计算机(无论intel/amd)浮点处理器大多支持开根号指令，基于硬件的，虽然可能硬件(arm/fpga)不同他们的实现或有区别，比如x86架构该指令就是基于牛顿迭代实现的，因为精度可控。不是所有的处理器都这么干。

不过intel系列最初是基于此提供FSQRT指令，后来提供了更快的指令，SQRTSS (对于双精度是 SQRTSD) 。如果算上支持4个浮点数那就更快了。

好了，写这些，就是希望记得有人问最快的开根号代码，切勿上来就撸一个卡马克算法。

不过也有硬件指令(RSQRTSS)是基于卡马克近似算法，性能据说比上述指令快，但是精度就略差了(小于六千分之一),而且仅限单精度。 RSQRTSS

sqrtss gives a correctly rounded result.  rsqrtss gives an approximation to the reciprocal, accurate to about 11 bits.

如果你对这方面了解更多，欢迎指教 mail: aXRob21hc2xhdUBxcS5jb20=

需要指出上面c代码测试结果会因平台而异，比如在Intel(R) Xeon某型号 + gcc (GCC) 4.4.7 20120313 (Red Hat 4.4.7-18)上就没有sqrt的性能优化，而是call指令调用，导致-O1编译没能够优化到常量(但是循环跑Invsqrt 则优化到常量时间了)。

同样的，上述测试的是比较的趋势，数值不作为参考，毕竟还可以有比ANSI C 的sqrt更快的选择，比如第三方(fdlibm)优化过的或者含有fpu判断的 __ieee754_sqrt之类函数。

但如果你是C和汇编语言高手的话，在Intel平台，你可以直接调用xmmintrin.h库提供的 _mm_rsqrt_ss/_mm_rsqrt sd 之类的函数实现模拟调用汇编开根号/ _m128d _mm_sqrt pd(__m128d a)，甚至借助 / _m128d 同时能给四个double数值开根号的特性实现四线程并行，这些就比ANSI C内置的sqrt更快了。

下面来看看Java里效果怎么样

// JAVA 1
public class MathFunc{
    public static double sqrt_(double t, Double precise){
        double x0 = t, x1, differ;
        double prec = 1e-7;
        while (true) {
            x1 = (x0/2 + t/(2*x0));
            differ = x1 * x1 - t;
            if (differ < prec && differ > -prec) {
                return x1;
            }
            x0 = x1;
        }
    }
   
    public static double sqrt_INT(double t, Double precise){
        double x0 = t, x1, differ;
        double prec = 1e-7;
        int last = (int)x0;
        while (true) {
            x1 = (x0/2 + t/(2*x0));
            differ = x1 * x1 - t;
            if ((differ <= prec && differ >= -prec)) {
                return x1;
            }
            x0 = x1;
            int cur = (int)x0;
            if (cur == last) {
                return cur;
            }
            last = cur;
        }
    }
    public static void main(String[] args){
        for (int i = 0; i < 1000; i++) {
            sqrt_(230, null);
            sqrt_INT(230, null);
            Math.sqrt(230000);
        }
        double MM = 2855;//230000
        double s1=0,s2=0,s3=0;
        long start = System.nanoTime();
        int NN = 1_0000_0000;
        for (int i = 0; i < NN; i++) {
            s1 += sqrt_(MM, null);
        }
        long start2 = System.nanoTime();
        for (int i = 0; i < NN; i++) {
            s2+= sqrt_INT(MM, null);
        }
        long start3 = System.nanoTime();
        for (int i = 0; i < NN; i++) {
            s3+= Math.sqrt(MM+i);
        }
        long start4 = System.nanoTime();
        System.out.println(String.format("%s-%s-%s", (start2-start)/1e3, (start3-start2)/1e3, (start4-start3)/1e3));
        System.out.println(String.format("%s-%s-%s", s1,s2,s3));
    }
}

上文问题3，这里解答下：

如果纯粹从理论上讲，通过判断减少循环次数，则耗时应该减少，但实际上并不这样，因为double强转int也是耗时的操作，当这个耗时足以弥补减少循环的耗时时，才会得到更少的耗时。可以简单修改 MM 来验证下，实际上MM=2500/100 则确实耗时更少，MM=2837耗时则几乎相等，MM=2839/2840/2850则耗时多。

因为接近平方数 INT版本收敛快。

上面验证较简单，就不展开。不过从这段代码，让我们看一点有趣的发现：

// JAVA 2
public class TwoSqrt{
    final static double MM = 2855;
    static int NN = 1_0000_0000;
    static double s3 = 0;
    static long start3 = 0, start4 = 0;
    static double j = 0;

    public static void main(String[] args){
        for (int i = 0; i < 100; i++) {
            s3 = 0;j=0;
            long d = D();
            System.out.println("D:" + s3 + "/t" + d / 1e3);
            s3 = 0;j=0;
            long dd = DD();
            System.out.println("DD:" + s3 + "/t" + dd / 1e3);
            s3 = 0;
            long di = DI();
            System.out.println("DI:" + s3 + "/t" + di / 1e3);
        }
       
    }
   
    static long D(){
        //int s3=0;
        long start3 = System.nanoTime();
        for (int i = 0; i < NN; i++) {
            s3 += Math.sqrt(MM);//0.1;//
        }
        long start4 = System.nanoTime();
        return start4 - start3;
    }

    static long DD(){
        long start3 = System.nanoTime();
        for (int i = 0; i < NN; i++) {
            s3 += Math.sqrt(MM + j++);
        }
        long start4 = System.nanoTime();
        return start4 - start3;
    }

    static long DI(){
        long start3 = System.nanoTime();
        for (int i = 0; i < NN; i++) {
            s3 += Math.sqrt(MM + i);
        }
        long start4 = System.nanoTime();
        return start4 - start3;
    }
}

注意上面main函数里，我把每个方法都跑了一百次，并且每个方法内部还有1亿次循环。

我在这里新写一个JAVA类，不仅仅是方便JAVA/JVM/汇编代码，还是为了引出下文发现的三个问题。

每个函数都赋值给静态变量(该方法作用域外)s3，是为了避免逃逸变量的编译优化，这个简单就不详述。

不过，让我们先看 D() 这个方法，含逃逸变量，即可能会被JIT优化的情况

static long D(){
    double s3=0;
    long start3 = System.nanoTime();
    for (int i = 0; i < NN; i++) {
        s3 += Math.sqrt(MM);//0.1;//i;//0.2;//
    }
    long start4 = System.nanoTime();
    return start4 - start3;
}

上述代码的 D() 在外层循环跑两轮之后就会发现 D() 耗时几乎为零了。

单单看JVM code可能看不出来

 0: iconst_0       
 1: istore_0       
 2: invokestatic    #20  // Method java/lang/System.nanoTime:()J
 5: lstore_1       
 6: iconst_0       
 7: istore_3       
 8: iload_3         
 9: getstatic       #21  // Field NN:I
12: if_icmpge       25   
15: iload_0         
16: iload_3         
17: iadd           
18: istore_0       
19: iinc            3, 1
22: goto            8   
25: invokestatic    #20  // Method java/lang/System.nanoTime:()J
28: lstore_3       
29: lload_3         
30: lload_1         
31: lsub           
32: lreturn

我们看看 C2优化之后的汇编代码

....
0x000000011408c5d6: callq *%r10  ;*invokestatic nanoTime
                                 ; - TwoSqrt::D@2 (line 25)
0x000000011408c5d9: mov %rax,%rbx
0x000000011408c5dc: movabs $0x10728ec10,%r10
0x000000011408c5e6: callq *%r10  ;*invokestatic nanoTime
                                 ; - TwoSqrt::D@25 (line 29)
0x000000011408c5e9: sub %rbx,%rax  ;*lsub
                                   ; - TwoSqrt::D@31 (line 30)
....

可以看到for循环其实优化掉了，所以上述改一下，return s3; 这样循环显然是不会消除优化掉的。

好了，我们再改一下：

static long D(){
    // double s3=0;
    long start3 = System.nanoTime();
    for (int i = 0; i < NN; i++) {
        s3 += i;
    }
    long start4 = System.nanoTime();
    return start4 - start3;
}

即使用全局变量 s3以及改为 s3+=i，那么耗时会是常量还是跟for循环次数有关？

更多的，如果这里改为 s3+=0.1 或者 s3+=1；又会怎样？

实际上三者耗时不一样，但都不是几纳秒，即这个时间不是常量，和循环次数有关，也就是说， 即便是 C2优化后，循环依旧。

0x000000011a2f1130: vaddsd %xmm1,%xmm0,%xmm0
0x000000011a2f1134: mov    %r11d,%r8d
0x000000011a2f1137: add    $0xd,%r8d
0x000000011a2f113b: mov    %r11d,%r9d
0x000000011a2f113e: add    $0xe,%r9d
0x000000011a2f1142: vcvtsi2sd %r8d,%xmm3,%xmm3
0x000000011a2f1147: vcvtsi2sd %r9d,%xmm4,%xmm4
0x000000011a2f114c: mov    %r11d,%r8d
0x000000011a2f114f: add    $0xc,%r8d
0x000000011a2f1153: mov    %r11d,%r9d
0x000000011a2f1156: add    $0xb,%r9d
0x000000011a2f115a: vcvtsi2sd %r8d,%xmm5,%xmm5
0x000000011a2f115f: vcvtsi2sd %r9d,%xmm6,%xmm6
0x000000011a2f1164: mov    %r11d,%r8d
0x000000011a2f1167: add    $0xa,%r8d
0x000000011a2f116b: mov    %r11d,%r9d
0x000000011a2f116e: add    $0x8,%r9d
0x000000011a2f1172: vcvtsi2sd %r8d,%xmm7,%xmm7
0x000000011a2f1177: vcvtsi2sd %r9d,%xmm8,%xmm8
0x000000011a2f117c: mov    %r11d,%r8d
0x000000011a2f117f: add    $0x7,%r8d
0x000000011a2f1183: mov    %r11d,%r9d
0x000000011a2f1186: add    $0x6,%r9d
0x000000011a2f118a: vcvtsi2sd %r8d,%xmm9,%xmm9
0x000000011a2f118f: vcvtsi2sd %r9d,%xmm10,%xmm10
0x000000011a2f1194: mov    %r11d,%r8d
0x000000011a2f1197: add    $0x5,%r8d
0x000000011a2f119b: mov    %r11d,%r9d
0x000000011a2f119e: add    $0x4,%r9d
0x000000011a2f11a2: vcvtsi2sd %r8d,%xmm11,%xmm11
0x000000011a2f11a7: vcvtsi2sd %r9d,%xmm12,%xmm12
0x000000011a2f11ac: mov    %r11d,%r8d
0x000000011a2f11af: add    $0x3,%r8d
0x000000011a2f11b3: mov    %r11d,%r9d
0x000000011a2f11b6: add    $0x2,%r9d
0x000000011a2f11ba: vcvtsi2sd %r8d,%xmm13,%xmm13
0x000000011a2f11bf: vcvtsi2sd %r9d,%xmm1,%xmm1
0x000000011a2f11c4: mov    %r11d,%r8d
0x000000011a2f11c7: add    $0xf,%r8d
0x000000011a2f11cb: mov    %r11d,%r9d
0x000000011a2f11ce: inc    %r9d
0x000000011a2f11d1: vcvtsi2sd %r8d,%xmm14,%xmm14
0x000000011a2f11d6: vcvtsi2sd %r9d,%xmm2,%xmm2
0x000000011a2f11db: vaddsd %xmm0,%xmm2,%xmm0
0x000000011a2f11df: vaddsd %xmm1,%xmm0,%xmm0
0x000000011a2f11e3: vaddsd %xmm13,%xmm0,%xmm0
0x000000011a2f11e8: vaddsd %xmm12,%xmm0,%xmm0
0x000000011a2f11ed: vaddsd %xmm11,%xmm0,%xmm0
0x000000011a2f11f2: vaddsd %xmm10,%xmm0,%xmm0
0x000000011a2f11f7: vaddsd %xmm9,%xmm0,%xmm0
0x000000011a2f11fc: vaddsd %xmm8,%xmm0,%xmm0
0x000000011a2f1201: mov    %r11d,%r8d
0x000000011a2f1204: add    $0x9,%r8d
0x000000011a2f1208: vcvtsi2sd %r8d,%xmm1,%xmm1
0x000000011a2f120d: vaddsd %xmm0,%xmm1,%xmm0
0x000000011a2f1211: vaddsd %xmm7,%xmm0,%xmm0
0x000000011a2f1215: vaddsd %xmm6,%xmm0,%xmm0
0x000000011a2f1219: vaddsd %xmm5,%xmm0,%xmm0
0x000000011a2f121d: vaddsd %xmm0,%xmm3,%xmm0
0x000000011a2f1221: vaddsd %xmm4,%xmm0,%xmm0  ;*dadd
                                              ; - TT3::D@18 (line 27)

0x000000011a2f1225: vmovsd %xmm0,0x70(%rdx)   ;*putstatic s3
                                              ; - TT3::D@19 (line 27)

0x000000011a2f122a: vaddsd %xmm0,%xmm14,%xmm1  ;*dadd
                                              ; - TT3::D@18 (line 27)

0x000000011a2f122e: vmovsd %xmm1,0x70(%rdx)   ;*putstatic s3
                                              ; - TT3::D@19 (line 27)

0x000000011a2f1233: add    $0x10,%r11d        ;*iinc
                                              ; - TT3::D@22 (line 26)

0x000000011a2f1237: vcvtsi2sd %r11d,%xmm0,%xmm0  ;*i2d
                                              ; - TT3::D@17 (line 27)

0x000000011a2f123c: cmp    %ecx,%r11d
0x000000011a2f123f: jl     0x000000011a2f1130  ;*if_icmpge
                                              ; - TT3::D@10 (line 26)

上面就是 s3 += i 版的 D() C2优化后的代码，即循环确实存在的。

GCC会怎么优化呢？实际上，C还是有点不同的地方：C这里优化后耗时为常量。

以上文C代码为例：

begin = clock();
for(int i = 0; i < loopcnts; i++)
  res += i;
end = clock();
printf("asm(%f) sumd = %ld, /t/tcost: %f/n", MM, res, (double)(end-begin));

比如在gcc -O3 时，当 res为long类型时[JAVA版long跟double都不会有这个优化]，耗时为恒定的时间约4ns，因为此时GCC编译器优化为:

leal  -1(%r15), %eax
leal  -2(%r15), %ebx
imulq %rax, %rbx
shrq  %rbx

即 n*M/2 直接得到结果了，而不用再跑完循环了。

但是当 res为double类型时，却是需要跑for循环了，即：

LBB0_17:                                ## =>This Inner Loop Header: Depth=1
  xorps %xmm0, %xmm0
  cvtsi2sdl %ecx, %xmm0
  leal  1(%rcx), %esi
  xorps %xmm1, %xmm1
  cvtsi2sdl %esi, %xmm1
  addsd %xmm3, %xmm0
  leal  2(%rcx), %esi
  xorps %xmm2, %xmm2
  cvtsi2sdl %esi, %xmm2
  addsd %xmm0, %xmm1
  leal  3(%rcx), %esi
  xorps %xmm0, %xmm0
  cvtsi2sdl %esi, %xmm0
  addsd %xmm1, %xmm2
  leal  4(%rcx), %esi
  xorps %xmm1, %xmm1
  cvtsi2sdl %esi, %xmm1
  addsd %xmm2, %xmm0
  leal  5(%rcx), %esi
  xorps %xmm2, %xmm2
  cvtsi2sdl %esi, %xmm2
  addsd %xmm0, %xmm1
  leal  6(%rcx), %esi
  xorps %xmm0, %xmm0
  cvtsi2sdl %esi, %xmm0
  addsd %xmm1, %xmm2
  leal  7(%rcx), %esi
  xorps %xmm3, %xmm3
  cvtsi2sdl %esi, %xmm3
  addsd %xmm2, %xmm0
  addsd %xmm0, %xmm3
  addl  $8, %ecxa
  cmpl  %ecx, %edx
  jne LBB0_17
## %bb.6:
  testl %eax, %eax
  je  LBB0_9

这里插一句，如果你怀疑上文测试是否足够准确，那可以验证下这段代码：

begin = clock();
for(int i = 0; i < loopcnts; i++)
  res += i;
end = clock();
printf("asm(%f) sumd = %f, /t/tcost: %f/n", MM, res, (double)(end-begin));
res=0;
begin = clock();
for(int i = 0; i < loopcnts; i++)
  res += sqrt(MM+i);
end = clock();
printf("asm(%f) sums = %f, /t/tcost: %f/n", MM, res, (double)(end-begin));

在我的Mac机器上前者耗时约为后者1/5以内，足以论证之前C的测试结果，说明sqrt指令比手写的卡马克开根号性能提升是大于40%。

java -XX:TieredStopAtLevel=4 TwoSqrt

使用上面指令可以指定JIT优化级别

使用下面命令看lebel级别：

java -XX:+PrintFlagsFinal -version | grep CompileThreshold

level 0 - interpreter

level 1 - C1 with full optimization (no profiling)

level 2 - C1 with invocation and backedge counters

level 3 - C1 with full profiling (level 2 + MDO)

level 4 - C2

java -XX:TieredStopAtLevel=1 -XX:+UnlockDiagnosticVMOptions -XX:CompileCommand=’print, DL.D ‘ DL

使用上面命令可以输出指定方法JIT后的ASM代码

如何查看java的JIT信息，网上可自行搜索教程。

好了，让我们切回刚才的JAVA代码。

再次对TwoSqrt.java做个改动，添加了DL，即跟DD一样只是double换成了 long，运行几次，我们可看到跑出来的数据：

DD :6.666951101979061e+11 176899.189
DI :6.666951101979061e+11 672762.318
DL :6.666951101979061e+11 211739.138
DIJ:6.666951101979061e+11 292870.325
II :6.666695084696677e+11 228337.677
------------
DD :6.666951101979061e+11 171914.117
DI :6.666951101979061e+11 732866.289
DL :6.666951101979061e+11 220012.218
DIJ:6.666951101979061e+11 405955.404
II :6.666695084696677e+11 223352.986
------------
DD :6.666951101979061e+11 168237.792
DI :6.666951101979061e+11 254137.578
DL :6.666951101979061e+11 450223.373
DIJ:6.666951101979061e+11 265043.075
II :6.666695084696677e+11 255281.007
------------
DD :6.666951101979061e+11 169216.409
DI :6.666951101979061e+11 258624.100
DL :6.666951101979061e+11 443366.215
DIJ:6.666951101979061e+11 265902.974
II :6.666695084696677e+11 255865.587
....

这个结果和我用JMH测试下来接近:JMHBnhSqrt.java

Benchmark                Mode  Cnt          Score           Error  Units
JMHBnhSqrt.benchLoopDD   avgt    3  170126363.572 ±  18897173.028  ns/op
JMHBnhSqrt.benchLoopDI   avgt    3  304366448.848 ±  10242902.502  ns/op
JMHBnhSqrt.benchLoopDIJ  avgt    3  380274911.025 ± 311415636.755  ns/op
JMHBnhSqrt.benchLoopII   avgt    3  270527877.936 ±  68992367.132  ns/op
JMHBnhSqrt.sqrt          avgt    3          1.482 ±         0.059  ns/op

其中：

benchLoopDD: sumkk += Math.sqrt(MM + kk++); // kk为double，0开始，++

benchLoopDI: sumkk += Math.sqrt(MM + i); // i为for循环的i，0开始，++

benchLoopDIJ:sumjj += Math.sqrt(MM + jj++); // jj为int，0开始，++

benchLoopII: sumii += Math.sqrt(285 + ii++); // ii为int，0开始，++

1) 需要指出的是，这里其实benchmark的是for循环+sqrt函数的代价，而不仅仅是sqrt的代价了，所以如果你注释掉循环来用JMH benchmark的话再除以循环次数得到的值不一样了(因为for循环下编译优化的缘故)。

2) 似乎 JMH得到的 DI/DIJ性能结果和上面2次后的数据有点出入，可能是对JMH使用掌握的还不够深入，也可能是JMH本身原因，立个flag后面再写篇JMH的文章讨论此问题。

3) 总之应该相信JMH测试标准的准确性，但是不要过分迷信JMH的优化，尤其在自己不太确定的影响性能的配置上。

源码中DNO和D因编译优化不再细说了，我们看看几个：

1) Java版的DD的耗时已经非常接近之前C代码测出来的sqrt指令耗时–167888000纳秒，差距千分之一以内了。

2) 我们看到 DI 的耗时经历了从600多毫秒到稳定在260毫秒左右，DIJ也类似–这是C2的优化。

3) 尽管有上面的优化，但是耗时还是大于 DD 很多，30%多了，而他们区别只不过是被加数是double/int的区别，这是什么原因呢？

4) DL的耗时，比较奇怪，C1优化耗时(前两轮循环)还在 211739138纳秒左右，但是第三轮之后竟然到了 450223373纳秒。

上述问题3，在C版本是无此区别的，也就是说，无论如何解释优化的原因，这里都可看出来，Java的C2对DI的优化(int类型的循环体)不如 GCC 编译器，后者可以做到 DI和DD同样的耗时。

上述问题4，我们是不是可以得到结论，C2的优化效果(性能) 未必比 C1 好？

后两个问题，我在下一篇技术文里再讨论。

其他：

DNO/D的常量消除，上文写了，这里不重复，但看下java相关代码：

// 如对于：
double ss=0;
for(int i = 0; i < 1000000; i++)
    ss += sqrt(230000);

实际上，GCC -C1就能把 sqrt(MM) 优化成常量了：

_main:                                  ## @main
...
    movl    $1000000, %eax          ## imm = 0xF4240
    movsd   LCPI0_0(%rip), %xmm1    ## xmm1 = mem[0],zero
.p2align        4, 0x90
LBB0_1:                                 ## =>This Inner Loop Header: Depth=1
    addsd   %xmm1, %xmm0
    decl    %eax
    jne     LBB0_1

LCPI0_0(%rip)就是 double 479.58315233127195

相比之下，JAVA的JVM code和C1都不会把sqrt(230000)优化成常量，C2则会优化。如java 版 C1优化的循环体内调用的还是下面这句：

0x0000000110f7c229: vsqrtsd %xmm1,%xmm1,%xmm1  ;*invokestatic sqrt ; - TwoSqrt::D@19 (line 46)

java 版 C2优化的循环体内调用的还是下面这句：

L0002: vaddsd -0x138(%rip),%xmm0,%xmm0  # 0x0000000110f84480 ;*dadd ; - TwoSqrt::D@22 (line 46) ;
0x0000000110f845b8: vmovsd %xmm0,0x70(%r10)  ;*putstatic s3 ; - TwoSqrt::D@23 (line 46)
0x0000000110f845be: inc %ecx  ;*iinc ; - TwoSqrt::D@26 (line 45)
0x0000000110f845c0: cmp %r8d,%ecx
0x0000000110f845c3: jl L0002  ;*if_icmpge ; - TwoSqrt::D@10 (line 45)

就是说循环的相加根号后的数值已经不存在sqrt指令了。即，C2优化和GCC达到同样效果了。

上述也可见，DD和DI耗时在优化前差距还是很明显的，C2优化后差距减少了些。

此外，对于for循环，编译器/Java也有优化，比如for循环步长由 1 变为16(add $0x10,%ecx)为一个批次，减少了跳转指令的使用。

JAVA C2 针对DD优化采用了和GCC一样的思路，代码在下面可以看到，不过似乎有点区别是JAVA这里16个寄存器都用上了，而C用了四个，看起来不如C优化的那么紧凑。

这里贴一下C2之后代码，看看能不能发现什么。

L0000: vaddsd -0xd8(%rip),%xmm1,%xmm2  # 0x000000010feed180
                                                    ; {section_word}
0x000000010feed258: vaddsd -0xe0(%rip),%xmm2,%xmm1  # 0x000000010feed180
                                                    ;*dadd
                                                    ; - TwoSqrt::DD@26 (line 64)
                                                    ; {section_word}
0x000000010feed260: vaddsd -0xe0(%rip),%xmm2,%xmm2  # 0x000000010feed188
                                                    ; {section_word}
0x000000010feed268: vaddsd -0xe8(%rip),%xmm1,%xmm5  # 0x000000010feed188
                                                    ; {section_word}
0x000000010feed270: vsqrtsd %xmm2,%xmm3,%xmm3
0x000000010feed274: vsqrtsd %xmm5,%xmm2,%xmm2  ;*invokestatic sqrt
                                               ; - TwoSqrt::DD@29 (line 64)
L0001: vaddsd %xmm4,%xmm0,%xmm0
0x000000010feed27c: vaddsd -0x104(%rip),%xmm1,%xmm1  # 0x000000010feed180
0x000000010feed284: vaddsd %xmm0,%xmm3,%xmm0
0x000000010feed288: vaddsd -0x108(%rip),%xmm1,%xmm3  # 0x000000010feed188
0x000000010feed290: vaddsd %xmm2,%xmm0,%xmm4
0x000000010feed294: vaddsd -0x11c(%rip),%xmm1,%xmm0  # 0x000000010feed180
0x000000010feed29c: vaddsd -0x124(%rip),%xmm0,%xmm1  # 0x000000010feed180
0x000000010feed2a4: vaddsd -0x124(%rip),%xmm0,%xmm2  # 0x000000010feed188
0x000000010feed2ac: vaddsd -0x12c(%rip),%xmm1,%xmm5  # 0x000000010feed188
0x000000010feed2b4: vaddsd -0x13c(%rip),%xmm1,%xmm0  # 0x000000010feed180
0x000000010feed2bc: vaddsd -0x144(%rip),%xmm0,%xmm1  # 0x000000010feed180
0x000000010feed2c4: vaddsd -0x144(%rip),%xmm0,%xmm6  # 0x000000010feed188
0x000000010feed2cc: vaddsd -0x14c(%rip),%xmm1,%xmm7  # 0x000000010feed188
0x000000010feed2d4: vaddsd -0x15c(%rip),%xmm1,%xmm0  # 0x000000010feed180
0x000000010feed2dc: vaddsd -0x164(%rip),%xmm0,%xmm1  # 0x000000010feed180
0x000000010feed2e4: vaddsd -0x164(%rip),%xmm0,%xmm8  # 0x000000010feed188
0x000000010feed2ec: vaddsd -0x16c(%rip),%xmm1,%xmm9  # 0x000000010feed188
0x000000010feed2f4: vaddsd -0x17c(%rip),%xmm1,%xmm0  # 0x000000010feed180
0x000000010feed2fc: vaddsd -0x184(%rip),%xmm0,%xmm1  # 0x000000010feed180
0x000000010feed304: vaddsd -0x184(%rip),%xmm0,%xmm10  # 0x000000010feed188
0x000000010feed30c: vaddsd -0x18c(%rip),%xmm1,%xmm11  # 0x000000010feed188
0x000000010feed314: vaddsd -0x19c(%rip),%xmm1,%xmm0  # 0x000000010feed180
0x000000010feed31c: vaddsd -0x1a4(%rip),%xmm0,%xmm1  # 0x000000010feed180
0x000000010feed324: vaddsd -0x1a4(%rip),%xmm0,%xmm12  # 0x000000010feed188
0x000000010feed32c: vaddsd -0x1ac(%rip),%xmm1,%xmm13  # 0x000000010feed188
0x000000010feed334: vaddsd -0x1bc(%rip),%xmm1,%xmm0  # 0x000000010feed180
0x000000010feed33c: vaddsd -0x1c4(%rip),%xmm0,%xmm1  # 0x000000010feed180
0x000000010feed344: vaddsd -0x1c4(%rip),%xmm0,%xmm0  # 0x000000010feed188
0x000000010feed34c: vaddsd -0x1cc(%rip),%xmm1,%xmm14  # 0x000000010feed188
0x000000010feed354: vaddsd -0x1dc(%rip),%xmm1,%xmm1  # 0x000000010feed180
                                                     ;*dadd
                                                     ; - TwoSqrt::DD@26 (line 64)
                                                     ; {section_word}
0x000000010feed35c: vaddsd -0x1dc(%rip),%xmm1,%xmm15  # 0x000000010feed188
                                                      ; {section_word}
0x000000010feed364: vsqrtsd %xmm3,%xmm3,%xmm3
0x000000010feed368: vaddsd %xmm4,%xmm3,%xmm3
0x000000010feed36c: vsqrtsd %xmm15,%xmm4,%xmm4  ;*invokestatic sqrt
                                                ; - TwoSqrt::DD@29 (line 64)
0x000000010feed371: vsqrtsd %xmm14,%xmm14,%xmm14
0x000000010feed376: vsqrtsd %xmm0,%xmm15,%xmm15
0x000000010feed37a: vsqrtsd %xmm13,%xmm13,%xmm13
0x000000010feed37f: vsqrtsd %xmm12,%xmm12,%xmm12
0x000000010feed384: vsqrtsd %xmm11,%xmm11,%xmm11
0x000000010feed389: vsqrtsd %xmm10,%xmm10,%xmm10
0x000000010feed38e: vsqrtsd %xmm9,%xmm9,%xmm9
0x000000010feed393: vsqrtsd %xmm8,%xmm8,%xmm8
0x000000010feed398: vsqrtsd %xmm7,%xmm7,%xmm7
0x000000010feed39c: vsqrtsd %xmm6,%xmm6,%xmm6
0x000000010feed3a0: vsqrtsd %xmm5,%xmm5,%xmm5
0x000000010feed3a4: vsqrtsd %xmm2,%xmm0,%xmm0
0x000000010feed3a8: vaddsd %xmm3,%xmm0,%xmm0
0x000000010feed3ac: vaddsd %xmm0,%xmm5,%xmm0
0x000000010feed3b0: vaddsd %xmm0,%xmm6,%xmm0
0x000000010feed3b4: vaddsd %xmm0,%xmm7,%xmm0
0x000000010feed3b8: vaddsd %xmm0,%xmm8,%xmm0
0x000000010feed3bc: vaddsd %xmm9,%xmm0,%xmm0
0x000000010feed3c1: vaddsd %xmm0,%xmm10,%xmm0
0x000000010feed3c5: vaddsd %xmm0,%xmm11,%xmm0
0x000000010feed3c9: vaddsd %xmm0,%xmm12,%xmm0
0x000000010feed3cd: vaddsd %xmm0,%xmm13,%xmm0
0x000000010feed3d1: vaddsd %xmm0,%xmm15,%xmm0  ;*dadd
                                               ; - TwoSqrt::DD@32 (line 64)
0x000000010feed3d5: vmovsd %xmm0,0x70(%r10)  ;*putstatic s3
                                             ; - TwoSqrt::DD@33 (line 64)
0x000000010feed3db: vaddsd %xmm0,%xmm14,%xmm0  ;*dadd
                                               ; - TwoSqrt::DD@32 (line 64)
0x000000010feed3df: vmovsd %xmm0,0x70(%r10)  ;*putstatic s3
                                             ; - TwoSqrt::DD@33 (line 64)
0x000000010feed3e5: add $0x10,%ecx  ;*iinc
                                    ; - TwoSqrt::DD@36 (line 63)
0x000000010feed3e8: cmp %r9d,%ecx
0x000000010feed3eb: jl L0000  ;*if_icmpge
                              ; - TwoSqrt::DD@14 (line 63)
L0002: cmp %r8d,%ecx
0x000000010feed3f4: jge L0004
0x000000010feed3f6: xchg %ax,%ax  ;*getstatic s3
                                  ; - TwoSqrt::DD@17 (line 64)

其它

硬件：MacBook Pro 2017

JAVA Version:

java version “1.8.0_221”

Java(TM) SE Runtime Environment (build 1.8.0_221-b11)

Java HotSpot(TM) 64-Bit Server VM (build 25.221-b11, mixed mode)

outro1：篇幅限制，后面两个问题留到下一篇技术文里解答了。

outro2： 8月初开始准备写，中间因太忙，放下一段时间，可能内容不是连贯的，欢迎反馈。

Ref

牛顿迭代法求开方根
如何通俗易懂地讲解牛顿迭代法
Best-Square-Root-Method
RSQRTSS — Compute Reciprocal of Square Root of Scalar Single-Precision Floating-Point Value
performance implications of using vmovups and vmovapd
Intel指令周期
To my father. [1955-2019.08.23]

原文 http://thomaslau.xyz/2019/09/07/2019-09-07-on_carmac_and_java_jit/

正文到此结束