FPU посвящается (часть 2): Fortran и Асмъ — Архив WASM.RU

FORTRAN мой друг,
Приятель старый.
И очень вежливый старик.
Перед расчётами приник
Чопорный ум его не малый.

Теперь уже он домосед,
Но рад всегда своим знакомым.
Зайду и я, вгляжуся снова
В моря его библиотек…

АС Edmond

1. От автора
2. Благодарности
3. Стек…
4. Первые шаги
5. Кодирование математических выражений на FPU
   1. Мышление FPU кодера
   2. Психофизиологические возможности программиста
   3. Восприятие регистров FPU как множества, или FXCH
   4. Элементарные преобразования
   5. Операции сравнения
6. Функции и математические абстракции
   1. Трансцендентные функции
   2. Получение основных функций
7. Вопросы точности и округления
   1. Основы вычислительной математики
   2. Основные принципы ситуаций потери точности
   3. Методы и способы увеличения точности
8. Функции, использующие FPU
9. Заключение

1 От автора

Ну, вот мы снова перед могуществом FPU. Ещё немного, и мы заглянем в его суть. Программирование FPU – это действительно отдельная тема, и отдельный удивительный подход, который имеет применение только здесь. Во всём пространстве кодов x86, пространство кодов FPU позволяет утверждать, что из всех вариантов, есть только один – оптимальный, и этот вариант достаточно предсказуем, и формализуем.

FORTRAN не только первый из ЯВУ, но и наиболее простой классический. Это та причина, по которой современные инженеры продолжают лелеять язык. Им невдомёк от сложных конструкций C++ и ненужной «красоты» кода. Им нужна простота.

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

2 Благодарности

Особой благодарности автора по праву заслуживает редактор WASM.RU CyberManiac, который с нечеловеческим терпением прочитал черновую версию статьи. Он внёс в статью не просто свои замечания. Он внёс свежий ветер перемен, который разбудил сонное сознание автора. Утро мудрее вечера – ибо мысли рождаются заново.
Пусть таких редакторов будет больше!!!

3 Стек…

Можно долго рассуждать, почему регистры FPU недоступны напрямую, а только в виде стека. Достаточно верным окажется предположение об экономии средств. Действительно, создание восьми новых регистров повлекло бы за собой увеличение затрат на формирование команд. Впрочем, стек FPU – очень интересный стек, и несколько отличается от стека столь привычного. Это небольшое отличие, о котором, наверное, любят рассуждать только такие снобы, как автор ?, может и не столь важно. И, тем не менее, фраза «доступ к регистрам FPU осуществляется как элементам стека» – не есть верной. То есть, в командах FPU операнды стека могут появляться непосредственно, и хотя это похоже на [ebp+xx] (ST(1), ST(2)), всё-таки придаёт эффект прямого доступа к содержимому регистра . Такое утверждение легко оспорить, и тем не менее автор рекомендует образовывать своё мышление именно так, а не иначе.

Для начала, следует уяснить чёткое правило модификации стека. Для уяснения правил хорошо запомнить следующий рисунок:

R0	<= TOP	ST(0)
R1		ST(1)
R2		ST(2)
R3		ST(3)
R4		ST(4)
R5		ST(5)
R6		ST(6)
R7		ST(7)

Как видите увеличение ST(x) происходит в сторону возрастания Rx. Если мы будем говорить «TOP увеличился», это значит, что он сместился «вниз». Иначе – «вверх». Например:

R0		ST(5)
R1		ST(6)
R2		ST(7)
R3	<= TOP = 3	ST(0)
R4		ST(1)
R5		ST(2)
R6		ST(3)
R7		ST(4)

Итак, существует две операции характерных для стека: push и pop.

PUSH

1. Выполнить уменьшение ТОР
2. Поместить содержимое в стек

 

*Во всех известных автору источниках эту операцию описывают:

1. Поместить содержимое в стек
2. Выполнить уменьшение ТОР

Причина, по которой, автор поменял местами две эти операции – достижение большей ясности, так как вершина стека всегда ST(0). Программист думает, что сперва ТОР указывает на новый ST(0), и только потом заполняет его. Это не корректно, однако удобно программисту.

POP

1. Скопировать операнд в вершине стека
2. Увеличить значение TOP и пометить прежний ST(0) пустым

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

4 Первые шаги

Для начала, дабы не переусердствовать рассмотрим выражение:

с = a/z ! Где все числа REAL(4)  

Вот его трансляция:

17: c = a/b
 00401013 fld dword ptr [esp+8]
 00401017 fdiv dword ptr [esp+0Ch]  

Здесь, как вы догадались, a и b – параметры какой-то функции, это видно по их методу адресации.
Вот так всё просто, и нечего более. Теперь усложним выражение.

z = a-b
 b = z**2
 c = a+b+z  

Где, a, b, c – это статические переменные, имеющие атрибут SAVE. Код, который мы получаем:

10: z = a-b
 00401013 fld dword ptr [esp+8]
 00401017 fsub dword ptr [esp+0Ch]
 0040101B fst st(1)
 11: b = z**2
 0040101D fmul st,st(1)
 12: c = a+b+z
 0040101F fadd dword ptr [esp+8]
 00401024 fadd st,st(1)
 ……………………………………………………………  

Этого примера достаточно, чтобы познакомиться с мышлением программиста FPU. Начиная с первой строчки программы, мы помещаем значение a в стек, после чего вычитаем значение b, которое храниться в стеке. То, что должно заинтересовать программиста – инструкция fst st(1). Давайте разберёмся, что происходит при её выполнении. При этом нелзя забывать, что мы рассматриваем код, генерируемый компилятором, а значит не обязательно код, служащий образцом.
Инструкция fst копирует операнд в память либо в незанятый регистр FPU. Это делается для того, чтобы после результат выполнения z умножить на себя, то есть возвести в квадрат. Однако этот приём преследует также ещё одну цель: сохранить значение в st(0), для последующих операций.
А теперь модифицируем пример:

z = a-b
 b = z**2
 c = a+b+z**2  

