Особенности вещественных чисел в Delphi - Примеры «неправильного» поведения вещественных типов

Written on . Posted in Delphi

ОГЛАВЛЕНИЕ

Страница 6 из 7

Примеры «неправильного» поведения вещественных типов

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

Все примеры построены одинаково: на форму надо кинуть два компонента – метку (TLabel) и кнопку (TButton). Так как это только примеры, я не стал придумывать имена для этих компонентов, пусть называются Button1 и Label1. Обработчик Button1Click содержит некоторый код, результаты работы которого выводятся на форму через Label1. Таким образом, нужно запустить программу, нажать на кнопку и посмотреть, что будет написано в метке. Я буду приводить только код обработчика Button1Click, так как всё остальное тривиально. Напомню, что в Паскале допускается не ставить точку с запятой перед end'ом (за исключением end'а в описании класса) и перед until'ом. Я предпочитаю пользоваться этой возможностью, так что не надо тыкать в меня пальцем, что я забываю ставить точки с запятой.

Пример первый – «неправильное значение»

Итак, напишем такой код: 
var R:Single;
 begin
  R:=0.1;
  Label1.Caption:=FloatToStr(R)
 end;
Что мы увидим, когда нажмём кнопку? Разумеется, не «0.1», иначе не было бы смысла писать этот пример. Мы увидим «0.100000001490116». То есть расхождение в девятой значащей цифре. Ну, из справки по Delphi мы знаем, что точность типа Single – 7-8 десятичных разрядов, так что нас, по крайней мере, никто не обманывает. В чём же причина? Просто число 0.1 не представимо в виде конечной двоичной дроби, оно равно 0.0(0011). И эта бесконечная двоичная дробь обрубается на 24-ёх знаках; мы получаем не 0.1, а некоторое приближённое число (какое именно – см. выше). А если мы присвоим переменной R не 0.1, а 0.5? Тогда мы получим на экране 0.5, потому что 0.5 представляется в виде конечной двоичной дроби. Немного поэкспериментировав с различными числами, мы заметим, что точно представляются те числа, которые выражаются в виде m/2n, где m, n – некоторые целые числа (разумеется, n не должно превышать 24, а то нам не хватит точности типа Single). В качестве упражнения предлагаю доказать, что любое целое число, для записи которого хватает 24-ёх двоичных разрядов, может быть точно передано типом Single.

Пример второй – сравнение

Теперь изменим код так: 
var R:Single;
 begin
  R:=0.1;
  if R=0.1 then
   Label1.Caption:='Равно'
  else
   Label1.Caption:='Не равно'
 end;

При нажатии кнопки мы увидим надпись «Не равно». На первый взгляд это кажется абсурдом. Действительно, мы уже знаем, что переменная R получает значение 0.100000001490116 вместо 0.1. Но ведь «0.1» в правой части равенства тоже должно преобразоваться по тем же законам, ведь в компьютере всё предопределено. Тут самое время вспомнить, что процессоры Intel работают только с 10-байтным типом Extended, поэтому и левая, и правая часть равенства сначала преобразуется в этот тип, и лишь потом производится сравнение. То корявое число, которое оказалось в переменной R вместо 0.1, хоть и выглядит страшно, но зато представляется в виде конечной двоичной дроби. Информация же о том, что это на самом деле должно означать «0.1», нигде не сохранилось. При преобразовании этого числа в Extended младшие, избыточные по сравнению с типом Single разряды мантиссы просто заполняются нулями, и мы снова получим то же самое число, только записанное в формате Extended. А «0.1» из правой части равенства преобразуется в Extended без промежуточного превращения в Single. А 0.1 – бесконечная в двоичном представлении дробь. Поэтому некоторые из младших разрядов мантиссы будут содержать единицы. Другими словами, мы получим хоть и не точное представление числа 0.1, но всё же более близкое к истине, чем 0.100000001490116. Из-за таких хитрых преобразований оказывается, что мы сравниваем два близких, но всё же не равных числа. Отсюда – закономерный результат в виде надписи «Не равно».

Тут уместна аналогия с десятичными дробями. Допустим, в одном случае мы делим 1 на три с точностью до трёх знаков, и получаем 0.333. Потом мы делим 1 на три с точностью то четырёх знаков, и получаем 0.3333. Теперь мы хотим сравнить эти два числа. Для этого приводим их к точности в четыре разряда. Получается, что мы сравниваем 0.3330 и 0.3333. Очевидно, что это разные числа.

