о сложении (программистское)

[avva]

Числа с плавающей точкой - вот еще одно в череде компьютерных задач, которыми я занимался мало и понимаю плохо. Сегодня я узнал, что числа с плавающей точкой надо складывать - если их много - не "наивно", а особенным образом. Мне неприятно, что я до сих пор этого не знал; или если знал, то забыл - кажется логичным, чтобы это упоминалось во время курса по численному анализу, который я брал на втором курсе много лет назад, но это был один из наиболее нелюбимых предметов, и после экзамена я с радостью забыл вообще все, что его касается.

Эта тема всплыла в одной внутренней рассылке на работе, где утверждалось, что средний выпускник CS-факультета этого не знает. Я уже давно не выпускник, но честно признаюсь: не знал. А вы, коллеги, знали?

Суть проблемы вот в чем. Предположим, вам надо суммировать миллион чисел с плавающей точкой, находящихся примерно в одном диапазоне - на расстоянии одного-двух порядков друг от друга. Тогда простой, как пробка, цикл от одного до миллиона с накоплением суммы - очень плохой способ суммировать эти числа. Дело в том, что вследствие того, что промежуточные суммы будут со временем намного больше среднего слагаемого, потери точности вследствие округления начнут играть заметную роль. Суммирование миллиона чисел таким образом легко может потерять 3-4 знака точности. Если, скажем, сложить таким образом миллион одинаковых чисел, а потом сравнить с результатом умножения исходного числа на миллион, это легко заметить.

В Перле:

DB<1> $a=0.1234567890123456
DB<2> $x = 1_000_000*$a;
DB<3> for (1..1_000_000) { $y += $a; }
DB<4> x $y;
0 123456.789009695
DB<5> x $x
0 123456.789012346
DB<6> x $y-$x;
0 '-2.65109702013433e-06'


Накопившаяся разница проявляется в шестом знаке после запятой, или, иными словами, четыре последних знака потеряны.

Что делать? В общем случае - помнить об ошибках накопления и учитывать их. В конкретном случае сложения большого количества примерно одинаковых чисел самым простым оказывается складывать их "равномерно". Вместо того, чтобы одна переменная копила сумму, сложим сначала пары чисел-соседей, и получим в два раза меньше промежуточных сумм. Их тоже сложим попарно, и так далее - log(N) шагов. Пример такого алгоритма можно увидеть на этой странице (написан на Питоне).

Comments

[vova-belkin.livejournal.com]

Это еще что! Числа с плавающей точкой и сравнивать нельзя
т.е.

(4.2 - 2.2)^2 не всегда равно (5.2-3.2)^2

Факт очевидный, но на практике всегда о нем зыбываешь.

[avva.livejournal.com]

Как раз вот этот факт я помню - и это, пожалуй, единственное, что помню о проблемах с плавающей точкой. Если бы я писал что-то нетривиальное обрабатывающее много таких чисел, до того, как прочитал об этой теме сегодня, то сравнивал бы я осторожно (сравнивая модуль разницы с подходящей маленькой константой), а вот складывал бы наивно, не задумываясь.

[mummy1.livejournal.com]

Я очень мало знаю о числах с плавающей точкой, и меня это сердит. И на численном анализе о них ничего не расказывали.

[spamsink.livejournal.com]

Прежде чем их складывать попарно, их еще и отсортировать полезно.

[eraplee.livejournal.com]

а можно чуть конкретней? что даст сортировка?

[eraplee.livejournal.com]

Мне кажется вопрос поставлен неправильно -- знал/не знал. То, что ошибка при арифметической операции будет (для очень маленьких чисел, не обязательно для всех с плавающей точкой) -- очевидно. И то, что она будет накапливаться очевидно тоже любому. Отсюда же можно легко сделать вывод, что наивный метод сложения даст довольно неточный результат. А вот как лучше сделать вместо наивного -- тут конечно уже не каждый сразу сообразит ;)

[plakhov.livejournal.com]

По-моему, проблема как раз обратная.
Очень легко написать цикл до миллиона не задумываясь.
Но как только тебе рассказали (или сам осознал), что в этом есть проблема, решение с равномерным сложением кажется почти очевидным. По крайней мере, это первое, что пришло в голову - еще до того, как я добежал глазами до конца записи.

[michk.livejournal.com]

