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

Когда я начинал этот проект в 2021 году, мне нужен был код на C++, который бы реализовывал алгоритмы на основе примитивов BCD и проверял результаты. Этот код превратился в подпроект Methods. Он заработал, но в нём возникли небольшие баги с шириной мантиссы, из-за которых тестовые результаты оказались ненадёжными в пределах одного-двух последних разрядов. Вместо того, чтобы патчить его, я в 2025 году переписал всё с нуля в рамках подпроекта Proto: более чистой архитектуры, правильного эталона и генератора тестовых векторов оборудования, относительно которого можно валидировать микрокод FPGA.

Внутренний числовой формат

Регистры этого научного калькулятора хранят:

  • 16-разрядную мантиссу BCD, нормализованную так, что первый разряд всегда равен 1–9 (за исключением истинного нуля, для которого все 16 разрядов равны нулю)

  • Бит знака для мантиссы

  • 2-разрядную экспоненту BCD, обеспечивающую диапазон от 10⁻⁹⁹ до 10⁹⁹

  • Бит знака для экспоненты

Эта структура аналогична той, которая использовалась в HP-35 и HP-41 (в этих машинах применены 10-разрядные мантиссы). Такой формат очень подходит для десятичного калькулятора: ввод пользователя сохраняется в том же виде, без погрешности двоичного представления. При вводе 1 ÷ 3 × 3 мы получим ровно 1, а не 0.999999999. Все промежуточные результаты остаются в области десятичных значений, как и ожидает пользователь.

Методы: исходная реализация

В проекте Methods я реализовал на C++ восемь функций: сложение/вычитание, умножение, деление, ln, exp, tan и atan. Каждая реализация включала в себя и алгоритм BCD, и тестовую обвязку, которая сравнивала результаты с double языка C++. Макрос MAX_MANT изначально должен был везде иметь значение 16, но в тот раз я реализовал это неправильно и получил смесь из 14- и 16-значной точности — это был самый первый прототипный код, и я одновременно экспериментировал со слишком большим количеством размеров буферов. Вызванные этим баги были малозаметными: результаты казались верными, но расходились на несколько последних десятичных разрядов. Я не мог отличить реальные погрешности от артефактов округления.

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

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

Поиск нужного эталона

Тип double с плавающей запятой не может полностью верифицировать 16-разрядную BCD-арифметику, потому что его предел — это 15 разрядов. Хотя 15 разрядов вроде бы достаточно для большинства сценариев, для моего это было неидеально, потому что мне требовалось 16. Кроме того, когда начинаешь выполнять цепочки операций, точность падает ещё сильнее.

Оказалось, что решение скрывается в компиляторе. На x86-64 GCC реализует long double в виде 80-битного типа повышенной точности при помощи FPU x87: полная 64-битная мантисса, которую Intel использовала в каждом процессоре x86, начиная с 8087 (1980 год). Это даёт нам приблизительно 18-19 значимых десятичных разрядов, то есть дополнительный ресурс в два-три разряда по сравнению с 16-разрядным BCD-результатом. Если BCD и long double не совпадают в разряде 16, то BCD ошибочна.

Проблема в том, что всё это работает только в Linux с GCC. В Windows MSVC отображает long double в простой 64-битный double, не обеспечивая никакой дополнительной точности. Хоть код Proto может работать и в Windows, для полной точности ему требуется Linux или WSL. Часть моей системы разработки в Windows использует WSL, так что всё сработало отлично.

Перед тем, как остановиться на long double, я рассмотрел множество альтернатив:

  • GNU MP / MPFR — это золотой стандарт арифметики произвольной точности на C, но его API слишком подробный, а система сборки довольно тяжела для моего проекта.

  • Boost.Multiprecision — это обёртка вокруг MPFR, обеспечивающая более чистый интерфейс на C++, но ей всё равно требуется библиотека.

  • Наиболее многообещающим был пакет для C++ Хенрика Вестермарка, обеспечивающий произвольную точность — только заголовок, портируемый, отсутствие зависимостей, однако в Linux GCC я столкнулся с багами пограничных случаев, приводившими к неправильным результатам в некоторых диапазонах входных данных; из-за этого замечательный пакет превратился в ещё одну спорную переменную, поэтому я решил отказаться от него.

  • Модуль decimal Python оказался полезным для точечной проверки конкретных значений.

  • Wolfram Alpha использовался для ручной проверки некоторых критичных результатов.

