5. Использование пифагоровых троек в вычислениях на компьютере

Часто в компьютерных задачах нужно вычислить таблицу тригоно­метрических чисел для равномерно изменяемых углов. Для этого удобно использовать формулы сложения углов (8). При этом следующую пару синуса X_i+1 и косинуса Y_i+1 вычисляют итерационно на основе предыдущей пары:

          X_i+1 = X_i·cosβ + Y_i·sinβ;                   (9)
          Y_i+1 = Y_i·cosβ − X_i·sinβ,

где X₁ = 0, Y₁ = 1, X_i+1 = sinіβ, Y_i+1 = cosіβ, і = 1,2,…,N.

Но при достаточно большом количестве итераций N в результатах накопляется слишком большая погрешность. Эта погрешность вызвана двумя причинами. Во-первых, на каждой итерации, чтобы не возникло переполнение результатов умножения, в них отбрасываются младшие разряды, значение которых равняется погрешности. Во-вторых, неточное представление коэффициентов sinβ, cosβ также вызывает погрешность.

Для исследования этих погрешностей была разработана программа Circle.PAS, которая строит окружность из точек с координатами (R·cosіβ, R·sinіβ) с радиусом R пикселов, в котором для примера N = 200 точек. Вычисления ведутся в целых числах, которые в компьютере представлены в диапазоне от −2³¹+1 до 2³¹−1. Для того, чтобы не возникало переполнение при умножении, числа должны быть представлены в диапазоне от −2¹⁵+1 до 2¹⁵−1, то есть от −32767 до 32767. Данные соответственно масштабируются, так что 2¹⁵ означает число 1,0. После вычислений в правых сторонах формул (10) результаты делятся на 2¹⁵, чтобы сохранить масштаб результата. Коэффициенты синуса и косинуса равняются 2¹⁵·sin(2π/N) ≈ 1029 и 2¹⁵·cos(2π/N) ≈ 32751.

Выполнение этой программы показало, что с каждой итерацией стремительно возрастает погрешность вычислений. На рис. 5,а показано положение рациональных точек в системе прямоугольных координат, сгенерированых этой программой на протяжении 5000 итераций. Видим, что из-за погрешности вычислений гипотенуза треугольника сошлась практически к нулю. Причем длина катета треугольника уменьшается наполовину за 1400 итераций. Если допустить, что погрешность накопляется линейно, то ее значение составляет 0,5 за N_i = 1400 итераций. То есть погрешность одной итерации равняется δ_i = 0,5/1400 ≈ 36·10⁻⁵. Это приблизительно 12 значений младшего разряда представления данных в программе.

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

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

Такой алгоритм был реализован в программе Circle_Pyth.PAS, которая выполняет то же самое, что и предыдущая программа, но с использованием пифагоровых троек. Результат работы этой программы после 5000 итераций показан на рис. 5,б. Как видим с рис. 5, погрешности вычислений стали значительно меньше, чем в предыдущей программе. После 50000 итераций для того же самого угла и угла 2π/13 получим результаты, как показано на рис. 6,а и б. При этом угол 2π/200 и 2π/13 получается в тройках (508,16125,16133) и (9940,18939,21389), которые задают эти углы с погрешностями 7,8·10⁻⁵ и 0,2·10⁻⁵ , соответственно.

                                    а                                                                      б

Рисунок 5. Положение рациональных точек в пространстве после сложения углов по обычному алгоритму (а) и с использованием пифагоровых троек (б)

                                    а                                                                      б

Рисунок 6. Положение рациональных точек в пространстве после 50000 сложений углов величиной π/100 (а) і 2π/13 (б)

Анализ результатов показывает следующее. Погрешность сложения углов в примерах оказалась минимальной — точки в пространстве формируют линии, которые приближаются к радиальным лучам. Существенную погрешность составляет пропорциональное уменьшение сторон треугольников δ_i = с_і+1 − с_і.

Допустим, что эта погрешность пропорциональна количеству итераций. Тогда, соответственно с рис. 6, за N_і = 50000 итераций гипотенуза треугольника уменьшается приблизительно вдвое. Это значит, что погрешность можно оценить как δ_i = 0,5/50000 = 10⁻⁵. Это приблизительно 1/3 значения младшего разряда представления данных в программе и в 36 раз меньше, чем погрешность предыдущей программы. Причем погрешность почти не зависит от величины угла. Таким образом, можно утверждать, что собственно сложение углов выполняется точно. Имеющаяся погрешность возникает по причине нормализации результата через деление на с, которое не может быть точным.