 |
Машинная арифметика и идиомы численного программирования
|
|
|
|
Научно-технические расчеты — простейшее применение вычислительной техники.
// Эдсгер Дейкстра
Большая часть ловушек, связанных с реализацией численных алгоритмов, сосредоточена в области вопросов точности представления вещественных чисел, машинных нулей и контроля сходимости итерационных процессов. Здесь мы рассмотрим примеры, демонстрирующие особенности представления вещественных чисел и некоторые идиоматические конструкции, используемые при численном программировании.
Представление вещественных чисел
В большинстве современных компьютеров аппаратно реализована вещественная арифметика, соответствующая стандарту IEEE. В соответствии со стандартом существует 2 основных представления вещественного числа — с одинарной (float) и двойной (double) точностью. Во всех представлениях один бит используется как знаковый. Есть еще одно представление вещественных чисел — 10-байтовое, используемое при вычислениях для хранения промежуточных результатов и называемое представлением с расширенной точностью.
тип длина мантисса показатель
float 4 байта 23 бита 8 битов
double 8 байтов 52 бита 11 битов
Машинное представление вещественного числа с плавающей точкой имеет следующий вид: +/- M * 2E, где M — вещественная мантисса, E — целочисленный показатель. Мантисса представлена в виде последовательности разрядов двоичной дроби:
M = M0 + M1/2 + M2/4 +... + Mn/2n
Для представления всех вещественных чисел кроме нуля используется нормализованная форма, когда нулевой разряд мантиссы M0 равен 1. При таком представлении нулевой разряд опускается и хранятся только разряды M1... Mn. Таким образом, мантисса M удовлетворяет неравенству 1 <= M < 2.
|
|
|
Машинное представление вещественных чисел содержит конечное число элементов. Только некоторые вещественные числа могут быть точно представлены описанным выше способом. Такие вещественные числа называются точно представимыми или просто представимыми в используемой машинной арифметике. Даже простые на вид рациональные числа такие как 0.1 оказываются непредставимыми. Непредставимое число заменяется ближайшим представимым числом.
Простые примеры
Непредставимые в машинной арифметике числа оказываются заданными приближенно и точность приближения зависит от того, какая используется арифметика — одинарной или двойной точности. Кроме этого следует помнить, что все промежуточные вычисления происходят в арифметике расширенной точности. Поэтому в следующем примере
float x4 = 0.1;
double x8 = 0.1;
не стоит надеяться, что x4 == 0.1 или x8 == 0.1. Дело в том, что при сравнении используется 10-байтовое представление для константы 0.1, не совпадающее с 8-байтовым в случае переменной x8 и 4-байтовым для x4. Ясно, что по аналогичной причине не будет выполнено равенство x4 == x8.
Последовательные операции с непредставимыми числами могут приводить к накоплению погрешности. Например, если просуммировать 1000 4-байтовых представлений (float) числа 0.1, то мы получим результат, отличающийся на 0.000953674 от точного. Это — вполне заметная погрешность.
Просуммируем вещественные числа вида 1.0/j при j = 1...1000 в разном порядке — по возрастанию и по убыванию. При использовании вещественной арифметики однократной точности (float) полученные результаты будут отличаться на 6.67572e-006.
Первые выводы
При работе с вещественными числами на точное равенство надеяться не стоит даже при проведении простейших вычислений.
Погрешности, связанные с неточным представлением вещественных чисел, имеют тенденцию накапливаться.
В машинной вещественной арифметике нарушаются даже простейшие свойства чисел (например коммутативность сложения), на которые мы привыкли надеяться.
|
|
|
Модельная арифметика
Для построения примеров удобно использовать следующую модельную арифметику. Пусть наш воображаемый компьютер использует десятичную систему счисления, причем мантисса имеет 4 десятичных разряда, а показатель — 2 десятичных разряда. Вот одно из вещественных чисел, представимых в такой арифметике: 9.876e-12. Ясно, что таких вещественных чисел имеется лишь конечное число, причем в нашей арифметике есть минимальное и максимальное числа:
min = -9.999e99,
max = 9.999e99.
Есть также минимальное положительное число:
min_pos = 1.000e-99.
Потеря значимости
Рассмотрим в описанной выше модельной арифметике операцию сложения чисел 1.234e00 и 5.678e-04:
1.234e00 + 5.678e-05 = 1.234 + 0.00005678 = 1.23405678
Это число непредставимо в нашей арифметике и будет заменено ближайшим представимым:
1.234e00 + 5.678e-05 = 1.234e00
Таким образом, при сложении чисел разного порядка результат может быть равен одному из слагаемых. Этот эффект называется потерей значимости.
Машинное эпсилон
Рассмотрим более подробно операцию сложения двух положительных чисел. Пусть x > y > 0. Тогда
x + y = x * (1 + y/x)
Если число y/x окажется настолько маленьким, что 1 + y/x = 1, то мы получим потерю значимости: x + y = x.
Машинным эпсилон называется минимальное положительное число, которое при сложении с числом 1.0 дает результат, больший, чем 1.0. Машинное эпсилон — наиболее важная характеристика вещественной арифметики.
В нашей модельной арифметике машинное эпсилон равно 0.001 = 1.000e-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) > 1.0.
|
|
|
Параметры машинной арифметики
Основные параметры машинной арифметики — максимальное представимое число, минимальное положительное число и машинное эпсилон. Максимальное представимое число и минимальное положительное число определяют диапазон чисел, представимых в машинной арифметике, называемый также динамическим диапазоном арифметики. Динамический диапазон определяется разрядностью показателя. Машинное эпсилон определяет точность представления чисел и зависит от разрядности мантиссы.
Несмотря на то, что для всех основных параметров машинной арифметики существуют переносимые способы их вычисления, обычно предпочитают пользоваться готовыми значениями соответствующих констант. Заголовочный файл float.h содержит определения макросов, дающих основные параметры машинной арифметики.
float double
максимальное представимое число FLT_MAX DBL_MAX
минимальное положительное число FLT_MIN DBL_MIN
машинное эпсилон FLT_EPSILON DBL_EPSILON
Идиома 1: проверка на малость
Пусть x — вычисляемое значение, а e — его погрешность. Тривиальная реализация проверки погрешности на малость оператором
fabs(e) < eps
часто приводит к проблемам в случае больших по абсолютной величине значений числа x. Дело в том, что вещественные числа расположены на вещественной оси неравномерно — гуще вблизи нуля и реже вдали (между любыми соседними целыми степенями двойки равномерно расположено одинаковое количество вещественных чисел). Поэтому даже два соседних вещественных числа могут находиться на расстоянии, большем, чем eps. В численных алгоритмах правильная постановка вопроса о малости ошибки звучит так: «является ли погрешность малой по сравнению с ответом». Такая формулировка приводит к следующей реализации:
fabs(e) < eps * fabs(x).
В этой реализации абсолютный порог малости числа зависит от значения величины, с которой она сравнивается — при больших по модулю значениях x даже относительно большие e будут признаны малыми. Действительно, число 1 мало по сравнению с 1e6 в той же степени, что и число 1e-6 по сравнению с 1. Оказывается, что и эта реализация может приводить к проблемам, но уже в случае малых по абсолютной величине чисел x. Если x — почти нуль, то и правая часть получается почти нулевой. В этом случае правая часть может оказаться настолько малой, что никакое e не сможет удовлетворить неравенству. Эту проблему решает такая реализация:
|
|
|
fabs(e) < eps * (1.0 + fabs(x)).
Защититься от возможности задания излишне малых значений eps можно следующим образом:
fabs(e) < (eps + 2.0 * FLT_EPSILON) * (1.0 + fabs(x)).
Самым важным соображением здесь является следующее: ***в большинстве случаев следует формулировать численные алгоритмы таким образом, чтобы они содержали только относительные проверки на малость***. Наиболее часто такая проверка встречается в следующей ситуации:
float x; // вычисляемая величина
float e; // ошибка ответа
x =...; // вычисляем первое приближение
e =...; // вычисляем ошибку
// пока ошибка велика по отношению к ответу
while(fabs(e) > (eps + 2.0 * FLT_EPSILON) * (1.0 + fabs(x))){
x =...; // вычисляем новое приближение
e =...; // вычисляем новое значение ошибки
}
Идиома 1': сравнение на равенство
Проблема со сравнением вещественных чисел на равенство решается просто — нужно всегда использовать только сравнение на приближенное равенство т.е. сравнивать числа с учетом точности их представления. При сравнении чисел x и y на равенство мы фактически проверяем разность x - y на малость. Остается только один вопрос — с каким числом следует сравнивать разность x - y. Если числа x и y равноправны, то обычно их считают равными, если разность x - y мала по сравнению с максимальным по абсолютной величине из чисел x и y:
fabs(x - y) < (eps + 2.0 * FLT_EPSILON) * (1.0 + fmax(fabs(x), fabs(y))),
где fmax() — функция, возвращающая значение максимального из своих аргументов.
Пример: метод бисекции
Метод бисекции — простейший метод нахождения корня уравнения f(x)=0, лежащего на отрезке [a,b]. Будем считать, что функция f(x) — непрерывна и принимает на концах отрезка значения разных знаков. Метод состоит в уменьшении отрезка неопределенности вдвое на каждом шаге с сохранением условия знакопеременности. Программа пишется сразу:
float fa = f(a);
while(b - a > 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 одновременно весьма малы) даже для чисел разного знака мы получим ноль. Следует проверять именно условие на знаки чисел, страхуясь таким образом от неприятностей.
|
|
|
Вторая проблема связана с проверкой на приближенное равенство. Для больших чисел a и b проверка b - a > eps может привести к бесконечному циклу. Эту проверку следует заменить на следующее условие:
b - a > (eps + 2.0 * FLT_EPSILON) * (1.0 + fmax(fabs(a), fabs(b)))
Третья проблема связана с тем, что число m, определенное формулой m = (a+b)/2.0, может оказаться лежащим за пределами отрезка [a,b]. В качестве примера рассмотрим модельную 4-разрядную арифметику и пусть a = 0.9882 и b = 0.9884. Тогда a + b = 1.9766, но это число непредставимо в нашей арифметике и реально получится a + b = 1.977. Следовательно, m = (a + b)/2.0 = 0.9885 > b. Избежать этой неприятности можно, изменив формулу для вычисления середины отрезка [a,b]:
m = a + (b-a)/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 *x, 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 *x, 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);
}
Более существенным здесь оказывается даже не возможность переполнения, а существенная потеря точности ответа, поскольку при возведении в квадрат и последующем извлечении корня мы фактически теряем половину значащих цифр результата. В этом случае часто говорят о потере значимости. Масштабирование решает и эту проблему.
Аналогичная проблема возникает при вычислении корня из суммы квадратов двух чисел (норма вектора на плоскости, норма комплексного числа), но эта проблема решается использованием стандартной функции hypot(x, y), которая проводит вычисления по формуле sqrt(x*x + y*y) с контролем возможности переполнения и потери значащих цифр. При этом функция hypot(x, y) выполняется быстрее, чем вычисления по формуле с использованием масштабирования.
Пример: деление комплексных чисел
В качестве примера относительно нетривиального использования процедуры масштабирования приведем код, реализующий деление комплексных чисел:
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-разрядной арифметике два непредставимых числа X и Y оказались представленными числами x = 1.001 и y = 1.002. В обоих случаях относительная погрешность представления принимает значения порядка машинного эпсилон (eps = 0.001). Вычислим относительную погрешность разности z = y - x = 0.001:
|(y-x)-(Y-X)|/|Y-X| = 2 * eps / 0.001 = 2000 * eps.
Таким образом, вычисленная разность не содержит ни одной верной цифры. Нахождение разности близких чисел — потенциально очень опасная операция и единственное лекарство состоит в реорганизации вычислений таким образом, чтобы избежать необходимости производить вычитание для близких чисел. Вероятно самый известный, но и самый бесполезный с практической точки зрения пример такой реорганизации вычислений — задача корректной реализации вычислений по формулам для корней квадратного уравнения.
Примеры безосновательной паранои
Операции с возможным делением на ноль — одна из ситуаций, в которой параноидальный стиль программирования приводит лишь к ухудшению качества кода. Начинающие в случае возможного деления на ноль склонны излишне осторожничать, ограничивая значения делителя снизу по модулю. В действительности перед делением в большинстве случаев достаточно произвести проверку на точное равенство делителя нулю — только в этом случае делить числа нельзя. Для всех остальных значений делителя результат представим в используемой арифметике, а любая проверка на малость является излишней и часто существенно сужает динамический диапазон. Иногда все же стоит отслеживать, на какое минимальное число происходило деление в процессе вычислений.
В известной книге Л.Аммерала «Принципы программирования в машинной графике» можно найти большой список бессмысленно параноидальных советов. В этой книге приведена следующая таблица:
x == a fabs(x-a) <= eps1 x == 0 fabs(x) <= eps
x!= a fabs(x-a) > eps1 x!= 0 fabs(x) > eps
x < a x < a - eps1 x < 0 x < -eps
x <= a x <= a + eps1 x <= 0 x <= eps
x > a x > a + eps1 x > 0 x > eps
x >= a x >= a - eps1 x >= 0 x >= -eps
В этой таблице рядом с каждым обычным сравнением показана его эквивалентная форма с учетом точности представления чисел, причем
eps1 = eps * (1.0 + fabs(a)).
Несмотря на все советы автора, здесь нет причин для страхов, а код, подобный приведенному в таблице демонстрирует иногда замечательно странное поведение.
|
|
|
