Антон Григорьев
дата публикации 13-03-2001 00:00 Неочевидные особенности вещественных чисел
Более полный вариант этой статьи вошёл в книгу "О чём не пишут в книгах по Delphi"
Взяться за эту статью меня побудили появляющиеся время от
времени вопросы на Круглом Столе, вызванные непониманием внутреннего
представления вещественных чисел. Когда-то описание внутреннего
представления таких чисел было неотъемлемой частью любой
сколь-нибудь серьёзной книги по программированию, но сейчас у
авторов появились более интересные предметы для обсуждения:
COM/DCOM, ActiveX, OLE и многое другое. На вещественные числа просто
не хватает места. И люди, начавшие программирование с Delphi и не
имеющие опыта работы в более старых средах, часто оказываются
совершенно беспомощными перед непонятным поведением программы,
содержащей дробные вычисления. Надеюсь, моя статья прольёт свет на
эти вопросы и сделает поведение дробей более предсказуемым.
Для начала – немного математики. В школе мы проходим два вида
дробей – простые и десятичные. Десятичные дроби, по сути дела,
представляют собой разложение числа по степеням десяти. Так, запись
13.6704 означает число, равное
1*101+3*100+6*10-1+7*10-2+0*10-3+4*10-4.
Но внутреннее представление всех чисел в компьютере, в том числе и
вещественных – не десятичное, а двоичное. Поэтому используются
двоичные дроби. Они во многом похожи на десятичные, но основание
степени у них двойка. Так, число
101.1101=1*22+0*21+1*20+1*2-1+1*2-2+0*2-3+1*2-4.
То есть в десятичном представлении это число равно 5.8125, в чём
нетрудно убедиться с помощью любого калькулятора.
Теперь вспомним научный формат записи десятичного числа.
Первым в этой записи идёт знак числа – плюс или минус. Дальше идёт
так называемая мантисса – число от 1 до 10. Затем идёт экспонента –
степень десяти, на которую надо умножить мантиссу, чтобы получить
нужное число. Итак, уже упоминавшееся число 13.6704 запишется в этом
формате как 1.36704*101 (или 1.36704E1 по принятым в
компьютере правилам). Если записываемое число меньше единицы,
экспонента будет отрицательной. Аналогичная запись существует и в
двоичной системе. Так, 101.1101 запишется в виде
1.011101*1010 (Везде использована двоичная форма
записи, так что 1010 означает 22). Именно
такое представление используется в компьютере. Двоичная точка в
такой записи не остаётся на одном месте, а сдвигается на величину,
указанную в экспоненте, поэтому такие числа называются числами с
плавающей точкой (floating point numbers).
В Delphi существует четыре вещественных типа: Single, Double,
Extended и Real. Их общий формат одинаков: Знак Экспонента Мантисса Знак – это всегда один бит. Он
равен нулю для положительных чисел и единице для отрицательных. Что
же касается размеров мантиссы и экспоненты, то именно в них и
заключается различие между типами.
Прежде чем перейти к конкретным цифрам, рассмотрим подробнее
тип Real, сделав для этого небольшой экскурс в историю. Real – это
стандартный тип языка Паскаль, присутствовавший там изначально.
Когда создавался Паскаль, процессоры ещё не имели встроенной
поддержки вещественных чисел, поэтому все операции над этим типом
сводились к операциям с целыми числами. Соответственно, размер полей
в типе Real был подобран так, чтобы оптимизировать эти операции.
Микропроцессор Intel 8086/88 и его улучшенные варианты –
80286 и 80386 – также не имели аппаратной поддержки вещественных
чисел. Но системы на базе этих процессоров имели возможность
подключения так называемого сопроцессора. Эта микросхема работала с
памятью через шины основного процессора и обеспечивала аппаратную
поддержку вещественных чисел. В системах средней руки гнездо
сопроцессора обычно было пустым, так как это удешевляло систему
(разумеется, вставить туда сопроцессор не было проблемой). Для
каждого центрального процессора выпускались свои сопроцессоры,
маркировавшиеся Intel 8087, 80287 и 80387 соответственно. Были даже
сопроцессоры, выпускаемые другими фирмами. Они работали быстрее, чем
Intel'овские, но появлялись на рынке позже. Тип вещественных чисел,
поддерживаемый сопроцессорами, не совпадает с Real.
Чтобы обеспечить в своих системах поддержку сопроцессорных
типов, Borland вводит в Turbo Pascal типы Single, Double и Extended.
Extended – это родной для сопроцессора тип, а типы Single и Double
получаются из него очень простым усечением. При загрузке числа типа
Single или Double во внутренний регистр сопроцессора последний
конвертирует их в Extended. Напротив, при выгрузке чисел этих типов
из регистра в память сопроцессор усекает их до нужного размера.
Внутренние же операции всегда выполняются с данными типа Extended
(впрочем, из этого правила есть исключение, на котором мы
остановимся позже, после детального рассмотрения формата различных
типов). Single и Double используются для экономии памяти. Ни один из
них также не совпадает с типом Real. В системах с сопроцессорами
новые типы обрабатываются заметно (в 2-3 раза) быстрее, чем Real
(это с учётом того, что тип Real после соответствующего
преобразования также обрабатывался сопроцессором; если же сравнивать
обработку типа Extended на машине с сопроцессором и Real на машине
без сопроцессора, то там на отдельных операциях достигалась разница
примерно в 100 раз). Чтобы программы с этими типами можно было
выполнять и в системах без сопроцессора, была возможность подключать
к ним программный эмулятор сопроцессора. Обработка этих типов
эмулятором была медленнее, чем обработка Real.
Начиная с 486-ого процессора Intel берёт курс на интеграцию
процессора и сопроцессора в одной микросхеме. Процент брака в
микросхемах слишком велик, поэтому Intel идёт на хитрость: если у
микросхемы брак только в сопроцессорной части, то на этой микросхеме
прожигаются перемычки, блокирующие сопроцессор, и микросхема
продаётся как процессор 80486SX, не имеющий встроенного сопроцессора
(в отличие от полноценной версии, которую назвали 80486DX). Бывали и
обратные ситуации, когда сопроцессор повреждений не имел, зато
процессор был неработоспособен. Такие микросхемы превращали в
«сопроцессор 80487». Но это уже из области экзотики, и, насколько
мне известно, до России этот сопроцессор не дошёл.
Процессор Pentium во всех своих вариантах имел встроенный
сопроцессор. Таким образом, с приходом этого процессора тип Real
стал как бы обузой, а на передний план вышли Single, Double и
Extended. Чтобы свести к минимуму необходимые переделки программ,
Borland ввела новую директиву компилятора: {$REALCOMPATIBILITY
ON/OFF}. По умолчанию стоит OFF, что означает отсутствие полной
совместимости. В этом случае тип Real в Delphi совпадает с типом
Double. Если же совместимость включена, тип Real совпадает со своим
прообразом из Паскаля. Существует ещё тип Real48, который всегда,
вне зависимости от настроек, совпадает со старым Real. Далее в этой
статье под словом “Real” я всегда буду подразумевать старый тип.
Отмечу, что всё это появилось только в Delphi 4, в более ранних
версиях тип Real48 отсутствовал, а тип Real был всегда старым,
шестибайтным.
Итак, теперь можно, наконец, добраться до размеров полей.
| Тип |
Размер типа, байт |
Размер мантиссы, бит |
Размер экспоненты, бит |
| Single |
4 |
23 |
8 |
| Double |
8 |
52 |
11 |
| Extended |
10 |
64 |
15 |
| Real |
6 |
40 |
7 | |
Другие параметры вещественных типов, такие как диапазон и
точность, есть в справке Delphi, в которую я рекомендую почаще
заглядывать.
| Внутренний формат
вещественных чисел |
Рассмотрим тип Single, так как он является самым коротким и,
следовательно, самым простым для понимания. Остальные типы
отличаются от него только количественно. В дальнейшем числа в
формате Single мы будем записывать как
s eeeeeeee mmmmmmmmmmmmmmmmmmmmmmm, где s означает
знаковый бит, e – бит экспоненты, m – бит мантиссы. Именно в таком
порядке эти биты хранятся в четырёхбайтном значении (здесь учтена
перестановка байтов; напоминаю, что в процессорах Intel байты в
многобайтных значениях переставляются так, что младший байт идёт
первым, а старший – последним). В мантиссе хранится двоичное число.
Чтобы получить истинное значение мантиссы, к ней надо мысленно
добавить слева единицу с точкой (то есть, например, мантисса
1010000000000000000000 означает двоичную дробь 1.101). Таким
образом, имея 23 двоичных разряда, мы записываем числа с точностью
до 24-ёх двоичных разрядов. Такая запись числа называется
нормализованной.
Экспонента по определению всегда целое число. Но способ
записи экспоненты в вещественных числах не совпадает с обычным
способом записи чисел со знаком. Ноль в этом представлении
записывается как 01111111. В обычном представлении это равно 127.
Соответственно, 10000000 (128 в обычном представлении) означает
единицу, а 01111110 (126) означает –1, и так далее (то есть из
обычного беззнакового числа надо вычесть 127, и получится число,
закодированное в экспоненте).
Из описанных выше правил есть исключения. Так, если все биты
экспоненты равны нулю (то есть там стоит число –127), то к мантиссе
перед её началом надо добавлять не “1.”, а “0.” (денормализованная
запись). Это позволяет увеличить диапазон вещественных чисел. Если
бы этого исключения не было бы, минимально возможное положительное
число типа Single было бы равно примерно 5.9*10-39. А так
появляется возможность использовать числа до 1.4*10-45.
Побочным эффектом этого является то, что числа, меньшие, чем
1.17*10-38, представляются с меньшей, чем 24 двоичных
разряда, точностью.
Если все биты в экспоненте равны единице, а в матрице – нулю,
то мы получаем комбинацию, известную как INF (от английского
Infinity – бесконечность). Эта комбинация используется тогда, когда
результат вычислений превышает максимально допустимое форматом
число. В зависимости от значения бита s бесконечность может быть
положительной или отрицательной. Если же при такой экспоненте в
мантиссе хоть один бит не равен нулю, такая комбинация называется
NAN (Not A Number – не число). Попытки использования комбинаций NAN
или INF приводят к ошибке времени выполнения.
Для задания нуля все биты мантиссы и экспоненты должны
быть равны нулю (формально это означает 0*10-127).
С учётом описанных выше правил если хотя бы один бит
экспоненты не будет равен нулю (т.е. экспонента будет больше -127),
запись будет считаться нормализованной, и нулевая мантисса будет
рассматриваться как единица. Поэтому никакие другие комбинации
значений мантиссы и экспоненты не могут дать ноль.
Тип Double устроен точно так же, разница только в количестве
разрядов и в том, какое значение экспоненты берётся за ноль. Итак,
мы имеем 11 разрядов для экспоненты. За ноль берётся значение 1023.
Несколько иначе устроен Extended. Кроме количественных
отличий добавляется ещё и одно качественное: в мантиссе явно
указывается первый разряд. То есть, мантисса 1010...
интерпретируется как 1.01, а не как 1.101, как это было в типах
Single и Float. Поэтому если 23-битная мантисса типа Single
обеспечивает 24-знаковую точность, а 52-битная мантисса Double –
53-битную, то 64-битная мантисса Extended обеспечивает 64-, а не
65-битную точность. Соответственно, при денормализованной форме
записи первый разряд мантиссы явно содержит 0. За ноль экспоненты
принимается значение 16383.
Тип Real, как уже упоминалось, стоит особняком. Во-первых, в
нём используется другой порядок следования битов, а, во-вторых, не
используется денормализованная форма. Я не стал детально разбираться
с типом Real, потому что сейчас это нужно разве что историкам, но
никак не программистам.
Напомню, чуть выше я сказал, что сопроцессор всегда выполняет
все операции в формате Extended, оговорившись при этом, что есть
исключение, к которому я вернусь чуть позже. Если вы прочитали всё,
что было изложено выше, теперь у вас достаточно знаний, чтобы
понять, что это за исключение.
У сопроцессора есть специальный двухбайтный регистр,
называемый управляющим словом. Установка отдельных битов этого
регистра диктует сопроцессору то или иное поведение. Прежде всего,
это связано с тем, какие исключения может возбуждать сопроцессор.
Другие биты этого регистра отвечают за то, как будут округляться
числа, как сопроцессор понимает бесконечность – всё это можно при
необходимости узнать из документации Intel. Нас же будут
интересовать только два бита из этого слова – восьмой и девятый.
Именно они определяют, как будут обрабатываться числа внутри
сопроцессора.
Если восьмой бит содержит единицу (так установлено по
умолчанию), то десять байт внутренних регистров сопроцессора будут
использоваться полностью, и мы получим «полноценный» Extended. Если
же этот бит равен нулю, то всё определяется значением бита 9. Если
он равен единице, то используются только 53 разряда мантиссы
(остальные всегда равны нулю). Если же этот бит равен нулю – только
24 разряда мантиссы. Это увеличивает скорость вычислений, но
уменьшает точность. Другими словами, точность работы сопроцессора
может быть понижена до типа Double или даже Single. Но это касается
только мантиссы, экспонента в любом случае будет содержать 15 бит,
так что диапазон типа Extended сохраняется в любом случае.
Для работы с управляющим словом сопроцессора в модуле System
описана переменная Default8087CW:Word и процедура
Set8087CW(CW:Word). При запуске программы в переменную
Default8087CW записывается то управляющее слово, которое
установила система при запуске программы. Функция Set8087CW
записывает новое значение в управляющее слово. Одновременно это
новое значение записывается в переменную Default8087CW.
Такое поведение этой функции не всегда удобно – иногда бывает
нужно сохранить старое значение переменной Default8087CW
(впрочем, это несложно сделать, заведя дополнительную переменную). С
другой стороны, если значение управляющего слова изменить, не
используя Set8087CW (а в дальнейшем мы увидим, что такие
изменения могут происходить помимо нашей воли), то с помощью функции Default8087CW
просто нет возможности узнать текущее значение управляющего слова.
В Delphi 6 и выше появилась функция Get8087CW, позволяющая узнать
значение именно контрольного слова, а не переменной Default8087CW.
В более ранних версиях единственный способ получить значение этого
слова – использование ассемблера, тем более что в Delphi
нет проблем с ассемблерными вставками.
Установить значение управляющего слова можно с помощью
команды FLDCW, прочитать – с помощью FNSTCW. Обе эти команды имеют
один аргумент – переменную типа Word. Чтобы, например, установить
53-значную точность, не изменив при этом другие биты управляющего
слова, надо выполнить такую последовательность команд: asm
FNSTCW MyCW
AND MyCW,0FCFFh
OR MyCW,200h
FLDCW MyCW
end; Современные сопроцессоры обрабатывают числа с такой
скоростью, что при обычных вычислениях вряд ли может возникнуть необходимость
в ускорении за счёт точности - выигрыш будет ничтожен. Эта возможность используется, в основном, в тех
случаях, когда вычисления с плавающей точкой составляют значительную часть
программы, а высокая точность не имеет принципиального значения (например, в
движках для 3D-игр). Однако забывать об этой особенности работы сопроцессора не стоит,
потому что она может преподнести один неприятный
сюрприз, о котором чуть позже.
Из школы мы все помним, что не каждое число может быть
записано конечной десятичной дробью. Бесконечные же дроби бывают
двух видов: периодичные и непериодичные. Примером непериодичной
дроби является число «пи», периодичной – число 1/3 или любая другая
простая дробь, не представимая в виде конечной десятичной дроби.
Для тех, кто забыл математику, напомню, что периодичные дроби
– это такие дроби, которые содержат бесконечно повторяющуюся
последовательность цифр. Например, 1/9=0.11111…, 1/12=0.08333333…,
1/7=0.142857142857… Для записи таких чисел используют скобки – в них
заключают повторяющуюся часть. Те же числа должны быть записаны так:
1/9=0.1(1), 1/12=0.08(3), 1/7=0.1(428571).
Впрочем, это было небольшое отступление. Нас сейчас не
интересует вопрос о периодичности или непериодичности числа, нам
достаточно знать, что не все числа можно представить в виде конечной
десятичной дроби. При работе с такими числами мы всегда используем
не точное, а приближённое значение, поэтому ответ получается тоже
приближённым. Это нужно учитывать в своих расчётах.
До сих пор мы говорили о только о десятичных бесконечных
дробях. Но двоичные дроби тоже могут быть бесконечными. Даже более
того, любое число, выражаемое конечной двоичной дробью, может быть
также выражено и десятичной конечной дробью. Но существуют числа
(например, 1/5) которые выражаются конечной десятичной дробью, но не
могут быть выражены конечной двоичной дробью. Это сильно усложняет
жизнь программистам.
| Примеры
«неправильного» поведения вещественных
типов |
Далее мы рассмотрим примеры, в которых вещественные типы
ведут себя необъяснимо для человека, не знакомого с внутренним
форматом их записи. К каждому такому примеру будет дано объяснение.
Надеюсь, что изучение этих примеров поможет понять, какие подводные
камни подстерегают программиста, использующего вещественные типы, и
как эти камни обойти.
Все примеры построены одинаково: на форму надо кинуть два
компонента – метку (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.
В этой статье я постарался объяснить внутреннее устройство
вещественных чисел с точки зрения процессоров Intel и упомянуть
некоторые проблемы, которые с ними связаны. На самом деле все
проблемы сводятся к двум: во-первых, не всякое вещественное число
может быть представлено точно, и, во-вторых, не всякое вещественное
число, представимое в виде конечной десятичной дроби, представимо в
виде конечной двоичной дроби. Вторая проблема, наверное, приносит
больше неприятностей начинающим пользователям, так как она менее
очевидна. Рецепты преодоления этих проблем я сознательно не излагаю,
так как оптимальный вариант очень сильно зависит от конкретной
задачи. Человеку же, понявшему причины появления проблем, не
составит труда в каждом конкретном случае подобрать наиболее
приемлемое решение. В этом, собственно, и заключается разница между
программистом и ламером: первый разбирается в задаче и находит для
неё решение, второй умеет только кидать на форму готовые компоненты
и передирать куски чужого кода. А эту статью я писал для начинающих
программистов, а не для начинающих ламеров, отсюда и такой стиль.
Огромное спасибо Елене Филипповой за помощь в поиске
информации.
Григорьев Антон
Специально для Королевства Delphi
[Вещественные числа]
 Обсуждение материала [ 01-07-2019 03:46 ] 77 сообщений |
