#include <cassert>

int main()
{
   double a = 0.014390783999999978;
   double b = 0.052996730719235739;
   double c = a - b; //-0.038605946719235763
   
   assert((b+c) == a); //b+c = 0.014390783999999976
}

Можно ли без округления вычислить c, чтобы гарантировать (b+c) == a? Спасибо.
http://coliru.stacked-crooked.com/a/eb63581b259878aa

long double

Здравствуйте, Kazmerchuk Pavel, Вы писали:

Забудь про «==« для плавающей точки. Только разницу с эпсилон сравнивай

KP>Можно ли без округления вычислить c, чтобы гарантировать (b+c) == a? Спасибо.
Нет. Даже (со)процессор делает округление.
Тип округления устанавливается в одном из регистров сопроцессора.
без округления выполняются операции со степенями двойки: 1/2 = 2^(-1), 1/4 = 2^(-2) и т.д.
Но и в этом случае выполняется нормализация, при которой младшие разряды результата могут исчезнуть.

Здравствуйте, LaptevVV, Вы писали:

KP>>Можно ли без округления вычислить c, чтобы гарантировать (b+c) == a? Спасибо.
LVV>Даже (со)процессор делает округление.

Да тут даже до процессора дело не всегда доходит.
Компилятор выкидывает ненужные вычисления и заменяет тело функции на прямой вызов assert_fail. Не нужны в машинном коде эти ваши сравнения, если и так всё понятно

Здравствуйте, LaptevVV, Вы писали:

LVV>Нет. Даже (со)процессор делает округление.
А вычислить дополнительную поправку к с, чтобы, например (b + c + correction) == a?

Здравствуйте, Kazmerchuk Pavel, Вы писали:

KP>

KP>   double a = 0.014390783999999978;
KP>   double b = 0.052996730719235739;
KP>   double c = a - b; //-0.038605946719235763
KP>

KP>Можно ли без округления вычислить c, чтобы гарантировать (b+c) == a? Спасибо.

Нет. double умеет не больше 15ти значащих десятичных знаков хранить, а у тебя — 17. А тут a и b сразу уже хранят не те значения, которые ты им задаёшь. По идее, компилятор должен предупредить об этом, но не факт

Здравствуйте, reversecode, Вы писали:

R>long double

Который в MSVC, внезапно, тот же double

Здравствуйте, Homunculus, Вы писали:

H>Забудь про «==« для плавающей точки. Только разницу с эпсилон сравнивай

У меня — можно — https://github.com/al-martyn1/marty_decimal

это или UB или не UB
а какое там представление он имеет в каком компилере до одного места
что то годболд msvс перестал у меня показывать
остальные gcc clang ok

Здравствуйте, Kazmerchuk Pavel, Вы писали:

KP>Можно ли без округления вычислить c, чтобы гарантировать (b+c) == a?
Разница двух double чисел не всегда представима в double без потери точности.
То есть для double ответ — нельзя.

Если необходимо считать точно подобные суммы, то тут поможет только длинная арифметика (с хранением всех значащих бит). Но понятно, что не всегда подходит и всё это довольно далеко от стандартного типа double.


А так, операции в IEEE754 c floating-point не ассоциативные. То есть в них принципиально не выполняется равенство (a + b) + c == a + (b + c).
Нужно с этим уметь жить и писать алгоритмы соответствующим образом.

Например, этот эффект можно учитывать и минимизировать как в алгоритме суммирования Кэхэна. Заметь, что там тоже повторяется похожий паттерн из действий c = a — b; t = b + c; почти как у тебя в примере.

Либо этот эффект можно даже эксплуатировать во благо, например, для быстрого округления чисел без преобразования типов:

double round_pos(double x) {
    const double magic = 0x1p+52;
    return x + magic - magic;
}

Демо: https://godbolt.org/z/d9qexTdWW



Ну и ещё можно упомянуть, что компилятору обычно можно сказать, что он должен считать, что floating-point ассоциативный (например, через ключ -ffast-math в gcc/clang). Тогда он будет генерировать машинный код исходя из предположений, что (a+b)-b можно заменить сначала на a+(b-b), и потом на просто a (да, случай с NaN тоже считается несуществующим). Но не советую идти этой дорогой: тут наоборот сломаются хорошие алгоритмы, которые написаны в рассчёте на гарантии и поведение ieee754. То есть они будут выдавать иногда чушь, но зато работать будет быстро

Здравствуйте, reversecode, Вы писали:

R>это или UB или не UB

С чего бы это UB?


R>а какое там представление он имеет в каком компилере до одного места