Вот его код:

11: z = a-b
 00401013 fld dword ptr [esp+8]
 00401017 fsub dword ptr [esp+0Ch]
 0040101B fstp st(1)
 12: b = z**2
 0040101D fmul st,st(0)
 13: c = a+b+z**2
 0040101F fld dword ptr [esp+8]
 00401024 fadd st,st(1)
 00401026 fadd st,st(1)  

Посмотрите внимательно на этот код, и попытайтесь самостоятельно разобраться.
Следует отметить основное правило операндов команд FPU.
Один из операндов FPU обязательно должен быть регистром ST(0).

Это незримое правило не найти в официальных руководствах, однако его полезно запомнить.
Теперь, если вы проанализируете оба примера, то заметите несколько приёмов, которые будут попадаться очень часто. Обратите внимание на команду fstp st(1). Она совершенно лишняя… (что это? Секреты оптимизации FPU кода? ? Вряд ли.) Теперь посмотрите на команду fmul st,st(0) – классика возведения операнда в степень двойки. Его следует запомнить. В первом примере компилятор поступил мудро, скопировав значение z в ST(1) так как оно понадобилось бы при последующих вычислениях, а в ST(0) формировался результат.

Ещё более усложним пример:

! Приведён готовый листинг с комментариями
 11: d = a+b
 00401013 fld dword ptr [esp+0Ch]
 00401017 fadd dword ptr [esp+8]
 0040101B fstp st(1) ! Не обращайте внимания

 12: z = a-b

 0040101D fld dword ptr [esp+8] ! ST(0) = a;
 ! st(1) = a+b

 00401021 fsub dword ptr [esp+0Ch]
 00401025 fst st(2) ! Запомним z в st(2)

 13: b = z**2

 00401027 fmul st,st(2) ! z**2
 00401029 fstp st(3) ! После этой команды st(2) = b

 ! А st(1) = z

 14: c = (a+b+z**2+d*(a-z*b))/(a+b+z**2)

 0040102B fld dword ptr [esp+8] ! ST(0) = a

 ! st(1) = a+b; st(2) = z; st(3) = z**2

 0040102F fadd st,st(3) ! +b
 00401031 fadd st,st(3) ! +z**2
 00401033 fstp st(4) ! После этой команды st(3) = a+b+z**2 так как она ещё нужна
 ! А st(0) = b
 00401035 fld st(1) ! Копирует st(1) так что, st(0) = b, st(1) = a+b+z**2
 ! Теперь st(3) = z 00401037 fmul st,st(3) ! z*b
 00401039 fsubr dword ptr [esp+8] ! a-z*b
 0040103E fmul st,st(1)
 00401040 fadd st,st(4)
 00401042 fdiv st,st(4)
 00401044 ffree st(4)
 00401046 ffree st(3)
 00401048 ffree st(2)
 0040104A ffree st(1) ! dword ptr [esp+8] это a
 ! а dword ptr [esp+0ch] это b  

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

5 Кодирование математических выражений на FPU

5.1 Мышление FPU кодера

Основная задача FPU кодера «видеть» математическое выражение в контексте FPU. Это вовсе не раздел об оптимизации. Хотя в своей сути это и есть самая настоящая оптимизация. Может показаться странным, что оптимизация FPU идёт раньше обучения основных приёмов вычислений. На то есть веская причина – особенность FPU. В отличие от основного набора инструкций x86, набор инструкций FPU предполагает достаточно формализованную базу кодирования.
Чтобы приступить к задаче оптимизации – как процессу «выбора из многих», следует овладеть навыком специального мышления FPU кодера. Можно сказать, что речь идёт об овладении так называемой «бессознательной оптимизацией». Пример, приведённый выше – это яркий пример способа мышления FPU кодера.

Как только что говорилось, набор команд FPU достаточно однозначен, а значит, подлежит формализации. Наша задача написать код на выражение:

d = a+b
 z = a-b
 b = z**2
 c = (a+b+z**2+d*(a-z*b))/(a+b+z**2)  

Ясно, что, решая эту задачу, программист должен «смотреть» на этот пример не так как математик.
Существуют несколько этапов «развёртывания выражений»:

1. Разбитие выражения на операции.
2. Определение параллельных операций, квантовой последовательности/ зависимости.
3. Смешивание параллельных, независимых.

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

Рассматривая наш пример, получим:

Квант 1.
Операция 1. d = a+b
Операция 2. z = a-b

Квант 2.
Операция 1. z**2
Операция 2. a+b+z**2
Операция 3. z*b
Квант 3.
a-z*b
Квант 4.
a+b+z**2+d*(a-z*b)
Квант 5.
(a+b+z**2+d*(a-z*b))/(a+b+z**2)

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

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

1. Я вижу две независимые операции: a-b, и a+b
2. Это простые операции я пропускаю их, смотрю на квант 2, поскольку в кванте 2 нужна z, значит, я мысленно ставлю плюсик, что в st(0), должно оказаться z, тогда я смогу возвести в квадрат z.
3. Я постоянно смотрю на последний связанный квант «третьим» глазом, то есть, пассивно обращая на него внимание.
4. Я замечаю простое сходство в выражениях a+b, и a+b+z**2.
5. Значит, это есть контекст, при чём a+b мы уже имеем в первом кванте.
   Цель этого анализа – формирование предварительного впечатления о контексте. Как правило, кодинг происходит «изнутри наружу», то есть от простых выражений выполняющихся первыми к сложным, и «от много зависимого к простому», то есть первыми кодируются выражения насыщенные контекстными элементами, а после – выражения, содержащие меньший относительный контекст.

В нашем случае, кодирование следует начинать с нижней формулы, перед этим схватив контекстные элементы. Схватывание контекстных элементов доводиться до автоматизма – это элементарный поиск «одинаковых частей в формуле», при этом нужно помнить, что, из всех выражений нужно искать такую последовательность вычислений при которой переход от одного вычисления к другому вызовет минимальные затраты на смену состояния потока:

 Например:
   a+b
   a**2
   2*a

