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

来源：互联网发布：haskell linux 编辑：程序博客网时间：2024/05/22 15:34

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

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

>测试代码

#include <stdio.h>#include <math.h>//MinGW GCC 4.7.2 32-bit Releaseint 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.0fld             // yfld             // x fyl2x           // y * log2(x)fadd            // add 1.0f2xm1           // 2 ^ (x-1)

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

阅读全文

1 0