Формальные подходы к тестированию математических функций Текст научной статьи по специальности «Математика»

Формальные подходы к тестированию математических функций

В. В. Кулямин

Институт системного программирования РАН (ИСП РАН),

Б. Коммунистическая, 25, Москва, Россия E-mail: kuliamin@ispras. ru

Аннотация

Данная статья рассматривает вопросы проверки корректности вычисления математических функций на числах с плавающей точкой, формат которых определяется стандартом IEEE 754. В ней описывается метод разработки тестов для реализаций таких функций, основанный на формальных спецификациях их поведения. Предлагаемый метод основан на технологии разработки тестов UniTESK и двух дополнительных методиках: методике формирования строгих требований к реализации математической функции и методике построения набора тестовых данных для ее тестирования. Описанные методики опираются на специфические свойства представления чисел с плавающей точкой и особенности поведения самой тестируемой функции.

1. Введение

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

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

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

Правильность их проверяется как соответствие этих результатов ожиданиям пользователей в соответствующей ситуации.

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

• Моделирование космогонических процессов: развития звезд, планетных систем и галактик. Моделирование космических катастроф, последствия которых могут угрожать существованию человечества.

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

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

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

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

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

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

В данной работе изучаются проблемы проверки корректности реализаций математических функций, работающих с числами с плавающей точкой, анализируются имеющиеся достижения в этой области и предлагается метод разработки тестов для таких функций на основе формальных спецификаций их поведения. Полученный метод базируется на технологии 1МТЕ8К [1-3], дополняя ее методиками формирования точных требований к реализациям математических функций и выбора тестовых данных для их тестирования.

Сперва рассмотрим проблемы, с которыми сталкиваются исследования в этой области.

2. Проблемы корректного вычисления функций

Основные проблемы корректного вычисления математических функций связаны с дискретностью представления действительных чисел в компьютерах. Чтобы иметь возможность эффективно выполнять операции с действительными числами, они представлены в виде так называемых чисел с плавающей точкой, формат и правила действий над которыми определены в стандартах ШЕЕ 754 [4] (он же — 1ЕС 60559 [5]) и ШЕЕ 854 [6].

ШЕЕ 754 определяет представление двоичных чисел с плавающей точкой, ШЕЕ 854 обобщает его, определяя и десятичные числа с плавающей точкой. Однако, поскольку в большинстве случае на практике используется двоичное представление чисел, мы будем рассматривать только его.

2.1. Числа с плавающей точкой

Двоичное число с плавающей точкой имеет следующую структуру [4,6,7].

• Число представлено в виде набора из п бит, из которых первый бит является знаковым битом числа, следующие к бит отданы под представление его экспоненты, а оставшиеся (п-к-1) бит представляют его мантиссу.

0 1 к к+1 п-1

1 1 1 1 1 1 1 1 1 1 1 1 1 1 II II II II II II II II II II II II

знак экспонента мантисса

Рисунок 1. Битовое представление чисел с плавающей точкой.

• Знаковый бит S, экспонента Е и мантисса М числа х определяют его значение по следующим правилам.

х = (-\)s-2e-m, где

о S — знаковый бит, 0 для положительных чисел, и 1 для отрицательных;

о еслиЯ>0,тое = Я-2(ы)+1; иначе, если Е = 0, е = -2(ы)+2;