В этих трёх независимых выражениях количество комбинаций вычисления достаточно велико. Однако человек умеет откидывать заранее неэффективные комбинации, пытаясь найти минимум «движений» для вычисления выражений. Это просто сделать, заметив глазом, что переменная a попадается намного чаще, чем b. А значит, я могу забыть про b, и строить вычисления так, как будто бы у меня есть только a.

 fld a
   fst st(1)
   fst st(2)
   fmul st(2),st(0)
   fadd st(1),st(0)
   fadd b

Отсюда первое незримое правило FPU кодера:
1. Ищите те места, в которых часто встречается переменная или несколько переменных.

Глядя на предыдущий пример, давайте попробуем повторить его кодинг.

1. Смотрим на последнее выражение, и пытаемся в первом выражении достичь последнего.
2. Интересный, например вопрос: «Будет ли лучше, если поменять местами операции d = a+b, z = a-b.»

Второе незримое правило FPU кодера:
2. Из двух-трёх независимых соседних выражений, ищите такое, чтобы с верху, к нему было легче «прийти», а снизу «пристыковаться».

Глядя на два независимых выражения, мы видим, что a+b – контекст последнего выражения, а z и z**2 так же контекстно зависимы. Смотрим на частоту использования переменной, так же учитываем и правило обращения к переменной.

На данном этапе мы пытаемся решить дилемму: какую из двух переменных положить в стек FPU, то есть, какая более используемей: a или b.

Этот вопрос решается просто, я вижу, что недалеко находиться выражение a-z*b, и при этом b уже было переопределено. Такое предположение может и не оправдаться, когда стек FPU будет переполнен, однако мы делаем его.

Третье незримое правило FPU кодера:
3. Кодер FPU считает, что регистров FPU число безграничное.

Таким образом, мы решает, что нужно загрузит a, а не b.

fld a
 ;; Делаем предположение, что может понадобиться и b
 fld b
 fst st(4)
 fadd st(0),st(1) ;; a+b
 ;; Запоминаем a+b
 fst st(3)
 ;; Внимание фокус!!!
 fsub st(0),st(4)
 fsub st(0),st(4)
 ;; Это и есть a – b
 ;; Запоминаем результат, так как это наше z
 fld st(0)
 ;; fld st(0) – это культовая команда,
 ;; которая делает копию вершины стека
 ;; Получаем z**2
 fmul st(0),st(1)
 ;; Следующая особенность в том, что мы смотрим, понадобиться ли
 ;; z**2 ещё где-нибудь далее, и видим, что понадобиться,
 ;; более того, что оно будет складываться само с собой
 fld st(0)
 fadd st(0),st(1)
 ;; и мы быстро получаем конструкцию a+b+z**2
 fadd st(0),st(3) ;; Вот и a Пригодилось…,
 ;; Теперь нам нужно z*b, заметим, что в st(2) у нас z, и как
 ;; ни странно в st(1) b. На самом деле это не совсем случайность
 ;; Кроме того, отметим, что z после следующей операции уже
 ;; не будет нужна.
 fxch st(2)
 fmul st(0),st(1) ;; z*b
 ;; Ну а теперь маленький штрих
 fsubr st(0),st(3)
 ;; Осталось только умножить это на a+b (d)
 fmul st(0),st(5)
 ;; И прибавить a+b+z**2, которое уже есть в st(1)
 fadd st(0),st(2)
 fdiv st(0),st(2)

 ;; УРА!!! ;; Для цельного восприятия выпишем ещё раз без комментариев
 fld a
 fld b
 fst st(4)
 fadd st(0),st(1) ;; a+b
 fst st(3)
 fsub st(0),st(4)
 fsub st(0),st(4)
 fld st(0)
 fmul st(0),st(1)
 fld st(0)
 fadd st(0),st(1)
 fadd st(0),st(3) ;; Вот и a Пригодилось…,
 fxch st(2)
 fmul st(0),st(1) ;; z*b
 fsubr st(0),st(3)
 fmul st(0),st(5)
 fadd st(0),st(2)
 fdiv st(0),st(2)

 ;; Для сравнения
 fld dword ptr [esp+0Ch]
 fadd dword ptr [esp+8]
 fstp st(1)
 fld dword ptr [esp+8]
 fsub dword ptr [esp+0Ch]
 fst st(2)
 fmul st,st(2)
 fstp st(3)
 fld dword ptr [esp+8]
 fadd st,st(3)
 fadd st,st(3)
 fstp st(4)
 fld st(1)
 fmul st,st(3)
 fsubr dword ptr [esp+8]
 fmul st,st(1)
 fadd st,st(4)
 fdiv st,st(4)  

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

А теперь разберём поближе тот или иной этап кодинга данного примера.

В начале, после долгих изъяснений, автор всё-таки решил, что можно поместить параметр b в стек FPU, при этом заметьте, что мы сохранили его значение в st(4). Почему не в st(3), или в st(2), или st(1), а почему не в st(6)? Причина такого поведения автора следующая:

Четвёртое незримое правило FPU кодера:
4. Размещайте данные в FPU регистрах, в зависимости от очерёдности их использования. Менее используемые – ближе к st(7), чаще используемые – поближе к st(0). А также, (по месту их использования в коде) далее используемые – ближе к st(7), а вот-вот используемые ближе к st(0).

Особенно часто используемые хорошо хранить в st(1), st(2), так как их всегда легко превратить в st(0), при помощи fxxp, одновременно выполнив операцию и превратив регистр st(1) в st(0).

Особенность мышления FPU кодера зависит и от психофизиологических особенностей человека:

5.2 Психофизиологические возможности программиста

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

1. Три-четыре команды – абсолютное внимание и распознавание
2. Шесть команд, десять – априорное внимание и осознание
3. До двадцати команд – подсознательное восприятие

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

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

