Отзывы к статье: IEEE 754 - стандарт двоичной арифметики с плавающей точкой ******************************************************************************************************************* От кого: Владимир Юровицкий Тема: Яшкардин Владимир от Владимира Юровицкого Дата: Thu, 12 Aug 2010 19:12:01 +0400 Кому: Уважаемый Владимир. Я случайно попал на Ваш сайт. И мне было приятно увидеть ссылку на мои работы и то, что вы разделяете мою обеспокоенность. В настоящее время я веду работу по переходу в новый формат чисел. См. на моем сайте. Ваша поддержка для меня имеет ценность. Всего наилучшего Владимир Юровицкий к.э.н., доцент МФТИ, с.н.с. НОВЦ РГСУ член Международной академии информатизации лауреат ТОП-100 "ТВОРЦОВ ПОСТСОВЕТСКОГО ПРОСТРАНСТВА - 2009" по версии экспертного сообщества "Global Intellect Monitoring" "за вклад в развитие теории денег" (http://www.polit.nnov.ru/2010/03/03/thinker100creat2009/) вместе с В.Путиным, Д.Медведевым, Н.Назарбаевым, М.Саакашвили, Глазьевым и другими. web www.yur.ru, vladyur.livejournal.com mail: vlad@yur.ru ******************************************************************************************************************* От кого: "Grigoriy A. Sitkarev" Отправитель: "Grigoriy A. Sitkarev" Тема: Статья о IEEE-754 - дополнения Дата: Sat, 13 Nov 2010 22:02:56 +0300 Кому: info@softelectro.ru Приветствую. Статья ваша обновилась, как я заметил недавно, и в общем, выглядит хорошо. Но я заметил ряд неточностей, в частности, ваш пример с фасовочной машиной и бункером. Дело в том что так вычитать/складывать совершенно точно нельзя. Я вам написал пример и тест чтобы его проверить, надеюсь что вы каким-то образом опубликуете этот код а также переделаете свой пример на Visual Basic чтобы всё встало на свои места. Кроме этого, я думаю что вам стоит привести алгоритм сравнения чисел с плавающей точкой, который я вам присылал (функция fcmp() там была). Потому что если мы сравниваем два числа, одно из которых известно заранее а второе - результат вычислений, последнее может быть верно лишь с некоторым приближением, выраженного через машинную точность (машинное эпсилон). Это как раз случай, который у вас в самом конце - сравнение с нулём. Сравнивать такие числа оператором `==', как в Си, нельзя - только через fcmp(). Пример компилируете компилятором Си (можно и с оптимизацией, она не меняет порядок операций в выражениях арифметических, это даже в стандарте жёстко закреплено). Я думаю что алгоритм вам будет ясен, и вы прокомментируете его у себя тоже. sitkarev@illo:~$ gcc -O2 -Wall -o subcomp subcomp.c sitkarev@illo:~$ ./subcomp Left in hopper: 200.00000 kg Held in bunker: 100.00000 kg Обратите также пожалуйста внимание что там используется одинарная точность, и этого достаточно. Здоровья и терпения вам! -- Гриша #include #include #include struct acc_comp { float value; float compens; }; void sub_compens(struct acc_comp *acc, float quantum) { float tmp, c; tmp = quantum - acc->compens; c = acc->value - tmp; acc->compens = acc->value - c - tmp; acc->value -= tmp; } void sum_compens(struct acc_comp *acc, float quantum) { float tmp, c; tmp = quantum - acc->compens; c = acc->value + tmp; acc->compens = c - acc->value - tmp; acc->value += tmp; } void sub_test() { struct acc_comp hopper; struct acc_comp bunker; float tablet; int n, i; n = 10000000; hopper.value = 300.0; hopper.compens = 0.0; bunker.value = 0.0; bunker.compens = 0.0; tablet = 0.00001; for (i = 0; i < n; i++) { sub_compens(&hopper, tablet); sum_compens(&bunker, tablet); } hopper.value -= hopper.compens; bunker.value += bunker.compens; printf("Left in hopper: %04.5f kg\n", hopper.value); printf("Held in bunker: %04.5f kg\n", bunker.value); } int main(int argc, char *argv[]) { sub_test(); return 0; } ******************************************************************************************************************* От кого: "Grigoriy A. Sitkarev" Отправитель: "Grigoriy A. Sitkarev" Тема: Статья - ошибка в программе при вычислениях с плавающей точкой Дата: Thu, 30 Sep 2010 15:25:19 +0400 Кому: info@softelectro.ru Приветствую. Дорогой коллега, к сожалению в вашей программе, во втором примере, есть ошибка. Я не программирую на Basic, пользуюсь обычно Си, т.к. этот язык особенно хорошо подходит для вычислений с плавающей точкой и системного программирования. Дело в том что в Basic также как и в Си, константы с плавающей точкой по умолчанию имеют тип double, т.е. двойную точность. В Си есть оператор приведения типов который позволяет к примеру привести число с двойной точностью к одинарной. Поэтому правильно бы ваша программа должна была выглядеть так: #include #include #include int main(int argc, char *argv[]) { float a, b, c, d; a = 1.0; b = 3.0; c = a / b; d = (c - 1.0f/3.0f) * 1.0e9f; printf("Result: %08.2f\n", d); return 0; } В моём примере не совпадают имена переменных `c' и `d' с вашим, я их именую последовательно по алфавиту. Здесь я константы 1.0 и 3.0, имеющие по умолчанию двойную точность, привёл к одинарной точности. Поэтому результат вычисления будет правильным. sitkarev@illo:~$ ./fp Result: 00000.00 Если бы я не приводил константы к типу float, то у меня также была бы ошибка, как и в вашем примере: #include #include #include int main(int argc, char *argv[]) { float a, b, c, d; a = 1.0; b = 3.0; c = a / b; d = (c - 1.0/3.0) * 1.0e9; printf("Result: %f\n", d); return 0; } sitkarev@illo:~$ ./fp Result: 9.934108 Если в Basic есть возможность приведения типов, то после того как вы приведёте константы к одинарной точности тоже получите верный результат. Если же привести типы нельзя, то получается что константа 1.0/3.0 будет посчитана с двойной точностью, а `a / b' с одинарной. Отсюда ошибка, которая потом при умножении на большое число оказывается значимой. В тексте статьи также есть неточности. Например, вы цитируете кого-то, кто высказал предположение что в принятии и обсуждении IEEE-754 "не принимали участия профессиональные математики", так вот именно профессиональные математики принимали участие в разработке этого стандарта и в первую очередь именно их пожелания там учитывались. В частности, профессор университета Торонто, Вильям Кахан. Рекомендую вам посетить его странцу и почитать его материлы, очень много интересного. http://www.cs.berkeley.edu/~wkahan/ А также страницу с материалами группы IEEE 754. http://grouper.ieee.org/groups/754/reading.html Чтобы для вас стали ясными некоторые нюансы при вычислениях с плавающей точкой, я бы рекомендовал вам документ "What Every Computer Scientist Should Know About Floating Point Arithmetic". http://dlc.sun.com/pdf/800-7895/800-7895.pdf В примерах аварий, которые вы привели с сайта В. Юровицких, проблемы были вызваны вовсе не ошибками в IEEE-754 арифметике, а в неправильном использовании и программировании а также переполнении в целочисленной арифметике. По этой ссылке изложены подробно причины аварий, прочитайте пожалуйста. http://www.ima.umn.edu/~arnold/disasters/disasters.html Да, конечно IEEE-754 не идеален, и требует соответствующего обучения специалистов, чтобы они могли анализировать источники ошибок в вычислениях и были аккуратны там где этого требуют особенности этого механизма. На то что сегодня, в 21-м веке, мы откатились на 50 лет в этом направлении, сетует и сам Кахан в своей записке "Back to the Future of Undebuggable Floating-Point Software in Science and Engineering". http://www.cs.berkeley.edu/~wkahan/BASCD08K.pdf Нам надо учиться лучше понимать то что сделали наши предшественники, и не боятся критиковать их с конструктивных позиций и учётом реалий современных вычислений. Большое спасибо вам за вашу статью. -- С уважением, Ситкарев Григорий ******************************************************************************************************************************************* От кого: info@softelectro.ru Тема: Re: Статья - ошибка в программе при вычислениях с плавающей точкой Дата: Thu, 30 Sep 2010 17:00:02 +0400 Кому: "Grigoriy A. Sitkarev" Здравствуйте Григорий! Спасибо за отзыв о статье ieee754 Я чуть позже отвечу на ваше письмо подробнее. Ответьте мне пожалуйста на вопрос по вашим примерам. Я вам прислал три скриншота с вашими примерами на Си и одним моим примером без использования констант. Почему результат не сходится с вашим описанием? С уважением Владимир. ******************************************************************************************************************************************** От кого: "Grigoriy A. Sitkarev" Отправитель: "Grigoriy A. Sitkarev" Тема: Re: Статья - ошибка в программе при вычислениях с плавающей точкой Дата: Thu, 30 Sep 2010 18:22:56 +0400 Кому: info@softelectro.ru Это надо задавать вопрос тем кто писал компилятор MS VC++. Я использую GNU C и для Windows не пишу, такой уж я счастливый человек. Боюсь что в MS VC++ могут быть какие-то нюансы, и они скорее всего не совпадают с ANSI рекомендациями, касательно приведения типов и арифметики с плавающей точкой. В ANSI С99 должно быть такое железное правило, что константы с плавающей точкой, определённые без суффика f, имеют тип double, т.е. двойной точности. Чтобы сделать её одинарной, нужно этот суффикс добавить, как в примере. Дальше, если в выражении хотя бы один из операндов имеет тип double, то все они должны быть "подтянуты" к этому типу, т.е. все станут double и код генерируется такой, который работает уже с double. При конверсии типов с меньшей точностью в тип с большей точностью потери не должно происходить. В ANSI C99 это называется `usual arithmetic conversion' и описано в пункте 6.3.1.8. Если у вас нет этого документа, то я могу попробовать его найти в интернете, потому что сам брал его оттуда. Могу предположить что в MS VC++ не реализованы рекомендации C99, в отличие от GCC, где они _почти_ все присутствуют так или иначе. -- Г.А. ******************************************************************************************************************************************* От кого: "Grigoriy A. Sitkarev" Отправитель: "Grigoriy A. Sitkarev" Тема: Re: Статья - ошибка в программе при вычислениях с плавающей точкой Дата: Thu, 30 Sep 2010 18:46:32 +0400 Кому: info@softelectro.ru Владимир, Я посмотрел что там с C99 у компиляторов MS, и вы знаете, этой поддержки там не было, нет, и боюсь что может и не быть в будущем. Тогда ничего удивительного в том что наши результаты не сошлись нет. ******************************************************************************************************************************************** От кого: "Grigoriy A. Sitkarev" Отправитель: "Grigoriy A. Sitkarev" Тема: Re: SoftElectro Дата: Fri, 01 Oct 2010 05:12:06 +0400 Кому: info@softelectro.ru Приветствую. Я думаю надо уточнить, дело в том что в настоящее время действует стандарт IEEE-754-2008. Стандарт 25-ти летней давности уже несколько устарел и появился можно сказать де-юре уже после того как большинство производителей вычислительной техники уже реализовали у себявычисления с плавающей точкой, правда реализации были не у всех одинаковые (в частности, округление). Стандарт приводил этот зоопарк в какую-то согласованность, чтобы программы вычислений работали одинаково на всех машинах. Когда арифметика с плавающей точкой, во многом похожая на IEEE-754, уже работала на машинах, то DOS-ом тогда ещё и не пахло, как вы понимаете, и были другие машины и другие ОС, и Unix был вместе с Си. В 50-х годах считалось что арифметика с плавающей точкой неизлечима, непредсказуема и к ней не было доверия ни у инженеров ни у учёных, которые её всё равно использовали. Случайные аномалии, которые возникали по причине ошибок округления, редко диагностировались корректно. Обычно в этом случае пытались манипулировать данными, переставлять порядок арифметических операций и т.д. В 1957-м году группа исследователей (А.М.Тюринг, Д.Х.Уилкинсон, Ф.Бауэр, В.Кахан) получила полезный метод анализа вычислительных ошибок, метод получил название "backward error analysis". До появления этого метода, для каждой из численных задач появлялись свои алгоритмы, каждый из которых хорошо работал только для одних данных, и плохо для других. В 60-х были разработаны сравнительно надёжные и гибкие алгоритмы, которые используются и по сей день. В 70-х много алгоритмов стало надёжными достаточно для того чтобы поместить их в библиотеки типа LINPACK и EISPACK, которыми могли пользоваться учёные, инженеры и т.д. которым было вовсе не обязательно понимать как они работают. В 80-х появились стандарты IEEE-754, практически реализованные к тому времени на всех (почти) машинах, поэтому библиотеки типа libm (математическая библиотека) из BSD Unix можно было без особого труда портировать на разные аппаратные архитектуры. Если кто-то изготавливал новую машину, достаточно было поставить к ней несколько небольших библиотек, оптимизированных под эту платформу, с функциями sqrt(), log(), acos() и т.д., пакетом основных подпрограмм линейной алгебры (BLAS) и затем после перекомпиляции сразу был доступ ко всему багажу относительно надёжного программного обеспечения, на котором можно было считать, а оно в свою очередь не требовало понимания КАК оно работает. В учебных программах инженеров нет курса по анализу ошибок в вычислениях с плавающей точкой. На сегодняшний день программистов, которые пишут работающие программы, и понимающих как работают вычисления с плавающей точкой стало гораздо меньше, по отношению к общему числу программистов, которые не имеют понятия об этом. Конструирование новых алгоритмов из комбинаций проверенных алгоритмов тоже не всегда спасает, т.к. точность вычислений с плавающей точкой не обладает гарантированным свойством транзитивности. Это действительно потенциальная опасность. И решать её надо, в первую очередь тем что необходимо соответствующим образом готовить программистов, во вторую очередь тем что необходимо разработать средства для отладки и диагностики таких аномалий, на сегодняшний момент это затруднительно сделать. По поводу нашего примера, мне кажется что вы меня не совсем поняли. Я говорю про то что ошибки нет, а она возникает у вас (как вам кажется) потому что число 1.0/3.0 в случае если 1.0 и 3.0 это числа с двойной точностью, будет отличаться от 1.0/3.0 в случае если 1.0 и 3.0 у нас числа с одинарной точностью. 1.0 / 3.0 != 1.0f / 3.0f Ваш пример написан так что заставляет подумать что IEEE-754 совершенно никуда негодный метод. Я же вам хотел показать что его "негодность" часто определяется нашим пониманием того как он работает. Конечно представить число 1/3 абсолютно точно не удастся в арифметике с конечной разрядностью, и мы получим лишь приближение к нему с точностью до 0.5 машинного эпсилон (ulp), а ulp для одинарной точности и двойной различны а значит и числа будут разные. Вот этот момент надо прояснить было, что я и сделал в своём примере. Можно было бы сказать и так, что приближения при одинарной и двойной точности к реальному числу 1/3 у нас *разные*. Это разные числа. Поэтому при вычитании одного из другого у нас будет небольшая ошибка округления, но это не ошибка IEEE-754. И я показал как сделать так чтобы её не было, чтобы 1/3 - 1/3 давало 0, как мы ожидаем. На вашем компиляторе примеры не работают потому что MS VC++ делаетчто-то по своему усмотрению, не так как он должен делать по С99. Я не могу объяснить что он там творит с константами и что делает в вычислениях, это надо задавать вопрос тем кто этот компилятор писал. Дальше, сравнивать два числа с плавающей точкой, если одно из них является результатом вычисления, оператором `==' вообще-то тоже нельзя, нельзя и сравнивать оператором `>',`<',`>=',`<=' и `!='. Почему? Потому что результат может быть точен с определённым допуском, из-за ошибок округления и т.д., как в нашем примере. Мы можем их сравнивать только с определённым допуском. В частном случаемы бы делали что-то вроде: if (fabs(x - y) <= eps) { /* x == y */ } if ((x - y) > eps) { /* x > y */ } if ((x - y) < eps) { /* x < y */ } Для общего случая, мы можем знать количество битов с ошибкой в мантиссе. Тогда можно всегда вычислить эпсилон: eps = (pow(2, biterr) - 1) * GSL_FLT_EPSILON; eps = (pow(2, biterr) - 1) * GSL_DBL_EPSILON; Здесь GSL_FLT_EPSILON и GSL_DBL_EPSILON соответственно машинное эпсилон для одинарной и двойной точности (можно взять из заголовка gsl_machine.h библиотеки GNU Scientific Library), biterr количество ошибочных битов в мантиссе. Д. Кнут когда-то предложил алгоритм аккуратного сравнения чисел с плавающей точкой, я сейчас его приведу ниже. int fcmp(double x, double y, double eps) { int exp; double delta, diff; if (fabs(x) > fabs(y)) frexp(x, &exp); else frexp(y, &exp); delta = ldexp(eps, exp); diff = x - y; if (diff > delta) return 1; else if (diff < -delta) return -1; return 0; } Тогда наши числа 1.0/3.0 и 1.0f/3.0f можно бы было сравнивать так: int res; double eps; eps = (pow(2.0, 29)-1) * GSL_DBL_EPSILON; res = fcmp(a/b, 1.0/3.0, eps); Или так как мы знаем что для нашего примера епсилон это ulp для арифметики с одинарной точностью можно сразу написать так: res = fcmp(a/b, 1.0/3.0, GSL_FLT_EPSILON); И eps из первого примера достаточно точно равен GSL_FLT_EPSILON. Теперь про ваш пример, вот что у меня получается. sitkarev@illo:~$ gcc -Wall -O0 -o fp fp.c -lm sitkarev@illo:~$ cat fp.c #include #include #include #include #include #include int main(int argc, char *argv[]) { float a, b, c, d, e; a = 1.0f; b = 3.0f; c = a / b; e=1.0e9f; d = (c - a/b) * e; printf("Result: %f\n", d); return 0; } sitkarev@illo:~$ ./fp Result: 0.000000 Так оно и должно быть, потому что все числа у нас одинарной точности. Поэтому приближения a/b и 1.0f/3.0f - одно и то же число. Здесь всё корректно. Если у вас получается код который генерирует отличное от нуля число - это лишний повод поругать MS, потому что код который генерирует их компилятор не соответствует принципу приведения типов, который был принят в Си испокон веков (K&R). Вы упомянули к месту, что вычисления проходят не в компиляторе а в сопроцессоре, это ясно как божий день. Но вы не забывайте что это компилятор генерирует машинный код, а здесь у него есть "простор" для творчества, что мы и наблюдаем в случае с MS VC++. Я не думаю что надо убирать одинарную точность из компиляторов. Проблема же не в том что IEEE-754 с одинарной точностью так работает, проблема в том что надо учить как ей корректно пользоваться и когда ей пользоваться нельзя. Вы же не будете предлагать убирать из Си указатели, потому что многие не умеют ими корректно пользоваться? Если мы что-то не понимаем, это повод задуматься над тем как работает наша система образования, чему, как и зачем мы учим наших специалистов. Я вам рекомендую посмотреть на библиотеку GNU MP. Это библиотека для вычислений с произвольной точностью. http://gmplib.org Нашу переписку можете публиковать только её нужно будет отредактировать наверное как-то чтобы стиль соответствовал. По поводу статьи, честно говоря, я пока не знаю что ответить, я подумаю. От себя ещё хочу сказать вот что. Мы упали очень сильно в уровне подготовки наших специалистов. Большинство промышленных автоматизаторов вряд-ли смогут выжить без Step7 и удобных оболочек. Редко встретишь специалиста который может взять контроллер с голой ОС, написать к нему драйвера и программировать на Си задачи автоматизации или же взяться за разработку своей среды исполнения. Мало кто из автоматизаторов понимает _на самом деле_ как работает тот же PID-регулятор или фильтр НЧ, потому что всю жизнь пользовались готовыми компонентами из сред разработки или коробочками-регуляторами. И это не только в промышленной автоматизации сейчас такой подход. Проблема в том что ошибки, заложенные в фундамент таких "удобных" компонентов будут тиражироваться а вероятность их обнаружения специалистами которые слабо представляют что это на самом деле бесконечно мала. И ещё скажу, что это всё не случайно, потому что кому-то очень хочется остановить технический прогресс, а для этого в первую очередь нужно лишить общество фундамента на котором базируются знания. Немного как-то сумбурно всё получилось, но оно и не удивительно, пять утра на часах. Надеюсь что смог чем-то помочь. -- Григорий *********************************************************************************************************************************************