Поиск
Тема: светлая

ИК1303. Выходящее за пределы опыта

Прокомментировать Просмотры: 3

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

Исследуя микрокод легендарного калькулятора МК-61, я не переставал поражаться: как инженеры прошлого умудрились вместить столь мощный функционал в крайне скромные аппаратные рамки, да еще при полном отсутствии аппаратных умножителей? Сегодня в мире целочисленных процессоров Cortex одна лишь операция с плавающей запятой способна раздуть размер прошивки до критических значений (пожалуй, я утрирую, но эмоциональный отклик именно таков).

Если обрисовать картину вкратце: перед нами ретро-калькулятор МК-61 образца 1983 года. Его архитектура базируется на пяти микросхемах, объединенных в однобитную последовательную кольцевую шину. Два чипа отведены под память, три — выполняют специализированные вычислительные задачи. Фактически они функционируют параллельно, синхронизируясь посредством этого коммуникационного канала. По сути, перед нами исторический прообраз современной парадигмы NoC (Network-on-Chip) применительно к сфере FPGA.

Теперь о самой эстетике решения. За математические вычисления отвечает микросхема ИК1303 (1980 года выпуска). В ее арсенале заложено восемнадцать операций: четыре базовых бинарных (+, -, *, /) и четырнадцать так называемых F-функций ($10^x, e^x, \lg, \ln, \arcsin, \arccos, \operatorname{arctg}, \sin, \cos, \operatorname{tg}, \sqrt{x}, x^2, x^y, 1/x$). Дополняют этот набор команды обмена регистрами $x \longleftrightarrow y$ и служебные рутины, скрытые от пользователя: нормализация, генерация констант, редукция аргумента.

Никаких алгоритмов CORDIC, никакой двоичной арифметики с плавающей точкой! Но каким образом достигнута подобная эффективность? Все невероятно изящно: подавляющее большинство расчетов завязано на одну-единственную функцию! Всего на ОДНУ!

В роли этого фундаментального ядра выступает аппаратная непрерывная (цепная) дробь Гауска-Ламберта.

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

y = GL(x, p, s):     # Ядро микрокода srom[0F | 81, A0, A1, B0, B1 | 71..7B | BC..BF]
    dy, y0 = 2.0/x, 1.0/x
    y = 10*dy - y0
    for k in (9,8,..,1):
        n = k**p
        y = k*dy - y0 + s*n*n/y
    return y

p = 1  ->  числители k^2   (дробь Гаусса)
p = 0  ->  числители 1     (дробь Ламберта)

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

Для наглядности сравним с аналогами: архитектура HP-35 на каждом шаге вычислений прогоняет цикл циклических вычитаний до смены знака, тогда как алгоритм CORDIC на каждой итерации вычисляет знак поворота $d_i = \operatorname{sign}(zi)$.

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

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

