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

10^00 и 5.678 * 10^-04:

1.234 * 10^00 + 5.678 * 10^-05 = 1.234 + 0.00005678 = 1.23405678

Это число непредставимо в нашей арифметике и будет заменено ближайшим представимым:

1.234 * 10^00 + 5.678 * 10^-04 = 1.234 * 10^00

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

Машинное эпсилон

Рассмотрим более подробно операцию сложения двух положительных чисел. Пусть х > у > 0. Тогда

х + у = х * (1 + у/х)

Если число у/х окажется настолько маленьким, что 1 + у/х = 1, то мы получим потерю значимости: х + у = х.

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

В нашей модельной арифметике машинное эпсилон равно 0.001 = 1.000 * 10^-03.

Вычисление машинного эпсилон

Величину машинного эпсилон можно вычислить как результат выполнения следующей функции:

typedef float real;

real calceps() {

static volatile real eps, eps2, eps21;

eps =1.0;

eps2 = eps * 0.5;

eps21 = eps2 + 1.0;

while( eps21 > 1.0 ) {

eps = eps2;

eps2 = eps * 0.5;

eps21 = eps2 + 1.0;

}

return eps;

}

Заметим, что в этой функции существенно использование переменных eps, eps2 и eps21. Спецификаторы static и volatile использованы для того, чтобы воспрепятствовать возможной оптимизации кода компилятором. Дело в том, что промежуточные вычисления проводятся с использованием 10-байтовых вещественных чисел. По этой же причине условие fabs(x) > calceps() нельзя заменить на следующее условие: 1.0 + fabs(x – y) > 1.0.

Параметры машинной арифметики

Основные параметры машинной арифметики - максимальное представимое число, минимальное положительное число и машинное эпсилон. Максимальное представимое число и минимальное положительное число определяют диапазон чисел, представимых в машинной арифметике, называемый также динамическим диапазоном арифметики. Динамический диапазон определяется разрядностью показателя. Машинное эпсилон определяет точность представления чисел и зависит от разрядности мантиссы.

Несмотря на то, что для всех основных параметров машинной арифметики существуют переносимые способы их вычисления, обычно предпочитают пользоваться готовыми значениями соответствующих констант. Заголовочный файл float.h содержит определения макросов, дающих основные параметры машинной арифметики.

	float	double
максимальное представимое число	FLT_MAX	DBL_MAX
минимальное положительное число	FLT_MIN	DBL_MIN
машинное эпсилон	FLT_EPSILON	DBL_EPSILON

Идиома 1: проверка на малость

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

fabs (е) < eps

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

fabs(e) < eps * fabs(x).

В этой реализации абсолютный порог малости числа зависит от значения величины, с которой она сравнивается - при больших по модулю значениях х даже относительно большие е будут признаны малыми. Действительно, число 1 мало по сравнению с 1е6 в той же степени, что и число 1е-6 по сравнению с 1. Оказывается, что и эта реализация может приводить к проблемам, но уже в случае малых по абсолютной величине чисел х. Если х - почти нуль, то и правая часть получается почти нулевой. В этом случае правая часть может оказаться настолько малой, что никакое е не сможет удовлетворить неравенству. Эту проблему решает такая реализация:

fabs(e) < eps * (1.0 + fabs(x))

Защититься от возможности задания излишне малых значений eps можно следующим образом:

fabs(e) < (eps + 2.0 * FLT_EPSILON) * (1.0 + fabs(x))

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