А это пока ты по граблям не походил

Здравствуйте, Kazmerchuk Pavel, Вы писали:

KP>

KP>{
KP>   double a = 0.014390783999999978;
KP>   double b = 0.052996730719235739;
KP>   double c = a - b; //-0.038605946719235763
   
KP>   assert((b+c) == a); //b+c = 0.014390783999999976
KP>}
KP>

можно вычислить ошибку (смотри D. E. Knuth, TAOCP, Vol 2, 4.2.2, Theorem B.):

#include <cassert>

double sum(double u, double v, double& t) {
    double s = u + v;
    double up = s - v;
    double vpp = s - up;
    up -= u;
    vpp -= v;
    t = -(up + vpp);

    return s;
}

int main()
{
    double a = 0.014390783999999978;
    double b = 0.052996730719235739;

    double err;
    const auto c = sum(a, -b, err);

    assert((b+c) + err == a); //b+c = 0.014390783999999976                                                                                                                    
}

Здравствуйте, Zhendos, Вы писали:

Z>Здравствуйте, Kazmerchuk Pavel, Вы писали:

KP>>

KP>>{
KP>>   double a = 0.014390783999999978;
KP>>   double b = 0.052996730719235739;
KP>>   double c = a - b; //-0.038605946719235763
   
KP>>   assert((b+c) == a); //b+c = 0.014390783999999976
KP>>}
KP>>

Z>можно вычислить ошибку (смотри D. E. Knuth, TAOCP, Vol 2, 4.2.2, Theorem B.):

Z>

Z>#include <cassert>

Z>double sum(double u, double v, double& t) {
Z>    double s = u + v;
Z>    double up = s - v;
Z>    double vpp = s - up;
Z>    up -= u;
Z>    vpp -= v;
Z>    t = -(up + vpp);

Z>    return s;
Z>}

Z>int main()
Z>{
Z>    double a = 0.014390783999999978;
Z>    double b = 0.052996730719235739;

Z>    double err;
Z>    const auto c = sum(a, -b, err);

Z>    assert((b+c) + err == a); //b+c = 0.014390783999999976                                                                                                                    
Z>}
Z>

Спасибо! То что нужно!

Здравствуйте, Kazmerchuk Pavel, Вы писали:

KP>Можно ли без округления вычислить c, чтобы гарантировать (b+c) == a? Спасибо.

Еще +0.5 копеек: в статье https://habr.com/ru/post/525090/ сравниваются алгоритмы суммирования Kahan и Rump–Ogita–Oishi

Здравствуйте, Zhendos, Вы писали:

Z>можно вычислить ошибку (смотри D. E. Knuth, TAOCP, Vol 2, 4.2.2, Theorem B.):

И даже это ничего не гарантирует

Здравствуйте, Kazmerchuk Pavel, Вы писали:

KP>Можно ли без округления вычислить c, чтобы гарантировать (b+c) == a? Спасибо.
Можно, но для начала тебе нужно открыть любой учебник по вычислениям на компьютерах и внимательно прочитать и всё поймешь, что можно, что нельзя. Раньше в школе основы оного преподавали. "Значашие цифры", "системы счисления", "периодические и не периодические дроби", "целые, рациональные и иррациональные числа" — это, например, из тех основ.

Здравствуйте, Homunculus, Вы писали:

H>Забудь про «==« для плавающей точки. Только разницу с эпсилон сравнивай
Чушь. Всё зависит от того, что именно ты хочешь сделать и где.

Здравствуйте, Marty, Вы писали:

M>Нет. double умеет не больше 15ти значащих десятичных знаков хранить, а у тебя — 17. А тут a и b сразу уже хранят не те значения, которые ты им задаёшь. По идее, компилятор должен предупредить об этом, но не факт
А еще в советской школе рассказывали про погрешности округления и что с ними происходит при вычислениях.

Здравствуйте, watchmaker, Вы писали:

W>Ну и ещё можно упомянуть, что компилятору обычно можно сказать, что он должен считать, что floating-point ассоциативный (например, через ключ -ffast-math в gcc/clang). Тогда он будет генерировать машинный код исходя из предположений, что (a+b)-b можно заменить сначала на a+(b-b), и потом на просто a (да, случай с NaN тоже считается несуществующим). Но не советую идти этой дорогой: тут наоборот сломаются хорошие алгоритмы, которые написаны в рассчёте на гарантии и поведение ieee754. То есть они будут выдавать иногда чушь, но зато работать будет быстро
Это нехорошие алгоритмы, хорошие обычно делают устойчивыми и им пофиг на твои ключики компилятору.