Публикации FreePascal

Косячная арифметика

15.12.2018 Вадим Исаев

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

Билл Клинтон, бывший президент США

1. Аннотация
2. Небольшое предисловие
3. Неправильные типы данных
4. Опять типы данных, но уже разные
5. Автоматическая типизация - друг или враг?
6. Чёрная дыра или откуда берутся вечные двигатели
7. Выводы
8. Ссылки

1. Аннотация

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

2. Небольшое предисловие

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

Статья написана на основе [1], но в отличие от оригинала, примеры приведены на основе языка Pascal. Опять же в, отличие от оригинала, стержень здесь - именно примеры, теории же совсем мало и приводится она только как объяснение возникающих ошибок в, казалось бы, правильном коде. Хотелось бы особо подчеркнуть, что подобные ошибки не зависят от языка программирования и связаны исключительно со стандартом IEEE 754. описывающим формы представления чисел с плавающей точкой.

С равным успехом подобные ошибки можно наловить в программах на языках C\C++, Basic и даже черезвычайно научном Fortran. В общем там, где используются типизированные данные с плавающей точкой. Некоторые ошибки очевидны и их можно спрогнозировать заранее, а некоторые совсем даже нет.

3. Неправильные типы данных

Рассмотрим следующий пример:

program primer1;

Var
  a, b, f: Single;

Begin
  a:=123456789;
  b:=123456788;     
  f:=a-b;	
  WriteLn('Результат: ', f);
End

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

Очень интересный результат, не правда ли? Впору, на американский манер, вскричать "ВАУ!". Но в данном конкретном случае проказник и весельчак "вау" здесь совершенно не при чём. В этом коде у нас две ошибки:

1. Видимая (но не главная) - связана с формой представления исходных данных в коде программы;
2. Невидимая (главная) - связана с формой представления чисел для процессора.

Начнём с первой ошибки. Нас вводит в заблуждение, что переменным "a"и "b" присваиваются целые числа. И они действительно были целыми, но до тех пор, пока не попали в ласковые объятия переменных. Тип "integer", как мы знаем, имеет величину 4 байта и может оперировать плюс-минус двухмиллиардными значениями. Значения переменных вроде бы входят в этот диапазон, но, вот ведь незадача, оперция даёт прямо таки чудовищную ошибку. И здесь пора перейти к ошибке № 2, которая всё нам и разъяснит.

Попав в нецелую переменную, целые числа сразу начинают иметь именно нецелую форму, которая способна сыграть злую шутку, если о ней не знать. Заглянув в WIKI [2] (а заодно и непременно в Free Pascal Reference guide [3]) мы увидим, что хотя самый скромный тип "Single" базируется на тех же 4-ёх байтах, но эти байты разделены на секции, каждая из которых определяет форму числа:

• У числа должен быть знак - "+" или "-";
• У числа должна быть мантисса. На рисунке 1 это цифры, которые расположены до буквы "Е";
• У числа должна быть экспонента. На рисунке 1 это цифры и знак после буквы "Е";
• У экспоненты, в свою очередь, тоже должен быть знак "+" или "-".

Таким образом казалось бы большой объём хранения числа становится вовсе и не таким уж большим.

В данном случае нас интересует именно мантисса. А у мантиссы для типа "Single" предусмотрено всего лишь 7 разрядов числа. Хотя в спецификации FreePascal написано 7-8, но надо ориентироваться именно на 7, которые определены в IEEE 754. А теперь давайте вернёмся к исходному коду и посчитаем, сколько у нас разрядов в исходных данных, ведь никакой экспоненты там нет и в помине, следовательно всё число у нас одна сплошная мантисса. Вот тут то и становится видна самая главная ошибочка, которая привела к тому, что после присваивания исходных данных переменной, в них оказались совсем не те числа, которые мы задавали - компилятор их просто урезал, чтобы втиснуть в 7-ми разрядную мантиссу. Мало то, даже и по тем исходным данным, что получились после обрезания и то результат неверный. Как говориться - приходи кума любоваться:

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

program primer1_2;

Var
  a, b, f: Double;

Begin
  a:=123456789;
  b:=123456788;     
  f:=a-b;
  WriteLn('Результат: ', f);
End.

Вот теперь всё правильно:

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

Несмотря на столь благостное разрешение проблемы, у думающего человека может возникнуть вопрос - а что будет, если размера мантиссы не хватит даже в самом большом типе данных - "Extended"? В таком случае либо нужно будет менять алгоритм вычисления, либо пользоваться программами\библиотеками символьного вычисления, где ограничение на длину числа накладывает только оперативная память компьютера.

