Я не вижу себя в качестве писателя. Но так случилось, что я столкнулся с прекрасным.

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

Если коротко: есть калькулятор МК-61 (1983). Внутри пять микросхем, которые соединены в однобитовую последовательную кольцевую шину: две памяти и три вычислителя со своей специализацией. Фактически они работают параллельно и синхронизируются через этот канал связи. Можно сказать, что это прообраз парадигмы NOC (network on chip) в современном железе применительно к FPGA.

Собственно о красоте. Чип ИК1303 (1980) отвечает за математические расчёты. В нём восемнадцать вычислительных операций: четыре бинарных (+ — * /) и четырнадцать F‑функций (10^x e^x lg ln asin acos atan sin cos tg sqrt x^2 x^y 1/x). Сверх них — обмен x<→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 = sign(z_i).

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

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

Дробь Гаусса [Гипергеометрическая функция, БРЭ]:

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

\tanh z=\cfrac{z}{1+\cfrac{z^{2}}{3+\cfrac{z^{2}}{5+\cfrac{z^{2}}{7+\dots}}}}

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

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

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

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

D_k=\frac{2k-1}{u}+s\,\frac{(k^p)^2}{D_{k+1}},\qquad 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 и знак. Фактически это два бита информации.
И ещё два — произвести небольшие преобразования до и после кофемолки.

Пример:

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 = post( GL(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» сидит арифметика вроде «поделить пополам», но там же и приведение аргумента, без которого дробь не сходится вовсе: тригонометрии нужен остаток по периоду π, логарифму — мантисса отдельно от десятичного порядка, экспоненте наоборот — целая часть уезжает в порядок результата, а дроби достаётся остаток. В post коде (void)x — это ровно те самые скрытые коррекции, a f это остаток x по ln10 f,n = x-n*ln10, floor(x / ln10). Все это я сознательно скрыл в таблице, дабы не засорять суть.

И вот что тут важно. Число уже лежит в кольце как мантисса и порядок по отдельности. Отделить их — способ чтения, а не работа.

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

Четыре клетки в микрокоде

Те же четыре клетки (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 заходит дважды, через log и через exp, а pow10 один раз через exp.

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

И привет CORDIC: ему нужна таблица atan(2^‑i), по константе на итерацию, и она растёт вместе с точностью. У цепной дроби таблицы нет вовсе. Все коэффициенты — 2k-1 и k² — арифметические, их порождает счётчик цикла.
Хранить нечего. Все компактненько.

Анализ

А почему 10 итераций? Тут интересный момент.
Большинство операций не требует 10 итераций и сходится за 6–8. Но есть atan(), которому нужно более 10, так что это число — компромисс по скорости и качеству.
И он тут жертва: atan(1) у нас равен 45.000002.

Посмотрим скорость сходимости для разных функций. Хорошо видно кто у нас хромой.

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

С другой стороны, есть другая ось ограничений — точность констант. Например, acos страдает от константы π/2, заданной как 1.5707963, а log10 — от ln10 (2.3025851).

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

Так что atan(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 не так ярко интересно, но кто знает…
Может, интересен анализ и картографирование функций? Граф связей?

Ссылки