Дробь Гаусса [<a href=»https://bre.ruwiki.ru/wiki/%D0%93%D0%B8%D0%BF%D0%B5%D1%80%D0%B3%D0%B5%D0%BE%D0%BC%D0%B5%D1%82%D1%80%D0%B8%D1%87%D0%B5%D1%81%D0%BA%D0%B0%D1%8F%D1%84%D1%83%D0%BD%D0%BA%D1%86%D0%B8%D1%8F» rel=»noopener nofollow»>Гипергеометрическая функция, БРЭ]:

₀F₁ (p=0, s=+1):

ИК1303. Выходящее за пределы опыта

₂F₁ (p=1, s=+1):

\arctan z=\cfrac{z}{1+\cfrac{(1z)^{2}}{3+\cfrac{(2z)^{2}}{5+\cfrac{(3z)^{2}}{7+\dots}}}}

Цепные дроби Ламберта формируют классические разложения прямых тригонометрических и гиперболических функций (тангенса), в то время как общая гауссова теория непрерывных дробей охватывает их через отношения гипергеометрических рядов, покрывая обратные функции вроде арктангенса и ареатангенса, а через последний — и натуральный логарифм. Частные аналитические случаи сводятся к двум семействам с числителями k^0 и k^2.

Четыре возможных конфигурации параметров $(p, s)$ задают четыре классические цепные дроби, исследованные еще Ламбертом (1761) и Гауссом (1812). Их канонический вид универсален:

D_k=\frac{2k-1}{u}+s\,\frac{(k^p)^2}{D_{10}}=\frac{19}{u}

Рекурсия обрывается на элементе $19/u$, что как раз соответствует тем самым 10 уровням, жестко закодированным в микропрограмме.

Параметр $p=0$ формирует дроби Ламберта для прямых отображений, а $p=1$ — дроби Гаусса для обратных. Переключатель знака $s$ позволяет переходить от тригонометрических функций к гиперболическим:

\operatorname{ctg}u=D_1\;(p{=}0,s{=}{-}1)\qquad \operatorname{th}w=\frac{1}{D_1}\;(p{=}0,s{=}{+}1)\operatorname{arctg}u=\frac{1}{D_1}\;(p{=}1,s{=}{+}1)\qquad \operatorname{arth}z=\frac{1}{D_1}\;(p{=}1,s{=}{-}1)

Дальнейшие тождества выводятся естественным образом:

\ln\frac{1+z}{1-z}=2\operatorname{arth}z=\frac{2}{D_1},\qquad e^{2w}=\frac{1+\operatorname{th}w}{1-\operatorname{th}w}=1+\frac{2}{D_1-1}

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

Вернёмся к имплементации

Практическое применение этого механизма интуитивно понятно и элегантно. Входной аргумент требует предварительной подготовки, а результат работы вычислительного ядра — постпроцессинга. Иными словами, специфика конкретной функции задается в четырех точках:
Две точки определяют константы ($p$ и знак $s$), что эквивалентно всего двум битам информации.
Еще две точки отвечают за простейшие преобразования до и после вызова расчетного ядра.

Пример:

atan(x) = 1/GL(x, 1, 1)

Если задействовать следующий набор модификаторов (включая функции post_ln и post_exp, выведенные из теории выше):

double pre_id(double x)    { return x; }
double pre_half(double x)  { return x / 2.0; }                 // exp
double pre_ln(double x)    { return (x - 1.0) / (x + 1.0); }   // ln
double pre_asin(double x)  { return x / sqrt(1.0 - x * x); }   // asin


double post_inv(double D, double x)  { (void)x; return 1.0 / D; }
double post_ln(double D, double x)   { (void)x; return 2.0 / D; }
double post_exp(double D, double x)  { (void)x; return 2.0 / (D - 1.0) + 1.0; }
double post_sin(double D, double x)  { (void)x; return 1.0 / sqrt(1.0 + D  D); }
double post_cos(double D, double x)  { (void)x; return   D / sqrt(1.0 + D  D); }

То через обобщенную схему $y = \text{post}(\text{GL}(\text{pre}(x)), x)$ мы получаем прямую или косвенную реализацию широкого спектра математических функций:

функция

pre(x)

p

s

post(D, x)

atan

x

1

+1

1/D

asin

x/sqrt(1-x²)

1

+1

1/D

acos

x/sqrt(1-x²)

1

+1

π/2 — 1/D

tan

rem(x, π)

0

-1

1/D

sin

rem(x, π)

0

-1

flip(x)·sgn(rem)/sqrt(1+D²)

cos

rem(x, π)

0

-1

flip(x)·abs(D)/sqrt(1+D²)

ln

(m-1)/(m+1)

1

-1

e·ln10 + 2/D

exp

f/2

0

+1

10ⁿ·(1 + 2/(D-1))

Приведение аргумента машине достаётся даром

На этапе предварительной подготовки (pre) выполняются простейшие операции вроде деления пополам, но здесь же происходит и жизненно необходимая редукция аргумента, без которой дробь попросту не сойдется. Тригонометрическим функциям требуется остаток от деления на период $\pi$, логарифмам — отделение мантиссы от десятичного порядка, а экспоненте, напротив, отделение целой части, уходящей в порядок результата, тогда как дробной части достается остаток. В коде постобработки конструкция (void)x как раз символизирует те самые скрытые коррекции, а f обозначает остаток $x$ по модулю $\ln 10$: $f, n = x — n \cdot \ln 10, \lfloor x / \ln 10 \rfloor$. Все эти детали я намеренно опустил в сводной таблице, чтобы не загромыхать общую картину.

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

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

Четыре базовые точки в микрокоде

Те же самые четыре конфигурационные точки $(p, s)$, представленные в виде конкретных адресов микропрограммы emu145:

ctg    srom[13]   p=0 s=-1   Ламберт, круговая прямая          -> sin cos tg
exp    srom[6C]   p=0 s=+1   Ламберт, гиперболическая прямая   -> e^x 10^x x^y
atan1  srom[B7]   p=1 s=+1   Гаусс,   круговая обратная        -> atan asin acos
log    srom[E7]   p=1 s=-1   Гаусс,   гиперболическая обратная -> ln lg

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

Четыре опорных блока:

ctg   - дробь с единичными числителями, s=-1
atan1 - дробь с числителями k^2,        s=+1
log   - дробь с числителями k^2,        s=-1
exp   - дробь с единичными числителями, s=+1

Через них реализуются еще 11 функций:

через ctg   (3):  sin, cos, tan
через atan1 (3):  atan, asin, acos
через log   (2):  log, log10
через exp   (3):  exp, pow10, pow

При этом операция возведения в степень (pow) задействуется дважды (через логарифм и экспоненту), а pow10 — единожды через экспоненту.

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

Стоит вспомнить и про алгоритм CORDIC: ему необходима громоздкая таблица констант $\operatorname{arctg}(2^{-i})$, размер которой масштабируется вместе с требуемой точностью. Цепная дробь вовсе не нуждается в таблицах. Все коэффициенты — последовательности $2k-1$ и $k^2$ — генерируются «на лету» арифметическим счетчиком цикла.
Хранить в памяти просто нечего. Всё исключительно компактно.

Анализ

Почему было выбрано именно 10 итераций? Здесь кроется любопытный нюанс.
Большинству функций вовсе не требуется 10 шагов — они успешно сходятся за 6–8 итераций. Исключением является арктангенс, которому десяти шагов впритык, так что выбранное число представляет собой компромисс между скоростью и точностью.
И именно арктангенс становится здесь «жертвой»: значение $\operatorname{arctg}(1)$ на выходе выдает 45.000002.

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

Сходимость цепной дроби по числу уровней
Сходимость цепной дроби по числу уровней

С другой стороны, существует иное ограничение — разрядность самих констант. К примеру, арккосинус страдает от точности константы $\pi/2$, записанной как 1.5707963, а десятичный логарифм — от точности $\ln 10$ (2.3025851).

То же, но с реальными восьмиразрядными константами машины
То же, но с реальными восьмиразрядными константами машины

Таким образом, результат $\operatorname{arctg}(1) = 45.000002$ — это не программный баг, а естественная плата за 8-разрядную точность железа.

Ложка дёгтя

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

The choice of algorithms for the HP-35 received considerable thought. Power series, polynomial expansions, continued fractions, and Chebyshev polynomials were all considered for the transcendental functions. All were too slow because of the number of multiplications and divisions required to maintain full ten‑digit accuracy.

Специалисты HP ориентировались на точность в десять знаков, и цепные дроби в этот лимит не вписывались. Советский МК-61 оперирует восемью знаками — и здесь этот подход работает отлично.
Имея дело с одной и той же математикой, разработчики приняли разные решения, поскольку перед ними стояли принципиально разные технические задания.

Резюме

  • Жестко заданные 10 итераций, хотя большинству функций с лихвой хватает пяти-шести.

  • Рекурсия обрывается на элементе $19/u$, причем последний уровень вычисляется целиком, без дробного продолжения.

  • Умножение реализовано посредством многократных сложений.

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

Послесловие

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

Выражаю искреннюю признательность Феликсу Лазареву за создание эмулятора emu145.
Отдельная благодарность участникам сообществ МК61 / МК52 / MK85 и лично Сугоняеву за терпение 😉
И, конечно же, мое глубочайшее восхищение адресовано Эйлеру, Ламберту и Гауссу. Без их фундаментального вклада ничего подобного бы не появилось. Не зря У. Джоунс и В. Трон прослеживают именно эту историческую преемственность:
Эйлер → Ламберт → Лагранж → Гаусс → гипергеометрические цепные дроби.

Надеюсь, эта тема покажется увлекательной не мне одному. Боюсь, что в чипах 1302/1306 всё устроено не столь элегантно, но кто знает…
Быть может, читателям будет интересен более детальный анализ, картографирование функций или построение графа связей?

Ссылки

 

Источник

Поделиться:

Похожие статьи

Comments - Cackle

Поиск по играм, новостям и статьям…

Введите не менее двух символов

Введите не менее двух символов