pow()在这里似乎有一个偏差。

Question

pow()在这里似乎有一个偏差。

cfloating-pointfloating-point-precisionpow

42

这里发生了什么：

#include <stdio.h>
#include <math.h>
int main(void) {
    printf("17^12 = %lf\n", pow(17, 12));
    printf("17^13 = %lf\n", pow(17, 13));
    printf("17^14 = %lf\n", pow(17, 14));
}

我得到了这个输出：

17^12 = 582622237229761.000000
17^13 = 9904578032905936.000000
17^14 = 168377826559400928.000000

13和14与Wolfram Alpha不匹配，参见:

12: 582622237229761.000000
    582622237229761

13: 9904578032905936.000000
    9904578032905937

14: 168377826559400928.000000
    168377826559400929

此外，它不是由于某种奇怪的分数而出错 - 它确切地仅出错了一个！

如果这是因为我已经达到了pow()能为我做的限制，那么是否有替代方法可以计算这个呢？我需要一个函数来计算x^y，其中x^y始终小于ULLONG_MAX。

- jsj

44

也许你需要了解《计算机科学家应该了解的浮点数算术》一文吗？该文链接为：http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html。 - Some programmer dude

9

有趣的事实：这个循环不会结束：for (float f = 0; f < INT_MAX; f++) { } - ntoskrnl

4

%lf 中的 l 是无效的。你正在打印 double，如果编译平台将 double 映射到 IEEE 754 的双精度格式，则不会以这种方式打印 9904578032905937.0，因为不存在这样的双精度数字。 - Pascal Cuoq

5

@ntoskrnl说，这个程序_might_（可能）不会停止。在Cray和int是16位的系统上，它确实会停止。 - MSalters

4

使用“double”类型，你总是可以依靠前15个小数位（有效数字）。对于17的12次方，这是小数点前的所有数字。对于17的13次方，你不能相信最后一位数字。没有格式的情况下，“double”通常会被写成“9.90457803290594E+15”。而17的14次方使用该表示法为“1.68377826559401E+17”。 - Jeppe Stig Nielsen

显示剩余6条评论

4个回答

14

如果您的输入参数是非负整数，则可以实现自己的pow。

递归方式：

unsigned long long pow(unsigned long long x,unsigned int y)
{
    if (y == 0)
        return 1;
    if (y == 1)
        return x;
    return pow(x,y/2)*pow(x,y-y/2);
}

迭代地：

unsigned long long pow(unsigned long long x,unsigned int y)
{
    unsigned long long res = 1;
    while (y--)
        res *= x;
    return res;
}

有效率地：

unsigned long long pow(unsigned long long x,unsigned int y)
{
    unsigned long long res = 1;
    while (y > 0)
    {
        if (y & 1)
            res *= x;
        y >>= 1;
        x *= x;
    }
    return res;
}

- barak manos

9

递归版本无用，因为它不会重复使用答案。只需进行一次递归调用即可：root = y ? pow(x, y / 2) : 1; return root * root * (y&1 ? x : 1); - tom

@chux：刚刚注意到缺少了 if (y == 1)...我简直不敢相信它在那么长时间里一直存在而我却没有注意到。非常感谢！！！！ :) - barak manos

@tom 你说得对。错过了“最终乘以pow(0,1)”这一步。 - chux - Reinstate Monica

14

你得到的数字过于庞大，无法精确地用 double 表示。一个双精度浮点数基本上具有53个有效二进制位，能够表示所有不超过 2^53 或 9,007,199,254,740,992 的整数。

对于更高的数字，最后几位会被截断，你的计算结果将四舍五入为下一个可以用 double 表示的数字。对于 17^13，它只是略高于限制，因此取最接近的偶数。对于大于 2^54 的数字，这是最接近的可以被四整除的数字，依此类推。

- M Oehm

2

啊..最接近的偶数有意义。 - jsj

2

在其他很好的答案中补充一点：在x86架构下通常有x87 80位扩展格式可用，大多数C编译器通过long double类型支持该格式。该格式允许操作最高达2^64的整数数字而不产生间隙。

在<math.h>中有一个类似于pow()的函数，专门用于操作long double数值- powl()。还应该注意到，long double数值的格式说明符与double数值的格式说明符不同- %Lf。因此，使用long double类型的正确程序如下：

#include <stdio.h>
#include <math.h>
int main(void) {
    printf("17^12 = %Lf\n", powl(17, 12));
    printf("17^13 = %Lf\n", powl(17, 13));
    printf("17^14 = %Lf\n", powl(17, 14));
}

正如Stephen Canon在评论中指出的那样，这个程序并不能保证给出精确的结果。

- Constructor

2

(a) long double 的长度修饰符是 L，而不是 ll。 (b) 不能保证 powl 会提供一个次于最后一位单位精度的准确结果（这对于给出“正确”答案是必要的），在某些平台上它并不能。 - Stephen Canon

@StephenCanon (a) llf 是非标准变体，这是我的错误，谢谢。 (b) 当然你是对的，但我没有写过这个程序会给出正确的结果。 - Constructor

我甚至指出powl在许多实现中会给出极其错误的结果。 - tmyklebu

@tmyklebu 即使CPU支持x87指令集？ - Constructor

是的，特别是如果它是x87。 - tmyklebu

1

谢谢你的建议。我有相当丰富的单元测试套件（这就是原始问题的来源），所以拥有另一个选项很好 - 如果它能通过我的测试，那么对于每一个可能的输入，它都能完美地工作，对吧？;) - jsj

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- user395760 · Accepted Answer

pow 使用 double 数字。这些数字表示形式为 s * 2^e，其中 s 是一个 53 位整数。因此，double 可以存储小于 2^53 的所有整数，但只有一些大于 2^53 的整数。特别地，它只能表示大于 2^53 的偶数，因为对于 e > 0，该值总是 2 的倍数。

17^13 需要 54 位来精确表示，因此将 e 设置为 1，因此计算出的值变成了偶数。正确的值是奇数，因此它比实际值少了一个。同样，17^14 需要 58 位来表示。它也巧合地少了一个（只要你不应用太多的数论），它恰好比 32 的倍数 少一个，这是该大小的 double 数字被舍入的粒度。

对于精确的整数幂运算，您应该一直使用整数。编写自己的无需double的幂运算程序。如果y很大，则使用平方幂运算，但我假设它始终小于64，因此这个问题不重要。