[Перевод] Реверс-инжиниринг алгоритма тангенса в процессоре Intel 8087

Это ещё одна моя статья о чипе арифметики с плавающей запятой Intel 8087¹. В 1980 году компания Intel представила процессор 8087, который сильно повысил скорость выполнения операций с плавающей запятой в IBM PC и других системах. В этой статье мы рассмотрим алгоритм, лежащий в основе команды вычисления тангенса. Один из распространённых подходов к реализации тригонометрических функций — это алгоритм CORDIC. Второй — это аппроксимация при помощи многочленов. Для обеспечения высокой точности и высокой производительности 8087 сочетал в себе оба этих подхода, создав огромный прирост скорости по сравнению с микропроцессором 8086: тангенс в нём вычислялся всего за 90 микросекунд вместо 13 тысяч микросекунд².

Изучив цепи и микрокод 8087, я смог объяснить алгоритм команды вычисления тангенса FPTAN. Для исследования электросхемы 8087 я при помощи стамески снял с чипа крышку и сделал микроскопом снимок высокого разрешения. ROM микрокода — это большая прямоугольная область в центре кристалла, содержащая 1 648 управляющих чипом микрокоманд. Нижняя часть чипа (в красном прямоугольнике) — это тракт данных: цепь, выполняющая вычисления с плавающей запятой над 80-битными значениями³.

На увеличенном снимке тракта данных видны соответствующие функциональные блоки. В exponent ROM хранятся значения экспонент, необходимые алгоритму. В constant ROM содержатся константы, в том числе и используемые алгоритмом CORDIC. Shifter (блок сдвига) — это большой компонент, выполняющий сдвиг 64-битного значения влево или вправо на произвольную величину. Adder (сумматор) — это ядро вычислений 8087; наряду с реализацией сложения и вычитания он применяется в цикле для умножения, деления и квадратных корней. В регистре B (B reg) хранится один из аргументов сумматора, а второй аргумент может браться из нескольких источников. В sum register хранится вывод сумматора. В восьми стековых регистрах и временных регистрах хранятся числа с плавающей запятой. Наконец, в сдвиговом регистре (shift reg) хранятся 16 бит состояния для вычислений CORDIC.

Алгоритм CORDIC

CORDIC — это хитрый алгоритм быстрого вычисления трансцендентных функций на простом оборудовании: в нём применяются команды сдвига и сложения, а также поиск по таблицам, но не используются умножение и деление. Этот алгоритм был разработан в 1956 году для B-58 Hustler, который стал первым бомбардировщиком, способным летать на скорости 2 маха. На борту самолёта имелся аналоговый навигационный компьютер, однако точность его компонентов была ограниченной. Инженеру Джек Волдеру поручили задачу проектирования цифрового компьютера на замену аналоговому⁷. Важнейшей проблемой при этом стало то, что аналоговый компьютер мог с лёгкостью генерировать синусы и косинусы при помощи электромеханического ресолвера, а в цифровом виде тригонометрические функции вычислять сложно, особенно учитывая низкую скорость транзисторов той эпохи.

Джек Волдер придумал быстрый способ вычисления тригонометрических функций на простом оборудовании. Он назвал этот алгоритм (а также реализующий его компьютер) CORDIC: «COordinate Rotation DIgital Computer». CORDIC преобразует угол в вектор, координаты которого соответствуют необходимым тригонометрическим функциям. Хитрость заключается в том, чтобы разбить угол на последовательность специальных углов, упрощающих вращение вектора. Эти специальные углы вычисляются заранее и сохраняются в таблицу, благодаря чему вычисления CORDIC можно выполнять быстро даже на оборудовании 1950-х. Каждая итерация CORDIC добавляет один бит точности, поэтому алгоритм быстро сходится. CORDIC получил популярность, в том числе и в научных калькуляторах, где вместо двоичного применялся десятичный CORDIC.

Немного тригонометрии

Я постараюсь не углубляться в математику, но в этом разделе мне всё-таки придётся вкратце объяснить, как работает CORDIC. На графике ниже показано, как тригонометрические функции соотносятся с координатами точки. Допустим, у нас есть угол θ, задающий точку (X, Y) на единичной окружности. Его базовые формулы: X=cos θ, Y=sin θ и Y/X = tan θ. Следовательно, если мы сможем определить координату (X, Y), то определим и значение тригонометрических функций. Если точка находится не на единичной окружности, например, (X', Y'), то мы всё равно с лёгкостью можем определить tan θ. (Спойлер: именно так и поступает 8087.) Однако с sin θ и cos θ всё становится сложнее⁴.

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

К сожалению, поскольку матрица вращения (см. 1 ниже) требует sin и cos, то она, похоже, не решает нашу проблему. Однако можно поделить матрицу на cos θ; кажется, что это ещё бесполезнее, ведь теперь матрице (2) требуется tan, который мы и хотим вычислить. (Кроме того, также увеличилась длина вектора.) Но главное в CORDIC заключается в использовании специальных углов α_n = arctan(2^‑n). Подставив один из таких специальных углов в матрицу, мы получим матрицу (3), которую легко вычислить аппаратно: умножение на степень двойки можно выполнить сдвигом бит.

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

Теперь мы можем начать разбираться, как работает CORDIC. Сначала мы разбиваем нужный нам входной угол на сочетание специальных углов, сумма которых равна нужному углу⁵. Затем применяем приведённую выше формулу вращения для каждого специального угла, начиная с единичного вектора (1, 0). Результатом будет точка (X, Y) в нужном угле, поэтому нужный тангенс просто будет равен Y/X⁶. Так как таблица специальных углов вычислена заранее, операции arctan не замедляют процесс. Стоит также отметить, что после первых членов специальные углы приближаются к 2^‑n, то есть на каждом шаге уменьшаются примерно вдвое.