Ни один из вариантов не давал особых преимуществ по сравнению с GCC, в котором всё это имелось «бесплатно».

Proto: золотой эталон

Проект Proto — это полностью переписанный код: он имеет согласованные 16-разрядные мантиссы BCD, верифицирован относительно 80-битного long double компилятора GCC в Linux, имеет готовый фреймворк допусков с тремя классами (Strict, Standard, Relaxed) и содержит генератор тестовых векторов для оборудования. Благодаря 18–19 разрядам эталонной точности для сравнения с 16 разрядами BCD какие-либо неопределённости отсутствуют: расхождение в любом разряде будет означать баг BCD.

Proto работает в трёх режимах, выбираемых флагами командной строки:

  • Режим разработчика (по умолчанию): прогоняет все тестовые векторы и выводит только расхождения (строки NEAR и MISS). При чистом прогоне сообщения не выводятся. Флаг -c включает вывод в цветах ANSI, позволяющий выделять расходящиеся разряды. Ну и выглядит это красиво.

  • Режим аппаратных векторов (-t): генерирует тестовые векторы для верификации оборудования; по одному вектору на строку в формате с фиксированными столбцами, который может напрямую обрабатывать тестовая обвязка Verilog. Таким образом генерируются векторы для каждого алгоритма (всех 25 операций).

  • Интерактивный калькулятор с RPN (-?): четырёхуровневый RPN-стек командной строки с полной арифметикой BCD. При вводе 3 enter 4 * получается результат 12. Он полезен для точечной проверки результатов и для сравнения с оборудованием. Этот режим я придумал уже на поздних этапах разработки: если уж я закодировал большинство функций, то почему бы не добавить строку ввода для выполнения RPN-выражений?

Тестирование как основная задача

Баг FDIV Intel Pentium 1994 года — живой пример недостаточного тестирования. Отсутствие пяти записей в таблице поиска из 1066 элементов заставляло процессор возвращать ошибочные результаты при делении. Миллионы процессоров были выпущены с этим багом, прежде чем профессор математики Томас Найсли не обнаружил расхождения при проверке вычислений простых чисел. Техническая первопричина была крошечной, но её последствия стоили Intel 475 миллионов долларов.

Proto предназначен для отлавливания эквивалента этого бага в BCD: погрешности округления a в 16-м разряде, смещения на единицу итерации CORDIC и пограничного случая в логике разряда защиты.

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

Методика проста: Proto генерирует тестовые векторы, выполняя свою эталонную реализацию BCD для большого сгенерированного им множества векторов ввода и выводя результаты. Тестовая обвязка FPGA передаёт эти векторы микрокоду и сравнивает результаты. Расхождение будет означать наличие бага или в микрокоде, или в алгоритме BCD; в любом случае он быстро проявится. В части 4 будет подробно рассмотрен весь конвейер работы с тестовыми векторами.

Разряд защиты, бит фиксации и округление

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

Точность зависит и от стратегии округления, и от алгоритма, используемого для вычисления каждой функции. В случае трансцендентных функций это алгоритм CORDIC. В моей реализации четыре основные операции (+,-,*,/) обеспечивают полную 16-разрядную точность (0,5 ULP), а трансцендентные функции теряют от одного до трёх разрядов, в зависимости от внутреннего количества итераций алгоритма.

Примечание: существует классический трюк для проверки внутренней точности калькулятора: ввести 9, а затем последовательно вычислить sin, cos, tan, arctan, arccos, arcsin (в режиме градусов). Математически arcsin(arccos(arctan(tan(cos(sin(9°)))))) должно вернуть ровно 9: мы просто применяем шесть функций и обратные им. На практике же, каждый калькулятор возвращает собственные результаты. HP-41C возвращает 9.000417403. HP-42S возвращает 8.99999864267. TI с чипом Mostek возвращает 8.99999614252. Отклонение от 9 — это признак семейства внутренних алгоритмов и количества разрядов защиты.

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

CORDIC: один алгоритм, чтобы править всеми

Добившись работы основных четырёх операций, можно реализовать трансцендентные функции. Я использовал CORDIC — то же семейство алгоритмов, которое применили в 1972 году в калькуляторах HP-35. Он применяется при проектировании BCD-калькуляторов, потому что требует только сложения, вычитания и сдвиги разрядов, без умножения во внутреннем цикле.

Легче всего понять пример с тригонометрией. Чтобы найти tan(32°), начнём с точки Q в координате (1, 0) на единичной окружности, и повернём её против часовой стрелки, чтобы суммарный угол поворота равнялся 32°. При этом тангенс будет равен Qy/Qx.

