О гарантированной точности решений задач вычислительной математики в арифметике с плавающей запятой и переменной длиной мантиссы Текст научной статьи по специальности «Математика»

УДК 537.322.2

А. Г. Бирюков1, А. И. Гриневич2,1

1 Московский физико-технический институт (государственный университет)

2000 «ГринМарк»

О о КМ

гарантированной точности решении задач вычислительной математики в арифметике с плавающей запятой и переменной длиной мантиссы

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

Ключевые слова: погрешность округления и точность решения задач вычислительной математики, машинное число с переменной длиной мантиссы, алгоритм типа Маркова вычисления значения функции, гарантированная точность решений задач ВМ.

1. Введение

В настоящее время в вычислительной практике при решении задач ВМ в арифметике с плавающей запятой преимущественно используются одинарная, двойная и четверная точности МЧ, что позволяет решать широкий круг практических задач. Однако существует множество задач, для решения которых четверной точности МЧ недостаточно по причине ошибок округления, и для их решения необходимо использование машинной арифметики с мантиссой числа большей длины. Классическим примером такой ситуации является задача решения систем линейных уравнений с плохо обусловленной матрицей. В настоящее время в бесплатном доступе получила распространение библиотека программ GNU GMP [4], реализующая стандарт IEEE 754 [1,9], в которой длина мантиссы в арифметике с плавающей запятой варьируется в широком диапазоне значений. Библиотека GNU GMP позволяет оперировать числами с длиной мантиссы от mm\n = 24 вплоть до mmax = 231 = 2 147 483 648 двоичных знаков, чему соответствует 8 и 646 456 993 десятичных знаков. Верхнее значение длины мантиссы МЧ mmax невообразимо огромно. В указанной библиотеке также реализована возможность динамического изменения длины мантиссы т в различных сегментах программы от rnmin до mmax. Появление в свободном доступе программного обеспечения с такими возможностями расширяет границы для получения решений широкого круга задач ВМ с гарантированной точностью высокого порядка.

Проблема анализа влияния погрешностей округления (ВПО) на решение задач ВМ актуальна со времени появления ЭВМ и остается таковой по сей день. Научные исследования над указанной темой ведутся в разных направлениях. Отметим классические работы по исследованию ВПО при решении задач линейной алгебры [2,3,7,11]; по исследованию ВПО в рамках интервального анализа [15,16,18]; по статистическому анализу ВПО [12,19]; исследованию новых моделей по выработке машинного числа [13]; алгоритмов с автоматической коррекцией ошибок округления первого порядка — метод CENA [14].

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

Для выделенного класса погрешностей получена оценка длины мантиссы МЧ, гарантирующей достижение требуемой точности решения.

Напомним определение машинного числа с плавающей запятой согласно стандарту IEEE 754 [1]:

Определение 1. Машинным числом будем называть число вида:

хт,р = ±0, sisis2...sm ■ be = ± (+ I2 + ... + ^) ■ Ье = ±р, ■ be,

где b £ {2, 3,...} — основание, у — мантисса числа, т — количество знаков в мантиссе (число знаков, длина, размер), Si £ {0, 1, ...Ъ — 1}, г £ [1, т] ^ значащие цифры, s1 = 0, е £ [emin, emax] — порядок чиела, р — размер порядка МЧ или количество знаков в представлении числа е. При заданных Ъ, т, р машинные числа образуют конечное множество, обозначаемое далее Мь>т>р. Очевидно, Mt,m,p С Q: где Q — множество рациональных чисел. Для вещественного числа х £ R ближайшее к нему машинное число хт,р есть машинное представление этого числа, т.е. число хт,р есть результат округления числа х. Число 5^ = Ъ~т назовем погрешностью (точностью) мантиссы или точностью машинных чисел (машинной арифметики). □

Машинное представление действительного числа х, полученное в результате его округления, имеет следующий вид [2,3]:

хт,р = Ж (1 + и) + V, (1)

где и ■ V = 0.

При V = 0 М ^ ^ и %т,р = 0; при U = 0 Хт,р = 0 и |ж| ^ So. Относительная погрешность представления числа ж при хт,р = 0 не бол ее $1 = 1 b1-m, абсолютная погрешность |ж — хт,р | ^ ^1 |ж|. Абсолютная погрешность при хт,р = 0 |ж — хт,р | = |ж| ^ So, So = 1 ■ bemin ~ b-bP— погрешность нуля, ближайшее к нулю число. На практике выбираются такие величины m и р, что можно всегда считать выполненным соотношение So << £i, а если учитывать что р = 64 [4], то 5о = 2-2 ми 5о = 10-5,55^10 .

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

Определение 2. Пусть р : Rn ^ R1 — некоторая функция, х £ Rn — вектор, хт — его машинное покомпонентное представление. Тогда: р(х), р(хт) — значения функции р в точках х и хт, рт(хт) — машинное представление значения ip (хт); рт (ж) = рт(хт)-Значения функций р (ж) и р (хт) в общем случае представляются бесконечным числом

□

Определение 3. Математические функции и операции, реализуемые в библиотеках программ стандарта С99, назовем базовыми или стандартными. К базовым операциям относятся: округление чисел, арифметические операции, логические операции, операции вычисления математических функций, таких как sin х, cos х, tan ж, и их обратные, ех, ах, ху, loga ж и т.д. □

Одной из библиотек, реализующих базовые функции стандарта С99, является GNU MPFR [4], в которой при заданной длине мантиссы т для значений базовых функций выполняется округление до последней значащей цифры, т.е.

Рт(хт) = Р(хт)(1 + и), |«| ^ 5l, рт(хт) = 0, (2)

где <р(хт) и <рт (хт) — значение одной из стандартных (базовых) функций, вычисленное в точкежт £ М5,т,р. Для случая, когда рт (хт) = 0 |<£ (хт) — рт (хт)1 = |<£ (хт)1 ^ fo.

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

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

Во-первых, появляется возможность постановки задачи о достижении заданной точности решения (или вычисления с гарантированной точностью): для заданного е > 0 найти решение задачи ВМ с погрешностью не более е.

Во-вторых, появляется возможность рассмотреть вопрос о «реабилитации» численных методов решения задач ВМ, признанных в вычислительной практике неустойчивыми по причине влияния ошибок округления.

В-третьих, возникает также вопрос о влиянии повышения точности МЧ на эффективность применяемых на практике методов. Например, можно ли, увеличивая длину мантиссы МЧ, уменьшить общее время вычислений? И если да, то в каких случаях?

В настоящей работе дается вариант ответа на первый вопрос.

2. Задачи ВМ и алгоритмы их решения

Существует несколько классификаций задач вычислительной (прикладной) математики [8]. Общим для них является то, что процесс численного решения задачи ВМ на ЭВМ представляет собой упорядоченную последовательность конечного числа вычислительных базовых операций, которая определяется алгоритмом решения данной задачи. Этот процесс можно представить как вычисление значения некоторой функции (вектор-функции, отображения) / £ Кк в точке х € С С Кп, где С — область определения этой функции. Определение 4. Задачей вычислительной математики будем называть совокупность понятий и условий Р(/, х, т, е), определяющих возможность вычисления значений вектор-функции £ : Кп ^ Кк,тр.е / (х) — решение данной задачи в точке х € С] е — точность решения, т — длина мантиссы. Условие достижения точности означает: найти значение вектор-функции /т(х) при заданном х € С такое, что

\\1ш(х) - /(ж)|| , _

II/(*)\| " ‘

□

Возможны и другие определения задачи ВМ, но для целей настоящей работы определения 4 достаточно. Здесь и далее под нормой ||-|| понимается евклидова норма

