Пример с плавающей точкой по Rump

1 min


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

Этот пример Зигфрида Румпа просит нас оценить

333,75 Y6 + Икс2 (11 Икс2Y2Y6 – 121 Y4 – 2) + 5,5 Y8 + Икс/ (2Y)

в Икс = 77617 и Y = 33096

Здесь мы оцениваем пример Rump с одинарной, двойной и четырехкратной точностью.

#include 
#include 

using namespace std;

template  T f(T x, T y) {
    T x2 = x*x;
    T y2 = y*y;
    T y4 = y2*y2;
    T y6 = y2*y4;
    T y8 = y2*y6;
    return 333.75*y6 + x2*(11*x2*y2 - y6 - 121*y4 - 2) + 5.5*y8 + x/(2*y);
}

int main() {  

    int x = 77617;
    int y = 33096;
    
    cout << f(x, y) << endl;
    cout << f(x, y) << endl;
    cout << (double) f<__float128>(x, y) << endl;

}

Это дает три ответа,

-1.85901e+030
-1.18059e+021
1.1726

ни один из которых не является отдаленно точным. Точный ответ -54767/66192 = -0,827 ...

Python дает тот же результат, что и C ++, с двойной точностью, что неудивительно, поскольку числа с плавающей запятой в Python являются C двойными числами.

Откуда проблема? Обычный подозреваемый: вычитание почти равных чисел. Давайте разделим выражение Румпа на две части.

s = 333,75 Y6 + Икс2(11Икс2Y2 - Y6 - 121Y4 - 2)

t = 5,5Y8 + Икс/ (2Y)

и посмотрите на них отдельно. Мы получили

s = -7917111340668961361101134701524942850.00         
t =  7917111340668961361101134701524942849.1726...

Ценности -s а также T согласен с 36 цифрами, но в четвёртой точности есть только 34 цифры (1). Вычитание этих двух чисел приводит к катастрофической потере точности.

Похожие посты

(1) Число с четверной точностью имеет 128 битов: 1 знаковый бит, 15 для показателей степени и 112 для дробной части. Поскольку ведущий ноль является неявным, это дает точность в 113 битов. С лога10(2113) = 34.01…, это означает, что квад имеет 34 десятичных знака точности.


0 Comments

Ваш e-mail не будет опубликован. Обязательные поля помечены *