Подведём итог: алгоритм CORDIC состоит в циклическом обходе таблицы сохранённых углов. Если сохранённый угол меньше нужного нам, сохранённый угол вычитается из нужного угла, давая новый нужный угол, а для получения нового вектора применяются приведённые выше уравнения (сдвиги, сложение и вычитание). В конце тангенс исходного угла определяется как Y/X.

Аппроксимация рациональным многочленом

Точность CORDIC зависит от количества используемых членов. При 16 членах точность примерно равна 2^-16, или 16 битам точности. Для получения 64 бит точности придётся вычислить 64 членов (и иметь таблицу 64 специальных углов). Чтобы получить ответ быстрее, 8087 использует 16 бит алгоритма CORDIC, а для оставшегося угла применяет другой алгоритм. (Оставшийся угол — это разность между суммой специальных углов CORDIC и исходного угла, поэтому он очень мал, примерно 2^-16.)

Для оставшегося угла 8087 использует аппроксимацию Паде, то есть соотношение двух многочленов. Существует целое семейство аппроксимаций Паде в зависимости от порядка многочленов. В 8087 применяется простая формула: 3 x/(3-x²)⁸. Несмотря на простоту этой аппроксимации, при малых значениях она очень точна; её погрешность пропорциональна x⁴. Так как x<2^-16, погрешность будет меньше 2^-64, что соответствует требуемой 64-битной точности 8087. Более того, 8087 не обязан выполнять деление в рациональном многочлене, поскольку FPTAN возвращает отдельные числитель и знаменатель. Таким образом, деление становится «бесплатным».

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

Соединяем всё вместе: алгоритм процессора 8087

Алгоритм вычисления тангенса 8087 состоит из трёх частей: определение битов выбора направления итераций CORDIC (называемое псевдоделением), вычисление рациональной аппроксимации и применение уравнений вращения в зависимости от битов выбора направления (называемое псевдоумножением)⁹.

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

На схеме ниже показан это процесс на примере исходного угла 0,95 радиан¹⁰. Процесс генерирует последовательность битов [1,0,0,1,0,1,0,1,0,0,1,0,0,1,1,1], где крайний левый бит обозначает arctan(2⁰) и так далее. На каждом шаге угол уменьшается примерно вдвое, поэтому остаточный угол очень мал.

Далее при помощи функции рациональной аппроксимации 3 x/(3-x²) вычисляется тангенс оставшегося угла. Результат используется в качестве начального вектора для следующего этапа. Деление здесь не выполняется; вместо этого числитель становится Y в начальном векторе, а знаменатель становится X. Таким образом мы избегаем затратного этапа деления, поскольку оно в неявной форме выполняется в ответе. Умножение на 3 выполняется легко (сдвиг влево и сложение), поэтому единственная затратная операция на этом этапе — возведение угла в квадрат, для чего требуется полное 64-битное умножение.

На последнем этапе применяется каждое вращение CORDIC, если соответствующий бит выбора направления из первого этапа был равен 1. Так как это чем‑то близко к двоичному умножению, которое прибавляет множимое на каждом этапе, если бит множителя равен 1, то это этап называют псевдоумножением. Как говорилось выше, каждое вращение вычисляется при помощи сдвигов, сложения и вычитания, поэтому вращение малозатратно. Каждое вращение на этом этапе соответствует уменьшению угла на первом этапе. Для снижения погрешности округления вращения применяются в обратном порядке, начиная с наименьшего. Для этого биты выбора направления сохраняются на первом этапе в 16-битный сдвиговый регистр, который на этом этапе сдвигается в обратном порядке.

На схеме выше показан процесс на примере исходного угла 0,95 радиан. Исходный вектор (зелёный) берётся из функции рациональной аппроксимации; он очень близок к (3, 0), но слегка повёрнут из-за крошечного остаточного угла. Каждый этап вращения генерирует новый вектор, соответствующий углу из соответствующего шага первой части; я показал две матрицы вращения. Обратите внимание, что векторы отклоняются от окружности: это следствие деления каждой матрицы на cos θ, но на тангенс это не влияет.

Команда FPTAN довольно любопытна тем, что не возвращает непосредственно тангенс. На самом деле, она возвращает частичный тангенс: два значения координат (X и Y), которые можно поделить одно на другое, чтобы получить тангенс. Можно задаться вопросом, почему чип не выполнял деление автоматически. Смысл в том, что в те времена деление выполнялось медленно, и отказ от него в некоторых случаях допускал возможности оптимизаций.

Аппаратные подробности о 8087

В этом разделе я объясню некоторые особенности 8087, важные для реализации FPTAN в микрокоде.

8087 поддерживает широкий спектр типов данных, но внутри него всё хранится в виде 80-битного числа с плавающей запятой, называемого «temporary real». Число состоит из трёх частей: бита знака, 15-битной экспоненты и 64-битной мантиссы (дробной части). В большинстве случаев число с плавающей запятой представляется в виде знак × мантисса × 2^{экспонента}. Преимущество чисел с плавающей запятой в том, что экспонента позволяет им охватывать широкий диапазон значений, от очень маленьких до очень больших. Мантисса — это 64-битное двоичное число вида 1.bbb…: единица в начале, за которой следует двоичная запятая (двоичный эквивалент десятичной запятой) и остальные биты¹¹. Важно то, что экспонента хранится с прибавленным к ней «перекосом» 16 383. Благодаря этому хранимая экспонента всегда положительна, даже если реальная экспонента отрицательна. Например, экспонента 1 хранится в виде 0x4000, а экспонента -16 — в виде 0x3fef.