число (2^"^-1) называется смещением экспоненты (bias)',

о если 0 < Е < 2k-1, то т имеет двоичное представление 1М, т.е. целая часть т равна 1, а последовательность цифр дробной части совпадает с последовательностью битМ;

если же Е = 0, то т имеет двоичное представление ОМ,

такие числа (с нулевой экспонентой) называются денормализоеанными.

• Максимальное возможное значение экспоненты Е = 2*-1 зарезервировано для представления положительной +со и отрицательной -со бесконечностей и специального значения NaN (not-a-number), которое возникает, если результат выполняемых действий нельзя корректно представить ни обычным числом, ни бесконечностью, например, 0/0 или (-со) + (+со).

+со имеет нулевой знаковый бит, максимальную экспоненту и нулевую мантиссу; -со отличается только единичным знаковым битом.

Любое число, имеющее максимальную экспоненту и ненулевую мантиссу, считается представлением NaN.

Все такие числа называются исключительными.

• Стандарты IEEE 754 и IEEE 854 определяют несколько возможных типов чисел с плавающей точкой, из которых чаще всего используются числа однократной точности (single precision), числа двойной точности (double precision) и числа расширенной двойной точности (double-extended precision).

Для чисел однократной точности п = 32 и к = 8. Соответственно, для мантиссы используется 23 бита и смещение экспоненты равно 127. Для чисел двойной точности п = 64ик=11. Для мантиссы используется 52 бита и смещение экспоненты равно 1023.

Для чисел расширенной двойной точности определенные значения кип

не фиксируются в стандартах, вводятся лишь ограничения 128 > п > 80 и к > 15. В процессорах Ые1 32-битной архитектуры используются значения п = 80 и к = 15. При этом для мантиссы используется 64 бита и смещение экспоненты равно 16383.

числа, имеющие одну экспоненту и разные мантиссы

интервалы между нормализованными числами с соседними экспонентами различаются в 2 раза

|111111111И| 1111 I I I |-1-1—I—Ч—---------1-------1_____----------------1—

денормализованные числа ^"

числа, имеющие нулевую мантиссу и различные экспоненты

Рисунок 2. Примерная картина распределения чисел с плавающей точкой.

Примеры чисел с плавающей точкой.

• Числа однократной точности

о 0 = (-1)°-2°"126-0.02 имеет представление

0 00000000 00000000000000000000000

0 127 127

о 1 = (— 1) -2 " -1.02 имеет представление

0 01111111 00000000000000000000000 о -17ю = (-1 )1-2131-127-1.0 00 12 имеет представление

1 10000011 00010000000000000000000 о 0.75ю = (-1)°-2126-127-1.12 имеет представление 0 01111110 10000000000000000000000 о самое маленькое положительное число, представимое в таком виде (денормализованное)

2-149 = (_1)».2°-126.0.000000000000000000000012 О 00000000 00000000000000000000001 о самое маленькое положительное число, представимое в нормальном виде 2"126 = (-1)°-21-127-1,02 0 00000001 00000000000000000000000 о самое большое положительное число, представимое в таком виде 2104-(224-1) = (-1 )°-2254-127-1.1111 1111 1111 1111 1111 1112 0 11111110 11111111111111111111111 о +со имеет представление

0 11111111 00000000000000000000000 о следующее за 1 число, представимое в таком виде 1+2"23 = (-1)°-2127_127-1.000000000000000000000012 О 01111111 00000000000000000000001

о предшествующее 1 число, представимое в таком виде

1-2"24 = 2"24-(224-1) = (-1)°-2126-127-1.111111111111111111111112

0 01111110 11111111111111111111111

• Числа двойной точности

о 0 = (-1)°-2°"1022 0.02 имеет представление

0 00000000000 0000000000000000000000000000000000000000000000000000 о 1 = (-і)°-21023“1023- 1 .Ог имеет представление

0 01111111111 0000000000000000000000000000000000000000000000000000 о -17ю = (-1)1-21027'1023-1.0 00 1 2 имеет представление

1 10000000011 0001000000000000000000000000000000000000000000000000 о 0.75ю = (-і)°-21022“1023- 1.12 имеет представление

0 01111111110 1000000000000000000000000000000000000000000000000000 о самое маленькое положительное число, представимое в таком виде (денормализованное)

2" 1074 _ j-jO 2»-1022 ^-52

0 00000000000 0000000000000000000000000000000000000000000000000001 о самое маленькое положительное число, представимое в нормальном виде 2'1022 = (-1)°-2Ы023-1.02 0 00000000001 0000000000000000000000000000000000000000000000000000 о самое большое положительное число, представимое в таком виде 2971-(253-1) = (-1)°-22046"1023-(253-1)-2"52 0 11111111110 1111111111111111111111111111111111111111111111111111 о +со имеет представление

0 11111111111 0000000000000000000000000000000000000000000000000000 о следующее за 1 число, представимое в таком виде 1+2"52 = (-1)°-2127_127-( 1 +2“52)

0 01111111111 0000000000000000000000000000000000000000000000000001 о предшествующее 1 число, представимое в таком виде 1-2"53 = 2"53-(253-1) = (-1)°-21022"1023-(253-1 )-2"52 0 01111111110 1111111111111111111111111111111111111111111111111111

В дальнейшем будем называть действительные числа, точно представимые в виде чисел с плавающей точкой, просто представимыми. Ясно, что не каждое действительное число представимо. Например, к — иррационально, а все представимые числа рациональны, точнее даже двоично-рациональны. 1/3 тоже не

тт ОІ0000000

представимо, поскольку не является двоично-рациональным числом. Числа 2 и

--Ю000000

2 не представимы, хотя и двоично-рациональны, потому что первое слишком

велико, а второе слишком мало по абсолютной величине, даже для представления с четырехкратной точностью. Числа, находящиеся слишком близко к представимым, тоже не представимы, например, для представления 1+2"1000 нужна мантисса из 1000 нулей и одной единицы.

Заметим, что существует число с плавающей точкой -0, отличающееся от 0. Стандарт IEEE 754 требует, однако, считать их равными. Кроме того, ни одна из операций над числами с плавающей точкой, описанных в этом стандарте, не должна давать

в результате -0, за исключением квадратного корня из -0. При сложении, вычитании, умножении, делении, вычислении остатка от деления и преобразованиях типов всегда в случае нулевого результата возвращается 0.

2.2 Требования стандартов к реализациям математических функций

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

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

Однако, ни стандарты ШЕЕ 754 и ШЕЕ 854, ни стандарт языка С 180/1ЕС 9899 [8], ни стандарт переносимого интерфейса операционной системы ШЕЕ 1003.1 [9] (известный как Р081Х), описывающие библиотеку математических функций языка С не фиксируют такого требования для большинства функций.

Стандарты ШЕЕ 754 и ШЕЕ 854 описывают работу только сложения, умножения, вычитания и деления чисел с плавающей точкой, а также вычисления остатка от деления, извлечения квадратного корня и преобразований между типами с плавающей точкой и между ними и целочисленными типами. Соответственно, только для этих операций требуется возвращать результат, полученный из точного приведением к ближайшему представимому числу согласно действующему режиму округления. Возможны 4 режима округления: просто к ближайшему, к 0, к +со и к -со. Кроме того, эти стандарты требуют аккуратного выставления флагов переполнения, слишком маленького результата или неточного результата при работе этих операций.

Стандарт С ссылается на требования ШЕЕ 754, добавляя только ограничения на значения результатов ряда функций для некоторых значений параметров (например, ехр(0) = соя(0) = 1, а яи(0) = ?£(0) = 0). Стандарт Р081Х, в свою очередь, ссылается на требования стандарта языка С, добавляя описание поведения реализаций математических функций в случае возникновения переполнения или слишком маленьких результатов и для тех значений параметров, где соответствующая функция не определена.

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

В последние 5-10 лет появились предложения стандартизовать необходимые для аккуратного моделирования требования к реализациям математических функций [10,11]. Многие инициаторы этой деятельности работают в проекте

Arenaire [12], совместно проводимом INRIA, CNRS и Высшей Нормальной школой Лиона, Франция. В результате активно разрабатывается набор стандартов ICO/IEC 10967 [13-15], формулирующий естественные ограничения на работу реализаций большинства математических функций. Эти ограничения касаются нескольких аспектов.

• Возможные погрешности вычисления функций выражены в терминах единиц последнего разряда (unit in the last place, ulp) [7]. Наилучшее приближение дало бы точность в 0.5 ulp, т.е. вычисленный результат отличался бы от точного значения функции не более чем на половину единицы последнего разряда мантиссы результата.

Однако природа математических функций такова, что такая точность не является практически обоснованной во многих случаях, хотя ее достижение потребовало бы значительных усилий от разработчиков библиотек. Это связано с тем, что число с плавающей точкой является приближенным представлением любого действительного числа, к которому оно является ближайшим. Таким образом, уже в значениях аргументов функции может иметься погрешность, которую ее вычисление не в силах исправить. Для широко используемых математических функций были проведены оценки возрастания погрешности при их вычислении, и на основании этих оценок были сформулированы более практичные требования к точности вычислений, ограничивающие погрешность результата величиной от 0.5 ulp до 2 ulp, в зависимости от функции [14].

• Стандарты серии ISO/IEC 10967 требует от реализации математической функции сохранения знака ее точного значения для данного значения параметра. Также требуется, чтобы во всех интервалах монотонности функции ее реализации были монотонны таким же образом, т.е. там, где сама функция убывает, ее численная реализация должна убывать, а там, где функция возрастает, — возрастать.

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

• ISO/IEC 10967 требует соблюдения специфических требований при вычислении функций в окрестностях точек, где они имеют известные представимые значения.

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

Это требование связано с тем, что плотность чисел с плавающей точкой в окрестности 0 гораздо больше, чем их плотность в окрестности 1 — между 1 и ближайшим к нему числом с плавающей точкой умещается много чисел, близких кО. Для двойной точности ближайшее к 1 число равно 1-2"53, а ближайшее к 0 — это 2"1074.

Рисунок 3. График функции, имеющей ненулевое значение в 0.

2.3. Дилемма составителя таблиц

Сформулированные требования к точности вычисления математических функций приводят к так называемой дилемме составителя таблиц (Table Maker's Dilemma) [16,17].

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

Для иллюстрации дилеммы составителя таблиц приведем следующий пример. Пусть вычисляется функция sin для двоичных чисел с плавающей точкой, имеющих 6 битов мантиссы. Синус числа II.IOIO2 = 3.625ю (выделены биты мантиссы) равен

0.011101101111110...2 = 0.063225984913... ю (снова выделены биты мантиссы). Приближенное вычисление 6-ти бит мантиссы результата может дать как 0.011101Ь, так и О.ОППОО2, поскольку точное значение очень близко к их среднему арифметическому. Только получив точный 14-й бит, мы сможем уверенно выбрать первое из них в качестве значения, ближайшего к точному результату.

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

соседние числа с плавающей точкой и их среднее арифметическое

Рисунок 4. Дилемма составителя таблиц.

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

Приведем более реалистичный пример для чисел с двойной точностью, имеющих 52 бита в мантиссе. Вычисляя значение натурального логарифма для

1.01100001001110010101010111011100100000000010111110002-2-35 получим

-10111.1111000000101111100110111010111101100000001101010 I 60 0011... 2

Обозначение I60 означает, что единица повторяется 60 раз. Таким образом, для получения корректно округленного значения нужно вычислять логарифм в этой точке с относительной погрешностью, не превосходящей 2"113.

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

Например, натуральный логарифм для числа 110101100.010100001011010000001001110010001010111011102 равен

110.00001111010100101111001101111010111011001111110011 161 0101...2.

Для практически всех часто используемых функций их значения в «обычных» (т.е., не равных 0, 1 или 2) двоично-рациональных числах не являются рациональными, и поэтому не могут ни быть представимыми, ни лежать в точности посередине между двумя представимыми числами. Из этого в силу конечности множества представимых чисел следует, что для каждой функции есть такое число е > 0, что вычисляя значения этой функции с погрешностью, не превосходящей е, можно всегда точно определить корректное округление, являющееся представимым числом. Однако вычислять функцию с такой точностью для всех значений аргумента может оказаться слишком неэффективно.

Например, для натурального логарифма на числах двойной точности при произвольном режиме округления такое е можно взять равным 2"118 [17]. Однако реально такая точность нужна только для единственного значения аргумента, во всех остальных случаях можно использовать меньшую. Для подавляющего же большинства представимых чисел двойной точности корректное округление их логарифма можно

„ __54

получить, вычисляя его с погрешностью, не превосходящей 2 .

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

3. Обзор работ по проверке корректности реализаций математических функций

Методам вычисления математических функций посвящено огромное количество работ. Одним из классических трудов на эту тему является сборник статей под редакцией Abramowitz и Stegun [18], хотя он был выпущен уже довольно давно и частично устарел. Более современное изложение методов вычисления элементарных функций (являющихся только подмножеством рассмотренных в [18]) можно найти в книге Muller [19].

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

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

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

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

• Работы по формальной верификации определенных алгоритмов.

• Работы по тестированию реализаций математических функций.

3.1. Работы по формальной верификации определенных алгоритмов

Формальная верификация корректности алгоритмов вычислений математических функций чаще всего проводится при проектировании блоков вычислений с плавающей точкой универсальных процессоров или при разработке специализированных вычислительных процессоров [28,32,34-52,54].

Основой для такой верификации всегда является точное знание реализованного алгоритма и формализация основных требований стандарта IEEE 754, необходимая для их строгой проверки. Формализация этих требований проводилась несколько раз в различных формализмах: в языках формальных спецификаций Z [20], Barrett в 1989 [21], и VDM [22], Wichmann в том же году [23], в формализмах инструментов автоматизации доказательств Nurpl [24], PVS [25], HOL [26] и ACL2 [27] —

O’Leary с соавторами для Nurpl в 1994 [28], Miner для PVS в 1995 [29,30], Саггепо для HOL в 1995 [30,31], Moore с соавторами для ACL2 в 1996 [32], Harrison для HOL в 1996 [33].

В работах Verkest с соавторами [34] и Cornea-Hagesan [35,36] все доказательства выполнялись вручную. Но в большей части исследований такого рода, ввиду значительной сложности верификации практически важных систем, применяются только формализации, которые могут использоваться инструментами для автоматизированного доказательства теорем.

Статьи [32,34,37-41] имеют дело только с алгоритмами умножения и/или деления чисел с плавающей точкой и представляют примеры верификации корректности таких алгоритмов по отношению к требованиям IEEE 754.

В работах [35,36,42-48] верифицируются также алгоритмы вычисления квадратного корня, иногда еще и вычисление остатка от деления для чисел с плавающей точкой и преобразования между различными такими типами и между ними и целыми числами.

В связи с тем, что в стандарте IEEE 754 определены только арифметические действия, вычисление остатка от деления и квадратного корня, проверке правильности вычисления других функций посвящено значительно меньше работ. В [49,50] описываются практические примеры формальной верификации алгоритмов, вычисляющих экспоненциальную функцию. В статье [51] верифицировалось вычисление синуса и косинуса, а в [52] — натурального логарифма. Все четыре работы были выполнены Harrison и его коллегами. В работе [53] нескольких исследователей из проекта Arenaire [12], представлен пример верификации вычисления натурального логарифма с помощью специализированного инструмента Gappa, позволяющего автоматизировать доказательства, касающиеся свойств элементарных функций. Других примеров верификации математических функции в доступной литературе найти не удалось.

В перечисленных работах верифицировались отдельные элементы и блоки процессоров компаний Intel (Pentium II, Pentium III, Pentium 4 и Itanium [35,36,41-43,46,47,49-52]) AMD (K5, Athlon [32,44,45]), IBM (Power4 [48]). При этом использовались следующие инструменты автоматизации доказательств: HOL ([43,49-52]) и близкий к нему Nurple ([29]), PVS ([38,40,54]), ACL2 ([32,44,45,48]) и его предыдущая версия Nqthm ([34]). Помимо этого применялись различные комбинации автоматизированного доказательства теорем с проверкой моделей (model checking) и символической проверкой эквивалентности моделей ([37,39,41,42,46,47]).

По-видимому, до сих пор не было удачных попыток полностью проверить соответствие блока вычислений с плавающей точкой стандарту IEEE 754. Во всех найденных работах либо проверяются не все операции, определяемые стандартом, либо не проверяется их работа для специальных чисел с плавающей точкой — денормализованных, -0, бесконечных и NaN. Работа [54], похоже, единственная, в которой систематически рассматривается поведение части описанных в IEEE 754 операций на таких исключительных значениях.

Можно также отметить, что в процессоре, верифицированном в проекте, описанном в [42] впоследствии была найдена ошибка в операции преобразования числа с плавающей точкой в целое (так называемый FIST bug [55]). Это показывает, что формальная верификация сложной системы, будучи тоже достаточно сложной деятельностью, из которой нельзя исключить участие людей, сама по себе подвержена ошибкам.

3.2. Работы по тестированию реализаций математических функций

В Интернет можно найти огромное количество различных программ для тестирования функций, работающих с числами с плавающей точкой, см. например, [56]. К сожалению, подавляющее большинство таких тестов крайне несистематично и проверяет какой-то один аспект вычислений, реже — два-три таких аспекта.

Как указывается в [56] (см. также иллюстрации из [57,58]), несмотря на то, что стандартизация вычислений с плавающей точкой началась около 20 лет назад, до сих пор многие поставщики библиотек и аппаратного обеспечения не придерживаются имеющихся стандартов достаточно строго, поэтому тесты на правильность поведения реализаций математических функций по-прежнему необходимы.

Среди наиболее систематичных работ по тестированию вычислений с плавающей точкой можно назвать следующие.

• Работы по тестированию на соответствие стандарту IEEE 754.

о В работе [59] описывается самый первый из известных систематических тестовых наборов для проверки корректности реализации операций над числами с плавающей точкой. Он появился еще до введения в действие стандарта IEEE 754, сделан в виде набора программ на Fortran и предназначен для тестирования только сложения, вычитания, умножения и деления.

о Специально для проверки на соответствие IEEE 754 был разработан тестовый набор, который описан в статьях [60,61] и может быть получен с сайта [62]. В этом наборе проверяются все требования стандарта к арифметическим операциям, вычислению квадратных корней и взятию остатков, а также преобразования между типами чисел с плавающей точкой и целыми.

о Программа PARANOIA [63,64] была создана одним из авторов стандарта IEEE 754 Кэханом (W. Kahan) и остается довольно популярным средством проверки на соответствие ему, хотя такая проверка менее тщательна, чем с помощью тестового набора, описанного выше. Она также проверяет только базовые арифметические операции.

о Другой подход к построению тестов для операций IEEE 754 используется в среде FPgen [65,66]. Здесь, помимо специальных значений, в качестве тестовых данных используются числа с плавающей точкой,

удовлетворяющие некоторым шаблонам — например, в которых нулевые и единичные биты мантиссы чередуются, или в которых мантисса содержит ровно 7 единиц.

• Работы по тестированию широкого набора математических функций.

о Тестовый набор UCBTEST [67] предназначен для тестирования базовых арифметических действий и достаточно широкого набора математических функций. Он оформлен как набор программ на разных языках, включая Fortran и С, и наборов предопределенных входных данных для разных функций. В каждом тесте проверяется, что для заданных значений параметров данная функция возвращает число с плавающей точкой, ближайшему к точному значению функции.

о Тестовый набор ELEFUNT [68], основанный на книге [69], также содержит тесты для многих математических функций в виде программ на С и Java и текстовых файлов с тестовыми данными и ожидаемыми результатами. Этот и предыдущий тестовый набор построены на основе проверки значений, возвращаемых реализациями функций для некоторых наборов аргументов. Методика выбора этих наборов аргументов, скорее всего, использовала несколько разных соображений.

■ Выделялись особые значения чисел с плавающей точкой: 0, -0, NaN +со, -со, минимальное положительное, максимальное положительное, числа, имеющие ровно один бит в мантиссе, и пр. (см. ниже).

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

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

из соображений, связанных со структурой известных алгоритмов для вычисления элементарных функций. Например, в ряде алгоритмов вычисления логарифма сначала значение аргумента при помощи умножения или деления на 2 приводится к интервалу (0.5, 1].

Соответственно, выбираются границы этого интервала и нескольких соседних с ним.

о Аналогичные подходы — использование ряда специальных значений, границ интервалов, определяемых часто используемыми алгоритмами вычисления данной функции, и случайных значений — применялись для построения более полных тестовых наборов, например, набора Беркли [70], а также в статье [71].

Отдельно стоит отметить работы в рамках проекта Arenaire [16,17,72], посвященные дилемме составителя таблиц и поиску чисел, для которых корректное вычисление функций с заданной точностью наиболее трудоемко. Эти числа можно использовать в

качестве «неудобных» тестовых значений для практически любой реализации соответствующей функции.

В рамках того же проекта был разработан инструмент МРСЬеск [73] для тестирования корректности реализации элементарных функций с точки зрения сохранения монотонности, симметрий, ограничений на область значений и корректности округления.

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

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

4. Предлагаемый подход

Подход к тестированию реализаций математических функций, предлагаемый в данном исследовании, основан на технологии 1МТЕ8К, использующей формальные спецификации требований к программному обеспечению для автоматизированного построения тестов на соответствие им.

Основные элементы технологии ишТЕвК [1-3], применительно к реализациям математических функций, таковы.

• Требования к поведению тестируемой системы представляются в виде формальных спецификаций, которые состоят из следующих частей:

о предусловия, описывающие области определения функций;

о постусловия, описывающие условия корректности возвращаемых функциями результатов;

о инварианты типов данных, описывающих условия целостности данных.

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

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

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

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

Такие ситуации также могут быть описаны в спецификациях с помощью специальных конструкций.

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

Возможности технологии 1МТЕ8К по автоматическому построению тестовой последовательности для тестирования математических функций не требуются. Единственное возможное исключение — наличие гипотез о зависимости работы реализаций математических функций от каких-то элементов внутреннего состояния тестируемой системы. В этом случае дополнительно должна быть построена обобщенная модель состояния, каким-то образом учитывающая те элементы, которые, как предполагается, могут влиять на работу функций.

Чтобы адекватно применять технологию ИшТЕвК для построения тестов для математических функций, необходимо ответить на два вопроса: какие требования должны предъявляться к поведению их реализации, т.е. что именно должно быть написано в спецификациях; и на каких значениях параметров функций должно проводиться тестирование.

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

4.1. Метод определения требований к математическим функциям

Данный метод определения требований к поведению реализаций математических функций заимствует часть идей из стандарта 1вО 10967 [13-15] и работ [10,11], посвященных вопросу разработки серии стандартов с повышенными требованиями к корректности вычисления математических функций и удовлетворяющих им библиотек. Некоторые элементы предлагаемого метода являются новыми и не встречаются в доступной литературе.

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

• Область определения функции и особые точки функции.

о Для всех значений аргументов, где математическая функция определена, ее реализация должна возвращать некоторый результат, который может быть равен +со или -со, если значение самой функции находится за пределами интервала чисел с плавающей точкой, но не должен быть №№

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

о Если функция имеет особенность в точке 0, не имеет там предела, равного +со или -со, нужно рассматривать односторонние пределы функции. Значение реализации функции в 0 нужно считать равным ее пределу при х —» +0, если он существует, а значение в -0 — пределу при х—»-0, если он есть. Примером такой функции служит котангенс.

о В остальных случаях должен возвращаться результат №№

Особо нужно рассматривать такие значения аргументов, по поводу которых нет однозначного мнения о принадлежности их к области определения функции или о возможном продолжении ее в эту точку по непрерывности. Примером служит значение 0°, которое иногда интерпретируется как 1, а иногда как №№

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

о Для функций, стремящихся к бесконечности при X —^ +СО ИЛИ X —» -СО, должны быть точно определены пределы представимости их значений. За этим пределами реализация также должна возвращать +со или -со в соответствии со знаком точного значения.

Например, для значений х > /и(2104-(224-1)) = 88,722839052... ехр(х)

не попадает в диапазон чисел однократной точности, поэтому реализация экспоненты для таких чисел должна возвращать +со.

Специальные значения, значения в 0, касательные и асимптоты.

о Нужно наиболее естественным образом определить значения функции для особых значений аргумента: -0, +со, -со. Обычно достаточно определять их как пределы, если те существуют, иначе как №№

о Значение функции для значения аргумента ИаИ должно быть равно №№

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

Например, ехр(0) = соя(0) = ск(0) = 1, яи(0) = ?ап(0) = агс$т(0) = 0 и т.п.

о Кроме этого, если в такой точке производная функции равна 0, то для любого аргумента из некоторой ее окрестности реализация должна возвращать то же самое значение.

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

Например, для экспоненты е = 1, при этом (е*-1) < 2'24 при X > 1 и(х-1) < /и(1+2"24) = 5.9604642999... -10"8, поэтому для всех таких х, которых довольно много, реализация экспоненты с однократной точностью должна возвращать

1.

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

Например, для значений х < /и(2'150) = -103.97207708399... значение е становится ближе к 0, чем к какому либо еще числу, представимому с однократной точностью. Поэтому реализация экспоненты для таких аргументов должна возвращать 0.

о Казалось бы, естественно предъявить аналогичные требования к функциям, имеющим негоризонтальные асимптоты или асимптотически близких к другим функциям. Например, с/г(х) = (ех+ех)!2 ~ е для достаточно больших значений аргумента, или яи(х) ~ х при х ~ 0.

Однако во многих случаях такое требование не может быть сформулировано достаточно аккуратно с учетом различных режимов округления, зафиксированных в ШЕЕ 754. Дело в том, что даже очень маленькая разность между двумя асимптотически близкими выражениями может дать отличие в

значимых битах мантиссы. Причина этого явления аналогична причине, порождающей дилемму составителя таблиц — слишком близкое расположение некоторых значений функции к представимым числам. Например, если рассматривать асимптотику е ~ 1+.т при х ~ 0 для чисел двойной точности, то при |.т| < 2"28 разница между е и 1+.т уже меньше 0.5 ulp, однако встречаются гораздо более близкие к 0 числа, представимые с двойной точностью, для которых мантиссы е и 1+.Т отличаются при выборе режимов округления К -СО ИЛИ К +СО.

Скажем, для х =

-1.100000000000000000000000000000000000000000000000000І2-2'52

ехр(х) =

1.111111111111111111111111111111111111111111111111110 1 053 Ю0...2-2Л

а 1 + х =

1.1111111111111111111111111111111111111111111111111100 1512-2_1.

Цветом здесь выделены биты, не помещающиеся в мантиссу числа двойной точности.

Похоже, что четкие требования можно предъявлять лишь по поводу соблюдения асимптотик вида fix) ~ х или fix) ~ -х (например, sin(x) ~ х и tg(x) ~ х при х ~ 0), поскольку при этом риск столкнуться с подобной проблемой невелик. Однако при этом нужно очень аккуратно вычислять границы действия подобных ограничений.

• Область значений функции.

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

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

— это 13176795/8388608 = 1.57079637050628662109375 > к/2. При этом tg(arctg( 230)) = -2.2877...-107, что противоречит основному свойству

обратных функций.

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

• Монотонность и сохранение знака.

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

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

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

о Знак значения реализации функции для некоторого аргумента должен совпадать со знаком значения самой функции. Если значение функции равно О для какого-то представимого числа, значение ее реализации для этого числа также должно быть равно 0.

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

Симметрии и периодичность.

о Если математическая функция является четной или нечетной, этим же свойством должна обладать ее реализация.

о В тех случаях, когда функция имеет симметрии относительно других представимых значений аргумента, например, выполняется правило fi 1 -х) = -fix), выполнение аналогичного свойства для реализации не всегда возможно, поскольку число (1-х) может быть представимым при непредставимом х. Следует особо рассматривать такие случаи и накладывать ограничения, быть может, касающиеся только представимых значений аргумента с обеих сторон такого равенства.

