[C/C++]C语言中math.h和cmath的pow()精度问题

参考链接： C++ pow()

帮小朋友们DEBUG的时候，他们有个题无论怎么提交OJ都不给过。 

我回来后想了想，估计是因为math.h库返回值转int时精度丢失的问题。 

>测试代码 

#include <stdio.h>

#include <math.h>

//MinGW GCC 4.7.2 32-bit Release

int main(){

    printf("math.h - double pow(double, double) 精度测试\n");

    int a=3;

    printf("%d\n",(int)pow(5,3));//1.输出125

    printf("%d\n",(int)pow(5,a));//2.输出124 这里丢精度了,结合下面的[3],我估计最后的结果是float->int 124.999999999999

    printf("%d\n",(int)round(pow(5,a)));//3.输出125 在[2]的基础上补上round()四舍五入函数，结果正常 

    printf("%lf\n",pow(5,a));//4.输出125.000000 显然，如果不转型成int，结果是没问题的

    return 0;

} 

>pow的精度问题研究 

math.h库里，pow函数是基于浮点运算的。 

试着找了一圈，没有找到源码，只在一些犄角旮旯里看到有人提到是在x86指令机上利用log和exp运算求出来的。也就是： 

double pow(double a, double b){

    //我这里简化了这个过程 大概写个伪码

    return x86_exp(a * x86_log(b));// 也就是：e^(a*ln(b))

}stackoverflow上也有关于这个问题的描述： 

https://stackoverflow.com/questions/14104711/what-algorithm-is-used-to-pow-function-in-c 

这个函数的具体实现过程则不得而知了。 

上面那段代码玄学的地方是这两行： 

    int a=3;

    printf("%d\n",(int)pow(5,3));//输出125

    printf("%d\n",(int)pow(5,a));//输出124

 这就是那种看上去都没错，但是一运行发现特么居然不对的神坑。 

你要说是幂参数a转到double b时转类型除了问题导致值不一样吧...那为啥单独给个数字3就没有任何计算问题。 

所以我想，要不然测试一下常量？ 

    const int b = 3;

    printf("%d\n",(int)pow(5,b));//输出125 没错，这次是125了orz 

  Emmmmmmmmm...const int 和 int 还有这种区别啊。难不成分配的长度不一样才导致转浮点时产生了什么误差？好吧我去看看汇编： 

看着好像没啥区别，就是存的地址不一样。那好吧，已知一个int四个字节，我多定义几个试试： 

int a=3,a2=3,a3=3,a4=3;

const int b=3,b2=3,b3=3,b4=3; 

于是： 

 每次增加4，至少在内存大小分配上const int 和 int 真的没啥区别。 

那好吧，我们再对比一下这几句的汇编: 

    int a=3;

    printf("%d\n",(int)pow(5,3));//输出125 I

    printf("%d\n",(int)pow(5,a));//输出124 II

    const int b = 3;

    printf("%d\n",(int)pow(5,b));//输出125 III 

int a的汇编是: 

   0x004013be <+14>:    movl   $0x3,0x2c(%esp) 

 首先是I:printf("%d\n",(int)pow(5,3)); 

   0x004013c6 <+22>:    movl   $0x7d,0x4(%esp)

   0x004013ce <+30>:    movl   $0x403064,(%esp)

   0x004013d5 <+37>:    call   0x401c80 <printf>然后是II:printf("%d\n",(int)pow(5,a));

   0x004013da <+42>:    fildl  0x2c(%esp)

   0x004013de <+46>:    fstpl  0x8(%esp)

   0x004013e2 <+50>:    mov    $0x0,%eax

   0x004013e7 <+55>:    mov    $0x40140000,%edx

   0x004013ec <+60>:    mov    %eax,(%esp)

   0x004013ef <+63>:    mov    %edx,0x4(%esp)

   0x004013f3 <+67>:    call   0x401c88 <pow>//注意这里是调用了pow函数

   0x004013f8 <+72>:    fnstcw 0x1e(%esp)// 开始转返回值（其实这个就是一般操作）

   0x004013fc <+76>:    mov    0x1e(%esp),%ax

   0x00401401 <+81>:    mov    $0xc,%ah

   0x00401403 <+83>:    mov    %ax,0x1c(%esp)

   0x00401408 <+88>:    fldcw  0x1c(%esp)

   0x0040140c <+92>:    fistpl 0x18(%esp)

   0x00401410 <+96>:    fldcw  0x1e(%esp)

   0x00401414 <+100>:    mov    0x18(%esp),%eax//最后存到0x18xxx去了

   0x00401418 <+104>:    mov    %eax,0x4(%esp)

   0x0040141c <+108>:    movl   $0x403064,(%esp)

   0x00401423 <+115>:    call   0x401c80 <printf> 

