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

[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(). Они имеются на всех платформах, работают везде одинаково и очень эффективно. Буквально несколько машинных инструкций. Часто помогают в подобных случаях.