о Если функция симметрична относительно непредставимого значения, как, например, синус — sin(n-x) = sin(x), его выполнение всегда может быть только приближенным.

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

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

представимы,— это свойства Д1+х) = х-Дх) и Д 1-х) = -х-Д-х) гамма-функции Дх). Иначе может нарушиться важное соотношение Г(п) = (и-1)! для целых положительных п.

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

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

• Корректное округление.

Помимо всех перечисленных ограничений, нужно потребовать, чтобы результат, возвращаемый реализацией, получался из точного результата функции для данного аргумента при помощи принятой в текущей конфигурации процедуры округления. Иногда практически бессмысленно требовать точности 0.5-1 ulp, но во всяком случае погрешность выше 1.5-2.0 ulp должна рассматриваться как неточность соответствующей реализации.

Стандарт IEEE 754 предписывает поддержку 4-х видов округления: к

ближайшему представимому числу, к +со, к - со и к 0.

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

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

• С помощью систем математических вычислений, например, Maple [74], Mathematica [75], MATLAB [76].

• С помощью библиотек корректно округляемых функций, таких, как разработанная на основе работ Ziv [77] IBM Accurate Portable MathLib [78], GNU MPFR [79], разрабатываемая Sun libmcr [80] или библиотек SCSLib [81] и CRlibm [82], разрабатываемых в рамка проекта Arenaire [16,83-86].

• Можно также разработать собственную реализацию функции на основе методов,

изложенных в книгах [18,19] и многочисленных статьях, или на основе методов интервальных вычислений [87-90].

4.2. Метод выбора тестовых данных

Предлагаемая методика выбора тестовых данных для тестирования математических функций основана на особенностях представление чисел с плавающей точкой, результатах работ [16,17,72] по вычислению «неудобных» значений аргументов и приведенном выше методе определения требований к реализациям таких функций, а также на технике построения тестов при помощи разбиения интервалов входных данных.

Основные шаги этого метода состоят в следующем.

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

• К исходному множеству добавляем числа 0, -0, +со, -со, минимальные и максимальные по абсолютной величине представимые числа, минимальное и максимальное денормализованные числа.