Для работы с 8087 программист хранит значения в восьми внутренних регистрах, упорядоченных в стек. В каждом регистре записано 80-битное число с плавающей запятой. Для оптимизации производительности с каждым значением в стеке регистров связано значение-«метка»: zero, valid, special или empty. Значение 0 помечается как zero. «Обычное» значение с плавающей запятой помечается как valid. Если значение представляет собой бесконечность, Not a Number (NaN) или денормализованное значение, то оно помечается как special. Метка empty означает, что в регистре не хранится значение; это позволяет распознавать выход за нижнюю границу стека (чтение пустого регистра) или переполнение стека (сохранение в непустой регистр).

Также у 8087 есть временные регистры для внутреннего использования: tmpA, tmpB и tmpC. (Эти регистры активно задействуются при вычислении тангенсов.) tmpA и tmpB — это 80-битные регистры с двумя битами меток. tmpC отличается от них тем, что не имеет знака, экспоненты и метки. Кроме того, его мантисса равна 68 битам; tmpC имеет дополнительный бит для предотвращения переполнения и три дополнительных бита низкого порядка для округления: guard, round и sticky.

В микрокоде каждая 16-битная микрокоманда выполняет одну низкоуровневую операцию. Многие операции копируют данные из одного регистра в другой. Другие обеспечивают условные переходы, вызовы подпрограмм и возвраты. Операция битового сдвига занимает две микрокоманды: первая копирует значение в блок сдвига и конфигурирует направление и величину сдвига. Вторая микрокоманда (которая может находиться на некотором расстоянии от первой) копирует результат из блока сдвига в целевой регистр.

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

Реализация алгоритма вычисления тангенса в микрокоде

В этом разделе мы поговорим о важных моментах микрокода вычисления тангенса. Ниже показан микрокод 8087 и мои комментарии. Каждая строка обозначает 16-битную микрокоманду. Поток управления довольно сложен, поэтому я добавил в примечания блок-схему¹².