Средний выпускник действительно не знает (т.е. может учит, но всё равно не знает). Но вот если с этим делом работаешь - быстро на такие вещи наталкиваешься. Самый простой пример - арифметическая прогрессия с маленьким шагом.

[e2k-4d-x-ussr.livejournal.com]

правильнее тогда выбирать 2 минимальных и их складвать? и так пока не останется только ОДИН. :)

[nine_k]

Imho, да.
Мой первоначальный вариант -- отсортировать и складывать от меньших к большим.
Если складывать сначала наименьшие пары, получится и сложений меньше, и потому потенциальная потеря точности будет ещё меньше.

[nine_k]

Да можно вообще оперировать числами произвольной точности. И даже местами абсолютной, применяя рациональные дроби.
Только медленно это и места много занимает.

[akater.livejournal.com]

Я вообще не выпускник факультета и не CS-ник. Однако я знаю об этом (читал в своё время древнюю, но всё равно полезную книгу под названием "Компьютерная алгебра", автора сейчас не вспомню), и считаю, что при моей специальности ("прикладная математика") не знать этого ещё хуже, чем при CS.

Более того, имхо, это проблема тех, кто преподаёт как раз математику, а не CS. Если курсы анализа и численных методов (в т.ч. компьютерных вычислений) не сопрягать, всё будет очень плохо.

[igorlord.livejournal.com]

"и так далее - log(N) шагов"

There are still N-1 steps. :)

It's just that the "outer loop" of the algorithm runs log(N) times.

[malaya-zemlya.livejournal.com]

Там еще сравнения...

[ygam.livejournal.com]

По-моему, это описывается в книжке по численному анализу в серии McGraw Hill International Series on Pure and Applied Mathematics; она 13 лет назад была в майкрософтовской корпоративной библиотеке.

[ygam.livejournal.com]

http://www.emis.de/journals/NNJ/DavIns-HM.html

Does anyone believe that the difference between the Lebesgue and Riemann integrals can have physical significance, and that whether, say, an airplane would or would not fly could depend on this difference? If such were claimed, I should not care to fly in that plane. --Richard W. Hamming

[squadette.livejournal.com]

WECSSKAFPA

[amarao-san.livejournal.com]

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

[vrml.livejournal.com]

На физтехе на численных методах рассказывали.

[http://users.livejournal.com/malfet_/]

+1
А кто был лектором, если не секрет? Не Чудов Л.А.?

[talash.livejournal.com]

Not a computer science graduate or student as of today, but the issue is not new to me. Programmers of the patriot missile defense system should have known better though (http://www.ima.umn.edu/~arnold/disasters/patriot.html)

[dr-tambowsky.livejournal.com]

Из чистого занудства, но что то там не то. Задача системы - отслеживать цель. Соответственно, смещение окна - это скорость цели, помноженная на (очень небольшое) время, прошедшее с момента последнего вычисления окна. Как они сами пишут: "The range gate's prediction of where the Scud will next appear is a function of the Scud's known velocity and the time of the last radar detection."

Насколько у системы отстали часы - совершенно не важно, важен только временной интервал. Как-то я не могу себе представить, чтобы ракету детектировали при запуске системы, а потом ждали 100 часов, чтобы вычислить новое положение окна - в этом случае да, 0.34 секунды, заработанные за 100 часов, сыграли бы свою роль... Там, видимо, другая проблема - с преобразованием: например, вместо того, чтобы вычислить разность показаний таймера в исходном точном представлении (целых числах), а потом конвертировать маленькую разность в плавающую точку, ребята сначала конвертировали близкие начальное и конечное показания, потом посчитали разность уже в плавающей арифметике. Когда система работает долго, и начальное и конечное показания - большие числа. Конвертировать и затем вычитать их на таком коротком регистре - действительно некрасиво, маленькая разность двух больших плавающих чисел утонет в погрешности. Из официального отчёта, вроде бы, следует, что так оно и было. Объяснение про 0.34 секунды - абсолютно неправильное.

Но так или иначе, история очень поучительная :)

[spamsink.livejournal.com]

Реальные примеры различия плавающих точек на разных архитектурах приводились в брошюре "Ошибки и ловушки при программировании на фортране" из серии "Библиотечка программиста". В качестве разных архитектур брались ЕС ЭВМ и БЭСМ-6. Первая проигрывала в длине мантиссы (24 против 40), а последняя - в диапазоне порядков (+/- 19 против +/- 38), так что нужно было заботиться и о потере значащих цифр - не только при сложении далеких чисел, но и при вычитании близких, и о пере-недо-полнениях в промежуточных результатах.

[dmpogo.livejournal.com]

+1

Dlja teh kto zanimaetsja chislennymi raschetami, skzazhem nauchnnymi, eto nachal'nye bazowye znanija. Ljuboj malo-malo prilichnyj chislennyj algoritm zabotitsja o wozmozhnoj poteri tochnosti.

[e2k-4d-x-ussr.livejournal.com]

как вариант 2:
написать математику сложения произвольных целых, записанных в виде строк. сконвертировать даблы в безразмерные целые (не более 1000 знаков) и поскладывать.

С умножением и делением сложнее.

[mivlad.livejournal.com]

Да есть такие библиотеки, даже и с умножением и делением и прочим.

Только вот милион чисел таким образом обрабатывать займёт довольно много времени.

[gaius-julius.livejournal.com]

я не помню чтобы меня этому учили, но несколько лет назад столкнулся с такой особенностью и поэтому знаю.

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

[azzo27.livejournal.com]

Про это еще в старинной блестящей книжке Мак-Кракена "Численные методы и программирование на ФОРТРАНе" было.
Надо сначала отсортировать по модулю, а потом складывать, начиная от меньших. Тогда точность была меньше и это было более актуально.

[itman.livejournal.com]

Я знал, потому что мне задали такой вопрос на собеседовании лет десять назад. Правда, там было в контексте вычисления среднего арифметического, а со средним арифметическим проще бороться.

[fantaseour.livejournal.com]

Со средним там и всякие медианы приходится делать :)

