Как стать автором
Обновить

Рекуррентные формулы для расчета ошибок итерационного суммирования двоичных чисел ограниченной длины

Время на прочтение8 мин
Количество просмотров6.1K
Всего голосов 18: ↑15 и ↓3+12
Комментарии19

Комментарии 19

Ещё в тему:
https://habrahabr.ru/post/266023/

У вас неплохой научный стиль изложения материала, но, будьте милосердны, не пишите на нём научно-популярные и прикладные статьи! Вы же на хабре, а не в Duke Mathematical Journal, будьте проще!

Кроме того, вам остро не хватает сорцов, Описываемая задача очень хорошо и приятно оборачивается в любом ООП-языке, позволяя спрятать в недра код для оценки погрешностей. Если бы вы написали обёртки над float/double/long double, позволяющие оценить погрешности вычислений, вы бы принесли сообществу реальную пользу и смогли бы измерить накопление погрешностей в популярных продуктах, например, в пакетах работы с графикой, особенно, форматах сохранения с потерей информации. А так полезность ваших статей, не смотря на их занимательность, остаётся на достаточно низком уровне.

P.S. Для получения точного значения арифметических и не только операций желательно использовать хотя бы MathCad (или, как минимум, Wolfram Alpha), который умеет выводить результаты арифметических операций с точностью до нужной цифры. Тот же Excel спокойно накапливает погрешность, хотя в нём и есть двоично-десятичные поля и numeric-типы. Ну и удобнее MathCad для таких задач.
Предлагаю посотрудничать. Как я вижу, вы грамотный программист, а я больше математик.
Пас. Не примите на свой счёт, просто в сутках всего 27 часов, и мне их и так категорически не хватает на все мои хотелки.
Как мы видим, абсолютная погрешность представления оказывается ничтожно малой по сравнению с погрешностью округления при итерационном суммировании.

И все? Целый пост ради того, чтобы сообщить очевидную вещь?

Имелась ввиду суммарная погрешность представления по сравнению с погрешностями округления. А это уже не так очевидно.

Может, и не очевидно — но уж точно общеизвестно. Я встречал упоминание этого факта так часто, что не могу вспомнить кто сказал мне об этом впервые.

На качественном уровне эта проблема обсуждается во многих учебниках. Но количественных оценок я не встречал.
Всё-таки для большинства практических целей алгоритм Кэхэна достаточно минимизирует погрешность. Если же погрешность настолько критична, то мне кажется, что стоит отказаться от представления чисел в таком формате и переходить на рациональные дроби. Конечно теряется «аппаратное ускорение» (FPU, SSE/AVX, GPU), да и проблема сокращения дробей и представления чисел условно бесконечной длины не то чтобы совсем тривиальна с точки зрения производительности, зато проблема погрешности устраняется как класс.
Когда производительность не критична, то можно использовать любые пакеты. Но бывают случаи, когда производительность и точность становятся во главе угла. См., например, статью пользователя MagisterLudi про проблемы американской ЗРК «Патриот», которая пропустила советскую ракету «керосинку»

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

Проверьте математику.


Мои результаты меньше в шесть раз:


  • js: Δ=53.105953216552734; δ=9.891748879763943e-8
  • c#:Δ=52,7059531211853; δ=9,81724309923268E-08

код

js:


"use strict";

var s = 0.0, N = 5368712233, M = N/4;

for (var i=0; i<M; i++) {
    s += 0.1;
    s += 0.1;
    s += 0.1;
    s += 0.1;
}

for (var i=0; i<N%4; i++) {
    s += 0.1;
}

console.log(Math.abs(0.1*N - s), Math.abs(0.1*N - s) / (0.1*N));

C#:


            double s = 0.0;
            long N = 5368712233;

            for (long i = 0; i < N; i++)
            {
                s += 0.1;
            }

            Console.WriteLine($"Δ={Math.Abs(0.1 * N - s)} δ={Math.Abs(0.1 * N - s) / (0.1 * N)}");
Очевидно потому, что ваш процессор, либо накапливает сумму до 80 байт, а потом округляет до double, либо округляет после каждого шага итерации. У меня рассмотрен случай без округления, т.е. все цифры, которые не помещаются в разрядную сетку машинного слова просто отбрасываются (усекаются)

Пересчитал на С++.


Использовал /fp:strict и _controlfp_s(nullptr, _RC_CHOP, _MCW_RC).


Код сложения:


fld         qword ptr [s]  
fadd        qword ptr [a]
fstp        qword ptr [s]
mfence // или lock add [s], 0

Вроде бы в таком режиме процессор ну никак не должен переиспользовать значения из 80ти разрядных регистров (хотя я не уверен).


Получилось Δ=110.431599 δ=2.05695e-07, то есть в два раза больше чем раньше (и в три раза меньше чем у вас). Причем из всех факторов повлиял лишь _RC_CHOP.

UPD: что-то я не в ту сторону полез. Строгий режим же подразумевает воспроизводимость результатов на других устройствах — а значит, хитрых оптимизаций быть не может, можно было не ставить барьеров. Добавил еще _controlfp_s(nullptr, _PC_53, _MCW_PC) для надежности. Ничего не изменилось.

Так что у меня неправильно?

Результаты ваших вычислений у вас неправильные. В три раза больше чем в действительности.

Результаты вычислений существенно зависят от выбранного способа округления, от компилятора, оптимизации.
Скажите мне пожалуйста, какой из возможных вариантов округления производится в вашем алгоритме и на каком итерационном шаге? Тогда можно будет сравнивать результаты.

Я же писал. Использовалось усечение на каждом шаге.

Вы правы. У меня приведен пример не совсем для double. Мантисса дабла предполагает 53 значащих цифры, а в моем примере только 51. Видимо поэтому и расхождения.
Зарегистрируйтесь на Хабре, чтобы оставить комментарий

Публикации

Истории