Приближенная рациональная арифметика с контролируемыми ошибками округления Текст научной статьи по специальности «Математика»

Вычислительные технологии Том 6, № 5, 2001

ПРИБЛИЖЕННАЯ РАЦИОНАЛЬНАЯ АРИФМЕТИКА С КОНТРОЛИРУЕМЫМИ ОШИБКАМИ ОКРУГЛЕНИЯ *

Г. Л. Литвинов, А. Я. Родионов, А. В. Чуркин Международный центр "Софус Ли", Москва, Россия e-mail: [email protected]

We describe an approximate rational arithmetic with round-off errors (both absolute and relative) controlled by users. The rounding procedure is based on the continued fraction expansion of real numbers. Results of computer experiments are given in order to compare efficiency and accuracy of different types of approximate arithmetics and rounding procedures.

Введение

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

Для повышения точности вычислений можно использовать различные варианты рациональной арифметики, которая оперирует с рациональными числами, представленными в форме отношения p/q целых чисел p и q (q > 0). Использование точной рациональной арифметики в принципе возможно (см., например, [1]), но, как правило, приводит к катастрофическому росту величин (и длин) числителей и знаменателей вычисляемых чисел, а также времени счета. Детально исследована (см., например, [2, 3]) приближенная рациональная арифметика, в которой фиксированы максимально возможные длины числителя и знаменателя (fixed-slash number systems) или сумма их длин (floating-slash number systems); при этом механизм округления использует разложение чисел в цепные дроби.

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, грант №99-01-01198. Материалы были доложены на XVI конференции по интервальной математике, Красноярск, 17-19 августа 1999 г.

© Г. Л. Литвинов, А. Я. Родионов, А. В. Чуркин, 2001.

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

Ниже описана приближенная рациональная арифметика с другим способом округления. Пользователь задает величины А и 8 абсолютной и относительной ошибок округления, где 0 < А < то, 0 < 8 < то. В частности, при А = 8 = 0 возникает точная рациональная арифметика. Если 8 = то, то фиксирована лишь абсолютная погрешность округления А. Если фиксирована лишь относительная погрешность 8, то возникает примерно такая же картина, как при вычислениях с плавающей запятой. Процедура округления применяется к дробям, у которых длина числителя или знаменателя превышает число M, задаваемое пользователем. При округлении число заменяется его наилучшим приближением, полученным путем разложения в цепную дробь, с заданными погрешностями А и 8. Результат округления всегда является несократимой дробью, но в отдельных случаях может совпадать с округляемым числом.

Арифметика этого типа была первоначально реализована средствами системы компьютерной алгебры REDUCE и использована для построения рациональных аппроксимаций (произвольной точности) к функциям от одной переменной [4]. В настоящее время такая арифметика реализована на языке программирования C++ (в форме, удобной для объектно ориентированного программирования в рамках проекта [5]) и средствами системы МАТЕМАТИКА 3.0. Приближенную рациональную арифметику произвольной точности предлагал использовать Ю. В. Матиясевич для своего апостериорного интервального анализа [6, 7].

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

1. Цепные дроби и приближенная рациональная арифметика

Напомним основные понятия теории цепных дробей. Обозначим через [а0; а\, а2,...] цепную дробь вида

1

ао +--1-. (1)

а1 +----

а2 + ...

Любое неотрицательное рациональное число ,p/q допускает единственное каноническое

представление в виде конечной цепной дроби

,

- = [ао; а1,а2, ...,а„], (2)

q

где все величины а^ (г = 0, ...,п) суть целые неотрицательные числа, причем а^ > 1 при г = 1,..., п — 1 и ап > 2 при п > 1. Аналогичное представление в виде бесконечной цепной

дроби вида (1) допускают иррациональные числа. Подходящая дробь pk/Qk порядка k для разложения (2) определяется равенством

— = [ао; аьа2 ,...,ßfc ], (3)

Qk

где k < n. Из (3) ясно, что подходящая дробь pn/qn совпадает с исходным числом p/q. Из теории цепных дробей [8] хорошо известно, что подходящая дробь вида (3) является наилучшим приближением к числу вида (2) в следующем смысле: для любой дроби r/s и подходящей дроби pk/Qk таких, что 0 < s < Qk и r/s = pk/Qk, справедливо неравенство

rp > Pk _ p

sQ Qk Q

Единственным и тривиальным исключением является случай р^ = а0 + 1/2, когда а0 и а0 + 1 одинаково хорошо приближают число р^. Для любой подходящей дроби рк/qk при к = 0 и к < п справедливо неравенство

Qk (Qk + Qk+i)

<

Р _ Pk

Q Qk

< — . (4)

QkQk+1

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