ц^ц = ^I г € Ка. Очевидно, что к вычислению значений вектор-функции сводятся, например, следующие задачи: решение системы линейных уравнений, вычисление производных функции разных порядков, задача интегрирования функции, задача Коши для дифференциальных уравнений, задачи поиска экстремума функции ф(г), г € С С Кп, И т.д.

Для задач ВМ разработаны численные методы их решения, которые в форме, адаптированной к ЭВМ, представляют собой алгоритмы решения этих задач. По своей природе алгоритмы бывают конечношаговыми (КША) и бесконечношаговыми (БША). К первым, например, относится метод Гаусса решения систем линейных уравнений, ко вторым — метод Ньютона решения систем нелинейных уравнений. Алгоритмы КША и БША вычисления значений функции / (х) в точной арифметике представим в следующем виде:

КША:

АЬС(/(х)) = ^1(аг) V ^2(а2) V ■ ■ ■ V -1(ам-1) V (ам), (3)

БША:

АЬС(/(ж)) = ^1(а1) V ^2(а2) V ■ ■ ■ V -1(ам-1) V (ам) V ■ ■ ■ N ^ ж, (4)

где символ V означает объединение операций, рг(аг), г € [1, N] — базовые операции, аг €

— аргумент г-й операции. Операции округления в (3) и (4) отсутствуют, а АЬС (/ (х))

представляет собой упорядоченную последовательность базовых операций, определенную алгоритмом решения задачи ВМ. Результатом реализации базовой операции (кроме логической) является действительное число. В (3) N — число шагов КША.

Пусть требуется вычислить значение функции /(ж), х € Кп, / € Кк при длине мантиссы т. Тогда некоторый КША вычисления значения /т(х) представляет собой упорядоченную последовательность N базовых операций <ргт(а^, г € [1,^], т.е.

АЬС(/т(х)) = <р]п(а{) V ^2т(а2) V ■ ■ ■ V ^-1(а^-1) V (ам)■ (5)

В отличие от (3), (4) в (5) присутствуют операции округления чисел.

БША для вычисления значения /т (ж) имеет вид

АЬС(/т(х)) = ^1т(а1) V V ■ ■ ■ V ^т0-1(а^о-1) V (^Ио), (6)

где аналогично (5) N0 — число шагов алгоритма, которое в (6) определяется по некоторому правилу его окончания. Однако при решении конкретной задачи ВМ для вычисления значений /т(х) могут существовать различные алгоритмы. Введем следующее определение. Определение 5. Алгоритм вычисления функции / € Кк (5) в точке ж € Вп при длине мантиссы т будет называться нормальным алгоритмом, для решения задач вычислительной математики (НАВМ), если вычисленное значение функции в точной арифметике по алгоритму (5), в котором логические операции не выполняются, дает точное значение функции /(ж), т.е. имеет место алгоритм

АЬС(/(ж)) = ^1(а1) V ^2(а2) V ■ ■ ■ V -1(ам-1) V (ам). (7)

□

Замечание. Будем считать, что при построении НАВМ использовались логические операции как из математической библиотеки, так и операции языка программирования, на котором реализуется алгоритм решения задачи ВМ. Структура АЬС (/т (ж)) (5), т.е. упорядоченная последовательность выполнения базовых операций, определяется базовыми логическими операциями и для данного т становится фиксированной. При выполнении АЬС (/ (ж)) базовые логические операции в (7) <^>г (а^) = 0, т.е. не выполняются, а остальные операции в точной арифметике выполняются в соответствии со структурой алгоритма (5). Отметим также, что в (7) операции округления из (5) представляют собой единичный оператор.

Аргументы аг € в (5) — это либо числа (векторы), являющиеся аргументами функции / (ж), либо чиела а*•, Ь € [1,з] — результаты выполнения операции ^>m(aj), ] < %, ] € [1,^] на ]-т итерации, полученные в процессе вычислений в (5).