Первая задача программиста – сконцентрировать внимание не на самих операциях, а на концентрировании одинаковых операций, которые не зависят друг от друга. Чтобы упорядочить процесс написания кода, он анализируется и раскрывается сверху вниз. Но при этом с постоянным подсознательным «взглядом» на предыдущие, написанные участки кода. Кодируя текущее выражение, программист должен охватить взглядом три следующих операции, чтобы приблизительно оценить контекст. Как правило, это не составляет особенного труда, учитывая, что основное время тратиться продумыванием текущего преобразования.
Хорошим подходом является итерационный. В этом случае, программист, не задумываясь, кодирует до определённого момента. Предположим 6-10 команд, после чего возвращается и корректирует код. К корректировке кода относиться вставка суффиксов p, r, инструкций FXCH, перестановка команд (см. дальше), и корректировка параметров.

Итак, давайте выделим два основных этапа:

1. Простое кодирование операций в соответствии с их последовательностью.
2. Корректирование типологии операций, их порядка и параметров.

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

Возвращаясь к нашему примеру, после того, как мы решаем поместить a и b в стек, пишем:

 fld a
 fld b
 fadd st(0),st(1) ;; a+b
 fst xx  

Здесь мы останавливается, так как далее должна пойти операция
fsub st(0),st(x), это есть сигнал, что следует остановиться и посмотреть на контекст, скорректировать код.

 fld a
 fld b
 fst st(4)
 fadd st(0),st(1) ;; a+b
 fst st(3)
 fsub st(0),st(4)
 fsub st(0),st(4)  

Красным выделена та часть, которая была скорректирована.

Для закрепления материала рассмотрим пример:

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

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

Другой объект вашего внимания – это максимальное использование контекста. Можно сформировать простое правило таких мест в выражении, где использование контекста будет максимальным:

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

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

Например:

жирным выделены те части выражения, на которые вы обязаны обратить внимания в первую очередь. Весь фокус в том, чтобы зрительно определить «большое автономное скопление операций». По моему часть, что левее «больше», чем та, что правее. Выбираем эту часть . И так далее, пока не доберёмся до наипростейшего выражения. Выбираем, в конце концов, самую элементарную скобку и только тогда начинаем сравнение. Время, затрачиваемое на эту операцию – около 2-4 секунд. Ещё большие выражения дробятся по строкам на листе, и рассматриваются аналогично. Вообще говоря, бессознательная часть способна обработать эту формулу , и найти «похожие части» в среднем за 5-10 секунд в зависимости от тренированности восприятия формул человеком. Взгляд расширяется от данной части выражения к другим, но не далеко (так как всё равно число регистров FPU ограниченно). Таким образом, рассматривается данная часть формулы, и обнаруживается идентичная часть, или части, которые имеют минимальное число операция для перехода друг к другу (имеют минимальное изменение контекста). В данном примере такие части это x, x^2, x^3, а так же…., x, 2x, и (a-b)^2, 2a-b. Теперь программист строит решение данной части, при этом стыки: начало и конец кода, он организует так, посматривая «одним глазом» на соседние выражения, чтобы значения, загруженные в FPU, находились в контексте стыковочных выражений. То есть, начало кода выражения должно содержать результаты выполнения предыдущего выражения и входные данные для следующего. Как вы поняли, этот принцип объясняется стремлением построить код с минимальным изменением контекста.

Кроме этого действует правило: чем чаще переменная встречается, тем больше она претендует быть загруженной в FPU. Поэтому, создавая код (в данном примере), программист считает, что переменные a и b должны быть в FPU, а вот переменная c, не обязательно. Кроме того, заметьте, что идентичными являются так же: 2a-b и a-b. Конечно, очевидным является, что: 2a-b = (a+a-b), правда это не всегда замечают.

Итак, как мы видим: есть два потока расчёта. Первый это с a+a-b, второй с x. Вопрос: «Какой выбрать?». Выбирается тот поток, в котором «завязано больше переменных и вычислений», то есть с x. Итак, мы начинаем писать кодинг выражения x^2+2x-c. Сперва вычисляем x и 2x, вычитаем c, потом +x2. Аналогично: a-b, (a-b)^2, и 2a-b, как a-b+a. И так далее, учитывая предыдущие рекомендации и особенности FPU.

5.2.1 Упражнения
Напишите порядок анализа и вычисления выражений
, .
При выполнении упражнений стоит помнить, что a^2-b^2 = (a+b)(a-b), а (a+b)^2 = a^2+2ab+b^2.

5.3 Восприятие регистров FPU как множества, или FXCH

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

Агнер Фог пишет, что аналогичные операции FPU способны спариваться. Да, способны, но только когда они не зависимы. С этой точки зрения код, написанный нами ранее, абсолютно не учитывал эту особенность. Основная проблема в том, что программист думает, st(0) – это всегда один регистр. Но ведь это не один регистр. Архитекторы Intel может специально, может нет, но, тем не менее, предоставляют возможность обращаться с FPU как с совокупностью регистров. В чём секрет? Секрет в маленькой команде FXCH. Эта команда меняет содержимое ST(0) и ST(x) местами. Однако инструкция реально не обменивает два регистра, она осуществляет аппаратное их переименование. Это означает, что вы можете в любое время превратить любой не st(0) регистр, и ещё при этом не тратить времени, так как инструкция FXCH обладает способностью спариваться с FPU инструкциями U-конвеера. Таким образом, особенности оптимизации должны учитываться вами при проектировании. В подробностях читайте 24, а так же обязательно 26-27 главы культового труда Агнера Фога «How to optimize for the Pentium family of microprocessors», который можно найти в русском переводе Aquilы/HI-TECH.

Мы же отметим, что легкие команды, такие как FADD, FSUB, идущие подряд можно распараллелить, если ввести между ними FXCH. Более тяжёлые команды FMUL, могут быть «разбавлены» двумя «лёгкими» инструкциями FLD, FXCH, как и FDIV. За более подробными рекомендациями обращайтесь в вышеупомянутое руководство.
Особенности процессора влияют на правила описанные выше, добавляя в них свои коррективы. Однако общий принцип сохраняется. Это значит, что при смешении различных независимых потоков операций, приоритет следует отдавать тем потокам, где имеются идентичные операции. Тяжёлые операции писать «разбавляя» служебным кодом типа FLD, либо независимым неFPU кодом.
Рассмотрим пример, деление комплексных чисел:

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