Для этой цели известен следующий алгоритм построения подходящей дроби [8]. Пусть задана дробь р^, где р, q — целые числа, р > 0, q > 1 (если р^ < 0, то приближается дробь |р^| и результат умножается на —1).

• Положим Ь-2 = р, Ь-1 = q, р-2 = 0, р-1 = 1, q-2 = 1, q-1 = 0.

• Для г = 1, 2, ... величины аг и Ьг последовательно вычисляются как частное и остаток от деления Ьг-2 на Ьг-1:

Ьг-2 = агЬг-1 + Ьг.

• Числитель и знаменатель подходящей дроби г-го порядка вычисляются по рекуррентным формулам

рг = агрг-1 + рг-2, qi = aiqi-l + qг-2.

• Если Ьг = 0, то приближение совпадает с приближаемым числом и процедура заканчивается.

• На каждом шаге (т. е. при каждом значении г = 0, 1, . . . ) проверяется, удовлетворяет ли приближение заданному критерию точности (см. ниже). При положительном ответе процедура заканчивается, в противном случае нужно сделать следующий шаг, увеличив г на 1.

В качестве критерия точности приближения выбирается любое из следующих условий:

1) абсолютная ошибка меньше А;

2) относительная ошибка меньше

3) одновременно выполняются эти оба условия.

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

1

2. Оценка ошибки округления

Из рекуррентной формулы для вычисления знаменателя подходящей дроби рк/Цк видно, что величины Цк минимальны при каждом фиксированном значении индекса к, если а = 1 при г = 0,1,..., к, т. е. ц = + Ц_2 = 1, = 1, Ц0 = 9-1 + Ц_2 = 1- Следовательно, для любой подходящей дроби рк /Цк справедливо неравенство

Цк > ^к+ъ

(5)

где /к — числа Фибоначчи, определяемые рекуррентными соотношениями /к = /к_1+/к_2 (к > 2) при /0 = 0, = 1. Хорошо известно, что числа Фибоначчи выражаются формулой

Я

„ = -1=(Фк - фк),

(6)

где Ф = 2(1 + л/5) ~ 1.618 ("золотое сечение"), ф = 2(1 — л/5) ~ -0.618. Также хорошо известно, что Ф — число, особенно плохо приближаемое своими подходящими дробями. Из (5) и (6) следует, что

ЦкЦк+1 > ^к+1^к+2 = ^(Фк+1 — фк+1)Фк+2 — Фк+2) =

5

1 /Ф2к+1 — (ФФ)к+1(Ф + ф) + ф2к+1^ = 1

5

2к+1 \ _ _ / ф2к+3

5

+ (—1)к + Ф

к , ф 2к+3

поскольку фф = —1 и ф + ф= 1. С учетом числовых значений Ф и ф получаем, что для любого к справедлива оценка

ЦкЦк+1 > ^Ф2к+2,

а для четных к — оценка

5

ЦкЦк+1 > ^Ф2к+3. 5

Из этих оценок и неравенства (4) следует, что абсолютная ошибка приближения будет

1п 5 1 , ,15 3

меньше Д, если к > — —— 1, а для нечетных к это верно при к > — —---.

2 Д 2 Д 2

3. Оценка числа итераций

Приведенные соотношения дают оценку сверху для количества итераций, требуемых для построения рационального приближения с абсолютной ошибкой, меньшей Д. Пусть Рк/Цк — подходящая дробь с требуемой ошибкой, а дробь рк_1/Цк_1 еще не дает требуемой точности. Для любого числа г обозначим через [г] целую часть числа г, т.е. наибольшее целое число, меньшее или равное г, а через [г] - наименьшее целое число, большее или равное г. Если при этом г целое, то [г] = [г], и [г] = [г] + 1 = [г + 1] — в противном случае. Тогда требуемая оценка сверху имеет вид

к

1] 5 ■ 2 §ф Д — 1

15

< 2 §ф Д

(7)

Если число

21о§ф А — 2

четно, то оценка (7) допускает усиление:

к

1] 5 21о§ф А

<

^ 5 21о§ф А

(8)

Для случая А = 10-М оценки (7) и (8) дают следующий результат.

Теорема. Если абсолютная ошибка, указанная в критерии точности округления, имеет вид А = 10-М, то справедлива оценка

к < [а + ЬЖ],

(9)

где а = ц 5 ~ 1.672 и Ь = 10 ~ 2.392. Если число

оценка (9) допускает усиление

а + ЬЖ — ^

четно, то

к

а — ^ + ЬЖ

(10)

Например, для N = 8 оценка (9) показывает, что к < 20, для N = 9 она дает к < 23, но в этом случае применима и оценка (10), так что к < 22.

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

