Машинная арифметика и идиомы численного программирования, страница 3

float x; // вычисляемая величина float e; // ошибка ответа

х = ...; // вычисляем первое приближение

е = ...; // вычисляем ошибку

// пока ошибка велика по отношению к ответу

while( fabs(e) > (eps + 2.0 * FLT_EPSILON) * (1.0 + fabs(x)) ){

x = ...; // вычисляем новое приближение

е = ...; // вычисляем новое значение ошибки }

Идиома 1': сравнение на равенство

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

fabs(x - у) < (eps + 2.0 * FLT_EPSILON) * (1.0 + fmax(fabs(x), fabs(у))),

где fmax() - функция, возвращающая значение максимального из своих аргументов.

Пример: метод бисекции

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

float fa = f(a);

while( b - а > eps ) {

float m = (a + b)/2.0;

float fm = f(m);

if( fa * fm <= 0.0 ) {

b = m;

} else {

a = m;

fa = fm;

}

}

Несмотря на простоту, эта программа содержит несколько принципиальных ошибок. В первую очередь отметим, что при вычислении произведения fa * fm може произойти как переполнение, так и потеря значимости. В случае потери значимости (если числа fa и fm одновременно весьма малы) даже для чисел разного знака мы получим ноль. Следует проверять именно условие на знаки чисел, страхуясь таким образом от неприятностей.

Вторая проблема связана с проверкой на приближенное равенство. Для больших чисел а и b проверка b - а > eps может привести к бесконечному циклу. Эту проверку следует заменить на следующее условие:

b - a > (eps + 2.0 * FLT_EPSILON) * (1.0 + fmax(fabs(a), fabs(b))),

Третья проблема связана с тем, что число m, определенное формулой m = (а + b)/2.0, может оказаться лежащим за пределами отрезка [а, b]. В качестве примера рассмотрим модельную

4-разрядную арифметику и пусть а = 0.9882 и b = 0.9884. Тогда а + b = 1.9766, но это число непредставимо в нашей арифметике и реально получится а + b = 1.977. Следовательно, m = (а + b)/2.0 = 0.9885 > b. Избежать этой неприятности можно, изменив формулу для вычисления середины отрезка [a, b]:

m = а + (b - а)/2.0

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

Исправленный вариант программы имеет следующий вид:

float abseps = (eps + 2.0 * FLT_EPSILON) * (1.0 + fmax(fabs(a), fabs(b)));

float fa = f(a);

while( b - a > abseps ) {

float m = a + (b - a)/2.0;

float fm = func(m);

if( fsgn(fa) * fsgn(fm) <= 0.0 ) {

b = m;

} else {

a = m;

fa = fm;

}

}

Здесь функция fsgn() возвращает знак аргумента т.е. +1.0 для положительных аргументов, -1.0 - для отрицательных и 0.0 для нулевых.

Идиома 2: переполнение и потеря значимости

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

float norm( float *х, int len ){

int j;

float cur, sum2;

sum2 = 0.0;

for( j = 0; j < len; j++ ){

cur = x[j];

sum2 += cur * cur;

}

return sqrt(sum2);

}

является некорректной. Простое масштабирование решает эту проблему, но код становится существенно более сложным:

float norm( float *х, int len ){

int j;

float cur, max, sum2;

max =0.0;

for( j  = 0; j < len; j++ ){

cur = fabs(x[j]);

if( cur > max )

max = cur;

}

if( max == 0.0 ) return 0.0;

sum2 = 0.0;

for( j =0; j < len; j++ ){

cur = x[j] / max;

sum2  += cur * cur;

}

return      max * sqrt(sum2);

}