Введенное понятие нормального алгоритма терминологически напоминает понятие нормальный алгоритм, Маркова (НАМ) [6]. НАМ содержит понятие алфавит, на основании которого строятся дальнейшие выводы грамматики. Алфавит может быть произвольным. В случае НАВМ алфавитом можно считать базовые функции библиотеки программ. На этом основывается аналогия между НАВМ и НАМ. Таким образом, его можно считать сужением НАМ на класс задач вычислительной математики.

3. Погрешности решений задач ВМ

В настоящем разделе изучаются погрешности решений задач ВМ, возникающие в итерационном процессе НАВМ.

Лемма 1. Пусть ^ — базовые вычислительные операции (кроме логических) из некоторой библиотеки программ, г € [1,^1] — номер базовой операции. Тогда значение ^гт(аг) можно представить в виде

<Рт(аг) = ^(аг) + 0.^1,

(8)

где — число базовых операций библиотеки, |ац| ^ |^>г(аг)|, = 2Ь1-т, аг € аг = ат,г

— вектор машинных чисел.

Доказательство.

Сначала рассмотрим операции по вычислению значений базовых функций, для которых справедливо равенство (2): <ргт(аг) = ^>г(аг)(1 + щ), где |иг| ^ 5^ 51 = 2Ь1-т. Преобразуем

его, и получим равенство (8):

¥т(аг) = <Рг(аг) + <£г (аг) ■ Пг = фг(йг) + ^(йг) ^ ^1 = фг(йг) + 0^1, Где а» = ^ ^(а^),

% € [-1, 1] и Ы < ИМ-

Для операции округления числа справедливо:

Рт(аг) = (аг)(1 + «г) + ^г, ГД6 «г^г = 0 Ы ^ 5Ъ | | ^ ^0, ^0 = Ь-е™п(р\

В случае, когда Ьг = 0, имеем <р>г(аг) = т.к. г — единичный оператор,

<Ргт (аг) = ат1 = аг (1 + щ) = аг + а^ь где а* = ^а,г, ^ € [-1,1], т.е. справедливо (8).

В случае, когда иг = 0 <ргт (аг) = 0 и ргт(аг) = рг(аг) + ьг = аг + щ, где |^| ^ £0, £0 = Ьеш[п(р\ Здесь значение аг из (8) равно аг = где |аг| ^ — очень малое число. Для одинар-

ной точности I0 = 2 ■ 10-29, для двойной — I0 = 2 ■ 10-308 [1]. Таким образом (8) также справедливо и для операции округления. □

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

Лемма 2. Пусть для чисел й, у, г и их приближенных значений йт, ут, хт выполнены условия Ай = йт - d = Ау = Ут - у = ау, Аг = гт - 2 = @г] Ут = VУгт = Ц-хь\ и — мантиссы чисел утш гт, £ — порядок числа; ц.у - ^ = ,цт-д6-^ 1 ^ д <т, ^т-я —

мантисса числа ц.у - ц,г, где т - д её длина; £ = {51, а = аё^ [3 = ^51. Тогда погрешность

А = У^ - У- = X = 5^ Н« - ^)•

Доказательство.

А

й + А^ й А(1 (у - г) - й (Ау - Аг)

у - г + Ау - Аг у - г (Ут - гт)(у - 2)

1 - ё (ау - ) = _^ (| - ^= ь<хь. <Э)

V У - 2 / ^т-д# у - -

Ут - У - ^ У ^т-дУ - ^

Значение числа % приведено в условии леммы. □

Отметим, что число % может быть большим, если у - £ мало по сравнению с ау -Число д характеризует порядок потери точности вычислений. Желательно, чтобы оно было существенно меньше т, например д ^ т

Замечание. Вычисления в НАВМ, в которых будут использованы в качестве промежуточных значений величины у , будут иметь оценку погрешности, в которой сомножитель Ья может сохраниться, возможны случаи, когда значение д возрастет, уменьшится или станет равным нулю. Следовательно, при анализе погрешности вычислений необходимо каким-то образом учитывать этот эффект потери точности вычислений. На практике данный эффект возникает при численном дифференцировании функций, при решении плохо обусловленных систем линейных уравнений и т.д.

Лемма 3. Пусть функция р : К™ ^ К1 удовлетворяет условию Липшица:

1р (х + А х) - <р (ж)| ^ Ь У А ж ||; х,х + А ж € С, (10)

где О С Кп — компакт. Тогда существуют такие числа 1г € [-Ь,Ь], что

<р(х + А ж) - р(х) = ^2Г=11гА хг.

Доказательство.

Из (10) следует, что для каждой пары точек х + А ж, х существует число Ь1 € [-Ь,Ь] такое, ЧТО

(х + Ах) - <р (х) = Ь1 ||АхЦ. (11)

Представим (11) в виде

(х + А х) - (х) = ||Аж||2 = Е?=1 А хг = ЕГ=1 ь1агА хг = £Г=1 ^Ажг,

где аг € [-1,1], т.е. 1г € [-£,£]. □

Рассмотрим свойства погрешностей конечношаговых алгоритмов.

Теорема 1. Пусть функция / (х) € ж € С С Дга; базовые функции рг(аг)

(кроме логических) либо являются функциями округления числа, либо непрерывны по Липшицу, т.е. в некоторой окрестности 0,г(аг) точки аг € ^удовлетворяют условию:

<рг(аг) - рг(Ьг) ^ Ьг1аг - Ьг1, г € [1,^], а алгоритм (5) вычисления функции /(ж) является нормальным для х € С. Тогда существует такой вектор С € Кк, что

1т(х) - /(х) = СЙ1, (12)

где $1 = 2Ь1-т, т - длина мантиссы.

Доказательство.

Для доказательства утверждения достаточно показать, что

ат,г = аг + а^1 и ргт(ат,г) = ^г(«г) + №1, V г € [1,Ж] , (13)

кроме логических операций рг, где аг и @г — некоторые константы, аг — точное значение аргумента.

Доказательство проводится индукцией по номеру г € [1,^] вычислительной операции алгоритма (5) (кроме логической). В качестве базы индукции можно взять операцию округления числа ^ — первого числа, с которого начинаются вычисления. По Лемме 1: ат,1 = ^т(а1) = ^1(°1) + а1^1 = а1 + |°1| ^ |°1|- Предположим, что для неко-

торого г € [1,^] выполнено (13). Тогда необходимо доказать, что (13) выполнено и для операции с номером г + 1. Компонентами аргумента ятг+1 могут быть либо числа, входящие в условие задачи (аргументы функции /(ж)), либо числа а3тг+1 = Ч>\П(аг), .1 € [1,«], ^ ^ г, для которых условие (13) выполнено по условию индукции. Если операция г + 1 есть разность двух близких чисел одинакового порядка, то по Лемме 2 справедливо (9). Если операция г + 1 является округлением, то справедливо (8). Теперь надо показать, что (13) выполнено для функции рг+1, удовлетворяющей условию Липшица. По Лемме 1 <^+>т,г+1) = {Рг+1(ат,г+1)(1 + Щ+1), ГДе |^г+1| ^ Учитывая, что аг+1 € № и Лемму 3, получим

(ат,г+і) = <р

— ,^+і

(аг+і + аі+і^і)(1 + щ+і) = (^г+1(аі+і) + ^ а>і+1 ^+^(1 + щ+і)

з=і

= рг+і(аг+і) +

иі+і _ ^г+і

5і

рг+і(аі+і) + (1 + щ+і)^ аІ+А

г+і

3=і

6і = рг+і(аг+і) + Рг+і^і,

(14)

где |і^+і| ^ Ьг+і, Є [—1,1], т.е. для операции і + 1 выполнено (13). Таким образом,

равенства (13) выполнены V і Є [1,^]• Так как компоненты вектор-функции /іт — это некоторые числа ц Є [1,^], і Є [1,к], значения которых по доказанному представимы в виде (р%_ = ірч(аг1) + Ріі5і, то, переобозначая Ргі в Сь, получим рт(ж) = /*(ж) + с^і и Іт (х) = !(х) + С5і, ГДЄ С = (сі, С2, ..., Ск). □

Таким образом, система уравнений (15) для г € [1, N] представляет собой систему условий для рекуррентного оценивания погрешностей вычисления значений функции / (ж), х € С.

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

Следствие.

1) Вектор С в формуле (12) назовем параметром погрешности (ПП) значения функции. Он является функцией точки ж и размера мантиссы т. Структура его сложна и может включать элементы погрешности вида (9). Тогда погрешность С&1 можно представить в виде С§1 = [А^111, А2ЪЯ2,..., Ак Ъ1к] $1. Возьмем компоненту вектора Аг0, для

которой

Ло Ъ1*

= тах

А№

Обозначим qг0 = д. Тогда существуют Сг такие, что

С51 = (Сь ..., Ск) Ъ15ц = (Сь ..., Ск) (5^)а = С (^)а,щеа = 1 - т $1 = 2Я

ц) = ^ (^м) > а = 1 т и1 = 2иц-Таким образом, теперь формулу погрешности (12) можно представить в виде

А/т (х) = !т (х) - / (х) = С (5ц)а , 0 < а ^ 1. (15)

Когда д = 0, то а = 1 и А/т (ж) = /т (ж) -/ (ж) = С5^, где С = 2С. Оценка погрешности

(15) дана для одной точки. Её обобщением является определение 9.

2) Из Теоремы 1 можно получить следующие выводы.

Во-первых,

С

< ж, т.е. компоненты сг, г € [1,к], для фиксированного достаточно большого т, ограниченные сверху и снизу числа. Этот вывод следует из того, что каждая компонента сг получена по рекуррентной формуле (14) за конечное число шагов.

С

ж при т ^ ж.

Во-вторых, значение параметра С может быть таким, что Например, если д = (1 - а) т при фиксированном значении а, т ^ ж. Известно [7], что погрешность значения функции / величина случайная, зависящая от длины мантиссы т. Но так как при увеличении длины мантиссы т модуль погрешности значений базовых функций уменьшается, то и величина колебаний значений параметра ||С||, при определенных условиях, также может уменьшаться и следует ожидать, что его значения будут иметь колебания около некоторой средней его величины. В определении (9) приведено понятие, удовлетворяющее этим условиям.

Рассмотрим теперь свойства бесконечношаговых алгоритмов вычисления / (ж). Для БША имеет место следующее определение:

Определение 6. Пусть в бесконечношаговом алгоритме выполнено N первых базовых операций и для / (ж) получено приближение значения вычисляв мой функции /П% (ж). БША вычисления функции / называется нормальным,, если для всех N алгоритм вычисления значения /П^(ж) будет нормальным. □

Определение 7. Бесконечношаговый алгоритм называется сходящимся, если /(ж) = Пт (ж). □

N

Теорема 1 может быть уточнена на случай сходящегося бесконечношагового алгоритма следующим образом:

Теорема 2. Пусть для вычисления значения функции fN (ж) БША является нормальным и выполнены условия Теоремы 1. Тогда существуют векторы С € Кк, 7 € Кк такие, что V е> 0: ||71| ^ е, имеет место представление

&(х) = / (х) + Сб1 + 7. (16)

Доказательство

Из сходимости бесконечношагового алгоритма следует, что V е > 0 3 N : ||fN(ж) - /(ж) Ц ^ е и fN (ж) - £ (ж) = 7, Ц7Ц ^ е. Из Теоремы 1 для данного N имеем равенство

/ттг (х) = fN(х) + С81, ГДе С € Дк. Тогда ДЛЯ значения ^ (ж) получим

/й (х) = / (х) - / (х) + fN (х) + ^1 = / (х) + ^1 +1.

□

4. О гарантированной точности решений задач ВМ

В этом разделе для случая, когда значение параметра погрешности ||С|| ограничено, получена оценка длины мантиссы МЧ, гарантирующая достижение требуемой точности решения задачи ВМ.

Определение 8. Будем говорить, что метод (алгоритм) вычисления значения функции / € йк называется корректным (КМ), если для любого е > 0 найдется такой размер мантиссы т, что

А = II/(ж) - /т(ж)|| < £, ИЛИ А < £, при ||/т(ж)|| = 0. (17)

11/т^Ц

Величину £ в (17) будем называть требуемой точностью. □

Обратимся к вопросу определения достаточной (гарантированной) точности, которую должно обеспечить ВУ — вычислительное устройство, т.е. ЭВМ, имеющая необходимое программное обеспечение.

Пусть ВУ имеет переменную длину мантиссы т, т.е. пользователь может выбрать размер мантиссы, необходимый для проведения вычислений. Пусть £ : О С Яп ^ а /т(ж)

— значение функции /, вычисленное в точке хт с помощью данного ВУ.

Определение 9. Будем говорить, что значение функции /т(х) имеет погрешность порядка, а, 0 < а ^ 1, относительно погрешности мантиссы 5^, если существуют константы С и Со такие, что V х € С:

А = II/(х) — fm(%)\\ ^ С (5ц)а , для абсолютной погрешности;

II/(х) — fm(x)\\ А . ла || . ...

---II /> / \||- = \\ Г ( \|| ^ U (др) прИ IIjm(%)II = 0, ДЛЯ ОТНОСИТбЛЪНОИ п0Грешн0СТИ,

\\}т(х) \1 \\jm\X) ||

(18)

где 5ц = Ь-т. □

Приведем без доказательства простое, но полезное утверждение, в котором гарантируется погрешность значения fm (ж) иорядка а.

Лемма 4. Пусть для функции f (ж) в некоторой точке х £ Q С G, где Q — компакт, выполнены условия Теоремы 1, а в равенстве (15) а удовлетворяет уеловию 0 < а ^ 1. Тогда значение функции fm (ж) имеет погрешность порядка ас, 0 < ас ^ 1, V х £ Q. □

Докажем теперь теорему о длине мантиссы МЧ, гарантирующей требуемую точность решения задачи ВМ.

Теорема 3. Пусть погрешности А значения функции /т(ж) или имеют порядок

а. Тогда для любого е > 0 данный метод вычисления функции будет корректным при т ^ |_1 — ^ logb -jj ил и т ^ 1 — ^ logb —■ , где LAI — целая часть ч исла А. Доказательство.

А

грешности аналогично.

Зададим число е > 0. Так как для достаточно больших т число (5^)а = Ь-ат может быть достаточно малым, то выберем т такое, что (5^)а = b-am ^ -j. Подставляя в (18) значение — вместо (5^)а , получим, что А = ||/(ж) — /т(ж)|| ^ е, т.е. выполнено (17) и метод вычисления будет корректным при выбранном т. Теперь из неравенства b-am ^ — найдем значение т: —am ^ logb —, откуда т ^ logb -j. Так как т — целое число, то очевидно: т ^ |_1 — ^ logj, . Значение мантиссы т для относительной погрешности определяется

по формуле т ^ 1 — ^ logb —■ . □

Замечание.

1. В важном частном случае, когда а = 1, оценки для т, из Теоремы 3 имеют вид

т ^ [1 — log6 и т ^ 1 — log6 f-

2. Константы С и Сс связаны следующими условиями:

< < ^с (6ц)а , (19)

А = \| f(x) — fm(x)\\ < С (5ц)С Ifm(x)\\ |fm(x)\\ < ||fm(x)\\ < ^ ^

ГДе °с = цЗм ™ °с = 7’ Т'6• \\f—x)\\ < f = °с> I = Ц/m^ll, 2 £ П CG

компакт. В случае, когда fm(x) = 0, относительная погрешность не рассматривается.

Ввиду особой важности Теоремы 3 для прикладных исследований переформулируем её в следующем виде:

Теорема 4. (правило гарантированной точности решений задач ВМ). Пусть погрешность (абсолютная или относительная) вычисленного значения функции f имеет порядок а, 0 < а < 1. Тогда для любой требуемой точности е > 0 существует такой размер мантиссы т, при котором достигается заданная точность е решения задачи. □

На практике Теоремы 3, 4 применимы до максимального значения мантиссы rnmax, которую обеспечивает данная библиотека программ. В частности, для библиотеки GNU GMP rnmax = 646456993 десятичных знаков.

5. Вывод

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

Литература

1. IEEE 754-2008: 754-2008 IEEE Standard for Floating-Point Arithmetic. - ISBN: 978-07381-5753-5.

2. Годунов С. К., Антонов А. Г., Кирилюк О. П., Костин В. И. Гарантированная точность решения систем линейных уравнений в евклидовых пространствах. — Новосибирск: Наука. Сиб. отд-ние, 1988. - 456 с. ISBN 5-02-028593-5.

3. Higham N. J. Accuracy and stability of numerical algorithms. — Philadelphia: Society for Industrial and Applied Mathematics, 1996.

4. Torbjorn Granlund \et al.\ GNU Multiple Precision Arithmetic Library 4.1.2 / http://swox.com/gmp/, 2002.

5. Fousse, Laurent and Hanrot, Guillaume and Lefevre, Vincent and Pelissier, Patrick and Zimmermann, Paul MPFR: A multiple-precision binary floating-point library with correct rounding. — ACM Trans. Math. Softw., 2007.

6. Марков А. А., Нагорный H. М. Теория алгорифмов. — М.: Наука, 1984.

7. Воеводин В. В. Вычислительные основы линейной алгебры. — М.: Наука, 1977. — 304 с.

8. Ива,нов В. В. Методы вычислений на ЭВМ: Справочное пособие. — Киев: Наукова думка, 1986. — 584 с.

9. IEEE standard for radix-independent floating-point arithmetic: ANSI/IEEE Std 854-1987, 1987. http://grouper.ieee.org/groups/754/

10. ISO/IEC 9899:1999 Standard for the С programming language (C99) — 1999.

11. Wilkinson J. Н. Rounding Errors in algebraic processes. — Englewood Cliffs, N.J.: Prentice-Hall, 1963. ISBN 0-486-67999-3.

12. Henrici P. Elements of Numerical Analysis. — John Wiley k, Sons Inc., New York, 1964.

13. Clenshaw C. W. and Olver F. W.J. Beyond floating point // J. Assoc. Comput. Mach. — 1984 - V. 31 - P. 319-328.

14. Langlois P. A Revised Presentation of the CENA Method. — ARENAIRE — Inria Grenoble Rheme-Alpes / LIP Laboratoire de l’lnformatique du Parallelisme

15. Шокин Ю.И. Интервальный анализ. — Новосибирск: Наука. Сиб. отд-ние, 1981.

16. Шарый С. П. Конечномерный интервальный анализ. — М., 2007.

17. Francoise Chaitin-Chatelin and Valdrie Fraysse. Lectures on Finite Precision Computations.

— Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 1996.

18. Алефельд Г., Херцбергер Ю. Введение в интервальные вычисления. — М.: Мир, 1987. — 356 с.

19. Воеводин В. В. Ошибки округления и устойчивость в прямых методах линейной алгебры. — М.: Изд-во МГУ, 1969. — 140 с.

Поступила в редакцию 05.06.2012.