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

Более существенным здесь оказывается даже не возможность переполнения, а существенная потеря точности ответа, поскольку при возведении в квадрат и последующем извлечении корня мы фактически теряем половину значащих цифр результата. В этом случае часто говорят о потере значимости. Масштабирование решает и эту проблему.

Аналогичная проблема возникает при вычислении корня из суммы квадратов двух чисел (норма вектора на плоскости, норма комплексного числа), но эта проблема решается использованием стандартной функции hypot(x, у), которая проводит вычисления по формуле sqrt(x*x + у*у) с контролем возможности переполнения и потери значащих цифр. При этом функция hypot(x,  у) выполняется быстрее, чем вычисления по формуле с использованием масштабирования.

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

Пример: деление комплексных чисел

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

typedef struct { float re, im; } cmplx;

cmplx cmplx_div( cmplx c, cmplx d ){

float r, p; cmplx res;

if( fabs(d.re) >= fabs(d.im) ){

r = d.im / d.re;

p = d.re + r * d.im;

res.re  = (c.re + r * c.im) / p;

res.im = (c.im - r * c.re) / p;

}else{

r = d.re / d.im;

p = d.im + r * d.re;

res.re = (c.re * r + c.im) / p;

res.im = (c.im * r – c.re) / p;

}

return res;

}

Потеря точности при вычитании

При нахождении разности близких по величине чисел может происходить катастрофическая потеря точности. Пусть в модельной 4-разрядной арифметике два непредставимых числа х и у оказались представленными числами х_* = l.001 и y_* = 1.002. В обоих случаях относительная погрешность представления принимает значения порядка машинного эпсилон (eps = 0.001). Вычислим относительную погрешность разности z* = у_* - х_* = 0.001:

| (y_* - x_*) - (y - x) | / | y – x | = 2 * eps / 0.001 = 2000 * eps.

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

Примеры безосновательной параной

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

В известной книге Л. Аммерала "Принципы программирования в машинной графике" можно найти большой список бессмысленно параноидальных советов. В этой книге приведена следующая таблица:

х == а	fabs(x-a) <= eps1	х == 0	fabs(x) <= eps
х != а	fabs(x-a) > eps1	х != 0	fabs(x) > eps
х < а	|х < а – eps1	х < 0	x < -eps
х <= а	х <= а + eps1	х <= 0	x <= eps
х > а	х > а + eps1	х > 0	x > eps
х >= а	х >= а – eps1	х >= 0	x >= -eps

В этой таблице рядом с каждым обычным сравнением показана его эквивалентная форма с учетом точности представления чисел, причем

eps1 = eps * (1.0 + fabs(a))

Несмотря на все советы автора, здесь нет причин для страхов, а код, подобный приведенному в таблице демонстрирует иногда замечательно странное поведение.