Если попробовать заменить число 0.1 на 0.5, то мы получим «Равно». Думаю, вы уже знаете почему, но для полноты текста объясню. 0.5 – это конечная двоичная дробь. При прямом приведении её к типу Extended в младших разрядах оказываются нули. Точно такие же нули оказываются в этих разрядах при превращении числа 0.5 типа Single в тип Extended. Поэтому в результате мы сравниваем два числа. Это похоже, как если бы мы делили 1 на 4 с точностью до трёх и до четырёх значащих цифр. В первом случае получили бы 0.250, во втором – 0.2500. Приведя их оба к точности в четыре знака, получим сравнение 0.2500 и 0.2500. Очевидно, что эти цифры равны.

Пример третий – сравнение разных типов

Немного усложним наш пример: 
var R1:Single;
    R2:Double;
 begin
  R1:=0.1;
  R2:=0.1;
  if R1=R2 then
   Label1.Caption:='Равно'
  else
   Label1.Caption:='Не равно'
 end;

Наученные горьким опытом, вы, наверное, ожидаете увидеть надпись «Не равно». Что ж, жизнь вас не разочарует, именно это вы и увидите. Тип Double точнее, чем Single (хотя его точности тоже не хватает для представления бесконечной дроби). В R2 мы получим не 0.100000001490116, а другое число, с точностью 15-16 десятичных знаков. Я не могу назвать точно это число, потому что FloatToStr воспринимает его как 0.1, так что, заменив в первом примере Single на Double, вы увидите 0.1 (только не надо обольщаться, всё равно это не 0.1, просто функция FloatToStr имеет такую особенность работы). Числа в обеих переменных приводятся к типу Extended, но при этом они не меняются и, как были не равны, так и остаются неравными. Это напоминает ситуацию, когда мы сравниваем 0.333 и 0.3333, приводя их к точности в пять знаков: числа 0.33300 и 0.33330 не равны.

Мне уже неловко надоедать вам такими очевидными замечаниями, но всё-таки: если в этом примере заменить 0.1 на 0.5, мы увидим «Равно».

Пример четвёртый – вычитание в цикле

Рассмотрим ещё один пример, иллюстрирующий ситуацию, которая часто озадачивает начинающего программиста 
var R:Single;
    I:Integer;
 begin
  R:=1;
  for I:=1 to 10 do
   R:=R-0.1;
  Label1.Caption:=FloatToStr(R)
 end;
Конечно, если бы в результате выполнения этого примера вы увидели бы ноль, я бы не стал тратить на него время. Но на экране появится -7.3015691270939E-8. Думаю, такой оборот дела уже никого не удивляет. Мы уже знаем про то, что число 0.1 не может быть передано точно ни в одном из вещественных типов, и про преобразования Single в Extended и обратно. При этом постоянно происходят округления, и эти округления приводят к тому, что мы получаем в результате не ноль, а «почти ноль».

Пример пятый – сюрпириз от Microsoft

Изменим в предыдущем примере тип переменной R с Single на Double. Значение, выводимое программой, станет 1.44327637948555E-16. Вполне логичный и предсказуемый результат, так как тип Double точнее, чем Single и, следовательно, все вычисления более точны, мы просто обязаны получить более точный результат. Хотя, разумеется, абсолютная точность (то есть ноль), для нас остаётся недостижимым идеалом.

А теперь – вопрос на засыпку. Изменится ли результат, если мы заменим Double на более точный Extended? Ответ не такой однозначный, каким его хотелось бы видеть. В принципе, после такой замены вы должны получить -6.7762635780344E-20. Но в некоторых случаях от замены Double на Extended результат не изменится, и вы снова получите 1.44327637948555E-16. Это зависит от операционной системы.

Всё дело в использовании «неполноценного» Extended. При запуске программы любая система устанавливает такое управляющее слово сопроцессора, чтобы Extended был полноценным. Но затем программа вызывает много разных функций Windows API. Какая-то (или какие-то) из этих многочисленных функций в некорректно работают с управляющим словом, меняя его значение и не восстанавливая при выходе. Такая проблема встречается, в основном, в Windows 95 и старых версиях Windows 98. Также имеются сведения о том, что управляющее слово может портиться и в Windows NT, причём эффект наблюдался не сразу после установки системы, а лишь через некоторое время, после доустановки других программ. Проблема именно в некорректности поведения системных функций; значение управляющего слова, устанавливаемое системой при запуске программы, всегда одинаково. Эта проблема известна: например, в исходных кодах VCL можно найти сохранение управляющего слова сопроцессра перед вызовом некоторых API-функций с последующим его восстановлением. Комментарии сообщают, что функция может изменить значение управляющего слова, поэтому необходимо его сохранение и восстановление.