где e,f – действительная и мнимая части результата.
Как видим, такое простое преобразование позволяет увеличить количество подобных частей в выражении, и уменьшить количество операций умножения. Такая оптимизация выражений, когда цель состоит не столько в упрощении выражения, сколько в увеличении количества подобных частей, и уменьшения числа тяжёлых операций, таких как деления, трансцендентных и геометрических функций – носит математический характер, но, тем не менее, связана с осознанием того, как сделать выражение «удобней» для вычислений.
Построим последовательность вычисления:

1. d/c
2. c+d(d/c)
3. (a+b*(d/c))
4. (b-a*(d/c))
5. (a+b*(d/c))/ c+d(d/c)
6. и так далее…

   fld d
    fst st(2)
    fdiv c
    fst st(1)
    ;; Поехали!!!
    fmul b ;; st(0) = b*(d/c)
    fxch st(1) ;; st(0) = d/c
    fmul st(2),st(0) ;; st(3) = d*(d/c)
    fmul a ;; st(0) = a*(d/c)
    fsubr b ;; st(0) = b-a*(d/c)
    fxch st(2) ;; st(0) = d*(d/c)
    fadd c ;; st(0) = c+d*(d/c)
    fxch st(1) ;; st(0) = b*(d/c)
    fadd a ;; st(0) = a+b*(d/c)
    fdiv st(0),st(1)
    fxch st(2) ;; st(0) = b-a*(d/c)
    fdiv st(0),st(1)
    fxch st(1)
    fstp st(0)
    ;; Всё!!! st(0) = мнимой части, st(1) = действительной
   

Формируя этот код, автор пытался сгруппировать одинаковые операции вместе, если это возможно, и, не задумываясь, использовал fxch. Если вы посмотрите на результат, то заметите, что написание кода при помощи fxch похоже на чехарду: нужный результат помещаем в st(x), а ненужный в st(0), который «по случайности» обычно оказывается ненужным потом, и затирается результатом следующей операции. Если после команды fxch, окажется, что вам будет нужно значение st(0), вы можете всегда вернуться и поставить fld st(0) – эта команда делает копию вершины стека.

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

5.4 Элементарные преобразования

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

• Имея вычисленный sin^2(a), легко получить cos^2(a): sin^2(x)+cos^2(x) = 1
• sin(2x) = 2*sin(x)*cos(x)
• cos^2(x/2) = (1+cos(x))/2
• выражения для быстрого умножения
• любая операция деления на константу – это операция умножения на константу в степени -1.

Единственной препоной перед использованием преобразований как альтернативе «чистым» операциям – потеря точности. Но про точность мы поговорим отдельно.

5.5 Операции сравнения

Операция сравнения чисел FPU определённо отличается от операции сравнения, к которой мы привыкли в целочисленном мире IA-32. Это отличие скрывается в большем количестве ситуаций, и возможных результатов сравнения двух чисел. Их можно разделить на 3 группы:

1. Сравнимы.
2. Несравнимы.
3. Логически сравнимы.

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

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

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

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

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

Обычной командой сравнения является команда FCOM x, которая выполняет сравнение st(0) с операндом x, и устанавливает коды условий (таб. № 4.5.1). При этом указатель стека ST не модифицируется, если это не команда FCOMP. Особой командой является FCOMPP, без операндов, которая извлекает из стека сразу два операнда st(0) и st(1). Кроме команд сравнения программисту также доступны несколько команд анализа. Это FTST и FXAM. Команда FTST сравнивает содержимое st(0) с нулём, и устанавливает коды условий (таб.4.5.2). Ещё больше информации о содержимом стека можно получить при помощи команды FXAM (таб. 4.5.3).

Таблица № 4.5.1 Коды условий сравнения (FCOM)
Условие	С0	С3	С2
ST(0) > x	0	0	0
ST(0) < x	1	0	0
ST(0) = x	0	1	0
ST(0) и x несравнимы	1	1	1

 

Таблица № 4.5.2 Коды условий проверки (FTST)
Условие	С0	С3	С2
ST(0) > 0	0	0	0
ST(0) < 0	1	0	0
ST(0) = 0	0	1	0
ST(0) и 0 несравнимы	1	1	1

 

Таблица № 4.5.3 Схема анализа (FXAM)
Если положительное С1 = 0, иначе С1 = 1
Если конечное С0 = 0	Число	С3	С2
	Ненормализованное	0	0
	Нормализованное	0	1
	Нуль	1	0
	Денормализованное	1	1
Если не конечное С0 = 1	Число	С3	С2
	Нечисло	0	0
	Бесконечность	0	1
	Пустой	1	0
	Пустой	1	1

Учитывая тот факт, что FPU разрабатывался как отдельный чип, существует потенциальная проблема в реализации условных переходов. Совершенно ясно, что команды условных переходов схемы x86 не работают с регистрами состояния FPU. Архитекторы Intel вышли из этой проблемы, создав особый «мостик» между регистром FLAGS и регистром состояния SR. Этим особым мостиком является команда FSTSW (FNSTSW), которая сохраняет регистр SR в 16 битной переменной, или регистре ax. А при помощи следующей за ней команды SAHF (которая помещает содержимое ah в FLAGS), мы получаем следующую картину:

Регистр состояния запоминается: FNSTSW ax	B	C3	TOP			C2	C1	C0
ah отображается на регистр FLAGS: SAHF ax	SF	ZF	-	AF	-	PF	-	CF

Как видим, флаг C1 теряется, но как раз он в обыденных сравнениях не участвует. Посмотрим, как выполняется сравнение чисел Фортраном:

9: if(a == b) then
 00401013 fld dword ptr [B]
 00401016 fcomp dword ptr [A]
 00401019 fnstsw ax
 0040101B and ah,40h
 0040101E je Fortran_Proc+1Ch (0040102c)

 9: if(a /= b) then
 00401013 fld dword ptr [B]
 00401016 fcomp dword ptr [A]
 00401019 fnstsw ax
 0040101B and ah,40h
 0040101E jne Fortran_Proc+1Ch (0040102c)

 9: if(a > b) then
 00401013 fld dword ptr [A]
 00401016 fcomp dword ptr [B]
 00401019 fnstsw ax
 0040101B sahf
 0040101C setnp cl
 0040101F setbe ah
 00401022 and cl,ah
 00401024 jne Fortran_Proc+22h (00401032)

 9: if(a < b) then
 00401013 fld dword ptr [A]
 00401016 fcomp dword ptr [B]
 00401019 fnstsw ax
 0040101B sahf
 0040101C setnp cl
 0040101F setb ah
 00401022 and cl,ah
 00401024 je Fortran_Proc+22h (00401032)

 9: if(a <= b) then
 00401013 fld dword ptr [A]
 00401016 fcomp dword ptr [B]
 00401019 fnstsw ax
 0040101B sahf
 0040101C setnp cl
 0040101F setbe ah
 00401022 and cl,ah
 00401024 je Fortran_Proc+22h (00401032)

 9: if(a >= b) then
 00401013 fld dword ptr [A]
 00401016 fcomp dword ptr [B]
 00401019 fnstsw ax
 0040101B sahf
 0040101C setnp cl
 0040101F setb ah
 00401022 and cl,ah
 00401024 jne Fortran_Proc+22h (00401032)
 

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

9: if(a > b) then
 00401013 fld dword ptr [A]
 00401016 fcomp dword ptr [B]
 00401019 fnstsw ax
 0040101B sahf
 00401024 ja Fortran_Proc+22h (00401032)

 ;; Команда ja это ZF = 0, CF = 0, или C3 = 0, C0 = 0

 9: if(a < b) then
 00401013 fld dword ptr [A]
 00401016 fcomp dword ptr [B]
 00401019 fnstsw ax
 0040101B sahf
 jnb @F
 jnc Fortran_Proc+22h (00401032)
 @@:
 ;; Если a >= b!!!
  

Ещё немного фантазии и вы без труда научитесь сравнивать с экзотическими видами чисел, учитывать операнды не-числа и другие случаи.
Хочется лишь напомнить, что, начиная с модели P6, вы можете писать код:

fcomp b
 fnstsw ax
 sahf
 

как

fcomip b  

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

6 Функции и математические абстракции

К элементарным трансцендентным функциям относятся:

1. Тригонометрические (sin, cos, tg…)
2. Обратные тригонометрические (arcsin, arccos, arctg…)
3. Логарифмические функции (log2, lg, ln)
4. Показательные функции (x^y, 2^y, 10^x, e^x)
5. Гиперболические функции (sh, ch, th, …)
6. Обратные гиперболические (arcsh, arch, arcth…)

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

6.1 Трансцендентные функции

Инструкция	Функция	Описание
FSIN	Sin(st(0))	Синус числа, находящегося в st(0). Операнд считается заданным в радианах, и имеет пределы 2^63 и -2^63. Если операнд выходит за эти пределы, флаг C2 = 1, а st(0) не изменяется.
FCOS	Cos(st(0))	Косинус числа. (Замечания аналогичны)
FSINCOS	st(1) = Sin(st(0)); st(0) = cos(st(0))	Косинус и синус числа (в одном флаконе )
FPTAN	Tan(st(0))	Тангенс числа в st(0). Вычисляет тангенс. Результат помещает в st(0), после чего помещает 1 в стек, так что: St(0) = 1, st(1) = tan(st(0)).
FPATAN	Atan(st(1)/st(0))	Вычисляет арктангенс числа st(1)/st(0). Результат записывается в st(1), а st(0) выталкивается. Результат имеет знак st(1), и меньше числа PI по абсолютной величине.
F2XM1	2^st(0)-1	Значение st(0) лежит интервале (-1,1)
FYL2X	y*log2(x)	Вычисляет st(1)*log2(st(0)). Результат размещается в st(1), а st(0) выталкивается из стека. При этом st(0)>0, если st(0) = 0, тогда результат (ZM = 1) равен бесконечности со знаком обратным st(1).
FYL2XP1	y*log2(x+1)	Вычисляет st(1)*log2(st(0)+1). St(1) равен результату, после чего st(0) выталкивается из стека. При этом первоначальный st(0) принадлежит (-a,a), где a = 1-1/SQRT(2). Наличие этой команды обусловлено тем, что она даёт большую точность, чем FYL2X для суммы st(0)+1 для чисел близких к нулю.

Вы, наверное, заметили, что все функции имеют свой особый диапазон вычислений. Для тригонометрических – это (2^64,-2^64). И естественно ставится закономерный вопрос о том, что всё это надо как-то учитывать. Хорошо, если переменная REAL(4), тогда, например можно не беспокоится о выходе за пределы. Однако во всех других случаях, когда такой уверенности нет, придётся использовать привидение значения st(0) к необходимому.

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

На то есть много хитростей, которые зависят от задачи. Например, если вы рассматриваете функцию y = f(x), где x – это угол, который «наматывается» чем-либо, для вычисления функции вы можете хранить число «намоток» и «чистый» угол без их учёта. Таким образом, для тригонометрических функций у вас есть гарантированная величина угла, а x – это всего лишь произведения «чистого» угла и количества «намоток».

Подобный приём раздельного учёта угла и количества полных кругов возможен в циклических задачах, однако, тем не менее, не употребляется в других случаях. Ну что ж, всё в ваших руках.
Другая особенность тригонометрических функций – это использования или градусов или радиан. Практически в большинстве задач, использование градусов обусловлено только удобством для человека. Поэтому чтобы убрать эту проблему, достаточно преобразовать при компиляции радианы в градусы. Однако всё-таки есть много случаев, когда содержимое переменной угол используется как в градусах, так в радианах. Здесь стоит думать, как поступить: хранить переменную в радианах, а переводить в градусы – или наоборот? Что легче? Это зависит от того, что чаще.

