Как сравнивать вещественные числа

[vak]

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

Возникает проблема: как проверять правильность вычислений? Можно игнорировать некоторое небольшое количество младших бит мантиссы. Вот пример кода.

bool nearly_equal(double a, double b, int ignore_nbits)
{
    if (a == b)
        return true;

    // Разбиваем числа на мантиссу и экспоненту.
    int a_exponent, b_exponent;
    double a_mantissa = frexp(a, &a_exponent);
    double b_mantissa = frexp(b, &b_exponent);

    // Экспоненты обязаны совпасть.
    if (a_exponent != b_exponent)
        return false;

    // Вычитаем мантиссы, образуем положительную дельту.
    double delta = fabs(a_mantissa - b_mantissa);

    // Определяем порог допустимой разницы.
    double limit = ldexp(1.0, ignore_nbits - 52);

    return delta < limit;
}

Запустим на нескольких тестовых значениях.

void try(double a, double b, int ignore_nbits)
{
    bool eq = nearly_equal(a, b, ignore_nbits);
    printf("%a ~= %a ? %s (ignore %d bits)\n", a, b, eq ? "True" : "False", ignore_nbits);
}

int main()
{
    try(1.000000001, 1.000000002, 22);
    try(1.0000000001, 1.0000000002, 18);
    try(1.000000000000001, 1.000000000000002, 2);
}

Получаем:

0x1.000000044b83p+0 ~= 0x1.000000089705fp+0 ? True (ignore 22 bits)
0x1.000000006df38p+0 ~= 0x1.00000000dbe7p+0 ? True (ignore 18 bits)
0x1.0000000000005p+0 ~= 0x1.0000000000009p+0 ? True (ignore 2 bits)

Вообще всем, имеющим дело с плавающей точкой, рекомендую функции frexp() и ldexp(). Они имеются на всех платформах, работают везде одинаково и очень эффективно. Буквально несколько машинных инструкций. Часто помогают в подобных случаях.

Comments

[ircicq]

Почему экспоненты обязаны совпадать?
Например у двоичных:
1.00000000001
0.11111111111

не совпадают, но мы хотим чтобы nearly_equal() позвращала true?

[vak]

Точно, надо будет учесть.

[juan_gandhi]

Здрасьте. А если try(0.9999999999999999, 1.000000000001, 7) ?

C экспонентой тщатильнее надо. (Программировал я когда-то FPP simulator).

[vak]

Да, с экспонентой я поторопился. Надо будет поправить.

[juan_gandhi]

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

[vak]

Да, и будет достаточно денормализовать всего на +/-1 порядок. Если разница больше, то понятно, что числа не равны.

[ircicq]

Около нуля это не так:

0.00001
0.000000001
разница экспонент: 4

[vak]

Как раз случай, когда числа заведомо не равны. Разница в несколько порядков это криминал.

[ircicq]

Контрпример:

A= 0.00010
B= 0.00001
C= 0.000000001

A и B различаются на 1 порядок, то есть при каком то значении ignore_nbits принимаются равными.

B и C - на 4 порядка

в то же время B ближе к C чем к A.

[vak]

Если значения различаются в разы, то понятно, что нет смысла говорить о количестве несовпавших битов мантиссы. ВСЕ биты не совпали.

Тут логично бы было ввести термин "относительная погрешность", но мне лениво.

[vit_r]

Сколько инструкций стоит обычное вычитание?

[vak]

Недорого, но как оно поможет?

[ichthuss]

Например, вычесть и посмотреть, насколько изменилась экспонента.

[vit_r]

Я бы просто из меньшего числа вычел бы большее (сравнение + вычитание), а потом прибавил бы дельту и посмотрел бы, получился ли положительный результат (вычитание + сравнение с нулём).

[dimorlus]

Ну а исправленная версия будет? А то я, когда надо, сравниваю как некий процент от большего числа, что, конечно, изрядно накладно.