Таким образом, приходим к неутешительному выводу: к тем проблемам с вещественными числами, которые обусловлены особенностями их аппаратной реализации, добавляются ещё и баги Windows. Правда, радует то, что в последнее время эти баги встречаются крайне редко - видимо, новые версии системы ведут себя более ответственно. Тем не менее, полностью исключать такую возможность нельзя, особенно если ваша программа будет использоваться на устаревшей технике с устаревшими системами (например, в образовательных учреждениях, финансирование которых оставляет желать лучшего). Чтобы наш пример всегда выдавал правильное значение -6.7762635780344E-20, достаточно поставить в начале нашей процедуры Set8087CW(Get8087CW or $0100), и программа в любой системе будет использовать сопроцессор в режиме максимальной точности. (Если вы используете старые версии Delphi, эту строку можно заменить на Set8087CW(Default8087CW), если, конечно, значения по умолчанию прочих флагов управляющего слова вас устраивают.)

Раз уж мы заговорили об управляющем слове, давайте немного поэкспериментируем с ним. Изменим первую строчку на Set8087CW(Get8087CW and $FCFF or $0200). Тем самым мы переведём сопроцессор в режим 53-ёхразрядной точности представления мантиссы. Теперь в любой системе мы увидим 1.44327637948555E-16, несмотря на использование Extended. Если же мы изменим первую строчку на Set8087CW(Get8087CW and $FCFF), то будем работать в режиме 24-ёхразрядной точности. Соответственно, в любой системе будет результат -7.3015691270939E-8.

Заметим, что при загрузке в 10-байтный регистр сопроцессора числа типа Extended в режиме пониженной точности «лишние» биты не обнуляются. Только результаты математических операций представляются с пониженной точностью. Кроме того, при сравнении двух чисел также учитываются все биты, независимо от точности. Поэтому код 

var R:Double; // или Single

 begin
  R:=0.1;
  if R=0.1 then
   Label1.Caption:='Равно'
  else
   Label1.Caption:='Не равно'
 end;
при выборе любой точности даст «Не равно».

Пример шестой – машинное эпсилон

Когда мы имеем дело с вычислениями с ограниченной точностью, возникает такой парадокс. Пусть, например, мы считаем с точностью до трёх значащих цифр. Прибавим к числу 1.00 число 1.00*10-4. Если бы всё было честно, мы получили бы 1.0001. Но у нас ограничена точность, поэтому мы вынуждены округлять до трёх значащих цифр. В результате получается 1.00. Другими словами, мы прибавляем к единице некоторое число, большее нуля, а в результате из-за ограниченной точности получаем снова единицу. Наименьшее положительное число, которое при добавлении его к единице даёт результат, не равный единице, называется машинным эпсилон.

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

Прежде чем искать машинное эпсилон программно, попытаемся найти его из теоретических соображений. Итак, мантисса типа Extended содержит 64 разряда. Чтобы закодировать единицу, старший бит мантиссы должен быть равен 1 (денормализованная запись), остальные биты - нулю. Очевидно, что при такой записи наименьшее из чисел, для которых вполняется условие x>1, получается, когда самый младший бит мантиссы тоже будет равен единице, т.е. x=1.00...001 (в двоичном представлении; между точкой и младшей единицей 62 нуля). Таким образом, машинное эпсилон равно x-1, т.е. 0.00...001. В более привычной десятичной форме записи это будет 2-63, т.е. примерно 1.084*10-19.

Теперь напишем программу для отыскания машинного эпсилон. 

var R:Extended;
 begin
  R:=1;
  while 1+R/2>1 do
   R:=R/2;
  Label1.Caption:=FloatToStr(R)
 end;

В результате на экране появится число 1.0842021724855E-19 в полном соответствии с теоретическими выкладками (если в вашей системе присутствует описанный выше баг с переводом процессора в режим пониженной точности, вместо этого числа вы получите 2.22044604925031E-16, т.е. 2-52. Чтобы этого не происходило, исправьте значение управляющего слова).

А теперь заменим тип Extended на Double. Результат не изменится. На Single – опять не изменится. Но такое поведение лишь на первый взгляд может показаться странным. Давайте подробнее рассмотрим выражение 1+R/2>1. Итак, все вычисления (в том числе и сравнение) сопроцессор выполняет с данными типа Extended. Последовательность действий такова: число R загружается в регистр сопроцессора, преобразуясь при этом к типу Extended. Дальше оно делится на 2, а затем к результату прибавляется 1, и всё это в Extended, никакого обратного преобразования в Single или Double не происходит. Затем это число сравнивается с единицей. Очевидно, что результат сравнения не должен зависеть от исходного типа R.

форум delphi