• Полученное исходное множество разбивает числа с плавающей точкой на набор интервалов.

Для выбранных числовых параметров п и к на каждом из этих интервалов возможные тестовые значения выбираются следующим образом. Сначала интервал разбивается на п более мелких интервалов, равных по количеству содержащихся в них чисел с плавающей точкой. При этом возникает (и+1) точка. Затем берутся все числа, лежащие в рассматриваемом интервале и отстоящие не более чем на к чисел с плавающей точкой от полученных точек. Получаемое множество точек назовем пробным множеством.

*iMtl 1111***»*! 1111 * it * ill * I I I I I ft*»** I I I-1-1-1-*—*—*-

Рисунок 5. Схема выбора тестовых данных на одном интервале с параметрами п=4, А 2.

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

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

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

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

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

4.3. Примеры применения предложенных методов

В этом разделе будут рассмотрены два примера использования предложенных методик для формулировки требований и выбора тестовых данных для тестирования двух математических функций — экспоненциальной функции и тангенса. Оба примера будут рассматриваться для случая чисел двойной точности.

4.3.1. Экспонента

В данном разделе термин «экспонента» везде обозначает экспоненциальную функцию ехр(х) = е .

• Область определения и особые точки.

Экспонента определена для всех действительных чисел, особых точек у нее нет.

Она стремится к +со при х —» +со, поэтому нужно найти соответствующие граничные значения — наибольшее число с плавающей точкой, для которого значение экспоненты еще может быть представлено конечным числом, и наименьшее такое, что экспонента от него должна быть равна +со.

Максимальное конечное представимое в виде числа с двойной точностью значение равно 2971-(253-1). Его натуральный логарифм равен 709.7827128933839967322233899106571455.„ю, ближайшее к этому значению представимое число есть

709.7827128933839730962063185870647430419921875ю = 1011000101.11001000010111111101111101000111001111011112 Экспонента от этого числа равна 1.1111111111111111111111111111111111111111111100101010 00011011...2-21023. Здесь и далее выделены биты, не помещающиеся в мантиссу числа двойной точности.

Экспонента от следующего представимого числа 709.782712893384086783044040203094482421875ю = 1011000101.11001000010111111101111101000111001111100002 уже не может быть представлена конечным числом двойной точности. Таким образом, для всех представимых чисел, больших или равных этому, реализация экспоненты должна возвращать результат +со.

Для числа

1011000101.11001000010111111101111101000111001111011112

при режимах округления к ближайшему, к 0 или к -со должно выдаваться

1.11111111111111111111111111111111111111111111001010Ю2-21023.

а при округлении к +со должен выдаваться результат

1.111111111111111111111111111111111111111111110010101Ь-21023.

Специальные значения, значения в 0, касательные и асимптоты.

о Специальные значения.

Естественно доопределить значения экспоненты в специальных точках следующим образом.

ехр(-О) = 1; схр( / ) = /:

ехр(- / ) = 0; е.гр(Ка>1) = КаК

о Точные значения и их окрестности, окрестность 0.

Известно только одно точное значение экспоненты: ехр(0) = 1. Учитывая большую плотность чисел с плавающей точкой в 0, нужно определить окрестность 0, в которой экспонента должна возвращать 1.

Для чисел, больших 0, имеем