На эвристическом уровне нетрудно оценить и среднее значение к параметра к при заданной абсолютной ошибке А. Рассматривая подходящие дроби рк^к к произвольному действительному числу х, А. Я. Хинчин доказал, что почти для всех х справедливо соотношение Ишк^те = 7, где 7 — константа (этот результат приведен в [10]); как показал

П. Леви ([9], с. 320), 1п 7

п

1.18657... , т. е. ^ ^^

3.27582 ... Грубо говоря, этот

121п2

результат означает, что при достаточно больших значениях к знаменатель qk подходящей дроби "в среднем" близок к 7к. Допуская вольность по отношению к математической строгости, подставим в неравенство (4) вместо qk и qk+1 соответственно 7к и 7к+1, чтобы оценить среднее значение к при заданной ошибке приближения А. В качестве верхней гра-

1п(1/А) 1

ницы получим число —----, а нижняя граница отличается от верхней на величину

21п 7

2

1п(1 + 1/7)

21п 7

1п(1/А) ^ близка к —--. Если А

21п 7

0.11. Таким образом, величина к (которая не обязана быть целым числом) 10-м,

то

- 1п(1/А) N 1п 10

к--= —-« 0, 97 ■ N - N.

21п 7 21п 7

(11)

Эта оценка становится реалистической лишь при больших значениях N. В противном случае к существенно меньше, чем N.

4. Примеры вычислений с использованием различных вариантов рациональной арифметики

Для сравнения различных вариантов рациональной арифметики рассмотрим классиче-

• п

ский модельный пример вычисления функции sin x в точках xm = —+ 2nm путем сумми-

6

рования ее ряда Тейлора. Ряд суммируется до тех пор, пока очередной его член (который отбрасывается) не становится по абсолютной величине меньше, чем 10-7. При этом число

355

п заменено его приближением с абсолютной ошибкой 2.7 • 10-7.

Для вычислений использовался процессор Pentium MMX с тактовой частотой 200 МГц. Различные варианты рациональных арифметик моделировались с использованием арифметики целых чисел произвольной длины, реализованных в системе аналитических вычислений МАТЕМАТИКА 3.0 (применение для этих целей системы REDUCE или языка C++ приводит к аналогичным результатам).

Рассматриваются следующие варианты десятичной арифметики рациональных чисел:

I. Точная рациональная арифметика.

II. Приближенная рациональная арифметика, описанная выше в разделе 2, при M = 9, А = 10-8, 6 = то (так что фиксируется только абсолютная погрешность округления А = 10-8); контроль ошибок округления ведется без использования неравенства (4).

III. Эта же арифметика при M = 9, А = 6 = 10-8.

IV. Эта же арифметика при M = 9, А = то, 6 = 10-8 (так что фиксируется только относительная погрешность округления).

V. Арифметика типа fixed slash [2, 3] с фиксированной максимальной длиной L числителя и знаменателя дроби при L = 6.

VI. Эта же арифметика при L = 9.

VII. Эта же арифметика при L = 12.

VIII. Арифметика типа floating slash [2, 3] с фиксированной максимальной суммой S длин числителя и знаменателя дроби при S =12.

IX. Эта же арифметика при S =15.

X. Эта же арифметика при S = 18.

Хи'Редукционная" арифметика, в которой округление состоит в отбрасывании младших разрядов с сохранением D десятичных знаков в числителе и знаменателе при D = 9. Результаты численного эксперимента представлены в таблице.

Для вычисления значений sin ^+ 2пт^ при m = 0,1, 2, 3, 5, 6, 7 указаны абсолютная

погрешность результата е (относительная погрешность совпадает с 2е), время вычисления t в секундах и сумма s длин числителя и знаменателя результата.

Значения погрешностей округлены до первой значащей цифры для того, чтобы результаты поместились в формат таблицы. С этой же целью пропущен столбец, соответствующий значению m = 4. Как видно из таблицы, для точной рациональной арифметики стремительный рост параметра s приводит к катастрофическому росту времени вычислений. Любопытно, что точная арифметика не всегда дает наибольшую точность результата вычислений. Отчасти это объясняется погрешностью представления числа п, а отчасти

связано с тем, что sin ( —+ 2пп) = — представляет собой рациональную дробь особо про-

62

стого вида. Природа этого эффекта обсуждается в [2]. Если фиксирована лишь относительная ошибка округления (вариант IV), то полностью повторяется картина, характерная

для арифметики с плавающей запятой: погрешность быстро растет и, начиная с т = 5, становится больше вычисляемой величины (ср. [10], гл. 3).

Сравнение различных рациональных арифметик

m 0 1 2 3 5 6 7

I е 4■10-8 5■10-7 10-6 3■10-6 7■10-6 10-5 4■10-5

s 62 214 372 504 820 995 1190

t 0.03 0.15 0.44 2.0 16.9 27.9 38.0

II е 2■10-8 5■10-7 10-6 10-6 2■10-6 3■10-6 3■10-6

