Описанный в предыдущих частях калькулятор уже работал. Реализация 2021 года запускалась на реальном оборудовании, обрабатывала нажатия клавиш, вычисляла результаты и отображала их. Арифметика была корректной в том смысле, что большинство результатов было точным до 12 значимых разрядов, а это больше, чем требуется обычному пользователю.
Но «большинство» это не «все», а «примерно 12» — это не 15-16 разрядов, которые может и должна обеспечивать 16-разрядная BCD-машина. Существовали пограничные случаи, в которых результаты оказывались совершенно неверными. Имелись итеративные алгоритмы с точностью приемлемой, но не такой, какой она могла быть. Кроме того, в процессе тестирования я обнаружил ошибки, при отладке которых обнаружились фундаментальные баги в коде прототипа на C++. Это привело меня в смятение, ведь для их устранения мне бы пришлось переделать заново код прототипа. В конечном итоге, так я и поступил. Старый код я оставил в репозитории (Pathfinding/Methods) и с нуля разработал совершенно новую версию (Pathfinding/Proof). Я пообещал себе, что занимаюсь этим последний раз в жизни, поэтому стремился делать всё идеально.
В итоге, версия 2025 года устранила найденные мной проблемы. Это был не патч, а почти полное переписывание арифметического движка, расширение набора команд CPU и существенное увеличение библиотеки функций. В этом посте я расскажу об изменениях и их причинах.
Более совершенный эталон
Proto 2025 года (описанный в части 3) был переписан с нуля с полными реализациями всех нужных калькулятору функций. Самое важное заключалось в том, что он генерировал тестовые векторы: тысячи пар входных и выходных данных, которые утилита calctest прогоняла для реальной Verilog-симуляции. При запуске ./proto -t -a > hw.txt выводились тысячи подобных строк:
ADD +1.234567890123456e+15 +9.876543210987654e+10 +1.234567890123456e+15 OK SQRT -1.000000000000000e+00 INVALID nan
В каждой строке указаны конкретные входные данные, ожидаемые выходные и код статуса. Утилита calctest считывает эти векторы и пропускает их через Verilog-симуляцию, сравнивая результаты. При возникновении расхождений можно точно узнать, в каком случае произошла ошибка. Таким образом я замкнул цикл между эталоном на C++ и оборудованием: больше мне не требовалось изучение вручную и не нужно было гадать. Qt-приложение тоже могло считывать эти тестовые векторы, которые можно использовать для отладки более сложных несоответствий.
Тестовые числа включали в себя прописанные вручную значения, случайные значения и все мыслимые пограничные случаи. Также в них входили значения, которые должны генерировать ошибки, и возникновение этих ошибок тоже проверялось. Особыми случаями были линейно изменяющиеся значения на границах перехода между диапазонами, которые позволяют обнаружить ошибку даже в одном младшем значащем разряде.
Переписывание арифметики: защитные разряды и биты фиксации
Фундаментальной проблемой арифметики в версии 2021 года было усечение. После каждой операции результаты просто обрезались до 16 разрядов. Это корректно в смысле получения 16-разрядного результата, но это неточное решение: отбрасываемый 17-й разряд часто содержал в себе информацию о правильном округлении 16-го разряда.
Арифметика в версии 2025 года отслеживает механизм защитных разрядов и битов фиксации для каждой из основных операций: сложения, вычитания умножения и деления. Они обеспечивают возможность банковского округления (часто называемого округлением до ближайшего чётного): если защитный разряд больше 5, округление выполняется наверх; если он меньше 5, то выполняется усечение; если он равен 5, то округление производится вверх только при установленном бите фиксации (означающем, что истинное значение было строго больше x.5) или если последний сохранённый разряд нечётный. Это устраняет систематический перекос в сторону округления наверх до точных половинных значений, что важно при комбинировании множества округлённых значений.
Название «банковское округление» ошибочно: есть мало исторических свидетельств того, что банки когда-либо использовали его в качестве стандарта, а Европейская комиссия отметила, что банковский сектор традиционно округлял половинные значения наверх. Этот способ также называется округлением по Гауссу или голландским округлением. На самом деле, широкое распространение он получил благодаря стандарту IEEE 754 1985 года, который назначал округление до ближайшего чётного режимом округления по умолчанию для всего отвечающего требованиям оборудования, работающего с числами с плавающей запятой. До внедрения IEEE 754 способы округления в разных машинах варьировались: часто использовалось усечение и округление половинного значения наверх; выбранный способ часто оказывался плохо задокументирован и специфичен для оборудования. Внедрение IEEE 754 стал первым случаем согласованного выбора понятия корректного округления для всех отрасли.
В случае вычитания бит фиксации выполняет другую, но столь же важную роль: он генерирует начальное заимствование в вычитании мантиссы. Это исключает класс ошибок, которому была подвержена реализация 2021 года. Рассмотрим пример:
1.000000000000001 × 10⁰ - 9.999999999999999 × 10⁻¹
При выполнении вычитания второй операнд сдвигается вправо на одну позицию, чтобы обеспечить одинаковость экспонент. В реализации 2021 года (15 разрядов, защитный разряд отсутствует) последняя 9 просто терялась. При вычитании мы получали 1 × 10⁻¹⁴, то есть значение примерно девять раз больше. В новой реализации (16 разрядов: 15 значимых + 1 защитный) сдвинутая 9 сохраняется в виде защитного разряда, а вычисляемый результат корректен (2 × 10⁻¹⁵, что приблизительно равно истинному значению 1.1 × 10⁻¹⁵). Бит фиксации применяется в случаях, когда разряды выходят за пределы даже защитного разряда, обеспечивая учёт потери дополнительной точности и позволяя применять соответствующие коррекции.
Итоговый эффект: основные арифметические операции (сложение, вычитание, умножение, деление) теперь обеспечивают 16 корректно округлённых значимых разрядов (в пределах 0,5 ULP). В случае трансцендентных функций улучшения меньше, но они тоже заметны: квадратный корень достигает точности примерно в 14,9 разряда, логарифм — примерно в 14,5, возведение в степень — примерно 13-14. Теперь все функции переписаны согласно эталонному Proto, и каждая переписанная функция предварительно валидировалась векторами calctest.
Новые команды CPU
Более сложная арифметика занимала больше пространства кода. Для отслеживания защитных разрядов и битов фиксации 16-ниббловых BCD-операций требовались дополнительные команды, а в наборе команд 2021 года не было хороших инструментов для создания сжатого кода.
Для решения этой проблемы я добавил четыре новые команды: CALLI (оптимизация «вызов с аргументами», описанная в части 6; само по себе это уже сэкономило 186 слов ROM), TBLCALL для диспетчеризации таблицы одиночных команд, BSHR для BCD-деления на 2 с объединением переносов в цепочки и DECA для декрементов счётчиков цикла, не влияющих на цепочку арифметических переносов. Также я избавился от двух команд: BRANC и TEST. В итоге получилось освободить достаточно пространства ROM для добавления всего набора гиперболических функций; при этом не понадобилось идти ни на какие компромиссы.
Набор тригонометрических функций
Благодаря тому, что новые команды освободили кодовое пространство, в версии 2025 года были улучшены следующие функции:
sin,cos,tanи для градусов, и для радианasin,acos,atanи для градусов, и для радианatan2(арктангенс с двумя аргументами, охватывающий все четыре квадранта)преобразование из полярных в прямоугольные и из прямоугольных в полярные координаты (в режимах градусов и радиан)
улучшенный код уменьшения диапазона
В тригонометрических алгоритмах применён тот же подход, что и в логарифмических со степенными: CORDIC, реализующий метод псевдоделения/умножения HP-35. При работе с градусами тригонометрические функции выполняют полное уменьшение диапазона, а при работе с радианами они сначала производят преобразование в градусы, а затем идут по тому же пути. Двойное преобразование (каждое добавляет не больше 0,5 ULP погрешности) — эта та цена, которой мы покупаем улучшенную точность благодаря отсутствию умножения на иррациональное π/2 при уменьшении диапазона.
Об уменьшении диапазона стоит поговорить подробнее, потому что при неправильной его реализации происходит быстрая потеря точности при увеличении входных значений.
Алгоритм CORDIC наиболее точен для малых углов (примерно от 0 до π/4 радиан, или от 0° до 45°). В случае бо́льших входных значений необходимо сначала уменьшить угол до эквивалентного значения в этом диапазоне. Обычно это делается при помощи x mod (π/2), вычисляемого как x - n*(π/2), где n = floor(x / (π/2)). Для небольших углов это работает нормально, но при их увеличении погрешности становятся всё серьёзнее. На эти погрешности влияют две независимые причины.
Возьмём для примера tan(1000) в радианах. 1000 / (π/2) ≈ 636.6198. Для уменьшения диапазона нужно вычислить 1000 - 636 × π/2 ≈ 0.9735. И 1000, и 636 × π/2 ≈ 999.027 имеют близкую (большую) величину, поэтому из-за требуемого вычитанию сдвига мы можем потерять до трёх младших разрядов мантиссы из шестнадцати. При бо́льших начальных значениях ситуация только усугубляется, потому что вычитаемые числа ближе друг к другу. Это катастрофическая потеря значащих цифр — первая причина погрешности.
Вторая заключается в том, что иррациональная константа π хранится со всего 16 разрядами мантиссы, что при выполнении операций с ней привносит систематическую погрешность. Погрешность в хранящемся значении π умножается на n, поэтому она растёт с величиной входных данных. Результат тригонометрической функции, вычисленный на основании таких входных данных, будет ошибочным вне зависимости от точности алгоритма CORDIC.
Предварительное преобразование значений в градусы с последующим уменьшением диапазона снижает влияние обеих погрешностей, поскольку деление на 180, 90, 45, то есть на точные целые числа в BCD, не вызывает на этом этапе аппроксимаций и потери точности. Первая причина (потеря значащих цифр) сохраняется (вычитание близких значений), зато вторая устраняется полностью: при уменьшении диапазона в градусах не задействуются иррациональные константы.
Только после уменьшения до диапазона [0°, 45°) угол преобразуется в радианы для алгоритма CORDIC. Радианные функции (sinRad, cosRad) выполняют на входе преобразование в градусы, после чего идут по пути вычислений с градусами. Исключением становится tanRad: она остаётся в формате радиан и использует уменьшение диапазона по границе π/4, что хорошо работает для малых углов, но увеличивает погрешность при больших входных радианах, поскольку снова появляются две описанные выше причины погрешностей.
Существуют методики уменьшения диапазона в радианах, работающие для более широких диапазонов. Алгоритм Коди–Уэйта разбивает константу на сумму членов (верхний, средний, нижний), чтобы каждое частное произведение n × c_i было точным; алгоритм Пэйна–Ханека использует предварительно вычисленные части 2/π с повышенной точностью. Оба алгоритма повышают сложность, а поскольку такие очень большие значения радиан на практике маловероятны, я не стал их реализовывать.
BCD-калькуляторы HP традиционно работали с градусами.
Синус и косинус вычисляются при помощи тождества половинных углов: sin(x) = 2t/(1 + t²), где t = tan(x/2). Благодаря этому почти вся реализация sin и cos общая с tan, а размер кода остаётся приемлемым. Для арксинуса и арккосинуса используется asin(x) = atan(1/sqrt(1/x² - 1)), аккуратно выраженный так, чтобы избежать конфликта регистров между вызовами sqrt и atan.
Большую долю получившихся функций реализует слой скриптинга. Многие варианты вычисления обратных значений и преобразований реализованы всего несколькими токенами скриптинга поверх базовых реализаций CORDIC.
Новые функции и клавиши
Наряду с тригонометрическими функциями, в версии 2025 года появилось множество других калькуляторных функций:
LASTX — это регистр X до последней операции. Он необходим для восстановления после ошибок: если пользователь нажмёт не ту кнопку, LASTX вернёт утерянное значение. Он реализован в виде отдельного регистра, который обёртки pre_calc сохраняют автоматически перед каждой унарной и бинарной операцией.
Впервые функция LASTX появилась в HP-45 в 1973 году и стала одной из множества постоянных фич линейки калькуляторов HP. Также в HP-45 впервые появились нумерованные регистры хранения (девять, наряду с четырёхуровневым стеком), выбираемые пользователем тригонометрические режимы (градусы, радианы, грады) и клавиша Shift общего назначения (жёлтый префикс «f» на основе более ограниченной клавиши «arc» калькулятора HP-35). В этой машине появилось всё то, что современные пользователи калькуляторов с обратной польской нотацией считают чем-то обыденным.
Регистры памяти STO/RCL: десять адресуемых пользователем регистров памяти (с 0 по 9), хранящиеся в ОЗУ. При нажатии STO и цифры значение X сохраняется в соответствующий регистр; нажатие RCL извлекает его. В первой версии использовалась одна ячейка памяти STO/RCL.
CLx/A: при первом нажатии сбрасывается X, а при втором — весь стек. Для двойного действия используется конечный автомат режима ввода, определяющий, сброшен ли X.
Преобразование HMS: →HMS и HMS→ выполняют преобразования между десятичными часами и форматом Ч.ММ.СС. Каждая из функций состоит из нескольких токенов скриптинга, использующих уже имеющиеся примитивы умножения и операций по модулю .
Случайное число: RND возвращает возвращает случайное число с равномерным распределением в диапазоне [0, 1). В качестве генератора случайных чисел используется аппаратно реализованный Galois LFSR (Verilog), доступ к которому выполняется через MMIO. Аппаратный LFSR с частотой 50 МГц генерирует новое значение на каждом такте; микрокод считывает валидные BCD-нибблы и нормализует результат.
GCD и LCM: наибольший общий делитель и наименьшее общее кратное, реализованные в виде алгоритмов Евклида в слое скриптинга. В коде LCM вызывает GCD. Это бинарные операции (берут значения из пользовательских регистров X и Y).
Комбинаторика: nPr (перестановки) и nCr (сочетания), реализованные в виде циклов слоя скриптинга, позволяющих избавиться от ограничений диапазона факториалов благодаря внутреннему сокращению множителей.
Режим DISP: интерактивный формат дисплея и селектор точности разрядов, заменивший предыдущие клавиши FScEn и Digits. При нажатии SHIFT+EXP включается режим DISP, отображающий текущую настройку (например, «FIX 4»). Далее клавишей EXP циклически переключается формат (FIX → SCI → ENG → RAW), цифровые клавиши 0–9 задают точность, точка выбирает полную точность, а все остальные клавиши работают выполняют выход и исполняются обычным образом. С изменением настроек регистр Y меняется интерактивно, сразу обеспечивая визуальную обратную связь.
Попробуйте сами:
Обе версии работают в виде WebAssembly непосредственно в браузере.
Прерывания (потому что… почему бы и нет?)
Разумеется, простому непрограммируемому калькулятору не нужны прерывания. У него нет фоновых задач, требующих предварительного планирования. Основной цикл опрашивает клавиши, обновляет дисплей и ожидает. Для калькулятора это совершенно адекватная архитектура. Однако я всё равно добавил прерывания, просто потому, что их проектирование и реализация станут интересной задачей!
Пока обработчик прерывания в микрокоде содержит одну команду: RETI.
Реализация обладает приятными архитектурными свойствами:
Использование незадействованной кодировки безусловной команды. Флагу 15 я присвоил функцию отключения прерываний FLAG_IRQ_DIS. Каждая условная команда CPU кодирует своё условие в 5-битном поле: один бит отрицания и четыре бита выбора флагов. Сочетание negation=1, flag=15 (все единицы) считается «всегда true» — именно так RET, JMP и CALL кодируются в виде безусловных вариантов своих условных противоположностей. Однако кодирование negation=0, flag=15 не использовалось. Оно означает «true, когда установлен флаг 15». Я превратил его в RETI: условный возврат, исполняемый только при отключенных прерываниях, то есть именно в том состоянии, в котором находится CPU при выполнении обработчика прерываний. Никакого нового пространства опкодов не понадобилось; RETI умещается в уже существующий паттерн RETC, не требуя внесения изменений в декодер. Это просто алиас.
Атомарный return-and-reenable. RETI выполняет в одной команде три действия: извлекает адрес возврата из стека команд, сбрасывает FLAG_IRQ_DIS (заново включая прерывания), и устанавливает внутренний триггер irq_defer. Вход IRQ чувствителен к уровню: пока линия содержит высокий сигнал, флаг pending защёлкнут. Если RETI только выполнит возврат и снова включит прерывания, то при следующем получении команды калькулятор увидит всё ещё ожидающее прерывание и снова войдёт в обработчик, а прерванный код так никогда и не сможет выполниться. Триггер irq_defer предотвращает это: он заставляет CPU пропустить ровно одну проверку прерываний после RETI, гарантируя, что как минимум одна команда прерванной программы выполнится перед следующим прерыванием, избегая таким образом исчерпания команд. Именно так выполнят возврат из прерываний Z80.
Нагрузочное тестирование, доказавшее это. Привязываем вход IRQ к логической 1: калькулятор загружается, запускается и принимает нажатия клавиш, параллельно обрабатывая прерывание после каждой команды основного кода.
Безопасность команд из нескольких операций. Команды наподобие PUSH и POP, перемещающие несколько регистров, выполняют внутренние циклы в течение множества тактов, не возвращаясь в состояние получения команды. Так как прерывания проверяются только при получении команды, подобные многотактовые операции завершаются атомарно. PUSH 5 , сохраняющая пять регистров, завершит все пять сохранений до проверки прерываний. Кроме того, коррекция указателя стека после выполнения операции, которую некоторые инструкции откладывают до следующего цикла выборки, выполняется с защитой: проверка на прерывание в явном виде требует, чтобы на этот момент не оставалось незавершённой корректировки.
Полный механизм включает в себя чувствительную к уровню защёлку, управление на основе флага, триггер отложенного выполнения, атомарную инструкцию RETI, защиту для многооперационных инструкций, примерно 15 строк кода на Verilog для процессора и три псевдоинструкции ассемблера (DI, EI, RETI). По данным Quartus, всё это занимает 24 логические ячейки, то есть около 0,5 % ёмкости Cyclone II. На вход IRQ подаётся уже существующий 100-миллисекундный сигнал CTC таймера , благодаря чему обработчик прерываний получает «сердцебиение» с частотой 10 Гц. Пока что это периодическое событие никак не используется. Однако вся необходимая инфраструктура уже реализована, протестирована и готова к работе, а при цене всего в 24 логических ячейки она определённо того стоила. Я надеюсь задействовать её в одной из будущих реализаций.
Тулчейн
Тулчейн, благодаря которому всё это становится возможным, выглядит так: ассемблер (casm.py) и скриптовый компилятор (cscript.py). Оба представляют собой двухпроходные скрипты на Python 3, считывающие файлы исходников и выводящие совместимый с $readmemh шестнадцатеричный код, который Verilator и Quartus загружают непосредственно в виде содержимого ROM. Ассемблер поддерживает прямые ссылки, условный ассемблерный код, многоуровневые include файлов, локальные метки внутри подпрограмм и вычисление выражений. Компилятор скриптинга использует те же стандарты препроцессора и синтаксиса, но вместо 12-битных команд генерирует 4-битные токены. Вместе они имеют размер меньше 1,2 тысячи строк на Python; достаточно мало для полного понимания кода и достаточная мощь для того, чтобы решать все задачи.
Ситуация с тестированием
Точность новой версии арифметики характеризуется в документации Proto следующими классами допусков:
Операция | Гарантированная точность |
|---|---|
Сложение/вычитание/умножение/деление | 16,0 разрядов (0,5 ULP)* |
Квадратный корень | ~14,9 разряда |
Натуральный логарифм | ~14,5 разряда |
Возведение в степень | ~13-14 разрядов |
Тангенс (в градусах) | ~14-14,5 разряда |
Синус/косинус | ~13,5-14 разрядов |
*При вычитании почти равных чисел точность ухудшается пропорционально числу сократившихся значащих цифр; это не ошибка, а неотъемлемое свойство самой операции..
Сложно переоценить важность систематического тестирования относительно эталона. Наиболее наглядной демонстрацией этого стал баг FDIV Pentium 1994 года: отсутствие записи в таблице поиска, используемого алгоритмом деления с плавающей запятой Intel привело к тому, что для определённых входных данных процессор возвращал неверные результаты, и погрешность встречалась уже на четвёртом значимом разряде. Баг прятался в продаваемом оборудовании несколько месяцев, пока профессор математики Линчбергского колледжа Томас Найсли не обнаружил его, проверяя расхождения в вычислениях простых чисел на разных машинах. Поначалу Intel пыталась преуменьшить его важность, но в конце концов вынуждена была отозвать проданные процессоры, что стоило ей 475 миллионов долларов. Технической причиной бага были пять отсутствующих записей в таблице из 1066 элементов: небольшое упущение привело к очень серьёзным последствиям. Этот эпизод стал непосредственным толчком к применению формальных способов верификации аппаратных блоков для операций с числами с плавающей запятой, которые теперь стали стандартной практикой для Intel и всех других компаний.
Сквозные тесты (вычисление tan(atan(x)) и проверка совпадения значений x) — полезная мера защиты от ошибок знаков. ошибок квадрантов и багов на границах областей определения, но они не могут служить заменой прямого сравнения с эталонными значениями. В моём тестовом наборе используются оба подхода.
Что дальше
В нашей серии статей остался один пост. В заключении мы подведём итог всему, чему научились: технологическому стеку, полученным знаниям и тому, что можно улучшить.
Эталонная реализация BCD Proto находится в github.com/gdevic/Proto. Документация по алгоритмам, в том числе анализ точности, подробности методики CORDIC и сравнение алгоритмов можно найти в папке docs/.
