Поплавок точнее двойного?

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

Я хочу показать ошибку округлениях (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;
}

, где после постоянного числа с плавающей точкой для вычисления x-массива добавляется 'f'.

Вывод кажетсяnormal:

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

Мой вопрос:

Почему результат типа данных с плавающей запятой равен точному результату в первой ситуации?

Что делает компилятор?

Answer 1

float не более точен, чем double, и ваши вычисления float не дали вам точного результата pow(4,-59)/3.

Что происходит, так это то, что ваше повторение рассчитано намаленькая ошибка округления и усиливать ее каждую итерацию.В точной математике, каждое значение должно составлять ровно одну четверть предыдущего значения, но если это не совсем четверть из-за ошибки округления, разница увеличивается на каждом шаге.

Поскольку четверть представимого значениявсегда представимый (до тех пор, пока вы не столкнетесь с ненормальными числами и проблемами с занижением), у повторения есть дополнительное свойство: если вычисления выполняются с точностью, достаточно превышающей точность, с которой сохраняются результаты, затем округляют результаты до более низкой точности для хранениябудет округлять ровно до четверти предыдущего значения.(Выбор коэффициентов 9./4 и 1./2 дает рекуррентности еще более сильную версию этого свойства, где результат составляет ровно четверть старого значения даже до округления для хранения.)

---

При двойных значениях, с используемыми настройками компилятора и компилятора, ошибка округления возникает и усиливается.В случае с плавающей точкой вычисления выполняются с двойной точностью, что исключает ошибку округления на этапах повторения благодаря описанным выше свойствам, поэтому нечего усиливать.Если бы вычисление для удвоений было выполнено с длинной двойной точностью, произошло бы то же самое.

---

Давайте внимательнее посмотрим на полученные точные значения, используя спецификатор формата %a для печатичисла с плавающей точкой в ​​ шестнадцатеричном обозначении.Это похоже на 0x1.5555555555558p-6, где часть между 0x и p является шестнадцатеричным числом, а часть после p является десятичным числом, представляющим степень двойки, на которую умножается шестнадцатеричное число.Здесь 0x1.5555555555558p-6 представляет 0x1,5555555555558 раз 2 ^ -6.Формат %a всегда печатает точное значение типа float или double, в отличие от округления %g.

Мы также покажем третье вычисление, сохраняя результаты как double, но выполняя математические операции для long 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 раза.

Answer 2

Мне кажется, вы хорошо знаете, что ошибки округления с плавающей запятой могут привести к совершенно неправильным результатам. На самом деле, кажется, что вас больше удивляет получение «правильного» результата в примере 1, чем вас удивляет неправильный результат в примере 2.

Что ж, ошибки округления могут привести к крайне неправильным результатам, но нельзя допустить, что ошибки округления всегда приведут к крайне неправильным результатам. Иногда ошибки округления вызовут только незначительную ошибку, а в других случаях ошибки округления приведут к нестабильности всего расчета и приведут к экстремальным ошибкам.

Ответ от @ user2357112 https://stackoverflow.com/a/55194247/4386427 дает хорошее описание вашего конкретного случая.

Но одна часть вашего вопроса все еще остается без ответа:

Что делает компилятор?

Я предполагаю, что вы спрашиваете, почему этот код

a) x[i] = 9./4*x[i-1]-1./2*x[i-2];

дает результаты, отличные от этого кода

b) x[i] = 9.f/4*x[i-1]-1.f/2*x[i-2];
            ^            ^

Ответ таков: случай a) требует, чтобы вычисления выполнялись с двойной точностью, поскольку 9. имеет тип double, тогда как случай b) позволяет выполнять вычисления с одинарной точностью, поскольку все типы являются числами с плавающей запятой.

Если компилятор решит использовать операции одинарной точности вместо двойных, ошибки округления будут разными в случае а) и случае б). Различные ошибки округления (могут) приводят к различным результатам, как обсуждалось выше.

Не существует единого ответа, объясняющего, что будет делать компилятор, поскольку разные компиляторы могут давать разные результаты. Ниже приведен один пример, сгенерированный с использованием https://godbolt.org/ и x86-64 gcc 8.3 и -O0.

Для простоты он охватывает только 9./4*x[i-1] против 9.f/4*x[i-1], и я скопировал только отличающиеся строки.

   Case 9./4*x[i-1]:

    movss   xmm0, DWORD PTR [rax]
    cvtss2sd        xmm1, xmm0
    movsd   xmm0, QWORD PTR .LC0[rip]
    mulsd   xmm0, xmm1

и

   Case 9.f/4*x[i-1]:

    movss   xmm1, DWORD PTR [rax]
    movss   xmm0, DWORD PTR .LC0[rip]
    mulss   xmm0, xmm1
    cvtss2sd        xmm0, xmm0

Как видно, разница заключается в использовании одинарной точности (mulss) и двойной точности (mulsd).

Итак, подытожим:

Что делает компилятор?

Он использует различную точность с плавающей точкой для расчета, что приводит к различным ошибкам округления, что снова приводит к разным результатам.

Answer 3

Это расчет, который не может быть выполнен с использованием чисел с плавающей запятой.Вы добавляете большие и маленькие числа, и ошибка округления слишком велика для такого рода вычислений.

А 1/3 и 1/12 - просто удачное начальное начало для вычисления с плавающей запятой.Для других начальных значений оба вычисления дают почти одинаковые результаты, и оба обычно неверны.