FPTAN:
#1039 st(0) -> tmpA        ;сохраняем аргумент в tmpA
#1040 stackPtr--           ;проверяем, есть ли место в стеке для добавления результата
#1041 stack overflow?
#1042 jmp #1022 if tmp empty/special/overflow/div ;выход при плохом аргументе или переполнении стека
#1043 jmp #1045 if tmpA:tag ZERO ;равен ли аргумент 0?
#1044 except:precision     ;исключение precision для всего, кроме 0
#1045 expconst 0x3ff0      ;проверяем экспоненту на >= -16
#1046 adder: tmpA:exp + 0 cin=1 ;экспонента аргумента + 1
#1047 expConst -> Breg
#1048 adder: sumreg:frac - Breg cin=1 ;вычитаем -15
#1049 jmp #1061 if adder pos ;CORDIC, если экспонента >= -16
#1050 expconst 0x3fff      ;путь без CORDIC:
#1051 tmpA:exp -> Breg     ;сравниваем исходную экспоненту
#1052 expConst -> tmpB:frac ;с экспонентой 0
#1053 adder: tmpB:frac - Breg cin=1 ;вычитаем
#1054 sumreg:frac -> expConv ;3fff - exp (то есть -exp без перекоса)
#1055 sumreg:frac -> shiftcount
#1056 jmp #1082 if exponent[6:14] != 0 ;если exp > -64, переходим к рациональной аппроксимации
#1057 high bit -> tmpB:frac ;в противном случае возвращаем исходный аргумент
#1058 expConst -> tmpB:sign,exp ;push 1 для X (делитель): результат - исходный угол
#1059 tmpB -> st(0)
#1060 RNI
#1061 expconst 0x0010      ;путь CORDIC
#1062 sumreg:frac -> loopcounter ;loopcounter = exp + 16
#1063 tmpA:frac -> tmpC    ;tmpC = исходный угол
#1064 trig const -> Breg/rnd
#1065 adder: tmpC - Breg cin=1, roundmode ;вычитаем первый угол CORDIC
#1066 adder sign -> cordic, shl ;сохраняем бит в сдвиговом регистре
#1067 jmp #1079 if not adder pos
#1068 sumreg:frac,rnd -> tmpC
#1069 adder: tmpA:exp + 0 cin=1 ;инкремент exp, если использован первый угол
#1070 sumreg:frac -> tmpA:exp
#1071 jmp #1079
#1072 shift tmpC L 0 bytes, 1 bits ;вершина цикла сканирования CORDIC
#1073 shift L -> tmpC      ;сдвиг угла влево
#1074 trig const -> Breg/rnd
#1075 adder: tmpC - Breg cin=1, roundmode ;вычитаем угол CORDIC
#1076 adder sign -> cordic, shl ;сохраняем бит в сдвиговый регистр
#1077 jmp #1079 if not adder pos
#1078 sumreg:frac,rnd -> tmpC
#1079 jmp #1072 if not const latch zero ;конец цикла сканирования CORDIC
#1080 tmpC -> tmpA:frac    ;обновляем tmpA, tmpB
#1081 expConst -> shiftcount ;16 -> shiftcount (в противном случае зависит от экспоненты)
#1082 tmpA:frac -> tmpB:frac ;аппроксимация Паде
#1083 tmpA:frac -> Breg    ;tmpA = ang
#1084 call SQUARE          ;возводим угол в квадрат
#1085 shift sumreg:frac R count byte bit
#1086 shift R -> sumreg:frac ;двойной сдвиг для масштабирования квадрата
#1087 shift sumreg:frac R count byte bit
#1088 shift R -> Breg;      ;Breg = ang^2
#1089 1.1 -> tmpB:frac     ;3 (с exp 1)
#1090 adder: tmpB:frac - Breg cin=0
#1091 sumreg:frac -> tmpB:frac ;tmpB (X) = 3-ang^2
#1092 shift tmpA:frac R 0 bytes, 1 bits
#1093 shift R -> Breg/rnd
#1094 adder: tmpA:frac + Breg cin=0, roundmode ;ang + ang/2
#1095 sumreg:sign,frac,rnd -> tmpC ;tmpC (Y) = 3×ang с соответствующей экспонентой
#1096 expconst 0x3ff0      ;экспонента -15
#1097 expConst -> Breg
#1098 adder: tmpA:exp - Breg cin=1
#1099 jmp #1118 if not adder pos ;пропускаем, если экспонента слишком мала
#1100 shift tmpC R 0 bytes, 1 bits ;путь псевдоумножения CORDIC
#1101 #15 -> loopcounter
#1102 shift R -> tmpC
#1103 jmp #1113 if not cordic[0] ;CORDIC: если сохранён бит решения...
#1104 shift tmpC R loopcount ;отправить tmpC (Y) в блок сдвига
#1105 tmpC -> Breg/rnd     ;отправить tmpC в сумматор
#1106 adder: tmpB:frac + Breg cin=0, roundmode ;tmpB (X) не сдвигается, потому что tmpC масштабирован
#1107 sumreg:sign,frac,rnd -> tmpC ;tmpC = tmpC + tmpB (косвенно >>n)
#1108 shift R -> sumreg:frac,rnd ;сдвинули старый tmpC в sumreg
#1109 shift sumreg:frac,rnd R loopcount ;снова сдвиг вправо, потому что tmpC масштабирован
#1110 shift R -> Breg/rnd  ;и отправляем в регистр B
#1111 adder: tmpB:frac - Breg cin=1, roundmode
#1112 sumreg:frac,rnd -> tmpB:frac ;tmpB = tmpB - tmpC>>n
#1113 cordic shr           ;сдвиг наружу бита сдвигового регистра
#1114 shift tmpC R 0 bytes, 1 bits ;начало сдвига tmpC (Y)
#1115 jmp #1118 if cordic==0 ;завершено, если сдвиговый регистр пуст
#1116 shift R -> tmpC      ;сдвиг вправо tmpC (Y)
#1117 jmp #1103 if not const latch zero ;завершено, если счётчик на нуле
#1118 expconst 0x3fff      ;CORDIC завершён
#1119 tmpB:frac -> sumreg:frac ;проверка верхнего бита tmpB (X)
#1120 jmp #1124 if reg bit 63 ;если он задан, используем экспоненту 0
#1121 expconst 0x3ffe      ;в противном случае используем экспоненту -1 и
#1122 shift sumreg:frac L 0 bytes, 1 bits ;сдвигаем влево на один бит для нормализации
#1123 shift L -> tmpB:frac
#1124 expConst -> tmpB:sign,exp ;сохраняем экспоненту в tmpB
#1125 tmpC -> tmpA:frac    ;сохраняем tmpC (Y) в tmpA
#1126 tmpC -> sumreg:frac,sign
#1127 jmp #1132 if not sumreg[64] ;проверяем бит переполнения
#1128 shift tmpC R 0 bytes, 1 bits ;если он задан, нормализуем сдвигом вправо
#1129 shift R -> tmpA:frac
#1130 adder: tmpA:exp + 0 cin=1 ;и выполняем инкремент экспоненты на 1
#1131 sumreg:frac -> tmpA:exp
#1132 tmpB -> st(0)        ;tmpB в вершину стека: X (знаменатель)
#1133 stackPtr++
#1134 tmpA -> st(0)        ;tmpA в st(1): Y (числитель)
#1135 stackPtr--
#1136 RNI                  Конец FPTAN

Микрокод FPTAN начинается по адресу #1039 (десятичному). Сначала микрокод копирует аргумент из вершины стека в регистр tmpA. Если верхнее значение пустое, это означает выход за нижнюю границу стека. Если следующий элемент в стеке (позиция, в которой содержится результат) не пустая, это обозначает переполнение стека. В обоих случаях микрокод сообщает об исключении invalid и завершает команду.

