при работе с целыми числами, округление результата - нетривиально. система команд процессоров персональных компьютеров для этого попросту не заточена
например, на x86 процессорах: перемножили два 32-битных числа, ответ получили в регистрах edx (старшие 32 бита) и eax (младшие 32). теперь, если мы хотим округлить результат до 32 бит, хорошо бы к нашему 64-битному числу прибавить половинку будущего “младшего разряда”, на что нужен ADD на младшие 32 бита и затем ADC, чтобы учесть перенос в старшие 32 бита. так как не хочется усложнять программу (чуть ли не удваивать её размер из-за этого), то чаще всего на это “забивают”, получая результат непосредственно из edx, а про eax забывая
.code64
.global _start
.text
_start:
# rounding to "the nearest whole number" algorithm
movb nums , %al # multiply two bytes
mulb nums + 1
mov %ax , %bx # save result in BX
xor %bl , %bl # drop the low-byte from result
clc
add %al , %al # to set CF in case of AL > 0.5*AH(0)
adc $0 , %bh # if CF is set
# then round result to upper
# else - leave result as it is
mov %rbx , result
mov $1, %rax # write
mov $1, %rdi # to stdout
lea result, %rsi # from buffer address
mov $2, %rdx # buffer length (in bytes)
syscall
mov %rax, %rdi
mov $60, %rax
syscall
.data
nums: .byte 200 , 100
result: .quad 0
с делением тоже тяжело, особенно, со ЗНАКОВЫМ делением. чтобы результат округлился до ближайшего целого, надо бы к делимому прибавить половинку последнего разряда делителя, если их знаки совпадают, или вычесть половинку последнего разряда делителя, если знаки разные
возникает мысль: а так ли это надо вообще? может, просто отбрасывать дробную часть (в смысле, ненужные нам биты числа), и дело с концом? фраза “точность снизится в два раза” звучит страшно. а фраза “точность упадёт на 1 бит” - совсем не страшно :) “Да ладна, у нас их и так 32, или даже 64, и вряд ли они все были значащими, и так сойдёт” :) а если у нас был напряг с точностью, то тем более надо как можно скорее перейти на числа большей разрядности!
в достаточно сложных вычислениях точность из-за отсутствия правильного округления может упасть не в 2 раза, а на ПОРЯДКИ…
плохое округление, “отсечение” лишних бит, приводит к систематической огрешности, которая имеет обыкновение складываться “по амплитуде”: 100 раз в разных местах ошиблись на 0.5 - получаем ошибку на 50 единиц. правильное же округление даёт случайную погрешность, которая складывается “по мощности”, то бишь складываются дисперсии, а среднеквадратичное значение ошибки - это КОРЕНЬ из дисперсии. поэтому 100 раз правильно округлили, в каждом конкретном случае получили ср.кв ошибку. но затем, если мы складываем эти 100 значений, то ошибка помножится на корень из ста, то есть мы получим ошибку 2.9 единиц. и вот она - разница более чем на порядок
простую демонстрацию можно устроить в Excel’е:
в ячейку A1 записали формулу:
=СЛЧИС()*256
в ячейку B1 - формулу для округления вниз, т.е. простого отсечения:
=ОКРУГЛВНИЗ (A1;0)
и в ячейку C1 - формулу для корректного округления:
=ОКРУГЛ (A1;0)
затем копируем это на 100 строчек, и суммируем по каждому столбцу. и можно нажимать F9 (полный пересчёт), чтобы получать новые выборки случайных чисел и смотреть на новые результаты. видно, что сумма правильно округлённых чисел действительно редко уходит от истинного значения более чем на 2-3 единицы, тогда как при отсечении мы всё время имеем существенно меньший результат. не говоря уже о том, что “случайные” ошибки в целом куда приятнее себя ведут, чем систематические. в следующий раз исходные данные немножко изменились - и это дало новую выборку ошибок округления. если мы убрали систематику, то ответы от раза к разу будут крутиться вокруг правильного значения. в противном же случае мы можем с досадным постоянством получать один и тот же ответ, и всякий раз НЕПРАВИЛЬНЫЙ!
практический пример
попробовал напрочь убрать правильное округление из своей процедуры решения системы линейных уравнений 6х6 для алгоритма сопровождения. пока оно было, всё работало в целом неплохо. к примеру, ошибка измерения дальности получалась на 23% больше, чем при использовании плавающей точки, и сопровождение шло уверенное - сколько раз не запускай - можно быть уверенным, что всю дистанцию мы будем получать данные, пусть и чуть зашумлённые. когда я поставил отсечение вместо округления, каждый второй раз происходил срыв сопровождения, причём повторный захват и выход на сопровождение продолжались недолго, как правило, катавасия со срывами продолжается до самого конца. и в целом, шумы существенно увеличились по сравнению с исходным вариантом. если раньше ошибки в целочисленном варианте были на 23% выше, чем в исполнении с плавающей точкой, то теперь ошибки выше на 70%. и ведь это не такие уж громоздкие вычисления - счёт операций пока ещё на сотни, а участков, где происходит округление - и вовсе по пальцам пересчитать! когда туда добавишь операции, ведущие к составлению матрицы (посчитать 12 частных производных помещённых в матрицу 2х6, потом помножать их на транспонированные, образуя матрицу 6х6, и сложение их всех в итоговую матрицу 6х6), к обновлению кватерниона и вектора, всё станет ещё хуже
не думай о младшем разряде свысока! правильная работа с ним позволяет избежать систематического накопления ошибок в процессе вычислений, иначе ошибки доберутся до следующих разрядов. так что и округлению банкира не надо так уж удивляться - это один из способов получить полностью несмещённое округление
это известно любому, кто занимался финматематикой. любой метод составления таблиц помесячных выплат по кредитам (cashflow) с разделением выплат по процентам и по погашению тела самого долга, даёт разные числа. со стороны расчёты кажутся совершенно тривиальными, цифры - проверяемыми на любом этапе.. но нет! 3-4 программы в банке (Excel , mainfraim с++, юникс с явой, функции SAP) - дают разные результаты. очумевший банковский математик залезает на сервер центробанка, находит там публичный веб-сервис, который считает тоже самое, но на компах центробанка, и получает ещё один вариант правильного результата… при покупке холодильника это не страшно, но, вот при финансировании покупки меткомбината, когда у всех сторон сделки есть и свои методы подсчёта кэшфлоу, начинается юридическая буря
>>> Excel , mainfraim с++, юникс с явой, функции SAP - дают разные результаты
даже одна и та же программа на x86-м, скомпилированная с разными флагами, не обязана совпадать сама с собой
даже когда речь о целых числах? хотя да, помню, что переполнение - это Undefined behaviour
там многие считают в Floating Point
там где-то ещё должны BCD обитать и отлично себя чувствовать?
напомню ещё “округление ленивого банкира”, которому лень считать копейки, но хочется быть честным, поэтому он поочерёдно округляет то вверх, то вниз:
IF (MOD(ROW(A1),2)=1,FLOOR (A1),CEILING (A1))
“Казаться, а не быть”. первое округление банкир точно сделает вверх. а сумма остальных в пределе меньше первого округления
понятно, что если мы умножаем два 16-битных числа и заносим результат в 32 бита, то ответ получается ТОЧНЫЙ: это попросту “умножение в столбик” и там неоткуда возникнуть ошибке! но что, если нам нужен только 16-битный результат, округлённый до ближайшего целого? можно ли для того не считать все 32 бита и по прежнему получить все 16 бит ТОЧНЫМИ (ровно такими же, как при получении 32-битного ответа и последующем округлении)?
ответ был для меня немного неожиданным, хотя задним числом - очевидным :)
сначала немного про округление
если мы просто “отрежем” старшие 16 бит от полного ответа, то получим округление “вниз”, что не есть хорошо
если предварительно прибавить к 32-битному ответу 0x 0000 8000 (т.е 1/2 от цены того разряда, до которого округляем), а только потом взять старшие 16 бит - то получим округление “до ближайшего целого”. оно не идеально - Кнут предложил бы “округление банкира”, когда в точности 0.5 округляется до ближайшего чётного целого числа - но мне кажется, что в нашем случае это перебор! такое происходит в среднем 1 раз из 65536, и мы каждый раз округлим вверх, тогда как по “округлению банкира” примерно в половине случаев произойдёт округление “вниз”, т.е лишь в 1 случае из 131072 (менее 1/1000 %) результат будет отличаться
но пока что будем считать, что для округления мы попросту инициализировали аккумулятор не нулями, а значением 0x 0000 8000 (это не усложняет модуль ВООБЩЕ, триггеру пофиг, как инициализироваться!), потом стали накапливать в нём результат и, наконец, взяли лишь 16 старших бит - и такой вариант будет для нас эталонным
итак: можно ли упростить себе задачу и не считать все 32 бита результата, прежде чем провести округление? да, можно обойтись 31 битом :)
действительно, младший бит равен единице, если младшие биты обоих операндов равны единице, и никакого переноса из него произойти не может при всём желании! другими словами, младший бит - это самый угол “параллелограмма”, который образуется, когда мы умножаем числа в столбик, поэтому на верхние 16 бит он повлиять не способен. это не так, если использовать “округление банкира” - тогда и его мы обязаны посчитать!
а вот если срезать аккумулятор хотя бы до 30 бит, то в 1 случае из 262144 мы будем ошибаться на единичку!
приведём пример подобного поведения для 8-битных чисел. помножим 255 на 127. ответ: 32385 (в hex оно выглядит так: FF * 7F = 7E81). Округлённое до 8 бит значение - 127. попробуем теперь помножить их “в столбик”, но без двух последних разрядов:
11111111
x 01111111
---------
11111111
11111111
11111111
11111111
11111111
11111111
11111111
------------------
0111111001111100
теперь прибавим 0x0080 для правильного округления:
01111110,01111100
+00000000,10000000
--------------------
01111110,11111100
теперь “обрезаем” дробную часть - и получается 126, а вовсе даже не 127, как в прошлый раз!
так что если мы страдаем перфекционизмом в плане “ни упустить ни доли процента точности, и ни в коем случае не получить статистически смещённый результат”, то упростить перемножитель не особенно получится!
попробуем более прагматичный подход: само по себе округление до 16 бит уже ограничивает нам точность примерно до 0.29 от цены младшего разряда (корень из 1/12 - таково среднекв.значение ошибки, когда мы идеально округлили числа, равномерно распределённые на числовой прямой). если из-за уменьшенной ширины аккумулятора будет возникать дополнительная ошибка в 0.132 от цены младшего разряда, то полная ошибка возрастёт не более чем на 10% (если они независимы друг от друга - кажется, что так оно и есть. в самом худшем случае 100% корреляции, будет не +10, а +46, тоже не так уж фатально, но думаю величина 10% куда ближе к истине!) - это вполне допустимо “для всех практических применений”. в конце концов, если хочется существенно повысить точность, то нужно увеличивать число бит в конечном результате, иначе всё равно далеко не уедешь!
чтобы среднеквадратичная ошибка составила не более 0.132 от ц.м.р, нам хватит аккумулятора шириной в 23 бита. при этом среднее значение ошибки (“постоянное смещение”) составляет –0.0117 от ц.м.р, т.е. всегда, если происходит ошибка, мы получаем значение на единицу меньше, чем надо, и происходит так 1 раз из 85
кажется, что всё хорошо, 23 бит вполне достаточно. но хочется проверить ещё одну вещь - а не нарушаются ли у нас фундаментальные законы, такие как коммутативность умножения? когда мы смотрим на умножение в столбик, и на то, как мы начинаем “обрубать” младшие биты, уверенность в коммутативности пропадает…
особенно после прочтения книги The End of Error: Unum Computing Джона Густафсона, где он приводит пример суперкомпьютера Cray-1, у которого для экономии “железа” упростили аппаратный перемножитель 64-битных чисел с плавающей точкой (48 бит мантиссы) настолько, что коммутативность умножения перестала соблюдаться, A×B≠B×A. надо отдать должное Сеймуру: данная особенность работы была хорошо документирована в руководстве по эксплуатации. что интересно, в дальнейшем этот суперкомпьютер всё-таки заставили соблюдать коммутативность, причём сделали это весьма своеобразным способом: прежде чем приступить к умножению, два операнда A и B сравнивали, и если A>B, то меняли их местами!
поковыряемся немножко “ручками” - попробуем найти произведение чисел 23 и 42, условно “отрезав” последние 4 разряда:
23 = 010111[2],
42 = 101010[2],
010111
* 101010
-------------
010111
010111
010111
--------------
001111000110,
101010
* 010111
-------------
101010
101010
101010
101010
--------------
001111000110
области, которые обрезаются, кажутся не особенно похожими друг на друга, но дают один и тот же результат в обоих случаях. а если записать умножение целиком, включая нулевые строки, то потихоньку становится понятно - цифры там одни и те же, только немножко в другом порядке:
010111
* 101010
-------------
000000
010111
000000
010111
000000
010111
--------------
001111000110,
101010
* 010111
-------------
101010
101010
101010
000000
101010
--------------
001111000110
и это не совпадение - так будет каждый раз, когда мы будем срезать “ровно”, по вертикали!
т.е, главной ошибкой в CRAY-1 было не то, что они “пожадничали”, а в том что сделали это слишком уж неаккуратно - вместо того, чтобы отпилить лишний кусок строго по уровню, они вырвали ломоть вообще непонятно как. до сих пор не вполне понимаю, чем это обусловлено, но видать были у них на то причины
получается, что приближённое умножение у нас всегда будет коммутировать, какую бы ширину аккумулятора мы не брали. так что я всерьёз рассматриваю возможность немножечко сэкономить и применить 23-битный аккумулятор для перемножения 16-битных чисел. впрочем, сделаем и наиболее “честную” реализацию, вплоть до округления банкира. пока места хватает - пусть живёт :) если же совсем всё заполним - начнём срезать углы
округление чисел в ПЛИС. как бы вы реализовывали округление?
“вниз” (простое отсекание ненужных бит)
“до ближайшего целого” (прибавить 0,5, затем отсечь)
“округление банкира” (рассмотреть случай 0,5 в дробной части)