Формула поворота выглядит так:

Qx_new = Qx - Qy × tan(θ)
Qy_new = Qy + Qx × tan(θ)

Хитрость заключается в том, на какие углы выполнять поворот. Если выбрать углы, тангенсы которых равны степеням десятки: 1, 0.1, 0.01, 0.001… то умножение на тангенс будет сдвигом десятичного разряда, а не умножением. Такими углами оказались 45°, 5.711°, 0.5729° и так далее; их значения заранее вычислены и сохранены в таблицу.

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

В проекте Methods для всех функций используются элементы таблицы CORDIC с K=8. Сырой CORDIC сходится примерно по одному разряду на итерацию, то есть 8 итераций обеспечивают примерно 8 разрядов сырой точности. Proto сохраняет K=8 для тригонометрических функций, но восстанавливает оставшуюся точность благодаря остаточному делению после цикла CORDIC и уточнения рядами Тейлора для небольших углов, достигая точности примерно в 14 разрядов для tan и atan. Для ln и exp Proto дополнил таблицу до 15 итераций, достигнув 15-разрядной точности.

Покажу это на конкретном примере: полной трассировке вычисления tan(32°). Алгоритм выполняется в две фазы.

Фаза 1 – извлечение разрядов (псевдоделение). Входной угол (32° = 0.5585 радиан после уменьшения диапазона) раскладывается на сумму заранее вычисленных значений atan. На каждом этапе j алгоритм многократно вычитает atan(10⁻ʲ) из оставшегося угла, пока он остаётся неотрицательным. Количество вычитаний становится извлечённым разрядом этого этапа:

Этап j

atan(10⁻ʲ)

Разряд

Оставшийся угол

0

0.7854 (= 45°)

0

5.585×10⁻¹

1

0.09967 (= 5.71°)

5

6.016×10⁻²

2

0.01000

6

1.641×10⁻⁴

3

0.001000

0

1.641×10⁻⁴

4

0.0001000

1

6.410×10⁻⁵

5

0.00001000

6

4.098×10⁻⁶

6

0.000001000

4

9.806×10⁻⁸

7

0.0000001000

0

9.806×10⁻⁸

На этапе 0 atan(1) = 0.7854 радиан (45°) больше нашего угла 0.5585, поэтому разряд равен 0 – мы не можем выполнить вычитание даже один раз. На этапе 1 atan(0.1) = 0.09967 радиан (5.71°) умещается 5 раз: 5 × 0.09967 = 0.4984, что оставляет нам 0.0602. На этапе 2 atan(0.01) умещается 6 раз, потратив почти весь наш остаток и оставив всего 1.641 × 10⁻⁴. На каждом этапе пропадает новый десятичный разряд угла. После 8 этапов оставшийся угол ничтожен (9.8 × 10⁻⁸), а извлечённые разряды равны [0, 5, 6, 0, 1, 6, 4, 0]. Это промежуточный результат.

Фаза 2 – поворот CORDIC (псевдоумножение). Теперь алгоритм воссоздаёт тангенс. Начав с единичного вектора (x, y) = (1, 0), он применяет повороты, используя извлечённые значения разрядов, двигаясь от самого мелкого масштаба (j=7) до самого крупного (j=0). Каждый отдельный поворот в масштабе j вычисляется так:

y_new = y + x × 10⁻ʲ
x_new = x − y × 10⁻ʲ

Умножение на 10⁻ʲ — это сдвиг на десятичный разряд, поэтому при каждом повороте используется только сложение и вычитание – во внутреннем цикле отсутствует умножение. Разряд, извлечённый в фазе 1, сообщает, сколько раз применять этот поворот в каждом из масштабов. Например, при j=1 разряд равен 5, поэтому поворот последовательно применяется 5 раз в масштабе 10⁻¹. После 8 раундов поворотов вектор (x, y) оказывается повёрнут на суммарный угол, равный исходному вводу.

Последний этап – деление. Тангенс — это соотношение получившихся компонент y и x:

y = 0.5434297703056166
x = 0.8696694255289085
tan(32°) = y / x

Это деление — единственная операция класса умножения во всём процессе вычислений; для достижения этого этапа ядро CORDIC не использовало ничего, кроме сдвигов и сложений.

Результат BCD: 0.6248693519093394
Ожидаемый:     0.6248693519093275
Погрешность:      1.2 × 10⁻¹⁴

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

Операция

