Идиома 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)).

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

 

 

Задание на лабораторную работу

1. Разработать программный продукт, реализующий методы машинной арифметики в соответствии с вариантом задания.

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

3. Ответить на контрольные вопросы.

 

Вариант	Параметры машинной арифметики	Машинные операции
1	1	9
2	4	10
3	3	10
4	3	8
5	2	8
6	1	10
7	2	6
8	1	8
9	4	7
10	2	5
11	3	6
12	3	7
13	4	5
14	2	6
15	2	7
16	2	9
17	4	6
18	1	5
19	2	7
20	4	6

 

Задания

1. Найти машинное эпсилон, машинную бесконечность, минимальный и максимальный диапазон для двух вещественных типов данных. Использовать любой язык программирования. Сравнить результаты.

2. Для модельной арифметики с десятичной системой счисления, где мантисса имеет 4 десятичных разряда, а показатель — 2 десятичных разряда, 1 разряд выделен для знака определить степень погрешности, минимальное и максимальное значение, машинное эпсилон.

3. Для модельной арифметики с десятичной системой счисления, где мантисса имеет 8 десятичных разряда, а показатель — 3 десятичных разряда, 1 разряд выделен для знака определить степень погрешности, минимальное и максимальное значение, машинное эпсилон.

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

5. Реализовать операции сравнения вещественных чисел разных типов данных. Показать результаты обычного сравнения, сравнить числа с учетом значения машинного эпсилон. Проанализировать результаты.

6. Реализовать метод бисекции для нахождения корней квадратного уравнения. Модифицировать метод таким образом, чтобы выполнялась проверка на приближенное равенство. Проанализировать результаты.

7. Реализовать операцию умножения чисел с переполнением данных. Предложить модификацию метода для выполнения проверки переполнения данных.

8. Реализовать операцию деления чисел таким образом, чтобы можно было продемонстрировать потерю точности при вычислениях. Проанализировать результаты.

9. Реализовать операцию вычитания чисел таким образом, чтобы можно было продемонстрировать потерю точности при вычислениях. Проанализировать результаты.

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

 

4. Контрольные вопросы

1. Принципы хранения и представления вещественных чисел в компьютере.

2. Понятие машинного эпсилон.

3. Понятие машинного нуля.

4. Нахождение минимального и максимального возможных значений для диапазона данных.

5. Причины возможных ошибок при выполнении операций сравнения.

6. Причины потери точности при машинных вычислениях.

7. Причины переполнения данных и возможные подходы для решения проблемы потери точности.