Сложение двух чисел с плавающей запятой без потери точности +123




Здравствуйте, друзья, как вы думаете, если мы напишем такой код:

s = a+b;
z = s-a;
t = b-z;

то не кажется ли вам, что в результате его выполнения получится, что t=0? С точки зрения привычной математики действительных чисел это и правда так, а вот с точки зрения двоичной арифметики с плавающей запятой в переменной t будет кое-что другое. Там будет то, что спасает нас от потери точности при сложении чисел $a$ и $b$. Кого интересует данная тема, прошу под кат.


По моей традиции текстовая статья дублируется на видео. По содержанию текст и видео идентичны.



Читателю статьи для её восприятия нужно понимать способ представления двоичных чисел с плавающей запятой в формате IEEE-754 и понимать, почему, например, 0,1+0,2?0,3, но если такого навыка у вас по каким-то причинам нет, но вы хотели бы его приобрести, то прошу обратить внимание на список источников конце статьи, на пункты [1] и [3-5].

Итак, у нас следующая проблема. Сумма двух чисел с плавающей запятой c=a+b может вычисляться с некоторой ошибкой. Знаменитый на весь интернет пример 0,3 ? 0,2+0,1 хорошо это показывает. Наша задача в том, чтобы всё-таки складывать числа без этой ошибки. То есть сделать так, чтобы мы могли как-то сложить 0,2+0,1 и хоть в каком-то виде знать точный результат. В каком смысле точный, если даже исходные числа 0,1 и 0,2 не имеют точного представления в формате IEEE-754? Вот сейчас и поясню.

Начнём с того, что чисел 0,1 и 0,2 в двоичной арифметике с плавающей запятой быть не может, а наиболее близкие к ним значения для типа данных double (число удвоенной точности binary64, так его называют в Стандарте IEEE-754) следующие:

$0{,}1 > 0{,}1000000000000000055511151231257827021181583404541015625 ={}\\ {}=\frac{3602879701896397}{36028797018963968}$


$0{,}2 > 0{,}200000000000000011102230246251565404236316680908203125 ={}\\ {}=\frac{3602879701896397}{18014398509481984}$


Десятичный результат получен с помощью правильного онлайн-конвертера, а дробь посчитана с помощью библиотеки MPIR и функции mpq_set_d, которая умеет переводить double к рациональной дроби.

К сожалению, это данность, от неё никуда не уйти, если хранить числа в типе данных double (или любом другом типе фиксированного размера из Стандарта IEEE-754). Но проблема, которую мы решаем, другая. Нам нужно эти два числа, наиболее близких к 0,1 и 0,2, сложить так, чтобы получить результат без погрешности. То есть чтобы в результате сложения иметь число

0,1000000000000000055511151231257827021181583404541015625 +
0,2000000000000000111022302462515654042363166809082031250 =
0,3000000000000000166533453693773481063544750213623046875.

К сожалению, если просто написать код s=0.1+0.2, то мы получим кое-что другое, а именно

s == 0.3000000000000000444089209850062616169452667236328125

что отличается от правильного ответа ровно на

$t = 2{,}77555756156289135105907917022705078125 \cdot 10^{-17}$


«Подумаешь, — скажете вы, — десять в минус семнадцатой! Мне чтобы пиксель на экран вывести такая погрешность не помеха!». И будете правы. Задача точного вычисления каких-либо выражений относится к очень узкой области Computer Science, связанной с разработкой математических библиотек для языков программирования. Когда вы вызываете какую-нибудь функцию из такой библиотеки, то можете и не подозревать, что в основе её работы лежит труд десятков и сотен человек, выполняемый на протяжении десятилетий, а работает такая функция правильно благодаря совершенно неочевидным алгоритмам… Вот я и хотел бы приоткрыть для вас этот удивительный мир, для чего и пишу такие научно-популярные статьи.

Зачем может понадобиться точное сложение двух чисел? Приложений много, но вот одно из них, которое будет понятно всем. Если вы хотите сложить все числа из массива X[N] и сделаете это вот так:

S = 0.0;
for (i=0; i<N; ++i)
  S = S + X[i];

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

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

У нас есть два числа $a$ и $b$ типа double. Нам не важно, что эти числа $a$ и $b$ уже содержат в себе какую-то погрешность, полученную, например, в результате конвертирования десятичной строки в формат IEEE-754. Важно, что они сейчас перед нами, а их предыстория нас не интересует. Нужно сложить эти два числа c=a+b так, чтобы в результате сложения не возникло погрешности.

Но ведь это невозможно!

Да, это невозможно, спасибо что дочитали, до новых встреч :)

Шучу, конечно. Это невозможно только если мы храним результат сложения в одной переменной того же типа данных. Но вспомните пример выше. У нас была переменная s, и переменная t. Причём s была округлённым результатом сложения a+b, а t была равна разности этого округлённого результата и точного значения суммы. При этом выполняется равенство s+t=a+b.

И вот хорошая новость состоит в том, друзья, что такое представление суммы a+b в виде суммы s+t можно выполнить всегда (если a+b??)! Если s=a+b оказывается точно-представимым значением в типе double, то очевидно, что t=0. Если это не так, то t будет равно некоторому очень маленькому (по сравнению с s) числу, которое является точно-представимым. Итак, вот оно, фундаментальное свойство суммы чисел с плавающей запятой: ошибка округления в результате суммирования чисел типа double всегда будет точно-представимым числом типа double! А это означает, что пара чисел (s, t) всегда даёт нам возможность сохранить сумму чисел a+b «как бы» без погрешности. Да, мы вынуждены хранить две переменные вместо одной, но прелесть в том, что их всего две, а не больше.

Теперь опишем это свойство на математическом языке. Введём обозначение RN(x) – это результат приведения произвольного вещественного числа x к типу данных double по правилу округления round-nearest-ties-to-even, то есть это округление к ближайшему, но в случае равного удаления от двух ближайших к тому, у которого последний бит мантиссы равен нулю (чётный). Именно этот режим работает почти везде по умолчанию, то есть если вы не знаете, о чём я сейчас говорю, то в вашем процессоре на 100% работает именно этот режим, так что не беспокойтесь.

Пусть a и b – числа типа double. Пусть |a|?|b| и RN(a+b) ? ?. Тогда следующий алгоритм

s = RN (a+b);
z = RN (s-a);
t = RN (b-z);
return (s, t);

вычисляет пару (s, t) таких чисел, для которых:
  1. s+t = a+b (в точности)
  2. s – число типа double, наиболее близкое к a+b.


Основная цель моих статей — объяснять сложные вещи простым языком, поэтому давайте я построю для вас образ, помогающий понять этот алгоритм. Пожалуйста, посмотрите на рисунок, а ниже даётся его описание.


Прямоугольник олицетворяет переменную с плавающей запятой фиксированного размера, поэтому прямоугольники, помеченные символами «а» и «b» имеют одинаковый размер, но $b$ смещён относительно $a$ вправо, потому что у нас есть условие: |a|?|b|. Третий прямоугольник «Точное а+b» — это некоторое промежуточное значение, которое может иметь больший размер, чем размер переменных $a$ и $b$, поэтому оно будет округлено, и «хвостик», вылезающий за границы нашего типа данных, будет отброшен. Таким образом, мы переходим к четвёртому прямоугольнику «Округлённое a+b», это и есть наше значение s, полученное в первой строке алгоритма. Если теперь из s снова вычесть $a$, то останется только «b без хвостика». Это делается во второй строке алгоритма. Теперь мы хотим получить «хвостик» от $b$, и это делается в третьей строке алгоритма, когда из исходной переменной $b$ мы вычитаем «b без хвостика». Остаётся «хвостик», это и есть переменная t, которая показывает ту самую погрешность, которая возникла в ходе округления. При этом из данных рассуждений видно, что такой «хвостик» всегда будет умещаться в одну переменную, потому что он не может превышать размера переменной $b$. В крайнем случае, когда смещение $b$ относительно $a$ будет настолько большим, что s=RN(a+b)=a, то в этом случае, очевидно, z=0 и t=b. Скучный рисунок на эту тему я вам предлагаю изобразить самостоятельно.

Если вы не находите картинки убедительными для себя, то в книге [1, раздел 4.3.1, теорема 4.3] есть строгое математическое доказательство, но оно не умещается в формат научно-популярной статьи. Его краткая суть в том, что в нём показано почему в строках 2 и 3 алгоритма не будет выполняться округления, то есть эти выражения вычисляются точно, а если так, значит t=b-z=b-(s-a) = (a+b)-s в точности, что нам как раз и нужно: t является разностью между реальной суммой a+b и её округлённым значением s.

Давайте рассмотрим несколько примеров, удобных для восприятия человека.

Пример 1.

$\begin{array}{rcl} a &{}={}& 9007199254740991 = 2^{53}-1, \\ b &{}={}& 2, \\ s &{}={}& 9007199254740992, \\ t &{}={}& 1. \end{array}$


Пояснение. Реальное значение a+b=253+1, однако в типе данных double младший бит, равный единице, не влезет в 52 бита дробной части мантиссы, по какой причине будет выброшен при округлении, но переменная t «подхватит» его и сохранит в качестве погрешности, а переменная s сумеет сохранить только 253 ровно. Легко видеть, что s+t будет равно 253+1.

Пример 2.

$\begin{array}{rcl} a &{}={}& 1152921504606846976 = 2^{60}, \\ b &{}={}& 1073741823 = 2^{30}-1, \\ s &{}={}& 1152921505680588800, \\ t &{}={}& -1. \end{array}$


Пояснение выполнено в двоичном коде.

  a = 1000000000000000000000000000000000000000000000000000000000000
  b =                                111111111111111111111111111111
a+b = 1000000000000000000000000000000111111111111111111111111111111
      -------------------------------------бит округления--^ 
  s = 1000000000000000000000000000001000000000000000000000000000000