Точность

Алгоритм

+, −, ×, ÷

16 разрядов

Разряд защиты + округление до ближайшего чётного

√x

~15 разрядов

Ньютон-Рафсон + коррекция ULP

ln

~14,5 разряда

CORDIC, 15 итераций

e^x

~13–14 разрядов

CORDIC (функция, обратная ln)

atan

~14,8 разряда

CORDIC, K=8 + остаточное деление

tan

~14–14,5 разряда

CORDIC + уменьшение диапазона

sin, cos

~13,5–14 разрядов

Формула половинного угла через tan

В статье Джона Вальтера 1971 года показано, что sin/cos/atan, exp/ln/sqrt и умножение/деление — это всё один алгоритм с единственным параметром режима, выбирающим, круговые, гиперболические или линейные координаты. В чипах HP использовалось именно это: один аппаратный цикл, три семейства функций.

Гиперболический режим работает аналогично, только с другим набором предварительно вычисленных констант (полученных из значений гиперболического тангенса) и немного отличающимися свойствами сходимости; для обеспечения сходимости некоторые итерации необходимо повторять. Линейный режим всё равно проще: он сводится к простому умножению и делению сдвигом и сложением. Изящество структуры цикла почти одинаково во всех трёх режимах; меняются лишь таблица констант и логика знаков.

Основа большинства десятичных алгоритмов CORDIC моего проекта заложена в работах французского инженера Жака Лапорта, потратившего годы на реверс-инжиниринг алгоритмов калькуляторов HP. Он разобрался, как HP-35 вычислял ln, exp, tan и atan, из опубликованных данных таймингов и поведенческого тестирования. Жак умер в 2015 году, однако его архивированный веб-сайт остаётся самым подробным публичным ресурсом об алгоритмах BCD-калькуляторов HP.

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

Функция

Алгоритм

Файл исходников

+, −, ×, ÷

Сдвиг и сложение/вычитание с разрядом защиты и округлением до ближайшего чётного

addsub.asm, mul.asm, div.asm

tan, atan

Круговой CORDIC (псевдоделение/псевдоумножение)

tan.asm, atan.asm

ln, e^x

Поразрядный метод Meggitt (гиперболический CORDIC)

log.asm, exp.asm

√x

Метод Ньютона-Рафсона с квадратичной коррекцией ULP

sqrt.asm

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

Функция

Получена из

sin(x)

2t / (1 + t²), где t = tan(x/2)

cos(x)

sin(x + 90°)

asin(x)

atan(1 / √(1/x² − 1))

acos(x)

90° − asin(x)

log₁₀(x)

ln(x) / ln(10)

10^x

e^(x × ln(10))

y^x

e^(x × ln(y))

ˣ√y

e^((1/x) × ln(y))

1/x

1 ÷ x

x²

x × x

%

(x / 100) × y

%CHG

100 × (x − y) / y

D→R

d × (π / 180)

R→D

r × (180 / π)

→HMS, HMS→

Преобразование формата десятичные часы ↔ H.MMSS

n!

Итеративное умножение

nPr, nCr

Итеративное произведение и факториал

GCD

Алгоритм Евклида

LCM

|a × b| / GCD(a, b)

sinh(x)

(e^x − e^(−x)) / 2

cosh(x)

(e^x + e^(−x)) / 2

tanh(x)

(e^(2x) − 1) / (e^(2x) + 1)

asinh(x)

ln(|x| + √(x² + 1))

acosh(x)

ln(x + √(x² − 1))

atanh(x)

ln((1 + |x|) / (1 − |x|)) / 2

x̄, Sx, σx

Скользящее среднее и стандартное отклонение (алгоритм Уэлфорда)

P→R

x = r·cos(θ), y = r·sin(θ)

R→P

r = √(x² + y²), θ = atan2(y, x)

36 операций по цене восьми. Разве это не чудесно? 

Что дальше

Мы получили работающий верифицированный алгоритмический движок BCD на C++: 16 разрядов, разряд защиты, бит фиксации, округление до ближайшего чётного; пройдены тесты, сравнивающие с 80-битным long double, сгенерированы тестовые векторы оборудования. Следующим этапом будет проектирование CPU для выполнения движка. Микрокод нельзя писать, пока нет набора команд для его написания, оборудования FPGA для его запуска и верификации. В части 4 я расскажу о фреймворке разработки, объединяющем все эти инструменты. В части 5 мы получим первый аппаратный прототип. Часть 6 посвящена проектированию CPU.