А = 10-8 s 16 13 12 12 12 11 11

t 0.05 0.33 0.66 0.83 1.65 1.98 2.41

III е 4■10-8 5■10-7 10-6 10-6 2■10-6 3■10-6 3■10-6

А = 5 = 10-8 s 15 13 12 12 12 11 11

t 0.06 0.44 0.77 0.86 1.75 2.12 2.56

IV е 4■10-8 3■10-7 3■10-4 0.2 0.8 1.17 0.12

5 = 10-8 s 15 13 9 9 8 8 9

t 0.05 0.38 0.61 0.66 1.37 1.94 2.34

V е 0 0 10-3 0.7 1.4 3.4 5.6

L = 6 s 2 2 12 12 11 12 12

t 0.07 0.44 0.6 0.6 0.77 0.94 1.04

VI е 4■10-8 5■10-7 9■10-7 0.01 0.3 0.5 1.0

L = 9 s 17 18 17 18 16 18 17

t 0.10 0.72 1.19 1.25 1.38 1.42 1.51

VII е 4■10-8 5■10-7 9■10-7 10-6 2■10-4 0.01 0.1

L = 12 s 24 24 23 23 23 24 24

t 0.17 1.16 1.98 2.14 2.14 2.19 2.42

VIII е 4■10-8 5■10-7 10-6 10-4 0.06 0.8 0.9

S = 12 s 11 11 10 11 10 11 11

t 0.06 0.61 1.82 1.1 1.23 1.27 1.38

IX е 4■10-8 5■10-7 10-6 7■10-6 0.02 0.5 0.6

S = 15 s 14 14 13 13 13 13 14

t 0.16 0.87 1.42 1.48 1.81 1.87 2.14

X е 4■10-8 5■10-7 10-6 10-6 3■10-5 0.01 0.2

S = 18 s 17 17 15 17 17 17 16

t 0.11 0.94 1.54 2.04 2.53 2.59 2.62

XI е 0 3■10-5 0.04 0.1 0.4 0.9 1.2

D = 9 s 18 18 18 18 18 18 18

t 0.04 0.12 0.41 0.58 1.07 1.35 1.58

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

Ниже приведены значения параметра k, оценивающего среднее число итераций (см. раздел 4) в зависимости от заданной абсолютной погрешности А = 10-N при N = 16, 18, 20, ..., 36:

N 16 18 20 22 24 26 28 30 32 34 36

к 13.9 16.4 18.9 20.9 22.7 24.7 26.5 28.4 30.9 33.0 34.9

Вычисления проводились для описанного модельного примера при m = 9. Видно, что соотношение (11) дает реалистическую оценку для k, если число N достаточно велико.

Заключение

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

Список литературы

[1] HOWELL J., GREGORY R. T. Solving linear equations using residue arithmetic. Pt II // Nordisk Tidskr. Inf. 1970. Vol. 10. P. 23-37.

[2] MATULA D.W., KORNERUP P. Approximate rational arithmetic systems: analysis of recovery of simple fractions during expression evaluation // Lecture Notes in Comp. Sci. 1979. Vol. 72. P. 383-397.

[3] MATULA D. W., KORNERUP P. Finite precision rational arithmetic: slash number systems // IEEE Transaction on Computers. 1985. Vol. C-34, No. 1. P. 3-18.

[4] KRYUKOV A. P., Litvinov G.L., RODIONOV A. Ya. Construction of rational approximations by means of REDUCE // Proc. of the 1986 Symp. on Symbolic and Algebraic Computation. Symsac'86 (July 21-23, 1986, Waterk«, Ontario), Univ. of Waterloo, 1986. P. 31-33.

[5] Litvinov G. L., Maslov V. P., RODIONOV A. Ya. Unifying Approach to Software and Hardware Design for Scientific Calculations. E-print LANL arch^, quant-ph/9904024. 1999. (http://arXiv.org).

[6] МАТИЯСЕВИЧ Ю. В. Вещественные числа и ЭВМ // Кибернетика и вычисл. техника. 1986. Вып. 2. С. 104-133.

[7] MATIJASEVICH Y. A posteriori version of interval analysis // Topics in the Theoretical Basis and Appl. of Comp. Sci. M. Arato, I. Katai, L. Varga (Eds.) / Proc. Fourth Hung. Computer Sci. Conf. Budapest: Acad. Kiado, 1986. P. 339-349.

[8] Хинчин А. Я. Цепные дроби. М.: Физматгиз, 1961.

[9] Levy P. Theorie de l'addition des variables aletoires. Paris, 1937.

[10] МАК-КРАКЕН Д., ДОРН У. Численные методы и программирование на Фортране. М.: Мир, 1977.

Поступила в редакцию 23 июля 2001 г.