Управление порядком вычисления выражений с плавающими точкой (ассоциативность)
Встретился, тут такой параноидальный код (scipy.special (cephes)):
/* Computes fl(a+b) and err(a+b). */
XSF_HOST_DEVICE inline double two_sum(double a, double b, double *err) {
volatile double s = a + b;
volatile double c = s - a;
volatile double d = b - c;
volatile double e = s - c;
*err = (a - e) + d;
return s;
}
Т.е. гарантировано корректно работает с любыми ключами компиляции, но столь же гарантировано медленно, ибо стандарты C/C++ требуют записи/чтения памяти.
Что собственно немного противоречит его изначальному предназначению: быстрая реализация четверной точности. Нет, конечно, программную реализацию IEEE bynary128 (_Float128 и float128_t) он и в таком виде обгоняет, но...
В библиотеке qd откуда предыдущий код был родом, он выглядел так:
/* Computes fl(a+b) and err(a+b). */
inline double two_sum(double a, double b, double &err) {
double s = a + b;
double bb = s - a;
err = (a - (s - bb)) + (b - bb);
return s;
}
Тут понятно, шустро, но при применении широко известного ключа -ffast-math или менее известного -fassociative-math тут же возвращаемый результат err становится тождественным 0.. На этот счёт Дэвид Бейли написал в документации: "не применяйте такие ключи".
Что немного обидно, пишешь рядом static_assert() и assert():
static inline constexpr
bool check_two_sum(double a, double b,
double _s, double _err) {
double err = 0.;
double s = two_sum(a, b, err);
return s == _s && err == _err;
}
...
static_assert(check_two_sum(1e-35, 1., 1., 1e-35));
assert(check_two_sum(1e-35, 1., 1., 1e-35));
static_assert() всегда проходит, а assert() срабатывает (хотя иногда нестабильно, вроде как не на всех компиляторах).
Немного более обидно, что если gcc выставляет два макроса: __FAST_MATH__ и/или __ASSOCIATIVE_MATH__, то clang выставляет только один __FAST_MATH__, что не позволяет детектировать проблему вызванную ключом -fassociative-math.
Что ещё более обидно:
- У
clangсоответствующая_Pragma(), исключающая только этот вид оптимизации, действует на охватывающий блок, т.е. её надо вставлять внутрь функции, и это единственный вариант; - У
gccесть два варианта: классический push_options/pop_options, и вариант с атрибутом функции, т.е. её надо вставлять перед описанием функции; - У
MSVCесть только классический вариант push/pop.
У icc (Intel), вроде как у gcc и т.п. В общем, из-за clang нет единой схемы.
Пока непонятно, как для разных компиляторов, но хотя бы в едином стиле детектировать и/или обходить проблему?
И основной вопрос. Правильно ли я понимаю, что, по крайней мере, в современных спецификациях C23/C++23 (C23, 6.5 Expressions, абзац 8, FP_CONTRACT) разрешается, в том или ином виде сокращение плавающих операций в рамках одного выражения, но не между операторами? Т.е. корректный по C23/C++23 код (в данном случае, корректность строго определена: a + b == s + err && s == s ⨁ err, где + - точная математическая сумма, а ⨁ - машинная сумма, Деккер, или Кнут, т.2, теорема 4.2.2A), должен быть таким:
/* Computes fl(a+b) and err(a+b). */
inline double two_sum(double a, double b, double *err) {
auto s = a + b;
auto c = s - a;
auto d = b - c;
auto e = s - c;
*err = (a - e) + d;
return s;
}
Т.е. сейчас в C/C++ не как в Fortran, на скобки надежды нет?