Взаимозаменяемость операций сложения и умножения с плавающей запятой IEEE 754

Взаимозаменяемо ли сложение x + x на умножение 2 * x в стандарте IEEE 754 (IEC 559) с плавающей запятой или, в более общем смысле, есть ли гарантия, что case_add и case_mul всегда дают точно тот же результат?

#include <limits>

template <typename T>
T case_add(T x, size_t n)
{
    static_assert(std::numeric_limits<T>::is_iec559, "invalid type");

    T result(x);

    for (size_t i = 1; i < n; ++i)
    {
        result += x;
    }

    return result;
}

template <typename T>
T case_mul(T x, size_t n)
{
    static_assert(std::numeric_limits<T>::is_iec559, "invalid type");

    return x * static_cast<T>(n);
}

c++ floating-point ieee-754 numerical-stability
person plasmacel    schedule 04.10.2016    source источник
comment
Обратите внимание, что, кажется, существует много способов суммирования n * x, но, на удивление, многие из них эквивалентны! Это как-то связано с stackoverflow.com/questions/21690585/is-3xx-always-exact и stackoverflow.com/questions/21676955/   -  person aka.nice    schedule 08.10.2016

Ответы (4)


arrow_upward
9
arrow_downward

Взаимозаменяемо ли сложение x + x умножением 2 * x в стандарте IEEE 754 (IEC 559) с плавающей запятой?

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

или, в более общем смысле, есть ли гарантия, что case_add и case_mul всегда дают точно такой же результат?

В общем, нет. Насколько я могу судить, похоже, что это сохраняется в течение n <= 5:

  • n=3: поскольку x+x является точным (т.е. не требует округления), поэтому (x+x)+x включает только одно округление на последнем этапе.

  • n=4 (и вы используете режим округления по умолчанию), тогда

    • if the last bit of x is 0, then x+x+x is exact, and so the results are equal by the same argument as n=3.
    • если последние 2 бита равны 01, то точное значение x+x+x будет иметь последние 2 бита 1|1 (где | указывает последний бит в формате), которые будут округлены до 0|0. Следующее добавление даст точный результат |01, поэтому результат будет округлен в меньшую сторону, аннулируя предыдущую ошибку.
    • если последние 2 бита равны 11, то точное значение x+x+x будет иметь последние 2 бита 0|1, которые будут округлены до 0|0. Следующее добавление даст точный результат |11, поэтому результат будет округлен в большую сторону, снова отменяя предыдущую ошибку.

  • n=5 (опять же, предполагая округление по умолчанию): поскольку x+x+x+x является точным, он сохраняется по той же причине, что и n=3.

Для n=6 это не удается, например принять x как 1.0000000000000002 (следующий double после 1.0), и в этом случае 6x равно 6.000000000000002, а x+x+x+x+x+x 6.000000000000001

person Simon Byrne    schedule 04.10.2016
comment
Ничего себе, ваше доказательство намного короче, чем мое доказательство путем анализа случаев, когда 3 * x находится в том же бинаде, что и 2 * x или 4 * x. - person Pascal Cuoq; 04.10.2016
comment
Из того, что я могу сказать на основе исчерпывающих тестов с IEEE-754 binary32 (= float), умножение на n идентично повторному сложению (n-1), когда n ≤ 5 для FE_TONEAREST, n ≤ 3 для FE_TOWARDZERO, n ≤ 3 для FE_DOWNWARD и n ≤ 3 для FE_UPWARD. Можно ли распространить доказательство на режимы направленного округления? - person njuffa; 04.10.2016
comment
@njuffa Аргумент для n ‹= 3 не зависел от выбора режима округления. Доказательство стало зависеть от режима округления на шаге n = 4. - person Patricia Shanahan; 04.10.2016
comment
@PatriciaShanahan Спасибо, что указали на это, я упустил из виду, что ответ уже решает проблему режима округления. - person njuffa; 04.10.2016
comment
Было бы здорово, если бы в ответ были включены случаи, когда режимы округления не используются по умолчанию, упомянутые @njuffa. - person plasmacel; 07.10.2016

arrow_upward
3
arrow_downward

Если n, например, pow(2, 54), тогда умножение будет работать нормально, но в пути сложения, когда значение результата будет достаточно большим, чем входное x, result += x даст result.

person Mark B    schedule 04.10.2016

arrow_upward
1
arrow_downward

Да, но в целом этого не происходит. Умножение на число больше 2 может не дать тех же результатов, так как вы изменили показатель степени и может немного упасть, если вы замените его сложением. Однако умножение на два не может немного упасть, если его заменить операциями сложения.

person Malcolm McLean    schedule 04.10.2016

arrow_upward
1
arrow_downward

Если аккумулятор result в case_add становится слишком большим, добавление x приведет к ошибкам округления. В какой-то момент добавление x вообще не даст никакого эффекта. Таким образом, функции не дадут такой же результат.

Например, если double x = 0x1.0000000000001p0 (шестнадцатеричная запись с плавающей запятой):

n  case_add              case_mul

1  0x1.0000000000001p+0  0x1.0000000000001p+0
2  0x1.0000000000001p+1  0x1.0000000000001p+1
3  0x1.8000000000002p+1  0x1.8000000000002p+1
4  0x1.0000000000001p+2  0x1.0000000000001p+2
5  0x1.4000000000001p+2  0x1.4000000000001p+2
6  0x1.8000000000001p+2  0x1.8000000000002p+2
person nwellnhof    schedule 04.10.2016