c - float 比 double 更精确？

Question

有一个序列:

x(n+2)=9/4*x(n+1)-1/2*x(n)

x(1)=1/3,x(2)=1/12

确切的结果是x(n)=4^(1-n)/3

我想在计算中显示 x(60) 的舍入误差。

我的代码是

#include <stdio.h>
#include <math.h>
int main(void)
{
    float x[60];
    x[0] = 1./3;
    x[1] = 1./12;
    for (int i = 2; i < 60; i++) {
        x[i] = 9./4*x[i-1]-1./2*x[i-2];
    }
    double y[60];
    y[0] = 1./3;
    y[1] = 1./12;
    for (int i = 2; i < 60; i++) {
        y[i] = 9./4*y[i-1]-1./2*y[i-2];
    }
    printf("single:%g, double:%g, exact:%g\n", x[59], y[59], pow(4,-59)/3);
    return 0;
}

我用 gcc 编译它:

gcc seq.c

输出是:

single:1.00309e-36, double:1.71429, exact:1.00309e-36

如果我像这样更改上面的代码:

#include <stdio.h>
#include <math.h>
int main(void)
{
    float x[60];
    x[0] = 1.f/3;
    x[1] = 1.f/12;
    for (int i = 2; i < 60; i++) {
        x[i] = 9.f/4*x[i-1]-1.f/2*x[i-2];
    }
    double y[60];
    y[0] = 1./3;
    y[1] = 1./12;
    for (int i = 2; i < 60; i++) {
        y[i] = 9./4*y[i-1]-1./2*y[i-2];
    }
    printf("single:%g, double:%g, exact:%g\n", x[59], y[59], pow(4,-59)/3);
    return 0;
}

其中'f'被添加在用于计算x-array的常量 float 之后。

输出看起来很正常:

single:-9.2035e+08, double:1.71429, exact:1.00309e-36

我的问题是:

为什么第一种情况float数据类型的结果等于exact结果？

编译器做什么？

Answer 1

float 并不比 double 更精确，您的 float 计算没有给您 pow(4 ,-59)/3。

发生的事情是，您的循环旨在采用微小的舍入误差并在每次迭代中放大它。在精确的数学中，每个值应该恰好是前一个值的四分之一，但如果由于舍入误差而不是四分之一，则差异在每一步都会被放大。

由于可表示值的四分之一始终是可表示的(直到遇到次正规数和下溢问题)，递归具有一个额外的属性:如果计算的精度足以超过结果的精度存储，然后将结果四舍五入到较低的存储精度将四舍五入到前一个值的四分之一。 (9./4 和 1./2 因子的选择为递归提供了此属性的更强版本，其中结果恰好是旧值的四分之一甚至在舍入存储之前。)

---

对于 double ，使用您正在使用的编译器和编译器设置，舍入误差会发生并被放大。使用 float ，计算以 double 执行，消除了由于上述属性而导致的递归步骤中的舍入误差，因此没有什么可以放大的。如果以长 double 执行 double 计算，也会发生同样的事情。

---

让我们通过使用 %a 格式说明符以十六进制 表示法打印 float 来仔细查看生成的确切值。看起来像0x1.5555555555558p-6，其中0x和p之间的部分是十六进制数，之后的部分p 是一个十进制数，表示与十六进制数相乘的 2 的幂。这里，0x1.5555555555558p-6 表示 0x1.5555555555558 乘以 2^-6。 %a 格式始终打印 float 或 double 的精确值，这与 %g 不同，它是四舍五入的。

我们还将展示第三种计算，将结果存储为 double ，但以长 double 进行数学计算。

我们修改后的程序如下所示:

#include <stdio.h>
#include <math.h>
int main(void)
{
    float x[60];
    x[0] = 1./3;
    x[1] = 1./12;
    for (int i = 2; i < 60; i++) {
        x[i] = 9./4*x[i-1]-1./2*x[i-2];
    }
    double y[60];
    y[0] = 1./3;
    y[1] = 1./12;
    for (int i = 2; i < 60; i++) {
        y[i] = 9./4*y[i-1]-1./2*y[i-2];
    }
    double z[60];
    z[0] = 1./3;
    z[1] = 1./12;
    for (int i = 2; i < 60; i++) {
        z[i] = (long double) 9./4*z[i-1] - (long double) 1./2*z[i-2];
    }
    printf("float:%a, double:%a, double2:%a, formula:%a\n", x[59], y[59], z[59], pow(4,-59)/3);
    for (int i = 0; i < 60; i++) {
        printf("%d %a %a %a\n", i, x[i], y[i], z[i]);
    }
    return 0;
}

这是输出。我本来打算删节这个，但事实证明，如果不掩盖模式中有趣的部分，很难做到这一点:

float:0x1.555556p-120, double:0x1.b6db6db6db6dap+0, double2:0x1.5555555555555p-120, formula:0x1.5555555555555p-120
0 0x1.555556p-2 0x1.5555555555555p-2 0x1.5555555555555p-2
1 0x1.555556p-4 0x1.5555555555555p-4 0x1.5555555555555p-4
2 0x1.555556p-6 0x1.5555555555558p-6 0x1.5555555555555p-6
3 0x1.555556p-8 0x1.555555555557p-8 0x1.5555555555555p-8
4 0x1.555556p-10 0x1.555555555563p-10 0x1.5555555555555p-10
5 0x1.555556p-12 0x1.5555555555c3p-12 0x1.5555555555555p-12
6 0x1.555556p-14 0x1.5555555558c3p-14 0x1.5555555555555p-14
7 0x1.555556p-16 0x1.5555555570c3p-16 0x1.5555555555555p-16
8 0x1.555556p-18 0x1.5555555630c3p-18 0x1.5555555555555p-18
9 0x1.555556p-20 0x1.5555555c30c3p-20 0x1.5555555555555p-20
10 0x1.555556p-22 0x1.5555558c30c3p-22 0x1.5555555555555p-22
11 0x1.555556p-24 0x1.5555570c30c3p-24 0x1.5555555555555p-24
12 0x1.555556p-26 0x1.5555630c30c3p-26 0x1.5555555555555p-26
13 0x1.555556p-28 0x1.5555c30c30c3p-28 0x1.5555555555555p-28
14 0x1.555556p-30 0x1.5558c30c30c3p-30 0x1.5555555555555p-30
15 0x1.555556p-32 0x1.5570c30c30c3p-32 0x1.5555555555555p-32
16 0x1.555556p-34 0x1.5630c30c30c3p-34 0x1.5555555555555p-34
17 0x1.555556p-36 0x1.5c30c30c30c3p-36 0x1.5555555555555p-36
18 0x1.555556p-38 0x1.8c30c30c30c3p-38 0x1.5555555555555p-38
19 0x1.555556p-40 0x1.8618618618618p-39 0x1.5555555555555p-40
20 0x1.555556p-42 0x1.e186186186186p-39 0x1.5555555555555p-42
21 0x1.555556p-44 0x1.bc30c30c30c3p-38 0x1.5555555555555p-44
22 0x1.555556p-46 0x1.b786186186185p-37 0x1.5555555555555p-46
23 0x1.555556p-48 0x1.b6f0c30c30c3p-36 0x1.5555555555555p-48
24 0x1.555556p-50 0x1.b6de186186185p-35 0x1.5555555555555p-50
25 0x1.555556p-52 0x1.b6dbc30c30c3p-34 0x1.5555555555555p-52
26 0x1.555556p-54 0x1.b6db786186185p-33 0x1.5555555555555p-54
27 0x1.555556p-56 0x1.b6db6f0c30c3p-32 0x1.5555555555555p-56
28 0x1.555556p-58 0x1.b6db6de186185p-31 0x1.5555555555555p-58
29 0x1.555556p-60 0x1.b6db6dbc30c3p-30 0x1.5555555555555p-60
30 0x1.555556p-62 0x1.b6db6db786185p-29 0x1.5555555555555p-62
31 0x1.555556p-64 0x1.b6db6db6f0c3p-28 0x1.5555555555555p-64
32 0x1.555556p-66 0x1.b6db6db6de185p-27 0x1.5555555555555p-66
33 0x1.555556p-68 0x1.b6db6db6dbc3p-26 0x1.5555555555555p-68
34 0x1.555556p-70 0x1.b6db6db6db785p-25 0x1.5555555555555p-70
35 0x1.555556p-72 0x1.b6db6db6db6fp-24 0x1.5555555555555p-72
36 0x1.555556p-74 0x1.b6db6db6db6ddp-23 0x1.5555555555555p-74
37 0x1.555556p-76 0x1.b6db6db6db6dbp-22 0x1.5555555555555p-76
38 0x1.555556p-78 0x1.b6db6db6db6dap-21 0x1.5555555555555p-78
39 0x1.555556p-80 0x1.b6db6db6db6dap-20 0x1.5555555555555p-80
40 0x1.555556p-82 0x1.b6db6db6db6dap-19 0x1.5555555555555p-82
41 0x1.555556p-84 0x1.b6db6db6db6dap-18 0x1.5555555555555p-84
42 0x1.555556p-86 0x1.b6db6db6db6dap-17 0x1.5555555555555p-86
43 0x1.555556p-88 0x1.b6db6db6db6dap-16 0x1.5555555555555p-88
44 0x1.555556p-90 0x1.b6db6db6db6dap-15 0x1.5555555555555p-90
45 0x1.555556p-92 0x1.b6db6db6db6dap-14 0x1.5555555555555p-92
46 0x1.555556p-94 0x1.b6db6db6db6dap-13 0x1.5555555555555p-94
47 0x1.555556p-96 0x1.b6db6db6db6dap-12 0x1.5555555555555p-96
48 0x1.555556p-98 0x1.b6db6db6db6dap-11 0x1.5555555555555p-98
49 0x1.555556p-100 0x1.b6db6db6db6dap-10 0x1.5555555555555p-100
50 0x1.555556p-102 0x1.b6db6db6db6dap-9 0x1.5555555555555p-102
51 0x1.555556p-104 0x1.b6db6db6db6dap-8 0x1.5555555555555p-104
52 0x1.555556p-106 0x1.b6db6db6db6dap-7 0x1.5555555555555p-106
53 0x1.555556p-108 0x1.b6db6db6db6dap-6 0x1.5555555555555p-108
54 0x1.555556p-110 0x1.b6db6db6db6dap-5 0x1.5555555555555p-110
55 0x1.555556p-112 0x1.b6db6db6db6dap-4 0x1.5555555555555p-112
56 0x1.555556p-114 0x1.b6db6db6db6dap-3 0x1.5555555555555p-114
57 0x1.555556p-116 0x1.b6db6db6db6dap-2 0x1.5555555555555p-116
58 0x1.555556p-118 0x1.b6db6db6db6dap-1 0x1.5555555555555p-118
59 0x1.555556p-120 0x1.b6db6db6db6dap+0 0x1.5555555555555p-120

在这里，我们首先看到 float 计算没有产生 pow 公式给出的精确值(它没有足够的精度)，但是它非常接近，差异被 %g 的四舍五入所隐藏。我们还看到，float 值每次都以正好 4 的倍数递减，更改后的double 计算的值也是如此。原始 double 版本的 double 值开始几乎这样做，然后一旦放大的错误超过计算量就会发散。这些值最终开始增加 2 倍，而不是减少 4 倍。