4. Опять типы данных, но уже разные

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

program primer2;

Var
  a : Double;
  b : Extended;
  c : Extended;

Begin
  a := 123456789.123457;
  b := 123456789.123457;
  c := a - b;
  WriteLn(c);
End.

и смотрим что получается:

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

program primer2_11;

Var
  a : Single;
  b : Single;
  c : Single;
    
Begin
  a := 123456789.123457;
  b := 123456789.123457;
  c := a - b;
  WriteLn(c);
End.

и вот что у нас получается:

Вот ведь ёшкин хвост - типы неправильные, а результат соответствует математическому прогнозу. Где же тут подвох? Не буду долго интриговать. Проблема здесь не в типах данных самих по себе, а в том, что они в предыдущем коде были разными, а в этом - одинаковыми. Кстати говоря, стандарт языка Си ISO/IEC 9899:1999 прямо и недвусмысленно говорит, что все данные в обязательном порядке должны быть приведены к одному типу, иначе сами видите (см. рис. 4) к чему это приводит. Давайте, как всегда, посмотрим, что у нас было на входе в предыдущем коде:

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

---

Информация к размышлению

Вот, кстати, по поводу этой циферки на хвосте "Extended". У типа "Double", между прочем, такого артефакта не наблюдается. А вот у "Extended" величина такого довеска на хвосте очень сильно зависит от величины экспоненты. Проверьте сами, что выдаст на экран переменная этого типа, подставляя в экспоненту величины от минимально возможного (1) до максимально возможного (4932). А потом задумайтесь, стоит ли верить этому типу и зачем он вообще в Pascal'е существует...

---

5. Автоматическая типизация - друг или враг?

Здесь ошибка будет такая же, как и в предыдущем примере, однако тут нам уже подкузьмил сам компилятор. Нужно сделать по этому поводу небольшое лирическое отступление. Дело в том, что если мы вставляем в свой код какие-то константы неудосужившись их затипизировать, то компилятор сам, по своему разумению, выбирает для них тип. И для констант с плавающей точкой этот тип будет "Double", но может быть и "Extended". Давайте посмотрим, какая медвежья услуга из-за этого получается:

program primer3;

Var
  a : Single;
  b : Single;
  c : Single;

Begin
  a := 1;
  b := 3;
  c := a / b;
  c := c - 1/3;
  WriteLn(c);
  WriteLn('Размер переменной с = ', sizeof(c));
  WriteLn('Размер 1/3 = ', sizeof(1/3));
End.

Как видите, размер переменной "с" соответствует типу, который мы ей назначили, а вот для операции с константами "1/3" компилятор выбрал тип "Double", отчего в ответе мы получили не "0", как это ожидалось. Исправить эту ситуацию можно двумя способами:

1. Может быть правильным, если мы поверим компилятору, что он всегда в подобных случаях будет присваивать тип "Double". Тогда все свои переменные мы тоже будем делать с типом "Double". Однако компилятор - не патриарх всея Руси, поэтому мы должны не верить, а знать, и для этого использовать способ
2. Правильный - т.е. собственноручно назначать, в данном случае операции над константами, тот тип, который нам нужен.

program primer3_1;
		
Var
  a : Single;
  b : Single;
  c : Single;
    
Begin
  a := 1;
  b := 3;
  c := a / b;
  c := c - Single(1/3);
  WriteLn(c);
End.

6. Чёрная дыра или откуда берутся вечные двигатели

Честно вам признаюсь, если бы мне показали написаный ниже код и попросили бы найти ошибку до того, как я прочитал статью [1], то наверняка потратил бы на поиски ошибки уйму времени. Но давайте рассмотрим эту отнюдь не лежащую на поверхности ошибку. Представьте себе, что есть некое гипотетическое фармакологическое производство. В бункер засыпают несметное количество таблеток и он выдаёт их по одной на упаковочный конвейер. Как только таблетки в бункере кончатся, нужно подать сигнал. То, что таблетки кончились определяется по весу - 0 грамм, подаётся сигнал. Работает такой код:

program primer4;
		
Var
  a: Single; // вес таблетки в кг
  c: Single; // продукции в бункере в кг
  n: LongInt;// количество циклов
Begin
  c := 300; // Вес таблеток в бункере
  a := 0.00001; // вес таблетки
  n := 10000000; // Количество упакованых таблеток
    
  WriteLn('Начальный вес бункера: ', c);

  For n := 1 To n Do
    c := c - a; //одна таблетка забирается упаковочной машиной

  WriteLn('Количество упакованых таблеток: ', n);  
  WriteLn('Измененный вес бункера: ', c); 