ехр{\.\ 1111111111111111111111111111111111111111111111111112-2'54) =

1.0000000000000000000000000000000000000000000000000000 0110бОО 1... 2;

е.тр(1.02-2"53) =

1.0000000000000000000000000000000000000000000000000000 1053100... 2;

ехр{\.\ 1111111111111111111111111111111111111111111111111112-2'53) =

1.0000000000000000000000000000000000000000000000000000 1105010... 2;

схр( 1,0:-2"-) =

1.0000000000000000000000000000000000000000000000000001 052100... 2.

Поэтому, для чисел от 0 до (253-1)-2‘106 =

1.1111111111111111111111111111111111111111111111111 llb-2'54 экспонента должна возвращать 1 при режимах округления к ближайшему, к О и к -со. При округлении к +со для всех таких чисел, больших 0, результатом должно быть 1+2"52.

Для чисел от 2'53 до (253-1)-2'105 =

1.1111111111111111111111111111111111111111111111111 llb-2'53 экспонента должна возвращать 1 при режимах округления к 0 и к -со, а при режимах округления к ближайшему или к +со, на этом интервале должно возвращаться следующее за единицей представимое число 1+2'52.

Для значения аргумента 2"52 должно возвращаться 1+2"52 при режимах округления к ближайшему, к 0 и к -со, а при округлении к +со — 1+2'51.

Аналогично, для чисел, меньших О

елт?(-1.00000000000000000000000000000000000000000000000000012-2'53) =

1.1111111111111111111111111111111111111111111111111110 1520 1 0... 2"2-1;

ехр(-1.02'2"53) =

1.1111 1111 1111 1111 1111 1111 1111 1111 1111 1111 1111 11111 111 0541 11... 2-2 1;

ел7?(-1.00000000000000000000000000000000000000000000000000012-2'54) =

1.1111111111111111111111111111111111111111111111111111 0152001...2-2'1;

exp{-1.02-Tii) =

1.1111 1111 1111 1111 1111 1111 1111 1111 1111 1111 1111 11111 111 1055111... 2'2-1.

Поэтому, для -(1+2'52)-2‘53 экспонента должна возвращать 1-2"52 при округлении к - со или 0, и 1-2"53 при округлении к ближайшему или +СО.

На интервале от -2'53 до -(1+2'52)-2‘54, включая его концы, нужно возвращать 1-2"53 при округлении к -со, 0 или ближайшему, и 1 при округлении к +со. Для чисел от -2"54 до 0 нужно возвращать 1-2"53 при округлении к -со или 0 (кроме 0, для которого результат всегда 1) и 1 при округлении к ближайшему или +со.

Итоговые требования к значениям экспоненты в окрестности 0 для всех режимов округления суммированы в Таблице 1.

Округление К -СО к 0 к ближайшему К +00

Интервал/ то чка

-(1+2'52)-2'53 1-2'52 1-2'52 1-2-5Э 1-2-5Э

от -2'53 до -(1+2'ы)-2'54 1_2-5э 1_2-5э 1-2-5Э 1

о-54 0-1074 от -2 до -2 1_2-5э 1_2-5э 1 1

0 1 1 1 1

от 2'1074до (253— 1 )-2"106 1 1 1 1+2'52

от 2'53 до (253-1)-2'105 1 1 1+2'52 1+2'52

2'52 1+2'52 1+2'52 1+2'52 1+2'51

Таблица 1. Требования к поведению экспоненты в окрестности 0.

Асимптоты.

Экспонента стремится к 0 при х —» -со, поэтому нужно определить значения аргумента, при которых значение экспоненты становится неотличимым от 0, а также важно знать, когда оно становится денормализованным. Минимальное положительное представимое число равно 2"1074, а минимальное положительное нормализованное число равно 2"1022. Их натуральные логарифмы равны, соответственно, -744.440071921381262314107298446.„ю и -708.396418532264106224411228130... ю.

Ближайшие к ним представимые числа —

-744.4400719213812180896638892590999603271484375ю =

-1011101000.01110000101010001000110110101110001110000112 и -708.396418532264078748994506895542144775390625ю =

-1011000100.01100101011110111010111101010111100110100102.

Использование различных режимов округления, однако, «размазывает» эти границы на целые интервалы. Чтобы вычислить эти интервалы, определим значения, экспоненты которых при округлении к ближайшему будут округляться вверх и давать 2"1074, а затем такие, экспоненты которых при этом будут округляться вниз и давать тот же результат. Имеем

ехр(-1011101001.00100010000110101010010110100110000010100102 = Хі) =

0.145011...2*2"1075

1011101001.00100010000110101010010110100110000010100012 = Х2) =

1.0431Ю...2*2"1075

ехр{-1011101000.01110000101010001000110110101110001110001002 = Хз) =

1.142011...2*2"1075

ехр{-1011101000.01110000101010001000110110101110001110000112 = х4) =

1.044110...2*2"1074

Аналогично, для минимального положительного нормализованного числа имеем

ехр(-Ш 1000100.01100101011110111010111101010111100110100112 = Х5) 0.1111111111111111111111111111111111111111111001111011 1011. °'1022

ел/7(-1011000100.01100101011110111010111101010111100110100102 = -Т6) =

т-1022

1.0000000000000000000000000000000000000000000001111011 1011...

2-2'1

Результирующие требования к значениям экспоненты, близким к минимальному положительному представимому числу и минимальному положительному нормализованному числу, приведены в Таблице 2.

Округление К -СО к 0 к ближайшему К +00

Интервал/ то чка

> д: 0 0 0 2-1074

х3 > х > х2 0 0 г4074 2-1074

х4 2-1074 2-1074 г'!074 2-1073

Х5 > X > х4 д енормализованно е

Х5 0.1430011110112*2'1022 0.1430011111002*2"1022

Хб 1.04511110112*2"1022 1.04511111002*2"1022

Таблица 2. Требования к поведению экспоненты при переходе ее значений от нуля к положительным и от денормализованных к нормализованным значениям.

Асимптотика ехр(х)~1+х с трудом поддается фиксации в виде четких требований по ее соблюдению (см. раздел 4.1), поэтому никаких связанных с ней ограничений не налагается.

• Область значений.

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

• Монотонность и сохранение знака.

В связи со сказанным в предыдущем абзаце тот факт, что экспонента имеет только положительные значения, можно не проверять отдельно.

Экспонента — возрастающая функция. Соответственно, ее реализация также должна возрастать.

• Симметрии и периодичность.

Единственное важное свойство экспоненты — ехр(х+у) = ехр(х)-ехр{у) — может выполняться лишь приблизительно для чисел с плавающей точкой, поскольку

значения экспоненты для всех таких чисел, кроме 0, трансцендентны согласно теореме Линдемана [91], следовательно, представимы лишь приблизительно. Если х или у равно 0, это свойство не добавляет ничего нового к условию ехр(0)= 1.

• Корректное округление.

Погрешность вычисления экспоненты согласно стандарту ISO/IEC 10967-2 [14] не должна превосходить 0.5-1.5ulp. Остановимся на таких требованиях. В качестве дополнительной информации можно выдавать предупреждение специального вида каждый раз, когда фактическая ошибка превышает 0.5 ulp. Для вычисления правильно округленного значения экспоненты можно воспользоваться методом, описанным в [84] или готовыми библиотеками [79,81,82].

• Тестовые значения.

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

о -со;

о [-1.1111 1111 1111 1111 1111 1111 1111 1111 1111 1111 1111 1111 11112-21023,

-1011101001.00100010000110101010010110100110000010100102];

о [-1011101001.00100010000110101010010110100110000010100012,

-1011101000.01110000101010001000110110101110001110001002];

о [-1011101000.01110000101010001000110110101110001110000112,

-1011000100.01100101011110111010111101010111100110100112];

о [-1011000100.01100101011110111010111101010111100110100102, -1.00000000000000000000000000000000000000000000000000012-2'53];

о [-1.02-2"53,

-1.00000000000000000000000000000000000000000000000000012-2'54]; о [-1.02-2"54, -1.02-2"1022];

о [-0.11111111111111111111111111111111111111111111111111112-2'1022, -1.02-2"1074];

о -0;

о 0;

Г1 л о-1074

о [1.02'2 ,

0.11111111111111111111111111111111111111111111111111112-2'1022]; г 1 л о-1022

о [1.02'2 ,

1.11111111111111111111111111111111111111111111111111112-2'54]; о [1.02-2"53,

1.11111111111111111111111111111111111111111111111111112-2'53];

о [1.02-2‘52.

1011000101.110010000101111111011111010001110011110111 lo];

о [1011000101.11001000010111111101111101000111001111 ЮОООо,

1.111111111111111111111111111111111111111111111111111l2-21023];

О +'СО.

На полученных интервалах нужно строить тестовые значения в соответствии с методом, описанным в разделе 4.2 — деля интервал на п частей и выбирая точки, отстоящие от полученных (л+1) точки на интервале не более чем на к чисел с плавающей точкой.

Кроме этого, к полученному множеству тестовых значений надо добавить нескольких представителей NaN и набор чисел, вычисление корректно округленного значения экспоненты для которых наиболее трудоемко. Приведем несколько таких примеров, взятых из [17].

ехр{-\. 11101101001100011000111011111011011000100111111010102-2'27) = 1.1111111111111111111111111000010010110011100111000100 16°000... 2-2 1

ел7?(-1.01000000000000000000000000000000000000000000001100102-2'46) = 1.1111111111111111111111111111111111111111111101100000 085101 ...2-2-1

ел7?(-1.00000000000000000000000000000000000000000000000000012-2'51) =

1.1111111111111111111111111111111111111111111111111100 0101Ю1... 2-2-1

ехр( 1.11111111111111111111111111111111111111111111111111112-2"53) =

1.0000000000000000000000000000000000000000000000000000 1105ОЮ... 2; это конец одного из выделенных интервалов;

едгр>( 1.11111111111111111111111111111111111111111111100000002-246) =

1.0000000000000000000000000000000000000000000001111111 184010... 2;

ехр{\.0001111111111111111111111111111111111111111110101llb-2"45) =

1.0000000000000000000000000000000000000000000010001111 184000... 2;

едр(110.000011110101001011110011011110101110110011111101002) = 110101100.01010000101101000000100111001000101011101110 058100...2

4.3.2. Тангенс

• Область определения и особые точки.

Тангенс определен для всех действительных чисел, кроме чисел вида ;t/2+jrn для целого п.

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

Любое число с плавающей точкой, большее 1, может быть представлено как т-2*, где m — натуральное число, не превосходящее 253-1, а к — натуральное число, не превосходящее 971. Если такое число расположено очень близко к

и/2+л;п = (2л+1)-лУ2, значит и/2 ~ яг-2*/(2л+1). Таким образом, можно искать подходящие числа как числители наилучших рациональных приближений к лУ2. Для этого разложим числа л-2к при к от 52 до -971 в непрерывные дроби и, обрывая их в каком-либо месте, вычислим соответствующие подходящие дроби. Эти дроби приближают исходные числа наилучшим образом среди всех рациональных чисел, имеющих меньшие или такие же знаменатели [92].

Выполнение соответствующих вычислений дает 920 дробей, которые однозначно определяют числа с плавающей точкой двойной точности, приближающие кратные лУ2 с относительной погрешностью не больше, чем 2"57. Вот несколько первых таких чисел.

29-п]2 = 101101.10001101100101111000100010111011100100011011101 113010... 2 14479-11/2 =

29327-И/2 =

1011001111110010.1011111001101110101001101110001010011

И несколько последних.

1”011

197618317288004252206364340334461523249938251852771804320737665131462 691696567271297312161164216527568595012444936741107384701949271514894 825117746798310262428988560159244225487630579194650628623442695547951 691080982624236830670934028824074039328312895391481654408301108215527 83846121469368255050296717288515-И/2 =

1.0110000110100011110110111000110010001101000100101001 01026111

258215239476700584226130900745995853560281748618927419956422039974202 286570698806549395630266453068250888197386507164887712080806671331638 295216484470333049129476424667965728976743334784463915624846487050650 888702720312158190710847644636354056948141468538860276324737125629625 98749697007058695124283117548543-и/2 =

1.1100111000010100001100101110101011000111000001000111 11020011...2-21021

Хотя относительная погрешность приближения кратного п12 может достигать 2"1079, абсолютное расстояние до него для большинства таких чисел находится в пределах 10"15-10"17. Только 25 из этих чисел лежат от полюсов на расстоянии, не превосходящем 10"17, и только два — на расстоянии, не превосходящем 10"18. В Таблице 3 перечислены значения 25 нечетных кратных лУ2, ближайших к числам с плавающей точкой и значения тангенса для соответствующих чисел с плавающей точкой.

29-И/2 =

101101.10001101100101111000100010111011100100011011101 113010...2

fg(101101.100011011001011110001000101110111001000110111102) = -1.0110011010111001111010111100010010000101000011000110 01Ю...2-260

9206271 -n/2 =

110111001010100011111000.10101011100101110101110100101 Q30111... 2

tg{\ 10111001010100011111000.10101011100101110101 llOlOOlOb) = 1.0000010101110101100001001100010000101001101100111010 0110...2-259

138049179777104367775-ЛУ2 = 1.0111100000101011011110100010000011011111011011010100 072101...2-267

rg( 1.01111000001010110111101000100000110111110110110101002-267) = 1.0101101001100011110010000010101101101101101011100010 0100...2-257

1104381301317041933816362171-лУ2 =

1.0110011010111101010101000010010011100101011001010101 О95101... 2'290

tg{\.0110011010111101010101000010010011100101011001010ЮІ2-290) = 1.0101101001011110001010111000100001101010001010101010 l3010...2-257

227664762687035133079286258488902728695\-id2 =

1.0101000001001100101011000101000111110001111010101111 0137100...2-2131

rg(l .01010000010011001010110001010001111100011110101011112-2131) = 1.1110011111000111110100001111010000111111100000011101 10511... 2'25S

961727228574135 5388977328677436746192449293-л;/2 =

l.oioiioioiioioioiioiooiioooioiiooioiioooiiiooiiooiooo i147oii...2-2143

tg{\.010110101101010110100110001011001011000111001100100І2-2143) = -1.1100011110101011011100111000100111000010010011010111 0110...2-256

3982124097056689019860931670115672371080901279195265295-H/2 =

1.0000010100111001101101001000110100010100110001010101 o186ioo...2-2182

tg{\.000001010011100110110100100011010001010011000101010І2-2182) = 1.1010011111001011011111001011010000101000000010110000 03110...2-256

243709684824005521714504726240306266173066571317974432621719-лУ2 =

1.1110011111100100010010100111100010101100000110001011 l20300 0...2-2197

fg(l.ll 10011111100100010010100111100010101100000110001IOO2-2197) = -1.0000100110111010011111100000011000011111001100000000 1001...2-258

10177187664528444025108553957230483888384395218254373185418290015777 4575619-л;/2 =

l.oiioiooiiiioioioioiiooooiooiiooooioioooioiiiiooiioii o251ioi...2-2246

tg{\.01101001111010101011000010011000010100010111100110112-2246) = 1.0111001001010110011111001011100100000100011110101000 l3001...2-257

15624927791716915161759873344092142005165621843671972171731789642962 257024455-лУ2 =

1.1011001000011001011000110110010011010111010100001010 1258011...2-2253

tg{\. 101100100001100101100011011001001101011101010000101 ь-2253) = -1.1101110100010110100010001100011010010110100110101100 13011...2-257

71610307166660314043397602176728270128783849452911391884618409783207 463478073649776722681-лУ2 =

1.1100010001011100110100010001000101010100110111111Ю1 о300іоо...2-2295

г^(1.110001000101110011010001000100010101010011011 111 110І2-2295) = 1.1111110110111010010001101111010101010011101000001000 1 0310...2-257

79531021962858779827512141402968504124268502489089591753854287169082 2757973300060931533046037509210740152847560365069332549-И/2 = 1.111000111100101010011011011011000110010101011100Ю11 о412іоо...2-2408

tg{\. 11100011110010101001101101101100011001010101110010112-2408) = 1.1001111110011110010101010001001001100110110000111100 01300...2-256

43776696795279815301055904915945700987301503939761311288501927384592 34187554875937911988865352801115519693682101780670312164186360396953 23195218767-11/2 =

1.1011100010001100101110110100111000110010010101110Ю1 і491оп...2-2487

tg{\. 10111000100011001011101101001110001100100101011101102-2487) = -1.1010001011000010110011001011110000011010011000110010 0101...2-256

11589330310381988626449129598171351352584091074116125961273394956149 19456683050991296076480489809564268564293673421768988424019510581475 74305293329168873163869152174973-лУ2 =

1.10001011001010000110011101101100110111001100010110Ю і560ооо...2-2555

tg{\. 10001011001010000110011101101100110111001100010110112-2555) = -1.0000101000010001011001101101100011001010011011100010 0101...2-257

40912440333305918481550440392603355393318151657874369033644916645778 80843430257120834658307295698261362420794202530508444695188784113250 59263543780254873744888889331055895843-11/2 =

1.0100110010010110110000010001000100110100110100110Ю1 і583оіо...2-2577

tg{\.01001100100101101100000100010001001101001101001101102-2577) = -1.0110010101011100111110011110001000111100010011010111 03101... 2-258

65429133524241263178237370719527517029251219555784828550385052067786 68733143053622851340353340621053795171428918994903139424630388205704 7519616430579116455387138549345777195516511603705-лУ2 =

1.1000001100000000100111100010111010011110001011101010 і619оіі...2-2614

tg{\. 10000011000000001001111000101110100111100010111010112-2б14) = -1.1110000101100010000111111011000100110011011010100100 12011...2-257

59287932212270746098210389073106720676473902052684288108187468991207 27939955976870793398691724674236683942417720739021351535557283241575 82755209798363213599827165644437783394500726283264184124810624461674 9-КІ2 =

1.1101101101000001111100111100101101110001110101111011 0686111...2-2680

tg{\. 11011011010000011111001111001011011100011101011110112-2б80) = 1.0001000010111010011011110011111101000110011111010110 12010...2-258

30636572260931163506332425320948767147146434821312896452845246587739 42999960543512781934069734059379292127565945745701814493759190791305 15095712133324352639295031679514263477241946753392456572308669671533 6885-лУ2 =

1.1101111110101000110100011000111100101011001111101101 1693011...2-2689

г^(1.11011111101010001101000110001111001010110011111011102-2689) = -1.1011100111111000000110011010110110001011010111100111 10511...2-256

94971440572386015595571166139691201230707691489980939180815478595748 85577578104219175989760881732308451273779191276410784862786933003137 63136525734828774533668148005758919334015282384829230967867833601055 56253079245416317281232917440947587-и/2 = 1.0110111010001101011101111000110010010100110101100110 0798100...2-2794

tg{\.01101110100011010111011110001100100101001101011001102-2™) = 1.1001111001110100010011011111100111010110100011111100 0101...2-256

33864178045159811206438920823311565991202393932998380352421215184285 37554064774221620930267583474709602068045686026362989271814411863708 49986972132271594662263430201169763297290792255889271083061603403854 1342154669787134871905353772776431251615694251273653-лУ2 =

1.0110101011000101101100100110001011001010000111111110 1857011...2-2849

г^(1.01101010110001011011001001100010110010100001111111112-2849) = -1.1101100110111010100110100111100101110101011000110101 1010...2-260

88685758294706272204189579187820280422409877285489361868254170927334 81411325935922524369633276190848778076190586339609015909756150954578 31282700320668157612019880948068064298959801704400266492369256967459 715372092514867668494308766581538954683 8546388224593403-л;/2 =

1.1100111111100100100000100010100001011111100011101101 0866111...2-2860

fg(l.110011111110010010000010001010000101111110001110110І2-28бо) = 1.0001101000100100011010110110111110010010010001010110 12000...2-258

37105172474923906782576160374346777164937755101941743463232913447611 26458820318029413443054532013548715976758198761455677008210344281793 85754829765785467791950462488227489950141733334740102736828752663860 04962920828212934138031374382514718179363737540556512625927642564615 90125145265-11/2 =

l.oioooooioooioiiioioioiiiooiiooiiiooioiiiiioioioooooi i944ooo...2-2939 fg(l.010000010001011101010111001100111001011 11101010000102-2939) = -1.0001110001110010100101001010111011100111101110100100 01401. . . 2-257

22795463807341151743454041517545705884702298848407846038370271485088 07014891141753237567442204768461160038124562964557539085826195999224 02145326867889285827944189828230765840337018227362357217203013966072 70621803425566565987094025434998287231727147847848195901676458199728 364279528014629-11/2 = 1.1110000110011000011100010010001010110111111000000110 0955101...2-2951

fg(l.ll 100 0 01100110 0 0011100 0 1 00 1 000 1 0 1 011011111100 0 0 0 01102-2951) = 1.0111000101110100000001111101001000000001111001100110 13011...2-256

24982829712749183769184949361370722806939675044978450598655040276224 22069326569285540059089298879071869587152103490382783870770007392787 73724242538703521095769835440603290627555914001656747466652577114915 40415860945893559496946376666124169735622128694086982311970782271594 626856299313449263532321973-11/2 =

1.1110000000001001110001010011000101001000101111100001 0997100...2-2991

г^(1.111000000000100111000101001100010100100010111110000І2-2991) = 1.1011100011001011111101111110110011111010111111111010 02111...2-258

19761831728800425220636434033446152324993825185277180432073766513146 26916965672712973121611642165275685950124449367411073847019492715148 94825117746798310262428988560159244225487630579194650628623442695547 95169108098262423683067093402882407403932831289539148165440830110821 552783846121469368255050296717288515-лУ2 =

1.0110000110100011110110111000110010001101000100101001 01026111...2-21021

fg(l.ll 1000000000100111000101001100010100100010111110000 Ь-21021) = 1.0001001010111100000100001001010101100001110011010100 1 0311...2-257

Таблица 3. Нечетные кратные їг/2, ближайшие к числам с плавающей точкой и соответствующие значения тангенса.

В 20-м ряду Таблицы 3 приведено число с плавающей точкой, ближайшее к нечетному кратному к12. Это же число упоминается в статье [93], там оно вычислено с помощью программы, созданной W. Kalian и S. McDonald [94].

Из этого следует, что приведенное в 20-м ряду значение тангенса является максимальным по абсолютной величине. Таким образом, для всех чисел с

плавающей точкой двойной точности при всех режимах округления значения тангенса ограничены по абсолютной величине 2.133485385753703936ю-1018. Более того, для всех чисел, не представленных в Таблице 3, значения тангенса должны по абсолютной величине быть меньше, чем минимальное значение тангенса в этой таблице, находящее в 23-м ряду и равное 1.0399184334316502410-1017.

• Специальные значения, значения в 0, касательные и асимптоты.

о Специальные значения.

Значение тангенса в точке -0 можно определить равным как 0, так и -0. Второй способ выглядит предпочтительнее, поскольку он более точно отражает асимптотику Г£(.г) ~ х при х ~ 0 и имеет аналогию в стандарте ШЕЕ 754, определяющем ,5дгГ(-0) = -0.

Поскольку тангенс не имеет пределов в +со и -со, его значением в этих точках можно считать только ИаМ. Таким образом, получаем следующее.

&-0) = -0;

tg(+ со) = №№

=

г^ИаМ) = КаК

о Окрестность 0, асимптотика.

Известно, ЧТО ?£(0) = 0. Кроме ТОГО, Г£(.г)~.г при .т~0. Определим окрестность 0, в которой значение тангенса должно совпадать со значением его аргумента.

Имеем

г^(1.00100101000010111111111000011011000010000010111101002-2-26 = *0 = 1.0010010100001011111111100001101100001000001011110100 015°01... 2-2'26 = к

г^(1.00100101000010111111111000011011000010000010111101012-2"26 = .т2) = 1.0010010100001011111111100001101100001000001011110101 105210... 2'2"26

= ы

г^(1.01110001001101110100010010010001001000111110111101012-2-26 = Л-3) = 1.0111000100110111010001001001000100100011111011110101 150011...2'2"26

= к

г^(1.01110001001101110100010010010001001000111110111101102-2-26 = .т4) = 1.0111000100110111010001001001000100100011111011110111 055100... 2-2-26

= и.

Следующие из этих соотношений требования к значениям тангенса в окрестности 0 для всех режимов округления сформулированы в Таблице 4. В этой таблице выражение .т++ для числа с плавающей точкой х обозначает непосредственно следующее за х число с плавающей точкой.

Анализ других нулей тангенса, л-п при п Ф 0, см. ниже, в пункте о монотонности и сохранении знака.

Другие точные значения.

Кроме tg(0) = 0 заслуживают внимания соотношения tg(кl 4+к-п) = 1 и tg(-'I^I А+к-п) = -1.

Округление К -СО к 0 к ближайшему К +00

Интервал/точка

X = -х4 -(*4 ++) -и

-х2 >Х> -Хз -((-*)++) л; -((-*)++) л:

0 > д: > -Хі -((-*)++) л:

х = 0 0

хі>х>0 л; д:++

х3>х>х2 л; Х+ +

х = х4 и ?4+ +

Таблица 4. Требования к поведению тангенса в окрестности 0.

Вычисление чисел, приближающих нечетные кратные лУ4 с погрешностью, не превосходящей 2"52, дает 726 значений, которые могут рассматриваться как тестовые данные. Однако для многих из них тангенс не всегда должен быть равен 1 или -1. Например, для числа, ближайшего к лУ4, значение тангенса равно

1.1111111111111111111111111111111111111111111111111111 01300... 2"2-1, что дает 1 только при режиме округления К +СО.

Таких чисел с плавающей точкой, для которых 1 является ближайшим числом к значению их тангенса 199. В Таблице 5 перечислены только 3 из них, которые отстоят от кратных и/4 не более чем на 10"19.

29-лУ4 =

10110.110001101100101111000100010111011100100011011101 113010... 2

&(10110.1100011011001011110001000101110111001000110111102) = 1.0000000000000000000000000000000000000000000000000000 08101. . . 2

9206271-71/4 =

1.1011100101010001111100010101011100101110101110100101 030111...2-222

£(1.101110010101000111110001010101110010111010111010010І 2-222) =

-1.0000000000000000000000000000000000000000000000000000 07111...2

33864178045159811206438920823311565991202393932998380352421215184 28537554064774221620930267583474709602068045686026362989271814411 86370849986972132271594662263430201169763297290792255889271083061 6034038541342154669787134871905353772776431251615694251273653-лУ4

1.0110101011000101101100100110001011001010000111111110

1857011...2-2848

rg(l.011010101100010110110010011000101100101000011111111І2-2848) =

1.0000000000000000000000000000000000000000000000000000 0s 100... 2

Таблица 5. Кратные їг/4, ближайшие к числам с плавающей точкой и соответствующие значения тангенса.

о Асимптоты.

У тангенса нет горизонтальных асимптот.

Асимптотика в 0 рассмотрена выше.

Область значений.

Область значений тангенса не имеет математических границ.

Для чисел с плавающей точкой, как мы выяснили выше, тангенс не должен превосходить по абсолютной величине 2.133485385753703936ю-1018. Однако это требование следует из ограничений на точность вычислений.

Из Таблицы 4 также следует, что значение тангенса денормализовано в точности тогда, когда денормализован его аргумент.

Монотонность и сохранение знака.

Тангенс возрастает на каждом отрезке от л!2+лп до 7і/2+7г(л+1). Монотонность имеет смысл проверять только для тех значений аргументов, которые попадают в один такой отрезок.

Тангенс отрицателен на отрезках от 7і/2-(2л-1) до л-п, равен 0 в точках л-п и положителен на отрезках от л-п до л/2-(2п+1).

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

Можно использовать те же 920 дробей, которые были найдены нами ранее для приближения л!2 — дроби для л отличаются только удвоенными числителями, а соответствующие числа с плавающей точкой имеют ту же мантиссу и на единицу большую экспоненту.

Ближайшие к нулям тангенса числа с плавающей точкой выглядят несколько иначе, чем ближайшие к полюсам, хотя соответствуют почти тем же числам, умноженным на л. Все такие числа, расстояние от которых до кратных л не превосходит 10"17, а также значения тангенса для них представлены в Таблице 6.

29-л =

1011011.0001101100101111000100010111011100100011011101 113010... 2

fg(1011011.000110110010111100010001011101110010001101 lllOo) = 1.0110110101100001101101011000110010011001110001000010 l4000...2-2'60

9206271 -л =

1101110010101000111110001.0101011100101110101110100101 030111... 2

tg{\ 101110010101000111110001.0101011100101110101ІІОІООЮІ2) = -1.1111010101001111010100100010011110100100111010000011 15011... 2'2“59

2276647626870351330792862584889027286951-л =

1.0101000001001100101011000101000111110001111010101111 0137100...2-2132

г^(1.01010000010011001010110001010001111100011110101011112-2132) = -1.0000110010110110000001001101001101001111001101000001 02100...2'2"58

243709684824005 521714504726240306266173066571317974432621719-л; =

1.1110011111100100010010100111100010101100000110001011 1203000... 2-2198

fg(l.ll 100111111001000100101001111000101011000001100011002-2198) = 1.1110110101000001010111110101010011001000110010111000 0101...2-2'58

10177187664528444025108553957230483888384395218254373185418290015777 4575619-л =

1.0110100111101010101100001001100001010001011110011011 0251101... 2'2247

tg{\.01101001111010101011000010011000010100010111100110112-2247) = -1.0110000111101100111011001001110001010111011111111101 02111... 2"2"57

15624927791716915161759873344092142005165621843671972171731789642962 257024455-л =

1.1011001000011001011000110110010011010111010100001010 1258011...2-2254

tg{\. 10110010000110010110001101100100110101110101000010112-2254) = 1.0001001010111011101111100000011111110010110111100001 02110... 2"2"57

71610307166660314043397602176728270128783849452911391884618409783207 463478073649776722681-л =

1.1100010001011100110100010001000101010100110111111101 о300юо...2-2296

г^(1.110001000101110011010001000100010101010011011 111 110І2-2296) = -1.0000000100100100001010000111011011010111110000011010 10410...2-2'57

40912440333305918481550440392603355393318151657874369033644916645778 80843430257120834658307295698261362420794202530508444695188784113250 59263543780254873744888889331055895843-л = 1.0100110010010110110000010001000100110100110100110101 1583010...2-2578

tg{\.01001100100101101100000100010001001101001101001101102-2578) = 1.0110111011000110011110111100111101110111010100100010 0100...2-2'58

65429133524241263178237370719527517029251219555784828550385052067786 68733143053622851340353340621053795171428918994903139424630388205704 7519616430579116455387138549345777195516511603705-и =

1.1000001100000000100111100010111010011110001011101010 1619011...2-2615

1ё{\. 10000011000000001001111000101110100111100010111010112-2б15) = 1.0001000001001000001100000100000000001001000100101001 1011...2-2'57

59287932212270746098210389073106720676473902052684288108187468991207 27939955976870793398691724674236683942417720739021351535557283241575 82755209798363213599827165644437783394500726283264184124810624461674 9-71 =

1.1101101101000001111100111100101101110001110101111011 0686111...2-2681

г^(1.110110110 1 00 0 0 01111100111100 1 01101110 001110 1 011110112-2681) = -1.1110000010011000011110001101001000101110001010011101 0-101... 2-2-58

33864178045159811206438920823311565991202393932998380352421215184285 37554064774221620930267583474709602068045686026362989271814411863708 49986972132271594662263430201169763297290792255889271083061603403854 1342154669787134871905353772776431251615694251273653-71 =

1.0110101011000101101100100110001011001010000111111110 1857011...2-2850

г^(1.01101010110001011011001001100010110010100001111111112-2850) = 1.0001010010101110011100101110011010111010001000101110 14010...2-2"60

88685758294706272204189579187820280422409877285489361868254170927334 81411325935922524369633276190848778076190586339609015909756150954578 31282700320668157612019880948068064298959801704400266492369256967459 7153720925148676684943087665815389546838546388224593403-71 =

1.1100111111100100100000100010100001011111100011101101 о86611 1...2-2861

г^(1.11001111111001001000001000101000010111111000111011012-2861) = -1.1101000010001111010110011100001100101011111110100010 12010 ...2-2'58

24982829712749183769184949361370722806939675044978450598655040276224 22069326569285540059089298879071869587152103490382783870770007392787 73724242538703521095769835440603290627555914001656747466652577114915 40415860945893559496946376666124169735622128694086982311970782271594 626856299313449263532321973-71 =

1.1110000000001001110001010011000101001000101111100001 0997Ю 0...2-2992

rga.ll 1000000000100111000101001100010100100010111110000 Ь-2992) = -1.0010100101011010001110110000101001100100101100011101 0100... 2-2'58

Таблица 6. Кратные п, ближайшие к числам с плавающей точкой и соответствующие значения тангенса.

Отметим, что ближайшее к числу с плавающей точкой кратное к число, находящееся в 11-м ряду Таблицы 6, соответствует точно такому же целому

числу, что и ближайшее к числу с плавающей точкой нечетное кратное л/2. Симметрии и периодичность.

Тангенс — нечетная функция, £(-х) = -£(х). Это свойство должно учитываться при построении всех тестов с отрицательными значениями аргумента. Периодичность с трансцендентным периодом л, как и симметрии вида tg(^^l2-x) = 1 /£(х), выполняются на числах с плавающей точкой лишь приблизительно и, поэтому не могут быть использованы для формулировки строгих требований.

Корректное округление.

Погрешность вычисления тангенса согласно стандарту 180/1ЕС 10967-2 [14] не должна превосходить 0.5-2.0 и1р. В качестве дополнительной информации можно выдавать предупреждение специального вида каждый раз, когда фактическая ошибка превышает 0.5 и1р.

Для вычисления правильно округленного значения тангенса можно использовать известные библиотеки вычислений с заданной точностью [79,81,82].

Тестовые значения.

Идеальный набор интервалов для дальнейшего построения тестовых данных для тестирования реализаций тангенса мог бы состоять из всех его периодов, в которые попадает хотя бы два числа с плавающей точкой. Однако, таких интервалов слишком много — около 5.7-1015. Поэтому необходимо выбрать только некоторые из таких периодов, руководствуясь каким-то правилом.

Предлагается в качестве основы для тестов выбрать периоды тангенса, содержащие хотя бы одну из упоминавшихся выше точек, близких к кратным л, л/2 и л/4 (включая 0). Таких точек меньше, чем 3000, поэтому получающийся набор тестовых данных будет одновременно достаточно представительным и обозримым. Кроме того, к ним можно добавить несколько периодов, близких к содержащему 0, например, от -29л до 29л (раз уж число 29 дает одно из самых близких кратных л к числу с плавающей точкой). Для каждого из выбранных интервалов нужно вычислить его левый и правый концы — точки, ближайшие к нечетным кратным л/2 на этом интервале, его точку, ближайшую к кратному л, а также две точки на этом интервале, ближайшие к нечетным кратным л/4. Эти точки разбивают каждый такой период на четыре под-интервала. Все их занесем в исходное множество.

Кроме ТОГО, В исходное множество попадут ТОЧКИ XI, -XI, Х2, —Х2, ХЗ, -ХЗ, Х4, -Х4, вычисленные при рассмотрении асимптотики тангенса в 0, сам 0, -0, +со и -со, а также минимальное и максимальное положительные денормализованные числа и противоположные им.

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

Кроме этого, к тестовым значениям надо добавить нескольких представителей ИаИ и набор чисел, вычисление корректно округленного значения тангенса для которых наиболее трудоемко. Вот примеры таких чисел из [17].

г^(1.11011111111111111111111111111111111111111111000111112-2'22) = 1.1110000000000000000000000000000000000000000101010001 017801...2'2"22

1ё{\.011001111111111111111111111111111 11110100001000101002-2'18) = 1.0110100000000000000000000000000000001000111001100001 158010... 2-218

г^(1.01010000010010000110101100101111100001110000000101002-2"5) = 1.0101000001111000110011101011111111111001110001110010 105710...2'2"5

1ё{\.01000110101011000011011100100010010000110101001101 \02-2л) = 1.0111101110100100100111110111001110011000001010011110 155001... 2-21.

5. Заключение

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

Несмотря на наличие серьезных исследований этого вопроса, например, в работах [16,17,72,73], ведущихся в рамках проекта Агепа1ге [12], работ, в которых для построения тестов таких функций использовались бы их формальные спецификации,

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

Методика формирования требований построена на основе анализа особенностей и специфических элементов поведения самой математической функции с учетом использования чисел с плавающей точкой для представления ее аргументов и значений.

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

описываемого одним набором условий, поведения тестируемой функции. Кроме того, используются полученные в работах [17,72] числа с плавающей точкой, для которых определение корректно округленного значения данной функции требует наиболее высокой точности вычислений.

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

Литература

[1] I. Bourdonov, A. Kossatchev, V. Kuliamin, and A. Petrenko. UniTesK Test Suite Architecture. In Proc. of FME 2002. LNCS 2391, pp. 77-88, Springer-Verlag, 2002.

[2] В. В. Кулямин, А. К. Петренко, А. С. Косачев, И. Б. Бурдонов. Подход UniTesK к разработке тестов. Программирование, 29(6):25-43, 2003.

[3] V. Kuliamin, A. Petrenko, N. Pakoulin. Practical Approach to Specification and Conformance Testing of Distributed Network Applications. In M. Malek, E. Nett, N. Suri, eds. Service Availability. LNCS 3694, pp. 68-83, Springer-Verlag, 2005.

[4] IEEE 754-1985. IEEE Standard for Binary Floating-Point Arithmetic. NY: IEEE, 1985.

[5] IEC 60559:1989. Binary Floating-Point Arithmetic for Microprocessor Systems. Geneve: ISO, 1989.

[6] IEEE 854-1987. IEEE Standard for Radix-Independent Floating-Point Arithmetic. NY: IEEE, 1987.

[7] D. Goldberg. What Every Computer Scientist Should Know about Floating-Point Arithmetic. ACM Computing Surveys, 23(l):5-48, 1991.

[8] ISO/IEC 9899:1999. Programming Languages — C. Geneve: ISO, 1999.

[9] IEEE 1003.1-2004. Information Technology — Portable Operating System Interface (POSIX). NY: IEEE, 2004.

[10] G. Hanrot, V. Lefevre, J.-M. Muller, N. Revol, and P. Zimmermann. Some Notes for a Proposal for Elementary Function Implementation in Floating-Point Arithmetic. Proc. of Workshop IEEE 754R and Arithmetic Standardization, in ARITH-15, June 2001.

[11] D. Defour, G. Hanrot, V. Lefevre, J.-M. Muller, N. Revol, and P. Zimmermann. Proposal for a standardization of mathematical function implementation in floating-point arithmetic. Numerical Algorithms, 37(l-4):367-375, December 2004.

[12] http://www.inria.fr/recherche/equipes/arenaire.en.html

[13] ISO/IEC 10967-1:1994. Information Technology — Language Independent Arithmetic

— Part 1: Integer and Floating Point Arithmetic. Geneve: ISO, 1994.

[14] ISO/IEC 10967-2:2002. Information Technology — Language Independent Arithmetic

— Part 2: Elementary Numerical Functions. Geneve: ISO, 2002.

[15] ISO/IEC 10967-3. Information Technology — Language Independent Arithmetic — Part 3: Complex Integer and Floating Arithmetic and Complex Elementary Numerical Functions. Draft. Geneve: ISO, 2002.

[16] V. Lefevre, J.-M. Muller, and A. Tisserand. Toward Correctly Rounded Transcendentals. IEEE Transactions on Computers, 47( 11): 1235—1243,November 1998.

[17] V. Lefevre, J.-M. Muller. Worst Cases for Correct Rounding of the Elementary Functions in Double Precision. Proc. of 15-th IEEE Symposium on Computer Arithmetic, Vail, Colorado, USA, June 2001.

[18] M. Abramowitz and I. A. Stegun, eds. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. Dover Publications, 1965.

[19] J.-M. Muller. Elementary Functions: Algorithms and Implementation. Second edition. Birkhauser, Boston, 2006.

[20] http://www.zuser.Org/z/

[21] G. Barrett. Formal Methods Applied to a Floating-Point Number System. IEEE Transactions on Software Engineering, 15(5):611—621, May 1989.

[22] http://www.csr.ncl.ac.uk/vdm/

[23] B. A. Wichmann. Towards a Formal Specification of Floating Point. The Computer Journal, 32:432-436, 1989.

[24] http://www.cs.cornell.edu/Info/Projects/NuPrl/nuprl.html

[25] http://pvs.csl.sri.com/

[26] http://www.cl.cam.ac.uk/Research/HVG/HOL/

[27] http://www.cs.utexas.edu/users/moore/acl2/acl2-doc.html

[28] J. O’Leary, M. Leeser, J. Hickey, and M. Aagaard. Non-restoring Integer Square Root: A Case Study in Design by Principled Optimization. In T. Kropf, R. Kumar, eds. Proc. of the 2-nd International Conference on Theorem Provers in Circuit Design (TPCD’94): Theory, Practice, and Experience. LNCS901, pp. 52-71, Springer-Verlag, 1994.

[29] P. S. Miner. Defining the IEEE-854 Floating-Point Standard in PVS. Technical report TM-110167, NASA Langley Research Center, 1995.

[30] V. A. Carreno and P. S. Miner. Specification of the IEEE-854 Floating-Point Standard in HOL and PVS. Proc. of 10-th International Workshop on Higher Order Logic Theorem Proving and its Applications, Aspen Grove, Utah, USA, September 1995.

[31] V. A. Carreno. Interpretation of IEEE-854 Floating-Point Standard and Definition in the HOL System. Technical report TM-110189, NASA Langley Research Center, 1995.

[32] J. Moore, T. Lynch, and M. Kaufmann. A Mechanically Checked Proof of the Correctness of the Kernel of the AMD5K86 Floating-Point Division Algorithm. IEEE Transactions on Computers, 47(9):913-926, September 1998.

[33] J. Harrison. Theorem Proving with the Real Numbers. Technical report UCAM-CL-TR-408, University of Cambridge Computer Laboratory. UK, November 1996.

[34] D. Verkest, L. Claesen, and H. De Man. A Proof on the Nonrestoring Division Algorithm and its Implementation on an ALU. Formal Methods in System Design, vol. 4, 1994.

[35] M. Comea-Hasegan. Proving the IEEE Correctness of Iterative Floating-Point Square Root, Divide, and Remainder Algorithms. Intel Technology Journal, 1998.

[36] M. Comea-Hasegan. IA-64 Floating Point Operations and the IEEE Standard for Binary Floating-Point Arithmetic. Intel Technology Journal, Q4, 1999.

[37] M. D. Aagaard and C.-J. H. Seger. The Formal Verification of a Pipelined Doubleprecision IEEE Floating-Point Multiplier. Proc. ICCAD, IEEE, Nov. 1995, pp. 7-10.

[38] P. S. Miner, and J. F. Leathrum. Verification of IEEE Compliant Subtractive Division Algorithms. Proc. FMCAD'96, November 1996.

[39] E. M. Clarke, S. M. German, and X. Zhao. Verifying the SRT Division Algorithm Using Theorem Proving Techniques. Proc. CAV'96, LNCS 1102, Springer-Verlag, 1996.

[40] H. Ruess, N. Shankar, and M. K. Srivas. Modular Verification of SRT Division. Proc. CAV'96, LNCS 1102, Springer-Verlag, 1996.

[41] R. Kaivola and N. Narasimhan. Formal Verification of the Pentium 4 Floating-Point Multiplier. Proc. of Design, Automation and Test in Europe Conference and Exposition (DATE). IEEE, 2002.

[42] Y.-A. Chen, E. M. Clarke, P.-H. Ho, Y. Hoskote, T. Kam, M. Khaira, J. W. O'Leary, and X. Zhao. Verification of All Circuits in a Floating-Point Unit Using Word-level Model Checking. In Formal Methods in Computer-Aided Design, LNCS 1166, pp. 19-33, Springer-Verlag, 1996.

[431 J. R. Harrison. Verifying the Accuracy of Polynomial Approximations in HOL. TPHOLs'97. August 1997.

[44] D. M. Russinoff. A Mechanically Checked Proof of IEEE Compliance of the Floating Point Multiplication, Division and Square Root Algorithms of the AMD-K7 Processor. LMS Journal of Computation and Mathematics, 1:148-200, 1998.

[45] D. M. Russinoff. A Mechanically Checked Proof of IEEE Compliance of the AAID K5 Floating-Point Square Root Microcode. Formal Methods in System Design, 14 (1), January 1999.

[46] J. O'Leary, X. Zhao, R. Gerth, and C. H. Seger. Formally Verifying IEEE Compliance of Floating-Point Hardware. Intel Technology Journal, 1999.

[47] M. D. Aagaard, R. B. Jones, and R. Kaivola. Formal Verification of Iterative Algorithms in Microprocessors. Proc. DAC 2000. ACM, 2000.

[48] J. Sawada. Formal Verification of Divide and Square Root Algorithms Using Series Calculations. Proc. of ACL2 Workshop, Grenoble, France, April 2002.

[49] J. R. Harrison. Floating-Point Verification in HOL Light: the Exponential Function. Technical Report LTCAM-CL-TR-428, LTniversity of Cambridge Computer Laboratory. UK. June 1997.

[50] A. T. Abdel-Hamid. A Hierarchical Verification of the IEEE-754 Table-driven Floating-point Exponential Function Using HOL. Master's thesis, Dept. Electrical and Computer Engineering, Concordia LTniversity, Montreal, Quebec, Canada, 2001.

[51] J. Harrison. Formal Verification of Floating Point Trigonometric Functions. Proc. of Formal Methods in Computer-Aided Design, FMCAD 2000, LNCS 1954, pp. 217-233, Springer-Verlag, 2000.

[52] J. Harrison. Floating Point Verification in HOL. E. T. Schubert, P. J. Windley, and J. Alves-Foss, eds. Proc. of 8-th International Workshop on Higher Logic Theorem Proving and its Applications, Aspen Grove, LIT, LISA, September 1995. LNCS 971, pp. 186-199, Springer-Verlag, 1995.

[53] F. de Dinechin, C. Lauter, and G. Melquiond. Assisted Verification of Elementary Functions. INRIA Research Report RR-5683, September 2005.

[54] C. Jacob, C. Berg. Formal Verification of the VAAIP Floating Point Unit. In Formal Methods in System Design, 26, pp. 227-266, Springer, 2005.

[55] http://support.intel.com/support/processors/flag/tech.htm, Discussion of Flag Erratum,

2002.

[56] http://www.math.utah.edu/~beebe/software/ieee/

[57] http://people.redhat.com/drepper/libm/

[58] W. Kahan. What Can You Learn about Floating-Point Arithmetic in One Hour'l http://http.cs.berkeley.edu/~wkahan/ieee754status, 1996.

[59] N. L. Schryer. A Test of Computer's Floating-Point Arithmetic Unit. Computer Science Technical Report 89, AT&T Bell Labs, 1981.

[60] B. Verdonk, A. Cuyt, and D. Verschaeren. A Precision- and Range-Independent Tool

for Testing Floating-Point Arithmetic I: Basic Operations, Square Root and Remainder.

ACM TOMS 27( 1 ):92—118, 2001.

[61] B. Verdonk, A. Cuyt, and D. Verschaeren. A Precision- and Range-Independent Tool

for Testing Floating-Point Arithmetic II: Conversions. ACM TOMS 27(1): 119—140, 2001.

[62] http://www.cant.ua.ac.be/ieeecc754.html

[63] R. Karpinski. PARANOIA: A Floating-Point Benchmark. Byte Magazine 10, 2 (Feb.), pp. 223-235. 1985.

[64] http://www.netlib.org/paranoia/

[65] A. Ziv, M. Aharoni, and S. Asaf. Solving Range Constraints for Binary Floating-Point Instructions. Proc. of 16-th IEEE Symposium on Computer Arithmetic (ARITH-16’03), pp. 158-163, 2003.

[66] M. Aharoni, S. Asaf, L. Fournier, A. Koifman, and R. Nagel. FPgen — A Test Generation Framework for Datapath Floating-Point Verification. Proc. IEEE International High Level Design Validation and Test Workshop (HLDVT’03), pp. 17-22,

2003.

[67] http://www.netlib.org/fp/ucbtest.tgz

[68] http://www.math.utah.edu/pub/elefiint/

[69] W. Cody and W. Waite. Software Manual for the Elementary Functions. Prentice-Hall, Englewood Cliffs, NJ, 1980.

[70] Z. A. Liu. Berkeley Elementary Function Test Suite. M.S. thesis, Computer Science Division, Dept, of Electrical Engineering and Computer Science, University of California at Berkeley, December 1987.

[71] P.-T. P. Tang. Accurate and Efficient Testing of the Exponential and Logarithm Functions. ACM Transactions on Mathematical Software, 16(3):185-200, September 1990.

[72] D. Stehle, V. Lefevre, P. Zimmermann. Searching Worst Cases of a One-Variable Function Using Lattice Reduction. IEEE Transactions on Computers, 54(3):340-346, March 2005.

[73] http://www.loria.fr/~zimmerma/mpcheck/

[74] http://www.maplesoft.com/

[75] http://www.wolfram.com/products/mathematica/index.html

[76] http://www.mathworks.com/products/matlab/

[77] A. Ziv. Fast Evaluation of Elementary Mathematical Functions with Correctly Rounded Last Bit. ACM Transactions on Mathematical Software, 17(3):410-423, September 1991.

[78] IBM Accurate Portable MathLib http://rpmfind.net/linux/rpm2html/search. php?query=libultim. so. 2

[79] http://www.mpfr.org/

[80] http://www.sun.com/do wnload/products.xml?id=41797765

[81] http://www.ens-lyon.fr/LIP/Arenaire/Ware/SCSLib/

[82] http://lipforge.ens-lyon.fr/projects/crlibm/

[83] F. de Dinechin, A. Ershov, and N. Gast. Towards the post-ultimate libm. Proc. of 17-th Symposium on Computer Arithmetic. IEEE Computer Society Press, June 2005.

[84] D. Defour, F. de Dinechin, J.-M. Muller. Correctly Rounded Exponential Function in Double Precision Arithmetic. INRIA Research report RR-2001-26, July 2001.

[85] F. de Dinechin, C. Lauter, J.-M. Muller. Fast and Correctly Rounded Logarithms in Double-Precision. INRIA Research report RR-2005-37, September 2005.

[86] S. Chevillard, N. Revol. Computation of the Error Functions erf and erfc in Arbitrary Precision with Correct Rounding. Proc. of 17-th IMACS Conf. on Scientific Computation, Applied Math, and Simulation, Paris, France, July 2005.

[87] W. Kramer. Multiple-Precision Computations with Result Verification. In E. Adams, U. Kulisch, eds. Scientific Computing with Automatic Result Verification, pp. 325-356, Academic Press, 1993.

[88] M. J. Schulte, E. E. Swartzlander. Software and Hardware Techniques for Accurate, Self-Validating Arithmetic. Applications of Interval Computations, pp. 381-404, 1996.

[89] N. Revol, F. Rouillier. Motivations for an Aarbitrary Precision Interval Arithmetic and the MPF1 Library. Reliable Computing, ll(4):275-290, 2005.

[90] MPFI Library http://perso.ens-lyon.fr/nathalie.revol/mpfi_toc.html

[91] F. Lindemann. Uber die Zahl k. Mathematische Annalen, vol. 20, pp. 213-225, 1882.

[92] А. Я. Хинчин. Цепные дроби. М: Наука, 1978.

[93] К. С. Ng. Arguments Reduction for Huge Arguments: Good to the Last Bit. 1992. Доступна как http://www.validlab.com/arg.pdf.

[94] W. Kahan. Minimizing q*m — n. 1983. Неопубликованные заметки, доступны по http://http.cs.berkeley.edU/~wkahan/testpi/nearpi.c.