Если выравнивание кода не поехало, то вы видите, где у суммы a+b будет бит округления, и что само округление будет выполнено вверх. Чтобы снова получить паровоз из 30 единиц, нужно вычесть единицу из s=260+230. Поэтому t=-1. Как видим, t может быть отрицательным, когда округление выполняется вверх.

Я не говорил этого явно, но алгоритм работает даже когда b<0, то есть фактически происходит вычитание.

Пример 3.

$\begin{array}{rcl} a &{}={}& 9007199254740991 = 2^{53}-1, \\ b &{}={}& -2251799813685247{,}75 = -\frac{2^{53}-1}{4}, \\ s &{}={}& 6755399441055743, \\ t &{}={}& 0{,}25. \end{array}$


В этом примере число $b$ является сдвигом вправо на 2 бита числа $a$, поэтому если $a$ влезает в double «впритык», ясно, что вычитание здесь точным получиться не может. Останется небольшой хвостик размером в два бита после запятой.

Из доказательства, приведённого в [1], следует один положительный момент данного алгоритма: в ходе его выполнения не может возникнуть переполнения во второй и третьей строках, если оно не возникло при вычислении s=RN(a+b).

Недостатком является то обстоятельство, что алгоритм работает только для |a|?|b| (однако если вы посмотрели доказательство, то знаете, что в действительности достаточно более слабого условия, чтобы экспонента $a$ была не меньше экспоненты $b$, но проверить это условие намного сложнее). То есть нам нужно либо делать лишнюю проверку, либо как-то иначе гарантировать правильный порядок чисел на входе алгоритма. Ветвление не всегда обрабатывается достаточно быстро, это зависит от процессора и от того, как на это ветвление посмотрел компилятор (иногда он может его убрать, иногда — нет, а иногда он его убрал, но оказалось так плохо, что лучше бы не трогал). Возникает вопрос: а можно сложить числа без ошибок округления, если их порядок друг относительно друга заранее неизвестен?