End.

Вот тебе, бабушка, и Юрьев день... 10 миллионов таблеток уплыли явно по какой-то коррупционной схеме. Вроде бы в коде учтены предыдущие проблемы - значения соответвствуют типу переменных,в операциях используются одинаковые типы данных. Здесь вполне уместно вспомнить стихотворение известного астрофизика Бориса Штерна:

«Ой, папа, ты сказал: вчера,
Открыта чёрная дыра.
Мне непонятно, что случилось,
И где же та дыра открылась?»

Конечно на счёт коррупции я немного погорячился, но то что в нашей программе присутствует чёрная дыра (или запрещённый когда-то французской академией наук вечный двигатель) - несомненный факт. Не буду вас интриговать, я не Агата Кристи, а вы не мисс Марпл. Проблема очень похожа на ту, что была в примере 1. Вот только представлена она немного по другому. Здесь общий вес таблеток (300 кг) и вес одной таблетки (10 мг) сильно далеки друг от друга. Для выявления проблемы их можно просуммировать:

300 + 0.00001 = 300.00001

и эта сумма явно не влезает в тип "Single". Хотя такую проверку вряд ли можно признать профессиональной, однако с практической точки зрения она вполне рабочая. И если вес одной таблетки взять не "10 мг", а "100 мг", то всё будет отлично считаться.

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

program primer4_1;
		
Var
  a: Double; // вес таблетки в кг
  c: Double; // продукции в бункере в кг
  n: Double;// количество циклов
Begin
  c := 300; // Вес таблеток в бункере
  a := 0.00001; // вес таблетки
  n := 10000000; // Количество упакованых таблеток
    
  WriteLn('Начальный вес бункера: ', c);

  For n := 1 To n Do
    c := c - a; //одна таблетка забирается упаковочной машиной

  WriteLn('Количество упакованых таблеток: ', n);  
  WriteLn('Измененный вес бункера: ', c); 
End.

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

program primer4_2;

Var
  bunker, bunker1, tablet, tablet1, compens: Single;
  n, i: longint;

Begin
  tablet := 0.00001; // вес таблетки
  tablet1 := 0.0;    // вес таблетки с учетом ошибки в предыдущих итерациях
  bunker := 300.0;   // исходный вес бункера 
  bunker1 := 0.0;    // вес бункера после очередной итерации 
  compens := 0.0;    // компенсация веса таблетки 

  n := 10000000;     // количество циклов 

  For i := 0 To n Do
  Begin
    // вес таблетки с компенсацией ошибки
    tablet1 := tablet - compens;
    // вес бункера после вычета скомпенсированной таблетки
    bunker1 := bunker - tablet1;
    // вычисление компенсации для следующей итерации
    compens := (bunker - bunker1) - tablet1;
    // новый вес бункера
    bunker := bunker - tablet1;
  End;
  WriteLn('Изменённый вес бункера: ', bunker);
End.

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

7. Выводы

Если иметь в виду наши домашние и рабочие персональные компьютеры, то можно посоветовать навсегда отказаться от четырёхбайтных типов вроде "Single". Ну, правда, одна морока с ними. Тем более, что в этом сегменте компьютеров уже днём с огнём не сыщешь 32-ух разрядных процессоров. Даже процессоры ARM-сегмента последние несколько лет упорно переходят на 64 разряда. Т.е. 8-ми байтные переменные никакого замедления в расчётах не вызовут.

Точность при использовании чисел с плавающей точкой до сих пор остаётся очень серьёзной проблемой. Чем больше мы выделяем байт под число, тем больше значащих цифр в мантиссе можем использовать, а значит тем больше точность результата вычислений. И здесь у меня большие претензии к компилятору FreePascal. Дело в том, что стандарт IEEE754 в редакции 2008 года ввёл дополнительный, 16-байтный тип данных. У компиляторов других уважающих себя языков такой тип тоже давно появился. К примеру в "GCC" это тип "long double", а в "gfortran" - "real(16)". Мы же продолжаем до сих пор пользоваться каким-то подозрительным типом "Extended" размером в 10 байт. Я, конечно, понимаю - традиции и всё такое. Однако стоит подумать и о новинках. Хотя прошло уже 10 лет и 16-ти байтный тип данных вряд ли можно назвать новинкой.

8. Ссылки

1. Яшкария Владимир. IEEE 754 - стандарт двоичной арифметики с плавающей точкой
2. Википедия. IEEE 754-2008
3. Free Pascal Reference guide