Округление результата при делении чисел с плавающей точкой

Question

Известно, что все целые числа из отрезка [-2**53, 2**53] могут быть точно представлены числом IEEE-754 двойной точности. Если мы складываем, вычитаем или перемножаем два целых числа IEEE-754 двойной точности из указанного выше отрезка и точный вещественный результат также принадлежит указанному выше отрезку, то получившееся значение будет совпадать с точным вещественным результатом.

Меня интересует операция деления. Пусть a и b — числа IEEE-754 двойной точности, оба строго положительные целые и принадлежат указанному выше отрезку. Пусть также точный вещественный результат деления a / b лежит между двумя целыми числами n-1 и n.

Будет ли результат деления a / b при вычислении в числах с плавающей точкой также строго меньше n? Если да, хотелось бы увидеть доказательство этого факта (можно ссылку на сайт/книгу, где можно почитать об этом). Если нет, то хотелось бы увидеть пример на котором в результате округления значение n достигается.

Если числа a и/или b не принадлежат отрезку [-2**53, 2**53], то число n может быть достигнуто при вычислениях в числах с плавающей точкой. Например, пусть

a = 45'035'996'273'704'952; // (2**53 - 2) * 5 + 2;
b = 9'007'199'254'740'991;  //  2**53 - 1;

Точный вещественный результат равен примерно (wolframalpha):

4.999999999999999666933092612453...

В данном случае точный вещественный результат строго меньше n = 5. Но при вычислении в числах IEEE-754 двойной точности число n будет достигнуто. Пример:

#include <iostream>
#include <iomanip>
#include <limits>
#include <cstdint>
using std::cout;
 
int main()
{
    typedef std::uint_fast64_t uint_type;
 
    double a, b, c;
    uint_type ta, tb, tc, ic;
 
    a = 45'035'996'273'704'952ULL; // (2**53 - 2) * 5 + 2;
    b = 9'007'199'254'740'991ULL;     //  2**53 - 1;
    c = a / b;
 
    ta = static_cast<uint_type>(a);
    tb = static_cast<uint_type>(b);
    tc = static_cast<uint_type>(c);
    ic = ta / tb;
 
    cout << std::setprecision(std::numeric_limits<double>::max_digits10) << std::scientific << std::boolalpha;
 
    cout << "a  == "  << a  << "\n";
    cout << "ta ==  " << ta << "\n\n";
 
    cout << "b  == "  << b  << "\n";
    cout << "tb ==  " << tb << "\n\n";
 
    cout << "c  = "   << c  << "\n";
    cout << "tc ==  " << tc << "\n";
    cout << "ic ==  " << ic << "\n";
 
    cout << "tc == ic: " << (tc == ic) << "\n\n";
 
    return 0;
}

Вывод:

a  == 4.50359962737049520e+16
ta ==  45035996273704952

b  == 9.00719925474099100e+15
tb ==  9007199254740991

c  = 5.00000000000000000e+00
tc ==  5
ic ==  4
tc == ic: false

Answer 1

IEEE-754 определяет несколько вариантов округления (для C/C++ смотри описания функций fegetround()/fesetround()), но по умолчанию истинный результат должен округлятся к ближайшему значению, представимому в данном формате (а в случае равноудалённости таковых — к нечётному значению мантиссы).

Будет ли результат деления a / b при вычислении в числах с плавающей точкой также строго меньше n?

Судя по всему так и есть, хотя доказательства и не мой конёк, но я всё же попробую кое-как строго это обосновать.

Итак, нам дано:

$a, b \in \mathbb{N} \quad a,b < 2^{t+1}$ , где — размерность мантиссы числа.
$n-1 < {{a} \over {b}} < n \quad n \in \mathbb{N}$
— функция деления чисел с плавающей точкой
Число с плавающей точкой представимо в виде $(1+{m \over 2^t}) \cdot 2^p, m \in \mathbb{Z}, 0 \le m < 2^t$
⇒ Интервал между числом и следующим за ним представимым числом $\varepsilon = 2^{p-t}$

Предположим, что существуют a и b удовлетворяющие условиям выше и такие, что .

Тогда, существует точное разложение вида:
$c={a \over b}=n - {k \over b}$
, где $k \in \mathbb{N}$ .

При этом $(1+{m \over 2^t}) \cdot 2^p \le c \le (1+{m+1 \over 2^t}) \cdot 2^p$ .
Заметим, что для нашего случая $p \ge 0$ и ⇒ $c \ge 2^p$ .

По условию из стандарта округление в большую сторону должно происходить только если истинный результат отличается от итогового меньше, чем на половину интервала представимых значений.

$div(a,b) - c = n - n + {k \over b} < {1 \over 2} \cdot \varepsilon$

Далее пара нехитрых алгебраических преобраpований:

${k \over b} < 2^{p-t-1}$

$b > k \cdot 2^{t+1-p}$

т.к. $k \in \mathbb{N}$ , то $b > 2^{t+1-p}$

$b \cdot 2^{p} > 2^{t+1}$

т.к. $c \ge 2^p$ , то $bc > 2^{t+1}$

т.о. $a > 2^{t+1}$

Таким образом, чтобы округлилось вверх до целого необходимым условием является, чтобы делимое было больше, чем $2^{t+1}$ . Что, очевидно, противоречит начальному условию.

Дисклеймер: В доказательстве могут быть ошибки. Я не написал пока ни одной монографии по программированию и не готов платить за обнаружение оных двоичными долларами; но исправления и уточнения приветствуются.

БЛОГ НА HUSL

Округление результата при делении чисел с плавающей точкой

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