Из всего набора инструкций, программисты чаще всего удивляются особенностям инструкций показательных и логарифмических функций. Совершенно непонятно зачем понадобилась такая функция 2st(0)-1 в виде инструкции F2XM1. Да ещё с такими жёсткими ограничениями для st(0) (-1,1). Ответ на первый вопрос кроется в повышении точности вычислений. Об этом мы ещё более подробно поговорим в разделе 7. А пока скажем лишь, что эта команда позволяет избежать большой потери точности, когда x близок к нулю. Ответ на второй вопрос кроется в разделе 6.2.

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

6.2 Получения основных функций

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

6.2.1 Обратные тригонометрические функции

Здесь везде sign(x) есть функция, возвращающая знак аргумента x.

6.2.2 Логарифмические функции

Log₂(x) = FYL2X(x);
Log_e(x) = log_e(2)*log₂(x) = FYL2X( log_e(2), x ) = FYL2X ( FLDLN₂, x )
Log₁₀(x) = log₁₀(2)*log₂(x) = FYL2X( log₁₀(2), x ) = FYL2X ( FLDLG₂, x )

6.2.3 Формулы для гиперболических и обратных гиперболических функций

; ;

Arsh(x) = [ sign(x) ]*[ log_e(2) ]*log₂(1+z) = FYL2XP1( sign(x)*FLDLN2,z )

Где

Arch(x) = [ log_e(2) ]*log₂(1+z) = FYL2XP1( FLDLN2,z )

Где

Arcth(x) = [ sign(x) ]*[ log_e(2) ]*log₂(1+z) = FYL2XP1(sign(x) * FLDLN2 ,z )

Где

6.2.4 Формулы для показательных функций

2^x = (2^x-1)+1 = F2XM1(x)+1

= 1+F2XM1(x*log₂(e)) = 1+ F2XM1(x*FLDL2E)

= 1+F2XM1(x*log₂(10)) = 1+ F2XM1(x*FLDL2T)

= 1+F2XM1(y*log₂(x)) = 1+ F2XM1(FYL2X(y,x))

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

7 Вопросы точности и округления

После выхода первой части данной статьи, один из читателей написал в редакцию замечание, что точности и так достаточно, и что если это не так – значит вы плохой программист. Точно? Достаточно? Инженеры Интел (один из них) мотивировали размер мантиссы числа REAL(10) как достаточного для соблюдения точности при операции возведения в степень yx. Об этом было написано в первой части. Точность!!! Так как точность зависит именно от размера мантиссы, а не размера порядка. Нужно сказать, что на адрес редакции пришёл и такой комментарий: «Какая разница между числами с фиксированной и плавающей запятой…». Разница принципиальная, так как в первом случае количество бит на целую и дробную часть фиксированы, а в другом случае число не делиться на части, и представляет собой мантиссу. Как мы видим, числа с фиксированной запятой не подвержены тем всем проблемам, которые будут описаны далее. Но для начала, посмотрим на мир глазами вычислительной математики.

7.1 Основы вычислительной математики

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

Относительная погрешность – это отношение абсолютной погрешности к приближённому значению числа:

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

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

Правила вычисления погрешностей при операциях:

1. При сложении/вычитании чисел относительная погрешность результата равна максимальной относительной погрешности из двух.
2. При умножении и делении чисел относительная погрешность равна сумме относительных погрешностей составляющих.
3. При возведении в степень относительная погрешность есть произведение погрешностей показателя и основания.

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

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

0.xxxx*10ⁿ

Сравните это определение в вычислительной математики с определением нормализованного числа в FPU. Принципиально они похожи. Рассмотрим примеры:

1234 = 0.1234*10⁴, 12.23 = 0.1223*10².

Таким образом, потеря точности, уже зависит не только от погрешности вычислений FPU, округления результатов, она напрямую зависит от порядка взаимодействующих чисел. Если мы говорим, что нам нужна точность до 3-его знака после запятой, то замете, что число 0.1234*10⁴, вообще не содержит чисел после запятой. Это значит, что любые операции с таким числом приведут к потери точности на несколько порядков превышающую заданную. Классическим примером, который демонстрирует данную ситуацию, является пример вычитания чисел близких по значению.
0.1234*10⁴ - 0.1233*10⁴ = 0.1*10¹
Даже при том условии, что эти числа получены с точностью до 3-его знака после запятой, их разность имеет точность до 1. Более подробный пример будет рассмотрен далее.
Здесь и далее будет действовать один единый принцип потери точности. Хотя он только подразумевается, его осознание важнее знания самих случаев.

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

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

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

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

7.2 Основные принципы ситуаций потери точности

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

1. Вычитание близких по значению чисел
2. Операции сложения/вычитания большого и маленького числа
3. Антипереполнение
4. Переполнение

Как уже было сказано, потеря точности обусловлена «съеданием» младших разрядов меньшего числа, большим числом. Особенно интересно то, что такому недугу не подвержена операция умножения/деления.

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

Всё вышесказанное заставляет обратить внимание на операции сложения/вычитания.
Нужно, например, отметить, что сам по себе первый случай вычитания близких по значению чисел не приносит урона точности. Его урон косвенный. В одном контексте – он велик, в другом – не больше погрешности вычислений. Так в каком случае? Например, если вы рассматриваете x и y, как числа с точностью до 3-его знака, а диапазон этих чисел велик, и происходит вычитание y и x имеющих большие значения, полученные в результате определённых преобразований, то результат может быть трактован как «точный» (в неком смысле) только в том случае, если далее он нигде не используется. Иначе мы имеем проблемы.

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

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

В случае переполнения речь идёт о превышении границы отображения числа в FPU. Читатель заметит: «Ну, раз переполнение, что ж тогда….». Дело в том, что многие математические выражения, результат которых оказывается в пределе пространства FPU чисел, имеют промежуточные результаты, которые выходят за этот интервал.

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

Если мы вспомним ГЛАВНЫЙ принцип потери точности, то поймём, что любая операция, повышающая порядок числа, вызывает потерю точности. И эта потеря будет тем больше, чем больше изменения порядков.

