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 i 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 разів менше, ніж похибка попередньої програми. Причому похибка майже не залежить від величини кута. Таким чином, можна стверджувати, що власне додавання кутів виконується точно. Наявна похибка виникає за причиною нормалізації результату через ділення на с, яке не може бути точним.