Особенности вещественных вычислений

В отличие от целочисленных выражений, которые всегда вычисляются точно, вещественные выражения дают приближенный результат и вещественные переменные содержат приближенные значения. Это обстоятельство определяется способом хранения вещественных чисел в памяти компьютера. Любое вещественное значение представляется как sm×2^p, где s - знак числа (плюс или минус), m - двоичная мантисса, которая принадлежит интервалу [1,2) для всех ненулевых чисел, p - двоичный порядок числа. И мантисса, и порядок хранятся в двоичной кодировке. Величины типа Real занимают в памяти 6 байтов, или 48 битов, старший бит определяет знак числа - если он единичный, то число отрицательное. В следующих 39 битах хранится дробная часть мантиссы, то есть число m-1. И в восьми младших битах хранится число p+129, то есть увеличенный на 129 двоичный порядок числа; такая кодировка выбрана для того, чтобы не хранить отрицательные порядки. Представления чисел типа Single, Double и Extended отличаются по расположению полей от представления типа Real - у них порядок хранится в старших битах, а мантисса - в младших. Способы кодировки этих типов таковы:

тип Single: s (1 бит), p+127 (8 битов),m-1 (23 бита),

тип Double: s (1 бит), p+1023 (11 битов), m-1 (52 бита),

тип Extended: s (1 бит), p+16383 (15 битов), m(63 бита).

Способ кодировки и определяет как допустимые для данного типа значения, так и количество верных десятичных цифр. Определим, например, наибольшее допустимое значение типа Real. Максимальное значение величины p+129, занимающей 8 битов, есть 255, следовательно, наибольший возможный порядок числа равен 126. Наибольшее возможное значение мантиссы чуть меньше двойки, а именно равно 2-2^-39, тогда наибольшее число типа Real равно 2¹²⁷ - 2⁸⁷, что примерно равно 1.7014×10³⁸. Это число называется машинной бесконечностью; конечно, машинная бесконечность своя у каждого вещественного типа. Теперь найдем машинный ноль - наименьшее ненулевое число типа Real. Наименьший возможный порядок числа равен, очевидно, -129 в случае, когда поле порядка содержит число 0, но все вещественные числа, независимо от их мантиссы, у которых поле порядка нулевое, считаются точными нулями, следовательно, машинный ноль имеет порядок -128. Наименьшее значение мантиссы равно единице и, таким образом, машинный ноль типа Real равен 2^-128, или примерно 2.9387×10^-39. Зная кодировки других вещественных типов, вы легко сможете найти для них машинные бесконечности и машинные нули. По длине поля мантиссы можно оценить количество верных десятичных цифр в значениях данного типа. В типе Real поле мантиссы содержит 39 битов, следовательно, погрешность составляет 2^-39, что примерно равно 10^-12, таким образом, в числе типа Real не более 12 верных десятичных цифр.

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

 

{$N+,E+}

{$IFOPT N+}

Function Abs(x:Extended):Extended;

{функция Abs для режима с сопроцессором}

Var b : Array[1..10] Of Byte Absolute x;

{совмещаем по памяти 10-байтовый массив с аргументом}

Begin

b[10]:=b[10] And $7F; {зануляем старший бит старшего байта}

Abs:=x;

End;

{$ELSE}

Function Abs(x:Real):Real; {функция Abs для режима без сопроцессора}

Var b : Array[1..6] Of Byte Absolute x;

Begin

b[6]:=b[6] And $7F;

Abs:=x;

End;

{$ENDIF}

Var

r:Real;

d:Double;

s:Single;

Begin

r:=-333;

WriteLn(Abs(r));

d:=-777;

WriteLn(Abs(d));

s:=-555;

WriteLn(Abs(s));

End.

 

Конструкции {$IFOPT N+}, {$ELSE} и {$ENDIF} в этой программе - это директивы условной компиляции. Если опция N была включена, то компилятор будет компилировать операторы от {$IFOPT N+} до {$ELSE} и проигнорирует операторы от {$ELSE} до {$ENDIF}, если же опция N выключена, то компилироваться будет фрагмент программы от {$ELSE} до {$ENDIF}. Таким образом, у нас есть два разных экземпляра функции, но при каждом запуске программы работать будет только один из них, причем именно тот, который нужен. Теперь проверим, верны ли сведения, изложенные в данном разделе, пока мы лишь убедились, что старший бит действительно отвечает за знак числа. Возьмем какое-нибудь вещественное число, например типа Single, и попробуем восстановить его значение по его внутреннему двоичному представлению.

 

Var

x :Single;

b :Array[1..4] Of Byte Absolute x;

p :Integer;

m,Pow :Real;

Mask,L,M_ :LongInt;

k :Byte;

s :Boolean;

Begin

x:=-0.375;

{определим знак числа - он кодируется старшим битом b[4]}

s:=b[4] And $80>0;

If s Then WriteLn('Число отрицательное')

Else WriteLn('Число неотрицательное');

{следующие 8 битов определяют двоичный порядок, увеличенный на 127 - это байт b[4] без старшего бита и старший бит байта b[3]}

p:=Byte(b[4] ShL 1)+b[3] ShR 7;

Dec(p,127);

WriteLn('Двоичный порядок числа=',p);{p=-2}

{оставшиеся 23 бита кодируют дробную часть мантиссы - это байты b[1],b[2] и байт b[3] без старшего бита, соберем их в переменной L}

L:=b[1]+LongInt(b[2])ShL 8+LongInt(b[3] And $7F)ShL 16;

{вычислим мантиссу}

Mask:=1;

M_:=0;

For k:=0 To 22 Do Begin

If L And 1>0 Then Inc(M_,Mask);

Mask:=Mask ShL 1;

L:=L ShR 1;