Совершенно ясно, что такую ситуацию следует избегать как:

или

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

А теперь проанализируем наши размышления и выявим один единственный принцип по борьбе с потерей точности:

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

Будем называть такой диапазон порядка – безопасным диапазоном порядка.
Безопасный диапазон порядка зависит от:

1. Размера мантиссы
2. Целевой точности, которая должна быть обеспеченна

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

Ну, как? Теперь Вам уже не кажется, что формат расширенной точности столь необъятный? Автор напугал вас? Думаю, что бояться не стоит:

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

И всё-таки следует быть осторожным, и помнить, что именно «съедание» бит приводит к необъяснимым результатам в одной и той же задаче, которая вроде бы написана правильно, но считает «когда как».

7.3 Методы и способы увеличения точности

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

Пусть b >> 4ac, тогда SQRT(D) по модулю приближённо равно b. Таким образом, мы можем иметь случай, когда –b+b = 0, значит x1 = 0, что не есть решение задачи. Преобразуем выражение для x1:

Где было сложно? Правда? Зато проблема решена.

7.3.2 Порядок сложения/вычитания чисел

Несмотря на то, что все мы учим с первого класса: «при смене мест слагаемых сумма не изменяется», теперь пришло время опровергнуть это утверждение.
Пусть мы имеем 4-значное представление числа в машине (см. выше). Пусть мы имеем такую сумму:
S = 0.2764 + 0.3944 + 1.475 + 26.46 + 1364.
Теперь просчитайте эту сумму первый раз. При этом, считая, что у вас есть только четыре разряда, и производя округление результатов.
И второй раз, но задом наперёд. Ну что? Помните как фокусе 2 = 0? .
Вывод: При операциях сложении чисел одного знака суммирование следует вести в порядке возрастания их величин.

Пусть мы имеем такое выражение: x–x+y. Пусть мы имеем дело с числами одинарной точности (max(10^(37))). Совершенно ясно, что это выражение равно y. Если x = 10^30, а y = 10^(-30), то в случае:
1. x – x + y = y; 10^30 – 10^3 + 10^(-30) = 10^(-30)
2. x – (x + y) = 0; 10^30 – (10^30 + 10^(-30)) = 10^30 – 10^30 = 0
Как видно этот пример так же попадает под наш общий принцип.

7.3.3 Случай переполнения (антипереполнения)

Очень сложно вообразить, что в технических расчётах могут встретиться такие числа как 10^4932 (для чисел расширенной точности), или даже 10^308 (для двойной). И, тем не менее, стоит помнить о случае переполнения. Пример.

Пусть нужно найти длину вектора

если имеется хоть один x_i , который близок к нулю, но ещё может быть представлен в формате нормализованного числа, то, будучи возведенным, в квадрат такое число, может «превратиться» в нуль. В крайнем случае, при всех x_i близких к нулю мы можем получить, что |x| = 0.

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

7.3.4 Случай близкого результата

Пусть у нам есть функция y(x) = B, где B – это какое-то выражение.

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

1. Масштабирование – умножить выражение B на масштабный коэффициент.
2. Написать логарифмическое уравнение.
   Log(1-y(x)) = log(1-B).

Замете особенность – 1-x!!! Что увеличивает точность вычислений. Теперь это должно быть ясным для вас, как и то, почему некоторые инструкции FPU выдают результат именно как 1-x!!!
Приятно осознавать, что нет ничего такого, что было бы случайностью.

8 Функции, использующие FPU

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

Функции, используемые FPU можно разложить на следующие виды:

1. Подпрограммы
2. Мини-функции
3. Рабочие функции-преобразователи
4. Рабочие функции
5. Контекстно-зависимые функции

По виду работы:

1. С конечным числом аргументов
2. Со структурами
3. С массивами
4. Контекстные

Немного по каждому виду:

• Подпрограммы предполагают большой участок кода, использующий FPU, поэтому имеет смысл выполнять передачу аргументов в стеке, и лишь при небольшом их числе в самом FPU. Как правило, подпрограмма, возвращающая одно значение, помещает ответ в st(0), предварительно очистив все регистры FPU.
• Мини-функции, как правило, принимают аргумент в стеке FPU st(0), и возвращают результат в нём же. При этом содержание FPU может меняться. (Функции затирают st(1) и st(2) и не освобождают их)
• Рабочие функции-преобразователи так же принимают, как правило, один аргумент, но часто в стеке, а результат помещается в st(0). Состояние FPU при этом сохраняется (кроме st(0))
• Рабочие функции ведут себя подобно функциям преобразователям, очищают FPU от своих промежуточных результатов, но принимают аргументы не только в стеке.
• Контекстно-зависимые функции (как и многие другие не возможные в ЯВУ), работаю с st(0) и текущим состоянием регистров FPU. Ничего не очищаю, и не заботятся о сохранении значений стека FPU и его состояния.
• Различные виды функций используются каждая для своих характерных задач. Способ вызова и работы функции выбирается в зависимости от:
  1. Частоты вызова функции
  2. Функциональной завершённости (вспомогательная или совершенная функция)
  3. Метода представления результатов и аргументов

 

• Совершенной называется такая функция, которая реализует независимый, самостоятельный участок работы. (Например, root2(a,b,c) – нахождения корней квадратного уравнения).
• Вспомогательной называется функция, которая может использоваться только в определённом контексте, и для определённой системы кода.

(Мы ещё вернёмся к Совершенным и Вспомогательным функциям в следующих выпусках)

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

9 Заключение

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

Автор надеться, что ему удалось не просто поговорить про FPU, а и предложить определённый тип мышления, сконцентрировать внимание на «истинно важных вещах», не тратить время на несущественные детали.

Я уверен, что умение «мыслить нужным образом» ценнее всех знаний мира, вместе взятых.
Если у вас возникли трудности, которых вы не ожидали, пишите мне EdmondXASM@mail.ru, и мы встретимся снова…. А пока… пока будем ждать продолжения.

© Edmond / HI-TECH