Управление порядком вычисления выражений с плавающими точкой (ассоциативность)

Встретился, тут такой параноидальный код (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, на скобки надежды нет?


Ответы (0 шт):