А мне помнится как достопамятный Златкис из УПЦ приводил пример, что нельзя считать C_pⁿ в лоб просто посчитав факториалы -- переполнение наступит очень быстро.

[os80.livejournal.com]

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

[(Anonymous)]

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

Это не очень эффективно. Есть способ, требующий O(log N) дополнительной памяти: в одном сумматоре копится сумма не более чем двух элементов последовательности, в следующем сумматоре - сумма не более чем двух предыдущих сумматоров и т. д. Работает это как-то так:

0. S1 = data[0];
1. S1 += data[1] (счетчик наполнился, переносим выше); S2 = S1; S1 = 0;
2. S1 = data[2];
3. S1 += data[3] (счетчик наполнился, переносим выше); S2 += S1 (и этот счетчик наполнился); S3 = S2; S2 = 0; S1 = 0;
и так далее.

Какие именно операции с какими счетчиками надо делать, легко определить по значению соответствующей цифры в двоичной записи номера шага.

[avva.livejournal.com]

Неплохо, да.

Кроме того, есть и простой линейный способ суммирования, не требующий дополнительной памяти, и эффективно снижающий ошибку - в статье, которая раньше тут упоминалась в комментах, о нем говорится (theorem 8).

[http://users.livejournal.com/_kulikov/]

У нас это на третьем курсе изучают в рамках численных методов.

[gt.livejournal.com]

Большинство программистов неправильно сравнивают числа с плавающей запятой - полно примеров в опен сорсе. Введение в правильные способы: http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm

[mopexod.livejournal.com]

Да, и скалярное произведение векторов деленное на произведение их модулей иногда >1, и arcsin падает.
И школьная формула квадратного уравнения на практике плоха или даже совсем не работает.
А уж что происходит при решении "неудобной" линейной системы уравнений с помощью numerical recipes - смотреть страшно.

[ltwood.livejournal.com]

Проблема с арксинусом относится к историческим курьезам (оно было так на одной из серий машины Cray), корректная реализация стандарта IEEE-754 гарантирует в этом случае правильное поведение и для бесхитростной реализации.

Наезд на NR соврешенно беспочвенный, там реализован т.н. алгоритм GEPP т.е. LU-факторизация Гаусса с частичным выбором ведущего элемента. Для матриц общего вида этот метод применяется везде и всюду, не только в NR. Проблемы возникают для слабо обусловленных систем, но слабая обусловленность является свойством задачи и проблемой того недальновидного вычислителя, который не проверил предварительно обусловленность системы.