嚯，真有够长的，再看看III吧：printf("%d\n",(int)pow(5,b)); 

   0x00401428 <+120>:    movl   $0x3,0x28(%esp)//这行做的是const int b = 3;

   0x00401430 <+128>:    movl   $0x7d,0x4(%esp)//接下来是printf

   0x00401438 <+136>:    movl   $0x403064,(%esp)

   0x0040143f <+143>:    call   0x401c80 <printf>III和I看起来一模一样，也难怪输出的都是125 了。

 那么问题又来了，为啥I和III都没有call <pow>。假如我把II注释掉，再去看汇编，是这样的: 

=> 0x004013be <+14>:    movl   $0x3,0x1c(%esp) // I

   0x004013c6 <+22>:    movl   $0x7d,0x4(%esp)

   0x004013ce <+30>:    movl   $0x403064,(%esp)

   0x004013d5 <+37>:    call   0x401c30 <printf>

=> 0x004013da <+42>:    movl   $0x3,0x18(%esp) // III

   0x004013e2 <+50>:    movl   $0x7d,0x4(%esp)

   0x004013ea <+58>:    movl   $0x403064,(%esp)

   0x004013f1 <+65>:    call   0x401c30 <printf> 

是的，依旧没有调用pow。为啥？？？这结果是怎么运算出来的？ 

本着科学的精神，我做了这个测试： 

果然依旧没有调用pow。好吧，先放过这个问题...毕竟我的专精不在C的编译和汇编上，也许是有什么我尙不了解的知识点我还没了解到，改天去问问写C的底层大佬。 

还是回归正题，我们去考虑一下II中调用了pow进行幂运算的误差问题，毕竟I II III中，只有调用了pow的II输出了124这个错误的值。出于好奇，我去看了一下机器实现浮点运算的方法。另外，对于pow原生的参数类型double做以下测试： 

    double s=5.000,sb=3.0;

    printf("%d",(int)pow(s,sb));有： 

   0x004013be <+14>:    mov    $0x0,%eax

   0x004013c3 <+19>:    mov    $0x40140000,%edx

   0x004013c8 <+24>:    mov    %eax,0x28(%esp)

   0x004013cc <+28>:    mov    %edx,0x2c(%esp) //double s=5.0

   0x004013d0 <+32>:    mov    $0x0,%eax

   0x004013d5 <+37>:    mov    $0x40080000,%edx

   0x004013da <+42>:    mov    %eax,0x20(%esp)

   0x004013de <+46>:    mov    %edx,0x24(%esp) //double sb=3.0

   0x004013e2 <+50>:    mov    0x20(%esp),%eax //s

   0x004013e6 <+54>:    mov    0x24(%esp),%edx

   0x004013ea <+58>:    mov    %eax,0x8(%esp)

   0x004013ee <+62>:    mov    %edx,0xc(%esp)

   0x004013f2 <+66>:    mov    0x28(%esp),%eax //sb

   0x004013f6 <+70>:    mov    0x2c(%esp),%edx

   0x004013fa <+74>:    mov    %eax,(%esp)

   0x004013fd <+77>:    mov    %edx,0x4(%esp)

   0x00401401 <+81>:    call   0x401c70 <pow>

   0x00401406 <+86>:    fnstcw 0x1e(%esp)

   0x0040140a <+90>:    mov    0x1e(%esp),%ax

   0x0040140f <+95>:    mov    $0xc,%ah

   0x00401411 <+97>:    mov    %ax,0x1c(%esp)

   0x00401416 <+102>:    fldcw  0x1c(%esp)

   0x0040141a <+106>:    fistpl 0x18(%esp)

   0x0040141e <+110>:    fldcw  0x1e(%esp)

   0x00401422 <+114>:    mov    0x18(%esp),%eax

   0x00401426 <+118>:    mov    %eax,0x4(%esp)

   0x0040142a <+122>:    movl   $0x403064,(%esp)

   0x00401431 <+129>:    call   0x401c78 <printf>对比一下，基本可以确定就是传参int a的时候的问题。传参的不同，在经过pow内部的运算会产生不同的结果。由于现在pow是一个黑盒（不知道源码），我们只好动手算算按int传参之后的后果了。我们假设stackoverflow上给出的pow内部运算方法是对的，按照IEEE754对单双精度的定义及刚刚stackoverflow里某人推测给出pow的运算方法: 

fld1            // 1.0

fld             // y

fld             // x 

fyl2x           // y * log2(x)

fadd            // add 1.0

f2xm1           // 2 ^ (x-1) 

 我大概算了一下，实际上在经过运算后原本的125这个逼近值最后被存为浮点时变成01000010111110011111111111111111(但愿我没算错Orz，总之大概思路对行了)，换句话说，124.999。