End;

m:=M_/(LongInt(1)ShL 23);

WriteLn('Дробная часть мантиссы равна ',m);{m=0.5}

m:=m+1;

WriteLn('Мантисса равна ',m);{m=1.5}

{вычислим 2 в степени p}

Pow:=1;

For k:=1 To Abs(p) Do Pow:=Pow*2;

If p<0 Then Pow:=1/Pow;

{теперь вычислим значение переменной x}

x:=m*Pow;

If s Then x:=-x;

WriteLn('Значение переменной равно ',x);{x=0.375}

End.

 

Следствием приближенности вещественных чисел является, например, тот факт, что уравнение x+a=a, которое в обычной математике имеет единственный корень x=0, в компьютерной математике может иметь много нетривиальных решений; если уравнение решается в типе Real, то любой x, удовлетворяющий неравенству x<10^-12×a, будет решением. Аналогично тождество (x/a)×a=x, как правило, не выполняется, даже если a не равно нулю. Вообще говоря, операции сравнения “равно” и “не равно” не имеют смысла для вещественных операндов, их следует применять в программах крайне осторожно. Пусть, например, мы хотим вывести значения некоторой функции вещественного аргумента на отрезке [0,1] с шагом 0.1, запишем этот алгоритм так :

x:=0; While x<>1 Do Begin WriteLn(f(x)); x:=x+0.1; End;

и программа зациклится - значение x никогда не станет равным единице. Лучше записать этот фрагмент следующим образом :

For i:=0 To 10 Do Begin x:=i*0.1; WriteLn(f(x)); End;

Однако приближенность представления вещественных чисел имеет не только свои минусы, но и свои плюсы. Решим реальную задачу, в которой используются приближенные вычисления: вычислить сумму ряда åxⁱ/i! , i=0,1,...,∞. Несмотря на то, что необходимо просуммировать бесконечное число слагаемых, эта задача легко решается за конечное время, так как общий член ряда быстро убывает и начиная с некоторого n прибавление очередного слагаемого уже не будет изменять сумму.

 

Var x,S1,S2,a : Real;

Const i : Word = 0;

Begin

Write('Введите x '); Read(x);

a:=1; {общий член ряда}

S2:=0; {предыдущая частичная сумма}

Repeat

S1:=S2;

S2:=S1+a;

Inc(i);

a:=a*x/i;

Until S1=S2;

WriteLn('Сумма ряда = ',S1);

End.

Решим еще одну полезную задачу: найти корень уравнения f(x)=0 методом бисекции. Метод бисекции заключается в следующем: пусть уравнение f(x)=0 имеет единственный корень на отрезке [a,b], это значит, что график функции один раз пересекает ось X на этом отрезке. Определим знак функции в точке a и в точкеx=(a+b)/2. Если эти знаки одинаковы, то очевидно, корень лежит на отрезке [x,b], в противном случае - на отрезке [a,x]. Таким образом, за один шаг метода мы ровно вдвое уменьшили наш отрезок. Будем повторять эти операции до тех пор, пока отрезок не станет очень маленьким, и в качестве корня возьмем середину этого маленького отрезка. Попробуем реализовать этот метод:

 

Const

a : Real=0;

b : Real=10;

Function f(x:Real):Real; Begin f:=exp(x)-2; End;

Const epsilon=1E-10; {"очень маленький" отрезок}

Var x,Fa,Fx : Real;

Begin

Fa:=f(a); {нет необходимости вычислять f(a) в цикле, т.к. ЗНАК функции на левом конце отрезка не изменится!}

While b-a>epsilon Do Begin

x:=(a+b)/2;

Fx:=f(x);

If Fa*Fx<0 Then b:=x Else a:=x;

End;

WriteLn(x);

End.

 

Программа вычислила 0.693147...; легко проверить, что это правильный результат (корень уравнения равен ln2). Казалось бы, мы написали хорошую программу, но давайте теперь вычислим корень уравнения ln(x)-50=0, a=1, b=1e30 : программа зациклится! Выведем внутри цикла значенияa и b : эти числа почти одинаковы (отличие в 12-й цифре) и не меняются, но поскольку их порядок 10²¹, b-a существенно превосходит наш “очень маленький” отрезок. Есть два способа, которыми мы можем исправить программу. Первый способ заключается в правильном подборе “очень маленького” отрезка, но надо понимать, что вам придется подбирать этот отрезок для каждого нового уравнения, то есть фактически для каждого уравнения писать свою программу. Очевидно, что это бесперспективный путь. Выведем в нашей зацикленной программе не только a и b, но иx, может быть, это поможет нам придумать второй способ: значение x, оказывается, в точности равно b.

Мы могли бы прийти к выводу, что рано или поздно x станет равным илиa, или b, рассуждая чисто теоретически. Действительно, на каждом шаге цикла мы уменьшаем отрезок в два раза, если бы мы работали на вещественной оси, то величины a и bстремились бы друг к другу бесконечно, но, поскольку множество вещественных чисел в компьютере дискретно (из-за конечного количества цифр мантиссы), настанет момент, когда между a и b больше не будет ни одного числа. После этого выражение (a+b)/2 будет давать либо a, либо b. Воспользуемся этим обстоятельством и напишем хорошую программу :

 

Const

a : Real = 1;

b: Real = 1E30;

Function f(x:Real):Real; Begin f:=Ln(x)-50; End;

Var x,Fa,Fx : Real;

Begin

Fa:=f(a); x:=(a+b)/2;

While (x<>a)And(x<>b) Do Begin

Fx:=f(x); If Fa*Fx<0 Then b:=x Else a:=x;

x:=(a+b)/2;

End;

WriteLn(x);

End.

Программа дала верный результат 5.184705...E21.

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