Можно. Для этого есть второй алгоритм, но в нём вместо 3-х операций их уже 6.

  s = RN (a+b);
  b' = RN (s-a);
  a' = RN (s-b');
  d_b = RN (b-b');
  d_a = RN (a-a');
  t = RN (d_a+d_b);
  return (s, t)

Как это работает? В случае когда |a|?|b|, алгоритм работает по сути так же как предыдущий, и строгое доказательство для этого случая, которое вы можете отыскать в [2, раздел 2.3, теорема 7], сводится к тому, чтобы показать, что «лишние» три операции не портят этого результата. А вот основная сложность доказательства в том, чтобы рассмотреть случай |a|<|b|, он значительно труднее и тоже выходит за рамки лёгкого научно-популярного изложения. Схематичная картинка общей ситуации при |a|<|b| даётся ниже с пояснением.


Первые два прямоугольника показывают соотношение между $a$ и $b$, на этот раз $a$ смещено относительно $b$ вправо. Следующие три прямоугольника показывают процедуру сложения и отбрасывания «хвостика». Дальше начинается основная сложность. Вычисляем b’=s-a, то есть мы вычитаем из округлённого значения s число $a$ и выброшенный «xвостик», а это значит, что b’ будет чуть меньше, чем $b$, и этот недостаток отмечен красным прямоугольничком. Другое дело, когда мы вычисляем a’=s-b’, то есть вычитаем из суммы величину «b с недостатком», а значит получим «а без хвостика с избытком», и этот избыток обозначается зелёным прямоугольничком. Если мы теперь вычтем из $b$ величину b’ (которая с недостатком), то получим избыток db. Далее вычитаем из $a$ величину a’ (которая с избытком), и получаем хвостик с недостатком da. Если теперь сложить хвостик с недостатком и избыток, мы получим «чистый хвостик».

У читателя может возникнуть вопрос: а в каких случаях можно с уверенностью сказать, что t=0? Иными словами, про какие пары чисел (a, b) точно известно, что их сумма будет точно-представимым (без округления) числом в том же типе данных?

Мне известны только три теоремы на этот счёт, которые также приводятся в [1] и [2].

Теорема 1 [книга 1, раздел 4.2.1, теорема 4.2]. Если сумма RN(a+b) оказалась денормализованным числом, то сумма a+b является точно-представимой.

Теорема 2 [статья 2, раздел 2.2, лемма 4] Если |a+b|?|a| и |a+b|?|b|, то число a+b является точно-представимым (аналогичное верно и для разности).

Теорема 3 [книга 1, раздел 4.2.1, лемма 4.1] Пусть a?0 и a/2 ? b ? 2a, тогда разность a-b будет точно представимой.

Бдительный читатель сразу вспомнит про так называемую ошибку критического сокращения старших значащих битов, называемую в литературе “Catastrophic Cancellation” и спросит: «Как же так, ведь давно известно, что когда мы вычитаем очень близкие друг к другу числа, там возникает такая погрешность, что ой-ой-ой! Как ты говоришь, что в этом случае разность будет точно-представимой».

Дело в том, что Catastrophic Cancellation – это не ошибка вычисления, а неизбежное свойство чисел с плавающей запятой; это свойство не порождает никакую новую погрешность, а как бы «высвечивает» ту, что уже изначально была заложена в уменьшаемом и вычитаемом. Если не смотря на уже имеющееся в интернете описание вы хотите подробного обсуждения данной особенности чисел с плавающей запятой именно в моём исполнении, то пишите об этом в комментариях, пойду навстречу.

На основе рассмотренных в этой статье алгоритмов в следующей мы будем более правильно складывать числа из массива.

Благодарю за внимание!

Список источников


[1] Jean-Michel Muller, “Handbook of floating-point arithmetic”, 2018.
[2] Jonathan Richard Shewchuk, “Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates”, Discrete & Computational Geometry 18(3), 1997, pp. 305–363.
[3] На Хабре: «Что нужно знать про арифметику с плавающей запятой».
[4] На Хабре: «Наглядное объяснение чисел с плавающей запятой».
[5] Учебный видео-курс для «самых маленьких», предельно наглядное разъяснение чисел с плавающей запятой в 8-ми уроках. Первый урок на на YouTube.




Комментарии (156):

  1. v1000
    /#22192320

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

  2. datacompboy
    /#22192644

    Пока непонятно что мы получили. Было два числа такие, что их сумма представляла некую сумму. Мы получили два других числа с той же суммой.
    Знаю, что дальше магический финт ушами будет, но пока польза не видна.

    • encyclopedist
      /#22192932 / +2

      См например Kahan summation algorithm

      • byman
        /#22195302

        Я попробовал этот алгоритм на сумме float массива из 1К элементов и действительно этот алгоритм дал в 10 раз меньшую ошибку по сравнению с ошибкой обычного суммирования. За эталон брал сумму в double формате. Возможно удачные попались массивы чисел :)

  3. nin-jin
    /#22192692 / +1

    Тема интересная, но подача всё же сложная для восприятия. "Физический смысл" прямоугольников я так и не понял. Скорее всего метафора не удачная или нужны дополнительные пояснения. В идеале, тут нужна такая визуализация, чтобы и без формул было всё понятно.


    Я являюсь организатором митапов PiterJS, где помогаю докладчикам улучшать свои выступления. А так же являюсь автором канала на близкую тематику — Core Dump.


    Поэтому, у меня к вам деловое предложение: когда будете готовить следующее видео, напишите нам в чат с предложением поревьюить черновой вариант, — постараемся помочь сделать ваше видео более доступным.


    Так же возможен вариант какой-либо коллаборации или даже объединения.

    • ArtemKaravaev
      /#22192698

      Благодарю, Дмитрий. Я приму к сведению ваше предложение.

    • avdx
      /#22193134

      По-моему с прямоугольниками очень даже понятно. Хотя возможно это индивидуально. Но конечно нужно уже иметь какое то представление, что такое floating point.

    • Mingun
      /#22194580 / +1

      По прямоугольникам следует сразу же пояснить, что они представляют размерность разрядной сетки. Поэтому и одинаковые. А сдвиг обозначает относительное положение десятичной запятой в этих двух числах. Тогда сразу всё становится понятно.

  4. ICELedyanoj
    /#22192772 / +4

    Ох уж эти плавающие запятые…
    Помню как в ВУЗе на каком-то предмете начали нас гонять по задачам, решаемым при помощи симплекс-таблиц. Расчётки, контрольные и т.д.
    На кафедре ИКТ у нас был преподаватель, у которого была своя софтина для расчёта симплекс-таблиц, и старшекурсники ходили к нему на поклон чтобы рассчитать таблицы для диплома. Понятное дело не за спасибо.
    Как-то раз я увидел одну из распечаток с результатами, и был неприятно удивлён, так как увидел промежуточные и конечные результаты в стиле 0.9999999999.
    В симплекс-таблицах вся методика построена на делении одних чисел на другие, причём итеративно. При решении в десятичных дробях гарантированно будет накапливаться погрешность. Даже без плавающей запятой — просто за счёт округления.
    Поэтому классический вариант предполагает использование простых дробей.
    В общем захотелось мне утереть нос преподу. Я тогда только-только начал перебираться из MS-DOS в Windows 95, и первое что мне попалось под руку из инструментов — VBA для Excel.
    В итоге мне удалось написать макрос, расчитывающий симплекс-таблицы прямо на листе Excel. При этом, дабы не терять в точности, я обучил макрос работе с простыми дробями.
    Помню как смеялась комиссия на олимпиаде, где я рассказывал о причинах, побудивших меня к написанию этого макроса. Сказал, что меня не устраивает 0.999999999 вместо единицы. Сказал без фамилий, но вся комиссия прекрасно знала о ком идёт речь, потому и развеселились.
    Насколько я знаю, в самом Excel встроенный метод расчёта симплексов появился только спустя несколько лет.

  5. tsukanov-as
    /#22192826 / +1

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

    Это про двоичные числа утверждение? Почему тогда в статье десятичные?
    Кажется единственное, что тут нужно понимать, что компилятор «портит» ваши десятичные числа, которые вы зачем-то даете ему. И наоборот, портит двоичные, когда вы их зачем-то хотите видеть как десятичные.
    Хотите десятичные вычисления? Не пользуйтесь двоичными.

    • aamonster
      /#22192938

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

    • netch80
      /#22193080

      Хотите десятичные вычисления? Не пользуйтесь двоичными.

      Да, сегодня это для соответствующих применений (по советской терминологии, "экономических") наилучший вариант.
      Но это таки (пока?) заметно дороже, при равных длинах и точности.
      А ещё во многих случаях там, где нам кажется, что значения десятичные, они на самом деле двоичные, но просто так пахнут :) например, датчики обычно отдают двоичные значения. Не зря в C99 ввели "%a" для printf.

      • tsukanov-as
        /#22193274

        Я не пробовал, но с некоторых пор мне кажется что «экономические» вычисления можно вполне выполнять в копейках. Другими словами, на целых числах. Где нужно какие-то проценты посчитать, конвертировать в double, а потом обратно в целое. Это будет быстрее любой десятичной реализации, с тем же результатом.
        double позволит использовать в промежуточных расчетах доли копеек, но в конечном результате они обычно не нужны. Ну а если доли нужны, то просто считаем в долях опять же на целых числах. Тут можно упереться только в максимальные значения, если оперировать триллионами рублей.

        • netch80
          /#22193342

          Другими словами, на целых числах.

          Очень часто именно так и делается. Тем не менее есть альтернативы. Есть десятичная плавучка (у IBM даже с аппаратной поддержкой). Есть Decimal из C#, у которого 28 цифр мантиссы. Есть decimal управляемой точности в Python и ещё много где.


          Тут можно упереться только в максимальные значения, если оперировать триллионами рублей.

          Ну тут уже больше триллионов надо :)
          2^63/100 ~= 90 квадриллионов (коротких).

          • Regis
            /#22194916

            Более того, в IEEE 754 от 2008 года добавлены типы decimal32, decimal64, decimal128. И они даже реализованы на аппаратном уровне в некоторых (немногих) процессорах Intel.

            • netch80
              /#22195160

              Более того, в IEEE 754 от 2008 года добавлены типы decimal32, decimal64, decimal128.

              Я об этом же :)


              И они даже реализованы на аппаратном уровне в некоторых (немногих) процессорах Intel.

              Вот этого не видел. Ни в "Software Developer’s Manual", ни в "Instruction Set Extensions and Future FeaturesProgramming Reference" такого не описывают. Это какая-то особая фишка, или это вообще не x86?

        • ArtemKaravaev
          /#22194926

          Тут можно ещё упереться в нарушение закона. Дело в том, что если где-то прописано, что процентная ставка 0.1%, то при временном переводе её в double, мы получаем другую ставку, и нарушаем закон. Так что здесь нужно сначала доказать, что такой перевод повлияет только на скорость, но не на ответ.

          • gecube
            /#22194980

            Дело в том, что если где-то прописано, что процентная ставка 0.1%, то при временном переводе её в double, мы получаем другую ставку, и нарушаем закон.

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

            • ArtemKaravaev
              /#22194990

              Мне известно только два способа: рациональные числа и арифметика с произвольной точностью (увеличиваем точность так, чтобы эти копейки, даже возведённые в степень, не сыграли бы роли). Если бы я разрабатывал финансовое ПО, то смотрел бы в сторону библиотеки длинной арифметики MPIR и там есть соответствующие типы данных mpq (рациональные) и mpf (с плавающей запятой). Слишком сложных вычислений там быть не может, поэтому каких-то особенных тормозов не будет.

              • gecube
                /#22195004

                Физика с Вами не согласна )
                В том плане — что мне из банка не могут выдать вклад (или его часть) дробными частями копеек.
                А раз — где-то в алгоритме должно быть учтено, что "мнимая" часть денег от вклада — должна в момент снятия вклада "исчезнуть", а то банки оставались б должны всем своим клиентам )

                • ArtemKaravaev
                  /#22195024

                  Я говорю о том, чтобы делать указанным способом именно промежуточные вычисления, а выдавать уже с округлением. Накопившиеся именно в ходе выдачи округления распределять по закону, но я законов этих не знаю. Так что по остальной части вопроса — это не ко мне :)

                • chapuza
                  /#22195040 / -1

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

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


                  В первый год вам дадут рубль, во второй — тоже, а в третий — два рубля.

                  • nikolayv81
                    /#22195610 / +1

                    Это так не работает…
                    Есть счёт накопленных процентов, на который происходит начисление в соответствии с правилами (по сути формула) у этого счёта (для РФ) точность до копейки, большей точности там просто не будет. С момент причисления средств к основному счёту(капитализации)/выплаты сумма просто просто переводится с одного счёта на другой.

                    • chapuza
                      /#22195856

                      для РФ

                      Волшебный прием: вы сужаете целевое множество утверждения «мне из банка не могут выдать вклад (или его часть) дробными частями» (все банковские институты) до знакомого (банки РФ), после чего следует волшебный пасс: «это так не работает».


                      Это так не работает в банках РФ потому что писали правила, а потом переносили их в код — не особо умные люди, не желающие подумать.

                      • nikolayv81
                        /#22195936 / +1

                        Это так не работает потому что это правила бух.учёта по РСБУ, детально по US GAAP расписать не могу(, но вы наверное знаете как это работает по правилам в США? Или хотя бы как это происходит в Испании в банковском секторе?

                        • chapuza
                          /#22197228 / -1

                          как это происходит в Испании в банковском секторе?

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


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


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

                • nikolayv81
                  /#22195594

                  Округление происходит по правилам, причём после каждого "периода" при внесении на внутренний счёт накопленных процентов, правила оговариваются (либо правилам математического округления, либо по другим (т.н. бухгалтерское округление), никаких промежуточных данных в банковском по, по сути, не сохраняется, только остатки на счетах, у которых есть объявленная точность...

            • netch80
              /#22195168 / +1

              Если мы считаем процентную ставку помесячно и округляем в меньшую сторону — на сроке год — она эффективно получается меньше (на доли копеек, но тем не менее). Как с этим бороться?

              Где-то должны быть просто записаны правила. Если не в законе, то в договоре. Или в каких-то нормативных документах.


              Таких проблем полно. Вот с которой я сталкивался когда-то: пусть НДС 20% (типовой для Украины), а цены записаны без НДС. Вот цена одной единицы товара 1 гривна 12 копеек, соответственно НДС будет 22 копейки. А вот две единицы товара — 2 гривны 24 копейки, тогда НДС уже будет 45 копеек, а не 44. Как правильно считать — по всей накладной, по каждому пункту или по каждой единице товара?
              Насколько я помню, несколько лет налоговая морочила голову, пока не постановила — по накладной. Но на это потребовалось особое разъяснение.

              • VolCh
                /#22195348

                Разъяснение регулятора по нечеткому закону — лучший способ уменьшить риски штрафных санкций. Недоплату может покрыть потребуется, но без санкций, если их опротестовать.

                • netch80
                  /#22195356

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

            • VolCh
              /#22195336

              Мы просто придумали алгоритм разбрасывания копеек, чтобы итоговая сумма процентов на дату очередного начисления была равна обычно (half up) округленной сумме за весь прошедший период. Запросили разъяснения у регулятора допустимо ли, что в отдельные периоды на копейку разница по сравнению с расчётом за один период — ответили что допустимо, если итоговая сумма на конец периода совпадает. То есть если на 6-й месяц клиент видит, что ему на копейку больше начислили чем за 5-й, но за 6 месяцев вместе ему начислили 1/2 годовой, то норм. С аннуитетом решили не рисковать не соблюдением термина "равные платежи" и доли копеек отбрасывали в пользу клиента.

            • BloodyEnterprise
              /#22196718

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

        • remzalp
          /#22196772

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

    • edo1h
      /#22193710 / +1

      И наоборот, портит двоичные, когда вы их зачем-то хотите видеть как десятичные.

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

      • netch80
        /#22195204

        Всё формально так, но дальше дело в длине записи и в обратном переводе (потому что если перевели в десятичные, то надо и вернуть в двоичные?)


        Разные системы конверсии тут имеют разные правила. Например, repr() в Python сейчас действует по принципу: показывает форму с минимальной длиной из тех, что конвертируются в данное число во внутреннем (двоичном) представлении при стандартном правиле округления. Поэтому, например, если вы ввели 0.1, оно подберёт на выходе текст "0.1", а если вы его умножили на 3, оно будет на одну йоту отличаться от того, что ближайшее к 0.3, и поэтому будет выведено как "0.30000000000000004". Всё из-за задачи предельно честного вывода. А вот какой-нибудь "{:g}".format(0.1*3) выведет 0.3, потому что у него умолчательное усечение до 6 значащих цифр.
        А если бы не было бы такого усечения — то всё вы выводилось максимально длинно:


        >>> '%.55g' % 0.1
        '0.1000000000000000055511151231257827021181583404541015625'
        >>> '%.55g' % 0.3
        '0.299999999999999988897769753748434595763683319091796875'
        >>> '%.55g' % (0.1*3)
        '0.3000000000000000444089209850062616169452667236328125'

        ну и кому такие подробности нужны?


        По отношению к этим конверсиям известны следующие цифры:


        • Сколько максимум десятичных цифр (на всех значениях порядка) могут быть сохранены при переводе в двоичную форму и обратно. Оно в C++. Для float32 и float64 (double) соответственно 6 и 15.
        • Сколько максимум десятичных цифр требуется для перевода любого числа из двоичного внутреннего представления в десятичное и обратно, чтобы получить то же число. Оно в C++. Для float32 и float64 (double) соответственно 9 и 17.

        И поэтому умолчательное (почти везде) форматирование в 6 значащих цифр это как-то смешно.

        • ArtemKaravaev
          /#22195224

          Я вам больше скажу, почти все компиляторы не умеют переводить десятичное число в двоичное правильно. Кроме GCC последних версий. Все остальные, что были доступны, были успешно нами завалены. Когда число очень близко к границе округления, компилятору не хватает внутренней точности чтобы округлить правильно. Например, вот это число:


          float x = 1.26765082690182057923967346278399
          99999999999999999999999999999999999999999999
          99999999999999999999999999999999999999999999
          99999999999999999999999999999999999999999999
          9999999999999999999e30f;

          Правильный ответ 0x71800001. Компиляторы Visual C++ и Intel C++ (любых версий) выдают на один бит больше. Таких примеров у нас сотни.

          • netch80
            /#22195272 / +1

            Кроме GCC последних версий.

            Clang в диапазоне от 3.8 до 10.0 включительно — тоже выдал с 01 в конце. Стандартная сборка Ubuntu, если это важно.


            Виндовых компиляторов у меня нет, а вот numpy.float32 выдал с 02. Голый Python — тоже, но у него нет float32, поэтому работает промежуточный double. Проверил на GCC/Clang — если явно константу задать в double и потом из неё уже делать float, то получается в конце 02. Может, у остальных проблема из-за той же промежуточной конверсии?


            За пример спасибо, сложил в копилку.

            • ArtemKaravaev
              /#22195354

              Я не знаю как работают компиляторы, исходников у меня нет, но знаю о том, почему так получается. Дело в том, что есть так называемые пограничные случаи, когда десятичное число очень-очень близко к середине между двух точно-представимых. Так близко, что разница проявляется только в каком-нибудь тысячном бите после запятой. Каким бы ни был алгоритм конвертирования, если он работает с фиксированным числом битов, он будет завален каким-нибудь таким тестом. Поэтому делаются эти тесты очень просто. Берём число, лежащее точно посередине. Прибавляем к нему (или вычитаем из него), например, 2 в степени минус тысяча (десять тысяч, миллион) — и всё, тест готов. Промежуточные округления тут роли не играют, причина ошибки — недостаток битов в процессе вычисления, алгоритм не может охватить все тысячу битов, чтобы увидеть направление округления. Нужна длинная арифметика.


              По поводу Clang, там я точно уже не помню, но он на long double у меня легко заваливался, сейчас уже не скажу, какие тесты были.

              • netch80
                /#22195820 / +1

                Я не знаю как работают компиляторы, исходников у меня нет

                GCC, Clang — они есть. MPFR, MPC — тоже. Код math.h в какой-нибудь glibc — тоже. Не знаю, может, у вас какие-то навязанные снаружи ограничения с использованием только MSVC хозяйства, но сейчас это становится тупо непрактичным: наука давно уже упорно переползает на open source именно потому, что оно открыто, проверяемо и в конце концов лучше развиваемо.


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

                Да, и ваша предыдущая статья была об этом. Это я помню (и замечал там, что в вопросе про ULP вы абсолютизируете влияние этой самой проблемы, там, где во многих местах вероятнее предположить просто более ранние версии, или намеренную жертву точности в пользу скорости работы). Но именно код рассматриваемых компиляторов, когда компилирует программу, использует неограниченную точность (GCC явно требует для своей сборки комплект из GMP, MPFR и MPC). Поэтому вот это


                если он работает с фиксированным числом битов

                там, скорее всего, не выполняется (не поручусь за совсем тонкие случаи, но у них цель и метод именно такие — компиляция делается с произвольной точностью).


                По поводу Clang, там я точно уже не помню, но он на long double у меня легко заваливался, сейчас уже не скажу, какие тесты были.

                Ну, может, старая версия, а, может, проблема прикрученности этого самого x87 long double ржавой проволокой сбоку (и то наверняка давно исправили). И опять же вопрос в платформе. Это на Unix обычно long double честный, а на Windows/64 он алиас для просто double.

                • ArtemKaravaev
                  /#22195858

                  Я имел в виду нету исходников Intel С++ И MS C++.


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


                  Наука на Open Source не совсем вся переползает. Например, многие используют Maple (закрытый) для своей работы и некоторые его функции настолько сильнее чем у других (бесплатных) систем компьютерной алгебры реализованы, что тем до Maple ещё как до Луны. И я не верю, что это можно сделать (перейти на OpenSource), пока в нашем обществе деньги в принципе как явление существуют. Но это уже философский вопрос. Сам я против Open Source, но ни коим образом других людей не осуждаю за то, что деньги зарабатывают, на бесплатном софте сидя. Право каждого, не моё дело.


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

    • ArtemKaravaev
      /#22195164

      Указанная в статье проблема, связанная с тем, что при сложении чисел с плавающей запятой возникает погрешность, характерна для любых целых систем счисления. В книге [1], на которую я ссылаюсь, теоремы доказываются именно для произвольной системы.

      • tsukanov-as
        /#22195230

        Да так то оно так, но имхо вы пример выбрали такой, из которого можно подумать будто погрешность есть всегда. Берем python:
        >>> x = 0.25
        >>> y = 0.125
        >>> z = x + y
        >>> z == 0.375
        True
        >>> z - y == x
        True

        Ваш пример не работает потому что эти числа непредставимы в двоичной системе, и получается что они там «на грани» точности.
        Начало статьи как-то сформулировано не очень, имхо. Причина казуса не прояснена. И складывается впечатление, что вы боретесь с конвертацией десятичных чисел в двоичные. Если принять что BIN() — это операция конвертации, то:
        BIN(0.1) + BIN(0.2) != BIN(0.3)
        С чего бы этому быть не так? Дело то не в погрешности тут. А вот если бы вы показали:
        >>> x = 0.1 # конвертировали в двоичное
        >>> y = 0.2 # конвертировали в двоичное
        >>> z = x + y # какая-то двоичная сумма
        >>> z - x == y # ожидаем так?
        False # внезапно!

        То тут уже явно видно что дело в погрешности, которая возникает из-за того что числа «на грани» в двоичной системе.

        • ArtemKaravaev
          /#22195252

          Пример на то и пример, чтобы проиллюстрировать проблему, но не описать всё её многообразие проявлений. Более того, я пояснил, что статья для тех кто в теме, кому не нужно пояснять что такое 0.1+0.2 и почему это не равно 0.3. В этом примере причина ошибки НЕ В ТОМ, что числа на грани, а в том, что в результате сложения происходит округление суммы. Разложите, пожалуйста, эти числа на двоичное представление, сложите руками на бумаге и увидите, почему происходит округление. В учебном курсе [5] даётся подробное разъяснение именно этого примера в каком-то уроке, не хотелось бы здесь делать учебную комнату, уж прошу меня простить :)

          • tsukanov-as
            /#22195296

            «на грани» взято в кавычки если что :).

            Разложите эти числа на двоичное представление

            Так их нельзя разложить без потерь, о том и речь. Формула 0.1+0.2=0.3 не обязана работать в двоичной системе безотносительно округлений при сложении.

            • ArtemKaravaev
              /#22195366

              Я имею в виду, возьмите соответствующие ТОЧНЫЕ значения (которые умещаются в double без округления), указанные в статье:


              0,1000000000000000055511151231257827021181583404541015625 +
              0,2000000000000000111022302462515654042363166809082031250 =
              0,3000000000000000166533453693773481063544750213623046875.

              Разложите в двоичный вид и увидите, почему будет округление. Именно ЭТО округление и есть причина описываемой в статье ситуации. Если бы его не было, то условие if (0.1+0.2==0.3) срабатывало бы всегда, даже если эти числа сами по себе точно не могут быть представлены в двоичной арифметике.

              • tsukanov-as
                /#22195450

                Если бы его не было, то условие if (0.1+0.2==0.3) срабатывало бы всегда, даже если эти числа сами по себе точно не могут быть представлены в двоичной арифметике.

                Так нет же. Это логическая ошибка. Еще раз:
                BIN(0.1) + BIN(0.2) != BIN(0.3)
                Нет никаких причин чтобы тут было равенство. Вы наделяете функцию BIN() каким то фантастическим свойством. Она работает с потерями, математические инварианты не обязаны сохраняться после этих потерь. Функция BIN() транслирует числа в другую систему координат.
                0.3 при этом не обязано попасть в сумму 0.1+0.2 вот и все.

                А описываемый вами казус тоже про невыполнение математических законов, вот только уже внутри новой системы координат:
                z = x + y
                z - x != y

                Надеюсь так понятно о чем я :)

                • ArtemKaravaev
                  /#22195518

                  Я понял о чём вы, но вы не понимаете о чём я. Я о том, что причина неравенства НЕ В ТОМ, что функция BIN работает с потерей, а в том, что СУММА работает с потерей. Думаю, я всё объяснил, вы уж простите.

                  • tsukanov-as
                    /#22195552

                    У меня и не было проблем с пониманием этого. В комментарии, на который вы отвечаете:
                    z = x + y
                    z - x != y

                    Вот тут сумма с потерей, да.
                    В общем мы друг друга поняли :)

    • maxzhurkin
      /#22196386

      А если встречается деление на 3, использовать троичные, шестиричные или девятиричные?

      • tsukanov-as
        /#22196832

        Если хочется чтобы работала математика соответствующих систем, то да.
        Обратите внимание, что в начале статьи в большой хайповой картинке фигурирует формула 0.1+0.2=0.3, которая к статье не имеет отношения. В статье другие числа, т.к. нету таких чисел в двоичной системе.
        И на вопрос статьи «Но как?» ответ простой: «Никак». А то, о чем на самом деле идет речь в статье — это уже другой вопрос.

        ps 0.30000000000000004.com

  6. defuz
    /#22192926 / +1

    Слишком много лирических отступлений, и шуток и повторяющихся объяснений. В итоге статья обрывается на середине, так с не объяснив, к чему это все было.


    Я рассчитывал по крайней мере узнать как складывать массивы чисел не теряя точность. В итоге узнал только как представить сумму одних чисел суммой других чисел. :/

    • aamonster
      /#22192944

      Ну, как бы есть алгоритм Кэхэна (см, википедию, там объяснение простое), но он не гарантирует абсолютную точность сложения.
      Интересно, будет ли в рамках второй части статьи обобщение алгоритма Кэхэна с хранением массива float (минимальной необходимой длины), а не пары?

    • ArtemKaravaev
      /#22192948 / +2

      Благодарю за отзыв. Однако я прошу принять во внимание два обстоятельства. Во-первых, статья в своём названии уже говорит, что складывать мы будем два числа, во-вторых, про сумму чисел из массива сказано, что это будет в другой статье, там я расскажу, как можно уменьшить погрешность (но не убрать её совсем). Мне трудно представить, как читатель может после этого ожидать, что я расскажу об этом именно здесь. Вы говорите о том, что я много повторяюсь, но вы сейчас показали, что всё-таки я повторяюсь мало. Нужно больше.


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

      • Barbaresk
        /#22193610

        А почему нужна «бесконечная точность»? Если мы не берём в расчёт числа с бесконечными порядками (а ограничимся порядками от -300 до +300, например). И не берём бесконечное число слагаемых, а ограничемся количеством в N (ну, например, 10^15, это несколько петабайт чисел). То для расчёта суммы без потери точности нам достаточно «складывателя» разрядности 300 + 300 + log10(10^15) ~= 615. Естественно, что система счисления «складывателя» должна совпадать с системой счисления вводимых чисем. Иначе, у нас вылезает еще и деление при конвертации. Ну и выходит, что разрядность нашего «складывателя» как O(logK(M)+logK(N)), где M — кол-во допустимых значений чисел, N — кол-во слагаемых, K — система счисления. Разве не так?

        • ArtemKaravaev
          /#22194502

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

      • defuz
        /#22194376

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

        Не очень-то полезный трюк, согласитесь. Самый простой способ представить сумму двух 32-битных чисел в виде двух 32-битных чисел, это не делать ничего вообще.

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

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

        • Xobotun
          /#22194446 / +2

          оказались… два числа той же размерности.

          Я так понимаю, что нет. Если складывать два числа разных порядков, то в результате мы получим два числа: одно близкое к наибольшему порядку, а другое – может быть сильно меньше порядков обоих чисел.


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


          С другой стороны, это может быть и к лучшему: иметь цикл статей, у которых можно прочитать только первую и понять "ага, такой алгоритм есть, он существует, но я им почти никогда не воспользуюсь. Но буду знать". А те, кому интересно изучить дальнейшее использование, могут читать следующие статьи, пока не узнают всё или не остановятся на том уровне детализации, которого им будет достаточно.


          Кстати, на курсере похожим образом иногда публикуются лекции: отрезками по 15-30 минут, затрагивая только одну подтему темы.

          • defuz
            /#22195362 / +1

            Я так понимаю, что нет.
            Как нет, если да? Сумма битов результата равна сумме битов аргументов. Это просто представление a + b в виде некоторого c + d той же разрядности. Такую операцию трудно назвать «суммированием».

            Если бы операндов было больше двух – тогда такой трюк сам по себе был бы полезным. А так я могу разве что посмотреть на сдачу чтобы узнать сколько я потерял.

            Опять же повторюсь, я не сомневаюсь что у приведенной техники есть практическая польза. Но в статье она не показана. Так что ждем следующей статьи.

            • netch80
              /#22195378 / +1

              Опять же повторюсь, я не сомневаюсь что у приведенной техники есть практическая польза. Но в статье она не показана. Так что ждем следующей статьи.

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


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


              Так что польза есть. Большинству, да, не приходится сталкиваться с её применением.

            • Xobotun
              /#22195808 / +2

              А так я могу разве что посмотреть на сдачу чтобы узнать сколько я потерял.

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


              Сумма битов результата равна сумме битов аргументов

              Полагаю, только для чисел в целочисленном представлении. Для float, емнип, суммирование идёт сложнее: fpu расширяет операнды до большей разрядной сетки (скажем, с 32 бит до 40, увеличивая только мантиссу), и складывает обе мантиссы, выровненные по запятой.


              И вот на моменте выравнивания и обратного обрезания до 23 бит происходит потеря точности. Впрочем, я полагаю, это и так всем понятно/известно. :)


              Это просто представление a + b в виде некоторого c + d той же разрядности

              В лучшем случае разрядность не меняется, да. Как в примере в начале статьи 0.1 + 0.2 == 0.3 ± 2.7e-17. Порядок погрешности в 17 раз меньше порядка операндов.
              А в случае сложения сильно разных по порядку чисел, вроде (2^53 - 1) + 2 == (2^53 + 1) ± 1, порядок погрешности сравним с порядком другого операнда, и вот тут бы хотелось увеличить разрядность, чтобы результат влез в сетку.


              Новое число d в этом случае увеличивает общую разрядность суммы два (или меньше?) раза. Да, при сложении двух float можно просто записать результат в double и не мучаться. Но при сложении двух double – в float128? А при сложении двух float128 – в float256? Полагаю, так оно и есть, просто мы вместо введения нового типа данных с удвоенной разрядной сеткой вводим вторую переменную, в которую складируем переполнение разрядной сетки.


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


              Такую операцию трудно назвать «суммированием».

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


              С другой стороны, у классического fadd тоже проблема с суммированием: он возвращает не сумму, а число, близкое к сумме. :D
              У метода из статьи с этой точки зрения больше прав на звание "суммирования", ибо точность результата сильно выше. А два возвращаемых результата – суть одно число, просто записанное в несколько непривычном виде.


              К слову, если я правильно помню, в каком-то x86-ассемблере деление двух целых чисел записывало в два регистра два результата: частное и остаток. Остаток получался побочно, но, поскольку он тоже часто используется, его решили не выбрасывать, а складировать в один из регистров.


              Так что, полагаю, технику из статьи можно обозвать «суммированием с побочным эффектом в виде значения погрешности».


              Самый простой способ представить сумму двух 32-битных чисел в виде двух 32-битных чисел, это не делать ничего вообще.

              Выходит, что так. Не ошибается тот, кто ничего не делает. Вот два числа, лежат в памяти рядом, никакой погрешности! :D

              • ArtemKaravaev
                /#22195862

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

  7. netch80
    /#22193062

    В каком смысле точный, если даже исходные числа 0,1 и 0,2 не имеют точного представления в формате IEEE-754?

    IEEE754-2008 устанавливает и десятичные форматы и методы работы с ними. Считаю важным уточнить, что в статье речь только о двоичных форматах — иначе может создаться ошибочное представление о стандарте в целом.

    • ArtemKaravaev
      /#22193088

      Благодарю за замечание, но когда в черновике я так сделал, то увидел, что получается слишком много лишних уточнений. Мне кажется, что по контексту всё должно быть понятно. Особенно если учесть, что я в некоторых местах оставил это уточнение:


      Начнём с того, что чисел 0,1 и 0,2 в двоичной арифметике с плавающей запятой...

      • Regis
        /#22194544

        У вас в начальной формулировке проблемы ни слова про двоичный формат нет.


        Так что само собой на статью напрашивается ответ "вместо binary64 возьмите decimal64 и терпите издержки, если вам нужны точные десятичные вычисления".

        • ArtemKaravaev
          /#22194546 / +1

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

          • Regis
            /#22194874 / +1

            Но а вообще в десятичном формате будут те же самые проблемы потери точности при сложении.

            Откуда? Первопричина проблемы в том, что над числами, заданными в одной системе счисления выполняются операции в другой системе счисления с внесением ошибки за счет конверсии на входе и выходе. Если весь процесс вычисления делается в той системе счисления, в которой заданы начальные числа, то описанный класс ошибок полностью исчезает.

            • ArtemKaravaev
              /#22194886 / +3

              Оттуда, что при сложении чисел с плавающей запятой в любой системе счисления вы вынуждены смещать одно число относительно другого, когда их экспоненты разные. Результирующая сумма будет занимать больше места, чем позволяет вместить переменная. Значит часть битов будет отброшена. В книге [1] эти вещи, что я описываю в статье, обсуждаются для произвольной абстрактной системы с плавающей запятой, то есть там произвольная система счисления и произвольная, но фиксированная длина мантиссы.

            • netch80
              /#22195212 / +2

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

              Нет, конечно. Вот у вас IEEE754 decimal32 с 7 значащими десятичными цифрами. Вот вы складываете 9999999 и 4. Результат 10000003 не лезет в 7 цифр и округляется до 10000000. Опаньки!


              И вот тут тот самый описанный автором Kahan summation (или его близкий родственник) помогает: вы подсчитали 10000000-9999999, получили 1, вычли его из 4 (второе слагаемое), получили 3 — вот это та самая ошибка округления, которую можно при желании накапливать и которая позволяет вторым таким компонентом фактически удвоить точность расчётов.

      • Regis
        /#22194880

        Кстати, вы не всё исправили.


        Вот тут фактическая ошибка:


        К сожалению, это данность, от неё никуда не уйти, если хранить числа в типе данных double (или любом другом типе фиксированного размера из Стандарта IEEE-754).

        decimal64 — тип фиксированного размера из Стандарта IEEE-754, но у него озвученной проблемы нет.

        • ArtemKaravaev
          /#22194894 / +2

          Есть, попробуйте в decimal64 ввести число, записанное в 13-тиричной системе счисления. В ряде случаев будет возникать бесконечная дробь. Данная проблема есть почти всегда, в книге [1] есть теорема, в которой написано при конвертировании из каких систем счисления с какие не будет такой проблемы.

  8. byman
    /#22193092

    Интересная для меня тема. Недавно пришлось писать DSP библиотеку для процессора, так совсем замучился с float/double числами. Как только используешь SIMD вычисления, сразу результат отличается от результата не SIMD алгоритма. И порой непонятно то-ли баг в программе, то-ли это такая разбежка в результатах :) И какой результат все-таки точнее. Боюсь, что такие точные вычисления сильно убавят производительность. Жду следующую статью :)

    • netch80
      /#22193130

      Как только используешь SIMD вычисления, сразу результат отличается от результата не SIMD алгоритма.

      Вы в 32-битной Windows или Unix пишете? Там по умолчанию плавучка на FPU, а у него может стоять режим, в котором мантисса 64 бита. Но при любом режиме порядок не усекается до конверсии в целевой формат.
      Там даже вот такие варианты бывают (казалось бы, пусть 0.1+0.2 считается неточно — но должно считаться всегда одинаково? ан нет, момент конверсии влияет).


      Уточните режим компиляции не-SIMD варианта, проверьте через MPFR (дорого, поэтому на мелком примере) — тут много есть о чём подумать...

      • ArtemKaravaev
        /#22193148 / +1

        Это математическая статья, описываемые явления будут во всех без исключения форматах с плавающей запятой фиксированного размера не зависимо от компилятора, и если формат допускает денормализованные числа. Смысл статьи не в том, чтобы бороться с врождённой проблемой формата и как-то складывать десятичные дроби, а в том, чтобы находить сумму двух чисел без потери точности в уже заданной математической абстракции под названием двоичная арифметика с плавающей запятой. Просто я привязал её к формату IEEE-754, так как он описывает как раз такие числа и всем известен.

        • netch80
          /#22193184

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

          Против этого никто не возражал. Но byman@ задал вопрос "почему у него с SIMD результаты другие, чем без SIMD", вот я и подсказал, куда копать.


          и если формат допускает денормализованные числа

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

          • ArtemKaravaev
            /#22193192

            А, прошу прощения, я не сообразил, что вы именно ему отвечаете :) Опять запутался.


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

      • byman
        /#22194560 / +1

        Я не типичный юзер :) у меня свой VLIW in-order DSP с SIMD возможностями. Никаких ОС. Ряд библиотечных функций приходится писать самому. При реализации различных стандартных функций есть возможность за те же такты вычислить точнее. Также порой нужен какой-то более точный (пусть и медленный) алгоритм, чтобы оценить погрешность более скоростного вычисления. Отсюда и интерес к подобным статьям :)

  9. NumLock
    /#22193178

    Задали студентам написать вычисление матрицы методом Гаусса в питоне. Одна студентка использовала библиотеку дробей и производила все операции в виде дробей. Это было самое точное вычисление матриц в классе.

    • netch80
      /#22193196

      Для учебного примера — отлично :)
      Если же элементы матрицы будут поступать из реальных измерений — станет бесполезно. Поэтому такой метод хорош для частной задачи, но грубо плох методически.
      Реальные библиотеки иногда делают для минимизации эффектов округлений итеративные коррекции результата.

  10. Antervis
    /#22193244

    а не практичнее например считать в long double а потом приводить к double?

    • ArtemKaravaev
      /#22193252

      Это не будет работать по нескольким причинам. Во-первых, long double уже де-факто отменили (в потоковых обработках), во-вторых, он даёт только 11 дополнительных битов мантиссы, а чтобы сохранить погрешность нужно 52 (для double), в-третьих, возникает ошибка двойного округления… короче, ещё 10 причин я просто не хочу писать, вы уж простите. Мне кажется, вы не вполне поняли проблему.

  11. AVI-crak
    /#22193266 / +1

    DECFLOAT спасёт ситуацию.
    А в реальности, первые калькуляторы считали в десятичной системе, результаты проверяли на бумажке (как умели) — и они совпадали.
    Сейчас все повально обленились, и даже приложения «калькулятор» для пк использует двойную точность — отчего результат не всегда совпадает с бумажкой.

    • VolCh
      /#22193798

      1/33 давало 1, а 0.333333333 3 давало 0.999999999?

      • 411
        /#22194358 / +2

        Судя по всему хабр съел ваши знаки умножения. По крайней мере в мобильной версии сайта у меня.

        • netch80
          /#22195218 / +1

          Можно определить по тому, где начинается и кончается курсив:


          1/3*3 давало 1, а 0.333333333*3 давало 0.999999999?

        • VolCh
          /#22195372

          Да, съел :( на мобиле сложнее корректно форматировать чем на сайте. Выше правильно рвзъяснено

  12. agalakhov
    /#22193292

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

  13. Ontaelio
    /#22193296 / +1

    Статья хороша и познавательна в плане математики (и, конечно, очень жду продолжения). Но в ней не хватает простого практического совета: если хочешь точности (а также скорости и компактности), никогда не пользуйся переменными с плавающей запятой. 0.1 == 1/10. Все эти float, double и иже с ними созданы для того, чтобы облегчить жизнь программиста, когда не требуется идеальная точность вычислений.

    • ArtemKaravaev
      /#22193310

      Благодарю за отзыв. Но вот в каком непростом положении я нахожусь. Дело в том, что у каждого читателя своя математическая культура, и подстроиться под каждого невозможно. Например, в моём научном коллективе и без слов понятно, когда какой тип данных применять (рациональные числа, с плавающей запятой, с конечной или бесконечной точностью, позиты Джона Густафсона, интервальную арифметику, числа половинной точности или увосьмерённой). У нас все и так знают, что к чему и почему. У других людей может не быть этих знаний, но есть другие знания. Как я могу угадать, какие у них знания? Напишу про точность, мне ответят, что я мог бы этого не писать. Если не напишу, ответят, что надо было написать. Если напишу категорично (никогда не пользуйся), ответят, что иногда можно и даже нужно, если не напишу, то скажут, что надо было написать. Поэтому я написал только по делу. Кто в теме, тому статья полезна будет, кто не в теме, тем я порекомендовал для начала ресурсы [3-5].

      • Ontaelio
        /#22193428

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


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


        Но с практической точки зрения (а Хабр читают практики) важно помнить, что это именно обертки, под которыми скрываются банальные и совершенно неделимые нули и единицы. И правильный (ИМХО) подход — оперировать нулями и единицами как они есть. Это я к тому, что в статье упоминается именно тип double, а не вообще дробные числа.


        Ну, как-то так. Я ни в коем случае не критикую статью, повторюсь, она очень хороша и познавательна.

        • ArtemKaravaev
          /#22193462

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


          Я взял тип данных double просто для единообразия описания, если вы откроете книгу [1], то увидите, что там знания даются с позиции вообще полностью абстрактной системы с произвольным основанием (не только двоичным или десятичным) и произвольной, но фиксированной точностью. Там есть вырвиглазные теоремы и я хочу их как-то обойти. Вот для этого взял и свёл всё к double, понятно, что для float то же самое будет, только битов меньше.


          А что касается причин появления разных обёрток, то она намного сложнее, чем вы сейчас пишите. Но я бы не стал дискутировать об этом в комментариях, с вашего позволения :) Это очень длинная тема, и мы уйдём в философию, утомив других читателей.

          • Ontaelio
            /#22193492

            А что касается причин появления разных обёрток, то она намного сложнее, чем вы сейчас пишите.

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

        • VolCh
          /#22193844

          Есть только нули и единицы, других чисел в принципе нет (замнём про троичные, десятичные и прочие ичные компьютеры для ясности). А есть различные интерпретации этих последовательностей. float и integer просто разные интерпретации, уместные в разных случаях. Есть и другие интерпретации, включая строки. Вон JS обходится (почти) без integer. JS — исключительный случай? Нет, самый популярный язык по некоторым метрикам.


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

      • defuz
        /#22195340

        позиты Джона Густафсона, интервальную арифметику
        Расскажите про свой опыт использования позитов и интервальной арифметики. Если у вас есть примеры из практики – вообще отлично будет.

        • ArtemKaravaev
          /#22195376 / +1

          Примеры не дам, исследования нашего научного центра не подлежат разглашению. Могу только сказать одно — они не работают правильно в задачах, связанных с исчислением бесконечно малых. Мы завалили их на простых тестах, где нужно было: интегрировать, решать дифференциальное уравнение численно, применять другие численные методы, где нужно было работать с бесконечно-малыми (имеется в виду, конечно, просто с очень маленькими числами). Алгоритмы, которые хорошо сходятся для чисел формата IEEE-754 на позитах не сходились и уходили в бесконечность, либо плясали в широком диапазоне около ожидаемого ответа.


          Интервальная арифметика часто выдает очень правильный, но бесполезный ответ, например, что число будет в интервале от 0 до бесконечности. Далее, возникают интересные ситуации, когда интервал [A, B] нужно поделить на интервал [A, B] (ответ, естественно, не единица), особенно когда A и B имеют разные знаки. Больше сказать ничего не могу, увы.

          • defuz
            /#22195414

            Мы завалили их на простых тестах
            Вы использовали только posit или quire аккумуляторы? Мне кажется, уходить в «бесконечно-малые» и «бесконечно-большие» на posit просто нельзя – они имеют крайне низкую точность на очень больших и очень маленьких экспонентах. Универсальное правило – все значения, которые представляются в posit должны быть нормализованы (иметь экспоненциальный порядок близкий к 0). А вот промежуточные значения должны хранится в quire регистрах с весьма высокой точностью.
            бесполезный ответ, например, что число будет в интервале от 0 до бесконечности
            Мне кажется, что это полезный ответ, в том смысле, что если бы те же вычисления делались с применением обычной арифметики аналогичной разрядности, то погрешность результата была бы «бесконечной». Или вы говорили конкретно о Густафсоновской арифметике?

            • ArtemKaravaev
              /#22195544 / +1

              Нет, quire аккумуляторы там не работают, ошибка была ведь не в том, что мы в результате каких-то накоплений потеряли точность, а именно в результате ошибок в единичных операциях с очень маленькими числами. Плотность маленьких чисел у позитов настолько низка по сравнению с плотностью денормализованных чисел в IEEE-754, что это и неудивительно. А если прямо совсем всё делать в quire, то тогда смысла нет было создавать операции сложения и вычитания для обычных позитов. Но они ведь зачем-то есть, и они работают плохо в задачах математического анализа. В этом проблема. Не знаю как для других исследователей, для нас это однозначный приговор позитам. По крайней мере, на сегодня точно. Если нужна хорошая точность, мы возьмём float128, например.


              Мне кажется, уходить в «бесконечно-малые» и «бесконечно-большие» на posit просто нельзя

              Да, нельзя. В этом и проблема. И инструмент quire задачу решает не лучше, чем если вместо double я возьму float128, например, или, если будет нужно, более крупный тип данных.


              Мне кажется, что это полезный ответ, в том смысле, что если бы те же вычисления делались с применением обычной арифметики аналогичной разрядности, то погрешность результата была бы «бесконечной».

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

      • AnthonyMikh
        /#22195586

        Например, в моём научном коллективе и без слов понятно, когда какой тип данных применять (<...>, позиты Джона Густафсона, <...>)

        А где реально позиты на практике используют?

        • ArtemKaravaev
          /#22195622

          В случае работы с числами, близкими к единице они хорошо себя показывают. Но вообще я их просто так в этом списке упомянул, чтобы пояснить собеседнику, что я в курсе существования разных типов данных, в том числе весьма экзотических.

    • netch80
      /#22193350

      если хочешь точности (а также скорости и компактности), никогда не пользуйся переменными с плавающей запятой. 0.1 == 1/10.

      Правильно говорить — если нужна абсолютная точность (которая отличается от относительной).

  14. Tim_23
    /#22193348 / +1

    Сначала вроде было интересно и интригующе, но далее смысл как-то растерялся.
    Какой практический вывод должен следовать из статьи?
    Я к примеру инженер. пишу какой-то код, также работаю с массивами и с десятичными дробными данными, циклами на сотни тысяч итераций. Думать мне об этих погрешностях или нет? Если они такие маленькие, то практического смысла от них в реальной жизни мало — погрешности при переходе от виртуального расчета в физический мир реальных объектов съедят эти несостыковки.
    Может быть для долговременных систем (типа космических аппаратов летящих 30 лет со второй космической) или для процессора это играет роль?
    Поэтому практической значимости не увидел, а если это академический интерес — то против ничего не имею.

    • ArtemKaravaev
      /#22193370

      В вашем случае думать об этом не надо. Я упомянул в статье, что это касается только тех, кто пишет библиотеки для языков программирования, то есть не при помощи языков программирования, а именно на стадии создании самого языка. Там в этой сфере много всяких таких алгоритмов, и они реально там нужны. По этим темам в книгах выделяют целые главы, поэтому практический смысл есть. Но если вы его не видите, значит именно вам он не нужен. Но вы можете не подозревать, что ваш любимый язык программирования работает правильно только потому, что где-то когда-то люди изобрели подобные алгоритмы, о работе которых вы никогда не узнаете, но работа которых миллионы раз выполняется когда вы нажимаете кнопку Build.

      • Tim_23
        /#22193374

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

        • ArtemKaravaev
          /#22193384

          В следующей статье я собирался показать как разные способы суммирования чисел из массива могут давать погрешности, отличающиеся на несколько порядков. Эти способы, в частности, работают с применением знаний из этой статьи. Если торопитесь изучить эту сферу, то в книге [1] всё есть, я лишь немного от неё отступлю в следующей статье, а по сути расскажу то же самое.

        • cyberly
          /#22193552

          Погрешности — это одно, а вот с равенствами все становится интереснее. Буквально на прошлой неделе тикет был. JavaScript, валидатор формы. Одна из проверок — значение в поле A должно быть суммой значений в полях B и С. Работало несколько лет, пока все числа имели ровно 2 знака после запятой.

          console.log(0.01 + 0.02); 
          console.log(0.1 + 0.2);
          

          получаем:
          0.03
          0.30000000000000004
          


      • 411
        /#22194366

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

        Возникает закономерный вопрос: а как проверить, что мой язык программирования работает в указанном вами случае (или в каком-то другом)?


        Я вот полагаю, возможно ошибочно, что из примера a+b=s+t в большинстве языков программирования при попытке вывода s я получу именно s, а не s+t

    • SmallSnowball
      /#22193744 / +2

      Ошибки имеют свойство быстро накапливаться до таких состояний, когда игнорировать их не получается. Например, это чуть ли не основная причина, почему в CAD системах на данный момент нет по настоящему робастных булевых операций, хотя внутренняя точность самих моделей может достигать пары десятков нанометров. Реальный пример из практики — флотовые ошибки в булевых операциях приводили к тому, что при точности исходных деталей в CAD 1.0e-10 метров (ну т.е. десятая доля нанометра, никто такой точности детали даже произвести не сможет) погрешности в расчетной сетке составляли уже несколько сантиметров, что было абсолютно неприемлимо.
      Тем более, если у вас вычисления на сотни тысяч операций, очень легко внезапно получить погрешность одного порядка с результатом.

      • Tim_23
        /#22193936

        В расчетной сетке чего? Имеется ввиду численный анализ? Если так, то построение точной расчетной сетки — отдельная проблема. Точность в данном случае условная. Масштабировать объект или его части для расчета можно как угодно, разве тут вылезет погрешность CAD системы?

        • SmallSnowball
          /#22193978

          Имеется ввиду CFD, там общий процесс выглядит как «собрать из деталей CAD-сборки b-rep модель для анализа -> построить сетку на модели -> численный анализ на сетке» Первый этап либо делается вручную, либо отдается автоматике, которая просто мержит все активные детали в сборке. И вот уже на этом первом этапе уже можно всю точность растерять. Потом на этой модели будет строиться сетка, но какой в этом толк, если исходная модель уже со значительной погрешностью построена. Да и в принципе строить сетку на моделях, где кривые пересечения поверхностей отстоят от самих поверхностей на сантиметры — само по себе занятие не из простых. Скалирование моделей не сильно помогает, т.к. погрешность будет обусловлена постепенной потерей значащих битов в мантиссе.

          • Tim_23
            /#22194080

            По своему опыту скажу, что перед CFD анализом мы сильно упрощаем исходные сборки и модели ( с учетом принятых допущений, и их значимости), чтобы в том же meshing или icem сделать более менее адекватную и не переразмеренную сетку. Строить сетку в лоб из сборок как-то неэкономично.
            Потом вопрос не в сантиметрах и метрах, в относительных размерах и числах подобия типа Рейнольдса. Какой смысл брать реальную сборку и считать сантиметры, если к примеру мы не моделируем погран слой, а смотрим общее нагружение и интергральную картину? Не вижу тут нерешаемых проблем.
            А если это погран слой, который нужно аккуратно моделировать с учетом всяких выступов, то всякое моделирование все равно будет не точным, так как реальный объект и виртуальная модель далеко не одно и тоже.

            • SmallSnowball
              /#22194128

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

              • Tim_23
                /#22196010

                Я вот из того что сразу вспоминаю, так это то что иногда в автокаде, якобы совпадающие границе при сильном приближении не совпадают. Причем этот глюк странный. Вроде по привязке все точка в точку, а приближаешь две границы стыкующихся объектов начинают сползать туда-сюда, причем тыкание мышкой на таком зуме уже не приводит к результату. И вот думаешь, это реально или погрешность автокада. Здесь правда есть хитрости. Например в автокаде же есть функция порога при котором эти зазоры не должны учитываться и к примеру стык будет учтен и тела можно слить в одно.
                Также, если брать пример CFD, например в ансис Spaceclaim есть всякие фишки для ремонта «плохой» геометрии. Правда результат лишь гарантирует, что тело будет и поверхности не пересекаются, про точность тут уже трудно что-то сказать.

                • SmallSnowball
                  /#22196250

                  Разъезжание границ довольно просто объясняется, если я правильно ситуацию понял. Для визуализации brep триангулируется, если 2 тела будут триангулироваться независимо друг от друга, то точки триангуляции будут немного отличаться, даже если контачащие поверхности абсолютно одинаковые. Происходит это потому что при триангуляции вершин и ребер нужно учитывать все соседние грани, а они на разных телах разные обычно. Плюс при визуализации этой триангуляции координаты точек обычно из даблов перегоняют во флоты, т.к. gpu гораздо быстрее с флотами работает, а над просмотрщиком моделей в координатах двойной точности не заморачиваются. А флоты начинают сами по себе «дрожать» при зуме, если зазумиться тысяч в 10 раз, т.к. при таком зуме матрица проекции уже «плохая», если в одинарной точности считать.
                  Залепухи и прочий ремонт плохой геометрии очень хорошо работает, если хорошо понятно, где его лепить. Иногда ошибки проявляются слишком далеко от места, которое их вызвало, и автоматика не справляется. Особенно много проблем получается, если рядом с плохим куском находится какая-то мелкая деталь, тогда «плохая» топология может сожрать всю мелкоту рядом. Иногда такой результат устраивает, иногда нет.

    • VolCh
      /#22193874

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

      • Tim_23
        /#22193938

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

        • ikle
          /#22194536 / +2

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


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


          Из такого месячного вычисляем реальный (сложный) годовой, из реального годового (сложный) подневный. И уже далее считаем: у меня получались очень близкие цифры к банковским и я на этом остановился.


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


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


          P.S. Смею предположить, что либо более точная методика закреплена в каких-либо подзаконных актах, либо можно запросить точный метод расчёта у банка. Знающие люди пусть подскажут.

          • edo1h
            /#22195086

            По-моему, банкиры тут явно лукавят, чтобы показывать меньший годовой процент

            всё определено законом же
            http://www.consultant.ru/document/cons_doc_LAW_155986/e52bee2d092465172cd750d5a23927f45bb3d017/

            • ikle
              /#22195162 / +1

              Спасибо за ссылку, добавил в закладки на будущее.


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

              Согласно вашей ссылке:


              ПСК = i x ЧБП x 100

              Ну вот только ПСК — это не годовая ставка, как нам преподносят в рекламе, годовая ставка (в процентах) будет ((1 + i)^ЧБХ — 1) ? 100. Простой пример, ЧБХ = 12, 1% в месяц, то есть i = 0.01:


              • ПСК = 0.01 ? 12 ? 100 = 12 — нам объявляют 12% годовых, но реально же
              • ((1 + 0.01) ^ 12 — 1) ? 100 ? 12.6825%.

              ПСК — это то, что в сумме мы отдадим за год (в процентах), но это не годовая процентная ставка, годовая равна ПСК при ЧБХ = 1 (или i = 0). Так что лукавят.


              Впрочем, на практике ПСК, скорее, полезнее: при расчёте годового бюджета нам нужен ПСК, а не годовая ставка.

  15. Rutel_Nsk
    /#22193392

    Может дробные числа хранить в виде двух целых чисел — целая и дробная часть + основание дроби в типе переменной?
    А работать с ними по принципам похожим на BCD числа (операция и далее коррекция результата)?

    • AllexIn
      /#22193578

      Я с вами согласен. Но где вы предлагаете тип переменной хранить?

      • Rutel_Nsk
        /#22193652 / +1

        А зачем хранить тип переменной, в современном программировании тип переменной
        задается только в тексте программы.
        BCD_Float A,B; // И больше ничего ненужно

        • AllexIn
          /#22193670

          То есть вы хотите поддержку на уровне компилятора? Это немного другая задача.

          • Rutel_Nsk
            /#22193684 / +1

            Создать свой тип данных и переопределить операции вполне позволяет обычный С++

            • AllexIn
              /#22193700

              Тогда уж шаблоны делать, тогда не придется городить вагон классов под каждый вариант основания.

              • Rutel_Nsk
                /#22193732 / +1

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

  16. firegurafiku
    /#22193506

    По-моему, вы изобрели «double-double arithmetic»: хранение результата сложения в виде пары вещественных чисел (s, t) по факту удваивает количество бит мантиссы результата, то есть от чисел двойной точности мы переходим к числам четырёхкратной точности.


    Из-за этого заголовок статьи вызывает у меня некоторое недоумение: «Сложение двух чисел с плавающей запятой без потери точности» невозможно, если и аргументы, и результат имеют одинаковое количество бит для хранения мантиссы. Лучшее, чего можно добиться в этом случае, — корректного округления результата к ближайшему представимому числу, ровно этого стандарт IEEE 754 и требует от реализаций основных арифметических операций.


    А ещё мне не понятно, как именно автор собирался сделать равенство из печально известного неравенства 0.2 + 0.1 != 0.3.

    • VolCh
      /#22193900 / +1

      Если слагаемых много, то, как я понимаю, нужна только одна переменная для накопления потерянной точности. И добавление её к итоговой сумме может дать меньшую ошибку в итоге чем лишь округление на каждом шагу.


      По сравнению, наверное, будет 0.2 + 0.1 — 0.3 = 0

    • netch80
      /#22195226

      По-моему, вы изобрели «double-double arithmetic»:

      Описывает, да.


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

      Если учесть, что представления чисел уже округлены до подходящего значения на разрядной сетке float или double — то формально всё честно.


      А ещё мне не понятно, как именно автор собирался сделать равенство из печально известного неравенства 0.2 + 0.1 != 0.3.

      Равенство и не делается, но ошибка накапливается отдельно.

  17. AllexIn
    /#22193570

    Если у нас два дабла чтобы хранить один дабл — то нам нафиг не нужна плавающая точка и даблы вообще.
    Берем две целочисленных переменных того же размера что и даблы и получаем fixed point без обозначенных проблем с огромной точностью.

    • AnthonyMikh
      /#22195612

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

      • AllexIn
        /#22195652

        16 битный fixed point никогда не проиграет 8 битному float по точности.
        И про «данные одного порядка» — у float ровно таже проблема. Вы либо храните большие числа с низкой точностью, либо маленькие с высокой.

        • netch80
          /#22195804

          16 битный fixed point никогда не проиграет 8 битному float по точности.

          Это если вы твёрдо знаете, как именно он должен быть fixed (позиция точки), и если вам нужна именно абсолютная точность.
          Да, полно контекстов, где это так и нужно (те же финансы). Но полно и таких, где это не проходит.


          Вы либо храните большие числа с низкой точностью, либо маленькие с высокой.

          Вы всё время циклитесь на абсолютной точности.
          Когда мы меряем расстояние до какой-нибудь Тау Кита, нам плюс-минус миллион километров до лампочки. Зато в позиционировании транзистора в процессоре играют роль нанометры.

  18. Busla
    /#22193576

    Очевидно, что t — не ошибка операции суммирования, а накопленная ошибка по всем трём операциям.

  19. edo1h
    /#22193728 / +1

    Десятичный результат получен с помощью правильного онлайн-конвертера, а дробь посчитана с помощью библиотеки MPIR и функции mpq_set_d, которая умеет переводить double к рациональной дроби.

    ИМХО не хватает такой же дроби для результата (и, возможно, краткого пояснения почему получили не 0.3000000000000000166533453693773481063544750213623046875)

    • ArtemKaravaev
      /#22194512

      Это есть в нашем учебном курсе [5]. Я в начале статьи предупредил, что читатель должен быть в курсе этого примера, чтобы мне не пояснять вещи, не относящиеся к теме.

  20. amarao
    /#22193818

    Остаток от сложения по основанию 2.


    /Я тролль? Или нет?

  21. VolCh
    /#22193920 / +1

    Давно такого не было: только в одном докладе на highload fwdays упомянули мельком о проблеме и подобном способе её решения, подумал завтра погуглить на свежую голову, а пока Хабр почитать и тут же первая статья на эту тему. Спасибо! Хабр торт!

  22. nikitazvonarev
    /#22194090 / +4

    Эта тема очень нужная, и имеет практические применения. Так называемые error-free transformations позволяют решить некоторые задачи с повышенной точностью без применения чисел с плавающей точкой большей разрядности, например, вычисление скалярных произведений или вычисление многочлена. Зачастую вычисление значения полинома возле его корня (особенно кратного корня) плохо обусловлено (т.е. ответ на такую задачу трудно искать с малой относительной ошибкой), а такие вещи порой нужно делать; у меня подобные знания нашли применение при обработке сигналов, если точнее — при обращении матрицы-циркулянта. Кому интересно — поищите статьи по Compensated Horner Scheme, автор Graillat et al.

    На современных процессорах x64 (Intel, AMD), кстати, есть инструкция FMA (fused multiply-add), которая ускоряет подобные алгоритмы. По скорости получается, конечно, медленнее, чем без compensated алгоритмов (с обычной точностью), но намного быстрее, чем с более длинными типами (длиннее double). Кстати, никогда не компилируйте программы с compensated алгоритмами с флагом -ffast-math – вам GCC «эквивалентными» преобразованиями всю точность убьёт =)

    • ArtemKaravaev
      /#22194516

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

  23. IronSwan
    /#22194520 / +1

    Ух, помню на втором курсе универа препод чуть ли не выгонял людей с пары, когда слышал неправильное произношение стандарта. Не «АЙ ИИИ 754», а «АЙ-трипл И». Навсегда запомнил)

    • ArtemKaravaev
      /#22194528

      Это ещё ладно, я слышал как-то "И ЙЕ ЙЕ ЙЕ — 754" :) Сразу выключил видео, на котором это было.

      • netch80
        /#22195290

        Значит, будем требовать "ИИЭР 754". Локализовать так локализовать :)

  24. citreraydo
    /#22194522

    Я понимаю, что исходная статья про стандарт IEEE, а не про конкретные языки и компиляторы, но всё же у меня возникает конкретный вопрос. Не объявит ли компилятор C++, что в программе происходит UB, если поймает программиста на том, что он не уверен, что t=0 (и не разнесёт ли всё в результате, ну а чё, стандарт говорит можно же)?

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

    • nikitazvonarev
      /#22194974

      UB тут ни при чём, но «разнести всё» компилятор может, если укажете флаг -ffast-math (об этом уже писал выше). Точность сразу упадёт. Поэтому да, программы с такими трюками и этот флаг жить вместе не могут :)

      • citreraydo
        /#22196442

        Да, спасибо, просто вашего комментария ещё не было, когда я писал свой, он долго проходил цензуру.

        Похоже, это -ffast-math (согласно man gcc) даже в -O3 не входит, так что всё вроде не так плохо.

    • Pand5461
      /#22195050 / +1

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


      Как сказали уже, можно включить флаг fastmath, позволяющий делать "небезопасные" оптимизации. Ряд SIMD оптимизаций, например, попадают в эту категорию.

  25. DarkWolf13
    /#22196418

    в свое время имея несколько ограниченное время, решил подобную задачу используя самописную библиотеку операций с большими числами, где числа приводились в не дробный вид для выполнения операций (сложение, вычитание, деление, умножение) и судя по конечному результату и просматривая аналогичные решения, решение оказалось достаточно жизнеспособным…

  26. citreraydo
    /#22196450

    Ещё один вопрос. Я правильно понимаю, что если я хочу точно сложить n чисел типа double, где n < (разница макс и мин двоичного порядка)/(число двоичных знаков в мантиссе), то мне в любом случае потребуется n переменных типа double для точного хранения результата, двумя я уже не обойдусь?

    • ArtemKaravaev
      /#22196704

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

  27. MinGW
    /#22196702

    Спасибо огромное!