Любопытная особенность 8087 заключается в том, что он помечает результат исключением precision, если результат неточный. (Это происходит очень часто, скажем, в одном из десяти случаев.) Единственный тангенс, который 8087 может вычислить точно — это tan(0), поэтому все иные входящие данные приводят к исключению precision (#1044).

Затем проверяется экспонента аргумента, что разбивает исполнение кода на три пути. Если экспонента равна -64 или меньше, аргумент настолько мал, что тангенс равен аргументу в пределах уровня точности системы¹³. В данном случае (#1057) исходный аргумент возвращается без изменений (с делителем 1), а код завершает исполнение. Второй случай — когда экспонента аргумента равна -17 или меньше. В этом случае этап CORDIC пропускается, а процедура переходит непосредственно к рациональной аппроксимации (#1082).

Псевдоделение CORDIC

Самый любопытный случай — это путь с CORDIC (#1061), по которому исполнение идёт при экспоненте от -1 до -16. Концептуально этот код выполняет 16 шагов CORDIC, используя 16 сохранённых углов. Однако код оптимизирован и пропускает большие углы в случае малых входящих данных. В частности, счётчик цикла инициализируется в зависимости экспоненты исходного угла (#1062). Цикл псевдоделения CORDIC (#1061—#1080) проверяет углы и вычитает те, которые не слишком велики¹⁴.

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

Код оптимизирован для применения 64-битной целочисленной арифметики, а не арифметики с плавающей запятой. Это усложняет понимание кода, поскольку значения следует рассматривать как числа с плавающей запятой, имеющие неявную экспоненту. Эти экспоненты нигде не хранятся, но их можно определить, проанализировав алгоритм. Даже угловые константы в ROM не имеют явно заданных экспонент¹⁵. Более того, для сохранения точности неявные экспоненты меняются на каждом шаге цикла. Грубо говоря, каждый такт цикла уменьшает значения вдвое, и угловые константы тоже уменьшаются соответствующе. Если значения имели бы фиксированную экспоненту, то в итоге они начинались бы 16 нулями, впустую тратя точность. Но благодаря масштабированию значений на каждом такте (сдвигом влево) числа продолжают использовать все 64 бита. Иными словами, оборудование выполняет быструю целочисленную арифметику, но математически это можно воспринимать как вычисления с фиксированной запятой, где экспоненты физически не хранятся в чипе. Подробности об изменении неявных экспонент в цикле см. в таблице в примечаниях¹⁶.

Аппроксимация рациональным многочленом

После цикла псевдоделения CORDIC микрокод вычисляет аппроксимацию рациональным многочленом (#1082). Здесь в поток исполнения снова вливается путь вычисления для малых углов. Микрокод для вычисления аппроксимации многочленом обладает любопытными свойствами. Самая затратная часть — это процедура SQUARE, которая умножает число с фиксированной запятой на само себя. Умножение — это сложный процесс, поэтому я приведу подробности в отдельном посте. Если вкратце, 8087 использует алгоритм Бута для умножения сразу на два бита за раз (основание 4), благодаря чему работает вдвое быстрее обычного двоичного умножения. Цикл, выполняющий сдвиги и сложения, реализован аппаратно, а не на микрокоде. Иными словами, одна микрокоманда выполняет 32 сложения: аппаратная логика проверяет биты, выполняет соответствующее сложение или вычитание, обновляет счётчик цикла и повторяет цикл. Для повышения производительности сдвиги выполняются специализированными блоками сдвига, а не универсальным блоком сдвига 8087.

За исключением возведения в квадрат, многочлены вычисляются с помощью сложений, а не умножений. Константа 3 берётся не из ROM констант, а формируется специальными транзисторами, которые устанавливают два старших бита мантиссы. Обычно этот механизм используется для формирования NaN. Если неявный показатель степени равен 1, это соответствует числу 3. Вычитание квадрата даёт числитель. Для получения знаменателя угол сдвигается вправо (то есть делится на 2), а затем складывается с самим собой. Это позволяет умножить его на 3 (точнее, на 1,5, но увеличение неявного показателя степени до 1 превращает результат в трёхкратное значение). Деление числителя на знаменатель было бы дорогостоящей операцией, поэтому вместо этого числитель и знаменатель используются в качестве начальных значений Y и X для последующих шагов алгоритма CORDIC.

Псевдоумножение CORDIC

Далее идёт цикл CORDIC, в котором вращения применяются к вектору (#1100). Вращения применяются в обратном порядке относительно их вычисления на первом этапе: наименьшие вращения применяются первыми для сохранения точности. Биты сдвигаются из регистра сдвига, указывая, следует ли выполнять вращение. Как только в регистре сдвига остаются только нули, цикл прекращается. Иными словами, для малых углов цикл выполняется всего несколько раз, а не все 16 итераций.

Вращения выполняются с помощью сдвигов и сложений (или вычитаний), но отслеживать неявные показатели степени несколько непросто. Поскольку вектор изначально почти горизонтален, значение X составляет приблизительно 3, а значение Y — приблизительно 2^-16. По мере вращения вектора X остаётся примерно равным 3, тогда как Y потенциально может удваиваться на каждой итерации. Чтобы обеспечить максимальную точность представления каждого числа, показателю степени для X (регистр tmpB) присваивается значение 1, а показатель степени для Y (регистр tmpC) изначально равен −14 и увеличивается на 1 с каждой итерацией, по мере того как Y сдвигается вправо. (Напомним, что эти показатели степени нигде не хранятся явно, а подразумеваются.) Поскольку у двух регистров разные неявные показатели степени, перед сложением значения необходимо сдвинуть. Величины сдвигов неочевидны; разобраться в них поможет таблица¹⁶ из примечаний.

После цикла CORDIC заключительная часть кода (#1118) нормализует значения X и Y, формируя числа с плавающей точкой, которые возвращаются в качестве результата инструкции. Напомним, что на этом этапе это ещё не «настоящие» числа с плавающей точкой, поэтому может потребоваться их корректировка. В частности, если у X старший бит не равен 1, значение сдвигается влево на одну позицию. Если же у Y два старших бита равны 1, значение сдвигается вправо на одну позицию. В обоих случаях показатель степени корректируется соответствующим образом. Значения X и Y помещаются в стек (#1132), завершая выполнение инструкции.

Заключение

Из-за своей сложности команда FPTAN довольно медленная по сравнению с другими командами 8087. В документации говорится, что обычно она занимает 450 тактов. (Время может сильно колебаться; в зависимости от значения команда может занимать от 30 до 540 тактов.) В случае выбранного мной значения (0,95), 33% времени занимало псевдоделение CORDIC, 15% — рациональный многочлен (в основном операция возведения в квадрат), 47% — псевдоумножение CORDIC, а 5% — оверхед.

CORDIC — популярный способ вычисления тригонометрических функций, но у него есть альтернативы, например, многочлены. Процессор 8087 необычен тем, что в нём сочетаются CORDIC и рациональный многочлен. В Pentium компания Intel перешла от CORDIC к аппроксимациям многочленами; благодаря схемам быстрого умножения Pentium использование многочленов стало практичным. Сегодня у Intel есть такие библиотеки, как MKL (Math Kernel Library) и Short Vector Math Library (SVML), обеспечивающие реализации, сильно оптимизированные под оборудование Intel. Говорят, что в этих библиотеках применяются аппроксимации многочленами, а для высокой производительности используются параллельные команды (SIMD), а не специализированное оборудование. Из-за этих библиотек операции x87 и 80-битные float практически полностью устарели. Тем не менее, я надеюсь, вам понравился этот анализ алгоритма 8087.

Я планирую продолжать заниматься реверс-инжинирингом микрокода 8087. Я работал над этим с участниками «Opcode Collective» и, в частности, с Smartest Blob и Gloriouscow, которые преобразовали образы ROM в данные микрокода и подробно проанализировали содержимое. Более подробную информацию см. в репозитории 8087 на GitHub.

Примечания и ссылки

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

3. Внутренне 8087 использует для хранения чисел с плавающей запятой 80 бит: один бит знака, 15 бит экспоненты и 64 бит мантиссы. Однако в 8087 задействованы три дополнительных бита низкого порядка для округления, которые называются Guard, Round и Sticky. Эти биты гарантируют, что значение всегда будет округляться в нужном направлении. В некоторых частях тракта данных есть дополнительные биты для знака или переполнения: блок сдвига имеет ширину 68, а сумматор — 69 бит. В основном я буду игнорировать эти дополнительные биты и говорить, что тракт данных имеет ширину 64 бит. 

4. Типичная формула для вычисления sin θ из значений X и Y 8087 имеет вид sin θ = Y/sqrt(X²+Y²); формула для cos θ похожа на неё.

5. Специальные углы (в градусах) приблизительно равны 45°, 26,6°, 14°, 7,1°, 3,6°, 1,8°, 0,9°, … Вы можете задаться вопросом, всегда ли сумма специальных углов сходится к нужному углу θ. Доказательством можно найти в статье. Стоит отметить, что arctan(2^‑n) быстро приближается к 2^‑n (в радианах), поэтому это похоже на выражение двоичной дроби, на каждом шаге уменьшающее разность приблизительно вполовину.

6. Реализация CORDIC в 8087 отличается от типичного применения CORDIC, при котором находятся sin и cos. При таком подходе используются все специальные углы, но на каждом шаге угол или прибавляется, или вычитается, что приводит к схождению к нужному углу. Напомню, что поскольку на каждом шаге матрица умножается на 1/sin α_i, конечный вектор перестаёт быть единичным вектором. В результате этого cos и sin не соответствуют окончательному вектору (X, Y).

   Хитрость заключается в том, что если используются все специальные углы, то суммарный коэффициент масштабирования остаётся тем же, просто произведением всех множителей cos α_i. (При этом направление вращения не важно, поскольку cos α_i = cos ‑α_i.) Таким образом, итоговый коэффициент масштабирования можно вычислить заранее, а исходный единичный вектор можно умножить на эту величину. В результате окончательный вектор (X, Y) непосредственно даёт нам (cos θ, sin θ), что и требовалось.

   Самое важно здесь то, что существует два подхода к CORDIC: использование положительных и отрицательных углов или использование или пропуск углов. (В первом случае каждый угол умножается на +1 или на -1, а во втором углы умножаются на 1 или на 0.) В большинстве реализаций CORDIC применяется первый подход, а в процессоре 8087 — второй.

7. Разработчик CORDIC Джек Волдер написал статью The Birth of CORDIC, в которой подробно описывает его историю, приводя в том числе фотографии навигационного компьютера B-58. Исходные статьи Волдера 1959 года: The CORDIC Computing Technique и The CORDIC Trigonometric Computing Technique.

   

   На фото показан компьютер CORDIC II, использовавшийся на борту бомбардировщика B-58. Это был 30-битный последовательный компьютер, работающий на частоте 198,4 кГц и весящий целых шестьдесят килограммов. Он был реализован на вентилях диодно-транзисторной логики (diode‑transistor logic, DTL), состоял из 2498 транзисторов и 4265 диодов и имел барабанную память (в верхней части снимка), хранящую эквивалент 39 килобайт. Любопытно, что операции CORDIC были реализованы в виде базовых команд набора команд, а не как подпрограммы; команда CORDIC занимала 5 мс, то есть столько же, сколько умножение или деление.

8. Эта аппроксимация Паде называется аппроксимацией [1,2], потому что числитель — это многочлен первого порядка, а знаменатель — второго. Аппроксимации Паде можно вывести из многократного разложения тангенса на дроби. См. Constructing Padé Approximates.

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

9. Алгоритм 8087 описан статье Implementation of transcendental functions on a numerics processor Рафи Наве, руководившего проектированием оборудования 8087.

10. Задокументированный допустимый диапазон аргументов FPTAN равен 0<θ<π/4, то есть использованное мной демонстрационное значение 0,95 находится за его границами. Мне этот диапазон кажется странным, ведь первая константа CORDIC — это arctan(1), π/4, которая полезна только при входных углах больше π/4. Таким образом, 8087 использует дефицитную ячейку ROM и дополнительный такт цикла для обработки значений за пределами допустимого диапазона. У меня на этот счёт есть две гипотезы. Возможно, 8087 был реализован для большего диапазона аргументов, но численная погрешность оказалась больше ожидаемой, поэтому диапазон ограничили. Или же команде арктангенса FPATAN могла требоваться первая константа CORDIC, и было проще задействовать цикл с одной и той же константой для FPTAN и FPATAN. (Я не изучал FPATAN, так что подтвердить этого не могу.) Как бы то ни было, для 0,95 код работает; я использовал это значение в качестве примера, потому что углы получаются больше и нагляднее.

    Ещё одна любопытная особенность заключается в том, что документация расходится в том, входил ли в диапазон 0. В книге The 8087 Primer, написанной архитекторами 8087 и 8086, говорится, что диапазон аргументов не включает 0, и что это решение было принято для упрощения микрокода. Однако при изучении микрокода (#1044) становится ясно, что он явным образом проверяет входные данные на ноль и обрабатывает этот случай.

    Похоже FPTAN работает для значений почти вплоть до π/2 с потерей точности, когда результат становится бесконечно близким к π/2. Однако значения от 1 до π/2 (находящиеся за пределами диапазона и имеющие экспоненту 0), вероятно, работают по совпадению, а не умышленно. Из-за большего значения экспоненты выполняется дополнительная итерация цикла; при обращении к таблице тригонометрических констант происходит переход через её границу, и алгоритм начинает работу с последней константы, atan(2^–15). Однако из-за масштабирования это значение в итоге интерпретируется как приблизительно 2. Так как этот угол CORDIC больше, чем π/2, он игнорируется. Затем алгоритм продолжает вычисления с оставшимися константами, выполняя расчёты должным образом.

    Примечание к примечанию: хотя математический смысл термина range («область значений») — это вывод функции, Intel использует его для описания ввода функции, который математики называют domain («область определения»). ↩

11. Одно из различий между форматом temporary real и внешними форматами в том, что в temporary real начальная 1 хранится в числе в явном виде. Во внешних форматах начальная 1 подразумевается косвенно и не хранится в числе, обеспечивая дополнительный бит точности, но усложняя арифметику.

12. На блок-схеме показан поток управления микрокода FPTAN. Основные блоки — это псевдоделение CORDIC, рациональная аппроксимация и псевдоумножение CORDIC.

    

13. Сравнения занимают несколько микрокоманд, например, выбор константы наподобие 0x3ff0 (напомню, что экспоненты смещены на 16 383, или 0x3fff, поэтому 0x3ff0 соответствует -15), копирование константы в сумматор, вычитание константы из экспоненты и ветвление в зависимости от результата.

14. Первый этап алгоритма CORDIC (#1061—#1068) обрабатывается отдельно. Мне потребовалось довольно много времени, чтобы разобраться, в чём же особенность первого этапа CORDIC. Причина в том, что углы CORDIC близки 2^‑n, но не равны этим значениям. Следовательно, количество циклов CORDIC не меняется в точности с изменением экспоненты исходного угла. И в этом проблема, ведь количество тактов псевдоделения зависит от экспоненты, но количество тактов псевдоумножения зависит от наибольшего использованного угла CORDIC. Так как числа масштабируются в обоих циклах, несоответствие количества циклов приведёт к погрешности результата в 2 раза. Первый этап CORDIC проверяет, используется ли первая константа CORDIC, и для компенсации этого выполняет инкремент экспоненты tmpA на 1. (Если вам ничего непонятно, то просто расслабьтесь.)

15. В ROM констант процессора 8087 хранятся мантиссы, но не экспоненты. Когда я впервые извлекал в 2020 году константы из ROM констант, мне не удалось найти экспоненты, поэтому я придумал правдоподобные экспоненты, дававшие логичные значения. В 8087 есть отдельный ROM экспонент. Теперь, после исследования микрокода, я понимаю, что экспоненты в ROM экспонент не связаны с мантиссами в ROM констант. Я рад, что тогда не потратил кучу на определение экспонент.

16. В таблице ниже показаны вычисления на каждом этапе алгоритма с 0,95 в качестве исходных данных. Имейте в виду, что вычисления выполняются в целочисленной арифметике; экспоненты не хранятся как части чисел, а подразумеваются алгоритмом. Экспоненты меняются на каждом этапе, поэтому при уменьшении значений на стадии псевдоделения экспоненты тоже уменьшаются, гарантируя, что значения смогут использовать максимум доступных битов. При псевдоумножении происходит обратное: изначально значения малы и увеличиваются с вращением вектора; чтобы компенсировать это, на каждом этапе экспонента увеличивается. В отличие от «настоящих» чисел с плавающей запятой, битовые представления не нормализованы; иногда у них есть два первых разряда-единицы слева от двоичного разделителя, а иногда в начале идут нули. Тем не менее, на каждом этапе экспоненты фиксированы; они не меняются в зависимости от значения.

    Операция	Биты	Значение	Флаг	Регистр
    load θ	1.1110011…×2-1	0.9500		tmpC
    cmp atan(20)	1.1001001…×2-1	0.7854	1	const
    sub angle	0.0101010…×2-1	0.1646		tmpC
    сдвиг влево	0.1010100…×2-2	0.1646		tmpC
    cmp atan(2-1)	1.1101101…×2-2	0.4636	0	const
    сдвиг влево	1.0101000…×2-3	0.1646		tmpC
    cmp atan(2-2)	1.1111010…×2-3	0.2450	0	const
    сдвиг влево	10.1010001…×2-4	0.1646		tmpC
    cmp atan(2-3)	1.1111110…×2-4	0.1244	1	const
    sub angle	0.1010010…×2-4	0.04025		tmpC
    сдвиг влево	1.0100100…×2-5	0.04025		tmpC
    cmp atan(2-4)	1.1111111…×2-5	0.06242	0	const
    сдвиг влево	10.1001001…×2-6	0.04025		tmpC
    cmp atan(2-5)	1.1111111…×2-6	0.03124	1	const
    sub angle	0.1001001…×2-6	0.009007		tmpC
    сдвиг влево	1.0010011…×2-7	0.009007		tmpC
    cmp atan(2-6)	1.1111111…×2-7	0.01562	0	const
    сдвиг влево	10.0100111…×2-8	0.009007		tmpC
    cmp atan(2-7)	1.1111111…×2-8	0.007812	1	const
    sub angle	0.0100111…×2-8	0.001195		tmpC
    сдвиг влево	0.1001110…×2-9	0.001195		tmpC
    cmp atan(2-8)	1.1111111…×2-9	0.003906	0	const
    сдвиг влево	1.0011100…×2-10	0.001195		tmpC
    cmp atan(2-9)	1.1111111…×2-10	0.001953	0	const
    сдвиг влево	10.0111001…×2-11	0.001195		tmpC
    cmp atan(2-10)	1.1111111…×2-11	0.0009766	1	const
    sub angle	0.0111001…×2-11	0.0002181		tmpC
    сдвиг влево	0.1110010…×2-12	0.0002181		tmpC
    cmp atan(2-11)	1.1111111…×2-12	0.0004883	0	const
    сдвиг влево	1.1100100…×2-13	0.0002181		tmpC
    cmp atan(2-12)	1.1111111…×2-13	0.0002441	0	const
    сдвиг влево	11.1001001…×2-14	0.0002181		tmpC
    cmp atan(2-13)	1.1111111…×2-14	0.0001221	1	const
    sub angle	1.1001001…×2-14	0.00009604		tmpC
    сдвиг влево	11.0010010…×2-15	0.00009604		tmpC
    cmp atan(2-14)	1.1111111…×2-15	0.00006104	1	const
    sub angle	1.0010010…×2-15	0.00003500		tmpC
    сдвиг влево	10.0100101…×2-16	0.00003500		tmpC
    cmp atan(2-15)	1.1111111…×2-16	0.00003052	1	const
    sub angle	0.0100101…×2-16	0.000004482		tmpC
    аппроксимация X	1.0111111…×21	3.0000		tmpB
    аппроксимация Y	0.0011100…×2-14	0.00001345		tmpC
    Y += 2-15X	1.1011100…×2-14	0.0001050	1	tmpC
    X ‑= 2-15Yold	1.0111111…×21	3.0000		tmpB
    сдвиг вправо	0.1101110…×2-13	0.0001050		tmpC
    Y += 2-14X	10.0101110…×2-13	0.0002881	1	tmpC
    X ‑= 2-14Yold	1.0111111…×21	3.0000		tmpB
    сдвиг вправо	1.0010111…×2-12	0.0002881		tmpC
    Y += 2-13X	10.1010111…×2-12	0.0006543	1	tmpC
    X ‑= 2-13Yold	1.0111111…×21	3.0000		tmpB
    сдвиг вправо	1.0101011…×2-11	0.0006543		tmpC
    сдвиг вправо	0.1010101…×2-10	0.0006543	0	tmpC
    сдвиг вправо	0.0101010…×2-9	0.0006543	0	tmpC
    Y += 2-10X	1.1101010…×2-9	0.003584	1	tmpC
    X ‑= 2-10Yold	1.0111111…×21	3.0000		tmpB
    сдвиг вправо	0.1110101…×2-8	0.003584		tmpC
    сдвиг вправо	0.0111010…×2-7	0.003584	0	tmpC
    сдвиг вправо	0.0011101…×2-6	0.003584	0	tmpC
    Y += 2-7X	1.1011101…×2-6	0.02702	1	tmpC
    X ‑= 2-7Yold	1.0111111…×21	3.0000		tmpB
    сдвиг вправо	0.1101110…×2-5	0.02702		tmpC
    сдвиг вправо	0.0110111…×2-4	0.02702	0	tmpC
    Y += 2-5X	1.1110111…×2-4	0.1208	1	tmpC
    X ‑= 2-5Yold	1.0111111…×21	2.9991		tmpB
    сдвиг вправо	0.1111011…×2-3	0.1208		tmpC
    сдвиг вправо	0.0111101…×2-2	0.1208	0	tmpC
    Y += 2-3X	1.1111101…×2-2	0.4957	1	tmpC
    X ‑= 2-3Yold	1.0111110…×21	2.9840		tmpB
    сдвиг вправо	0.1111110…×2-1	0.4957		tmpC
    сдвиг вправо	0.0111111…×20	0.4957	0	tmpC
    сдвиг вправо	0.0011111…×21	0.4957	0	tmpC
    Y += 20X	1.1011110…×21	3.4797	1	tmpC
    X ‑= 20Yold	1.0011111…×21	2.4884		tmpB
    окончательный X	1.0011111…×21	2.4884		tmpB
    окончательный Y	1.1011110…×21	3.4797		tmpA

Если эта публикация вас вдохновила и вы хотите поддержать автора — не стесняйтесь нажать на кнопку