Ускоряем программу для 50-летнего процессора на 180000% / forpes.ru

Главная
Ускоряем программу для 50-летнего процессора на 180000%

Ускоряем программу для 50-летнего процессора на 180000% +245

14.11.2023 08:28

mark_ablov 34 35000 Источник

Введение

В прошлом году я написал программу, вычисляющую 255 цифр числа π на самом первом микропроцессоре от Intel - 4004. В той статье я упоминал рекорд ENIAC'a - 2035 цифр [^1], но побить его не смог. Настало время закрыть гештальт. В этот раз возьмём одного из преемников от Intel - 4040. Определим рамки, в которых будем разрабатывать наш код:

Основная цель это вычислить как минимум 2035 первых цифр числа π быстрее, чем это сделал ENIAC. В исходной работе с описанием данного результата [^2] озвучены 70 часов машинного времени. Значит нам нужно уложиться в 69 часов, 59 минут и 59 секунд.
Можем использовать только память, адресуемую напрямую процессором. То есть никаких внешних переключателей банков RAM/ROM.
Избегаем каких-либо предвычисленных значений, сохраненных в ПЗУ. Все вычисления должны производиться i4040. К примеру, мы не можем рассчитать таблицу логарифмов на ПК и включить её как константную секцию в код, а должны построить её на лету.
Никакой подстройки под конкретное целевое значение числа вычисляемых цифр π. Да, мы хотим вывести 2035+ цифр, но код должен корректно работать и для меньших значений.

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

Технические ограничения

Несмотря на то, что 4040 вышел на несколько лет позже чем предшественник, почти все ограничения 4004 сохранились:

Тактовая частота по-прежнему 740 килогерц. Одна инструкция исполняется за 8 или 16 тактов, поэтому у нас не более 92500 инструкций в секунду.
Набор инструкций (ISA) несколько расширился, но именно что "несколько". Всё еще нет ни умножения, ни деления, ни эффективных операций с памятью. Чего уж там, даже сдвиговые инструкции отсутствуют.
Размер ОЗУ не увеличился - 1280 байт. Да-да, чуть больше килобайта. Во многих источниках указана неверная информация про меньший объём, но существовали реальные системы [^3] с 1280 байтами на борту (адресуемых напрямую, как нам и надо).

Тогда что же нам принёс полезного 4040 [^11]? Поддержку внешних прерываний, пошаговую отладку, меньшее энергопотребление в некоторых случаях, но для нашей практической задачи более интересен следующий список:

Добавились инструкции AND/OR. И это именно нововведение, потому что 4004 не имел такой роскоши. Но не всё так радужно - бинарные операции не имеют аргументов, а выполняют побитовую дизъюнкцию/конъюнкцию над конкретными регистрами. К примеру, AN6 совершает побитовое ИЛИ над регистром rr6 и аккумулятором.
Стек вызова вырос с 3-х уровней до 7. Теперь переполнение стека менее вероятно, хотя я с этим всё же столкнулся.
Регистров стало больше - 24 вместо 16. Но ISA не поменялась и 4040 имеет бинарную совместимость машинного кода с 4004. То есть индекс регистра всё еще закодирован 4-мя битами. Каким же образом мы можем адресовать дополнительные регистры? Не самым эффективным образом - переключая регистровые файлы (банки) туда-обратно с помощью инструкций SB0 / SB1
К счастью, поддерживаемый объем ПЗУ увеличился. Добавился второй банк ROM с дополнительными 4 кибибайтами. Тем не менее переключение между банками непрозрачное и требует аккуратного тайминга в коде и использования инструкций DB0 / DB1.

Теоретическая подготовка

Выбор алгоритма

Давайте рассмотрим какие варианты у нас есть:

Краниковый алгоритм [^4]. Я описывал его в прошлой статье. Если вкратце, то алгоритм работает с массивом цифр, хранимых в ОЗУ. Размер этого массива пропорционален количеству вычисляемых цифр. Нам просто не хватит RAM для 2000+ цифр.
Формулы Мэчина. С вычислительной точки зрения вполне нас устраивают, но даже для таких простых формул нам не хватит памяти. Дело в том, что нам нужно хранить как минимум два длинных числа (разрядностью не менее чем целевое количество вычисляемых цифр π). Никак не упаковать 4000+ цифр в 1280 байт.
Улучшенный краниковый алгоритм [^5]. В нём не используется какой-либо массив, но, так же как и в формулах Мэчина, мы работаем с весьма длинными числами. Даже для 255 цифр необходимый объем памяти уже 7 кибибайт.
Формула от Bailey–Borwein–Plouffe [^6]. Весьма интересная идея, но по этой формуле мы получаем число π в 16-ричной или 2-чной системах счисления. Я всё же ориентирован на 10-ричную.
Алгоритм Плуффа [^7]. Очень медленный. С такой асимптотикой будет очень сложно уложиться в заданное время (скорее всего и невозможно).
Алгоритм от Фабриса Беллара, вдохновенный работой Плуффа [^8]. Не самый быстрый, но относительно простой и не требующий сложной математики.
Алгоритм Ксавье Гурдона [^9]. Очень быстрый, но весьма сложный и объемный. Для его реализации скорее всего не хватит ПЗУ, ведь у нас вся программа не может превышать объема в 8000 инструкций, а нам придётся написать все базовые математические функции.

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

Исходная формула для π

Плуфф (а затем и Беллар) используют весьма занятную формулу для π:

$\pi+3=\sum_{k=1}^{+\infty}\frac{k2^{k}}{\left(\begin{array}{c}2k\\ k\end{array}\right)}$

Вывод этой формулы (как и несколько других увлекательных равенств) можно найти в статье Лемера. [^10]

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

$\left(\begin{array}{c}2n\\ n\end{array}\right)=\frac{2^n(2n-1)!!}{n!},\:n!!=n \cdot (n-2) \cdots 3 \cdot 1$

Получаем окончательную форму числового ряда, с которым будем работать:

$\pi+3=\sum_{k=1}^{+\infty}\frac{k2^{k}}{\left(\begin{array}{c}2k\\ k\end{array}\right)}=\sum_{k=1}^{+\infty}\frac{k \cdot k!}{(2k-1)!!}$

Вычисление n-ой цифры дроби

Начнём с упрощеннго примера. Пусть у нас есть дробь $s = \frac{b}{A}$ и мы хотим узнать цифру на позиции этого числа. Как мы это можем сделать?

Для начала запишем факторизацию :

$A=\prod_{i=1}^{m}a_{i},\:i \neq j \Rightarrow a_i \perp a_j$

Затем определим несколько переменных:

$b_i = b \bmod a_i;\ A_i = \frac{A}{a_i};\ A_i^{-1} = \frac{1}{A_i} \bmod a_i$

После применения китайской теоремы об остатках мы получим:

$b' = \left( \sum_{i=1}^{m} b_i A_i A_i^{-1} \right) \bmod{A}$

Введём еще одну переменную $f_i = b_i \cdot A_i^{-1} \bmod a_i$ :

$b' = \left( \sum_{i=1}^{m} f_i A_i \right) \bmod{A}$

Так как $b' = b \bmod A$ , то мы можем использовать эту подстановку для вычисления дробной части :

$\{s\} = \frac{b'}{A} = \frac{\left( \sum_{i=1}^{m} \frac{f_i}{a_i} A \right) \bmod{A}}{A} \stackrel{(1)}{=} \left\{ \frac{\sum_{i=1}^{m} \frac{f_i}{a_i} A}{A} \right\} = \left\{ \sum_{i=1}^{m} \frac{f_i}{a_i} \right\}$

Возможно, тождество (1) очевидно, но для полноты статьи приведу небольшое доказательство:

$\frac{x \bmod N}{N} = \frac{x - N \lfloor \frac{x}{N} \rfloor}{N} = \frac{x}{N} - \lfloor \frac{x}{N} \rfloor = \{\frac{x}{N}\}$

Теперь мы можем записать формулу для вычисления цифр числа с позиции:

$D_n = \{s \cdot 10^n\} = \left\{ \sum_{i=1}^m \frac {f_i \cdot 10^n}{a_i} \right\} \stackrel{(2)}{=} \left\{ \sum_{i=1}^m \frac {f_i \cdot 10^n \bmod {a_i}}{a_i} \right\}$

Переход (2) нужен только для упрощения вычислений: мы стараемся избежать больших чисел (а может быть весьма большим числом). И так как мы ищем только дробную часть числа, то можем перейти к сложению остатков и модульной арифметике.

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

Вычисление n-ой цифры числового ряда

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

Пусть будет множеством простых чисел, которые встречаются в разложении для какого-либо. При этом для пары кратность простого числа в факторизации может быть 0 (это конкретное число не является множителем конкретного ). Также объявим $v_i^{max}$ - максимальная кратность для среди всех разложений . Суммируя все вводные, можем определить числовой ряд:

$S_N=\sum_{k=1}^N\frac{b_k}{A_k}, A_k=\prod_{i=1}^{m}a_i^{v_{i,k}}, 0 \leq v_{i,k} \leq v_i^{max},\:i \neq j \Rightarrow a_i \perp a_j, v_{i,k} > 0 \Rightarrow a_i \nmid b_k$

По схеме из прошлого раздела, объявим несколько переменных:

$A_{i,k} = \frac {A_k}{a_i^{v_{i,k}}}\\ A_{i,k}^{-1} = \frac{1}{A_{i,k}} \bmod a_i^{v_{i,k}}\\ f_{i,k} = \begin{cases} 0 &\text{if }v_{i,k}=0\\ b_k \cdot A_{i,k}^{-1} \bmod a_i^{v_{i,k}} & \text{if }v_{i,k} > 0\\ \end{cases}$

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

$\frac{b_k}{A_k} = \left\{ \sum_{i=1}^{m} \frac{f_{i,k}}{a_i^{v_{i,k}}} \right\} \\[15px]S_N = \left\{ \sum_{k=1}^N \sum_{i=1}^{m} \frac{f_{i,k}}{a_i^{v_{i,k}}} \right\} = \left\{ \sum_{i=1}^{m} \sum_{k=1}^N \frac{f_{i,k}}{a_i^{v_{i,k}}} \right\}$

Следующим шагом перейдём к использованию общего делителя $a_i^{v_i^{max}}$ , который не зависит от :

$f'_{i,k} = f_{i,k} \cdot a_i^{v_i^{max}-v_{i,k}} \\[15px] S_N = \left\{ \sum_{i=1}^{m} \sum_{k=1}^N \frac{f'_{i,k}}{a_i^{v_{i,k}} \cdot a_i^{v_i^{max}-v_{i,k}}} \right\} = \left\{ \sum_{i=1}^{m} \sum_{k=1}^N \frac{f'_{i,k}}{a_i^{v_i^{max}}} \right\} \\[15px] S_N = \left\{ \sum_{i=1}^m \frac{f_i}{a_i^{v_i^{max}}}\right\}, f_i=\left( \sum_{k=1}^N {f'_{i,k}} \right) \bmod a_i^{v_i^{max}}$

Применяем тот же трюк, что и для дробей - так как целая часть нас не интересует, то мы можем использовать не такие большие числа из класса вычетов по модулю $a_i^{v_i^{max}}$ .

Важное наблюдение - мы можем заметить, что $f'_{i,k}$ так же принадлежит к тому же классу вычетов по модулю $a_i^{v_i^{max}}$ . Это не столь очевидно, но сейчас мы это покажем, используя упрощенный пример:

$(x \bmod n^p) \cdot n^{q-p} = (x - n^p \cdot \lfloor \frac{x}{n^p} \rfloor) \cdot n^{q-p} = x \cdot n^{q-p} - n^q \cdot \lfloor \frac{x}{n^p} \rfloor \\(x \cdot n^{q-p}) \bmod n^q = x \cdot n^{q-p} - n^q \cdot \lfloor \frac{x \cdot n^{q-p}}{n^q} \rfloor = x \cdot n^{q-p} - n^q \cdot \lfloor \frac{x}{n^p} \rfloor$

Видно, что два выражения эквивалентны, и значит мы можем записать финальное выражение для :

$f'_{i,k} = b_k \cdot A_{i,k}^{-1} \bmod a_i^{v_{i,k}} \cdot a_i^{v_i^{max}-v_{i,k}} = b_k \cdot A_{i,k}^{-1} \cdot a_i^{v_i^{max}-v_{i,k}} \bmod a_i^{v_i^{max}}\\[15px]f_i=\sum_{k=1}^Nb_k\left(\frac{A_k}{a_i^{v_{i,k}}}\right)^{-1}a_i^{v_i^{max}-v_{i,k}}\bmod a_i^{v_i^{max}}$

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

$D_n=\left\{ \sum_{i=1}^m\frac{f_i \cdot 10^n \bmod a_i^{v_i^{max}}}{a_i^{v_i^{max}}} \right\}$

Мы можем записать значения , и для формулы π:

$b_k = k! \\A_k = (2k-1)!!\\v_{i,k} = v_{a_i}({A_k})\\u_{i,k} = v_{a_i}({b_k})\\f_i=\sum_{k=1}^N k\cdot \left( \frac{b_k}{a_i^{u_{i,k}}} \right) \cdot \left( \frac{A_k}{a_i^{v_{i,k}}} \right)^{-1} \cdot a_i^{v_i^{max}-v_{i,k}+u_{i,k}} \bmod a_i^{v_i^{max}}$

Функция $v_{a_i}(x)$ - степень вхождения простого числа в разложение .

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

Наибольшая кратность простого числа в факторизации делителя

Мы определили почти все переменные, которые нам нужны, за исключением $v_i^{max}$ . У нас была вводная, что не является множителем , но это не так для нашей формулы, ибо очевидно имеет в множителях . К счастью, это имеет влияние только на значание $v_i^{max}$ . Для нашей формулы мы можем записать степень вхождения как:

$v_{a_i}\left(\frac{A_k}{b_k} \right) = v_{a_i}\left(\frac{(2N - 1)!!}{N!}\right)$

Раскроем двойной факториал и используем формулу Лежандра:

$\begin{align*} v_{a_i}\left(\frac{(2N - 1)!!}{N!}\right) & = v_{a_i}\left( \frac{(2N)!}{2^N \cdot N! \cdot N!} \right) \\ & \stackrel{(3)}{=} v_{a_i}\left( \frac{(2N)!}{N! \cdot N!} \right) \\ & = \sum_{t=1}^L \lfloor \frac{2N}{a_i^t} \rfloor - 2\sum_{t=1}^L \lfloor \frac{N}{a_i^t} \rfloor \\ & = \sum_{t=1}^L \left( \lfloor \frac{2N}{a_i^t} \rfloor - 2 \lfloor \frac{N}{a_i^t} \rfloor \right) \end{align*}$

В (3) мы просто исключили потому что этот член делится только на 2, а среди у нас нет 2 (по определнию двойного факториала имеет только нечётные множители), а значит он никак не влияет на функцию кратности. Теперь взглянем на каждый член суммы и отметим, что он равен либо 0, либо 1:

$\lfloor \frac{2N}{a_i^t} \rfloor - 2 \lfloor \frac{N}{a_i^t} \rfloor = \lfloor \frac{N}{a_i^t} \rfloor + \lfloor \frac{N}{a_i^t} + \frac{1}{2} \rfloor - 2 \lfloor \frac{N}{a_i^t} \rfloor = \lfloor \frac{N}{a_i^t} + \frac{1}{2} \rfloor - \lfloor \frac{N}{a_i^t} \rfloor \leq 1$

Так как $L = \lfloor log_{a_i}{2N} \rfloor$ и, даже если каждый член суммы будет равен 1, то максимальное значение суммы будет огранично сверху $v_i^{max} \leq \lfloor log_{a_i}{2N} \rfloor$ .

Определение достаточной верхней границы для частичной суммы числового ряда π

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

$\frac{b_M}{A_M} \cdot 10^n < 1$

Умножая на , мы сдвигаем десятичную точку на позиций вправо. Цифры слева от точки уже более-менее корректны, и значит сильное увеличение не будет оказывать влияния на них, и можно принять их за финальные. Поэтому находим с какого целая часть произведения примерно определена, и добавляем некоторый $\epsilon$ в качестве запаса на точность.

Чтобы решить это неравенство, используем другую формулу центрального биномиального коэффициента:

$\left(\begin{array}{c}2n\\ n\end{array}\right)=\frac{4^n(2n-1)!!}{(2n)!!}$

Трансформируем числовой ряд π, используя формулу выше:

$\frac{M \cdot 2^M}{\left(\begin{array}{c}2M\\ M\end{array}\right)} \cdot 10^n = \frac{M \cdot 2^M}{\left( \frac{4^M \cdot (2M - 1)!!}{2M!!} \right)} \cdot 10^n = \frac{M \cdot 2M!!}{2^M \cdot (2M - 1)!!} \cdot 10^n$

Исследуем часть этого выражения и определим его границы:

$1 < \frac{2M!!}{(2M - 1)!!} \leq M$

Исходя из определения двойного факториала, достаточно очевидно почему у нас есть верхняя граница в виде , когда $M \geq 4$ . Теперь вернёмся к основному неравенству (заменив кусочек сверху на $M^t, 0 < t \leq 1$ ):

$1 < r \leq 2 \\[15px]\frac{M^{r} \cdot 10^n}{2^M} < 1$

Проведя серию нехитрых импликаций, получим искомую оценку для :

$\begin{align*} \frac{M^r \cdot 10^n}{2^M} < 1 & \Rightarrow M^r \cdot 10^n < 2^M \\ & \Rightarrow log_{2}{M^r} + log_{2}{10^n} < M \\ & \Rightarrow r \cdot log_{2}{M} + n \cdot log_{2}{10} < M \end{align*}$

Можно, конечно, продожить, но так как это всё равно только приблизительное число и мы в любом случае будем добавлять какой-то $\epsilon$ для обеспечения точности вычислений, то можно просто докинуть $r \cdot log_{2}{M}$ к этому $\epsilon$ . Для наших не таких грандиозных целей ( $n \leq 2048$ ), $\epsilon$ может быть весьма небольшим. Я остановился на 70.

Алгоритм

Крутим два цикла: внешний цикл по простым числам (), а внутренний служит для вычисления . Сразу же на уровне алгоритма закладываем небольшую оптимизацию - вместо вычисления факториалов на каждом шаге цикла, мы можем использовать реккурентные формулы: $A_k = (A_{k-1} \cdot (2k - 1)) \bmod m$ и $b_k = (b_{k-1} \cdot k) \bmod m$ . Кроме этого, мы сразу делим и на текущее простое число и одновременно обновляем значения $v_{i, k}$ и $u_{i, k}$ , которые будут использоваться для возведения в нужную степень.

$\begin{align*}&\begin{aligned}&N&&=\lfloor(n+\epsilon) \cdot log_{2}{10}\rfloor\\&D_n&&=0\\\end{aligned}\\&for\:(a \in \mathbb{P} | 2 < a <2N)\:\{\\&\quad\begin{aligned}&vmax &&=\lfloor \log_{a}{2N} \rfloor \\&m&&=a^{vmax}\\&v&&=0\\&u&&=0\\&f&&=0\\&A&&=1\\&b&&=1\\\end{aligned}\\&\quad for\:(k=2..N)\:\{\\&\quad\quad\begin{aligned}&b &&= \left( \frac{k}{a^{v(a,k)}} \cdot b \right) \bmod m\\&A &&= \left( \frac{2k-1}{a^{v(a,2k-1)}} \cdot A \right) \bmod m\\&u &&= u + v(a, k)\\&v &&= v + v(a, 2k-1)\\\end{aligned}\\&\quad\quad if\:(v - u > 0)\:\{\\&\quad\quad\quad f = (f + k \cdot b \cdot A^{-1} \cdot a^{vmax-v+u}) \bmod m\\&\quad\quad \}\\&\quad \}\\&\quad\begin{aligned}&d&&=(f \cdot 10^n) \bmod m\\&D&&=\left\{ D + \frac{d}{m} \right\} \\\end{aligned}\\&\}\end{align*}$

Прежде чем начать писать код на ассемблере, я разработал референсную версию алгоритма на JS:

const getNinePiDigits = (n) => {
  const N = Math.trunc((n + 20) * Math.log(10) / Math.log(2));

  let digits = 0;
  for (const a of primes.filter((prime) => prime > 2 && prime < 2 * N)) {
    const vmax = Math.trunc(Math.log(2 * N) / Math.log(a));
    const m = a ** vmax;

    let f = 0, A = 1, b = 1, v = 0, u = 0;
    for (let k = 2; k <= N; k++) {
      b = (b * (k / (a ** multiplicity(a, k)))) % m;
      A = (A * ((2 * k - 1) / (a ** multiplicity(a, 2 * k - 1)))) % m;
      u += multiplicity(a, k);
      v += multiplicity(a, 2 * k - 1);
      if (v - u > 0) {
        f = (f + (k * b * modularInverse(A, m) * (a ** (vmax - v + u)))) % m;
      }
    }

    d = (f * modularPower(10, n, m)) % m;
    digits = (digits + (d / m)) % 1;
  }

  return Math.trunc(digits * 1e9).toString().padStart(9, '0');
};

Инструментарий для разработки

Для более-менее простых програм (таких как прошлый проект) мне не требовалось что-то большее, чем простой компилятор ассемблера в машинный код. Но разработка текущего проекта становилась всё сложнее. Так что мне пришлось расширить набор инструментов.

Препроцессор

От препроцессора я хотел получить две возможности:

поддержка модульности - до этого весь код содержался в одном файле
примитивные макросы

К примеру, есть у нас модуль, который реализует функцию foo (foo.i4040):

%define FOO_PARAM1_VAL1 0x5
%define FOO_PARAM1_VAL2 0xA
%define FOO_PARAM2_VAL1 D

foo:
  ...
  BBL 0

А в другом модуле мы хотим использовать эту функцию:

%include "foo.i4040"

bar:
  FIM r0, $FOO_PARAM1_VAL1 . $FOO_PARAM2_VAL1
  JMS foo
  FIM r0, $FOO_PARAM1_VAL2 . 0
  JMS foo
  BBL 0

Макросы оказались очень удобны для работы с переменными в памяти. Я несколько раз пересматривал структуру размещения переменных в RAM (скажем, переменная, в которой хранилось значение , могла сменить регистр RAM с #5 на #9). Без этой фичи, мне бы пришлось искать в исходном коде конкретное число и думать о том, относится ли оно к операции с памятью или просто константа.

Компилятор и компоновщик

Первая версия компилятора была бесхитростной - берём инструкции одну за одной, транслируем их в машинный код и записываем в итоговый образ ПЗУ. Заключительным этапом подставляем конкретные адреса в инструкции перехода по меткам. Мы не могли сделать этого при первом проходе, потому что метка могла быть объявлена позже, чем инструкция перехода, которая ссылается на эту метку. Однако архитектура 4004/4040 имеет несколько особенностей связанных с переходами и структурой машинного кода:

Некоторые команды перехода могут работать с адресами, только внутри той же страницы ПЗУ (размером 0x100 байт), где и размещена сама инструкция. Например, инструкция по адресу 0x2AA может совершить переход на 0x255, но не на 0x655.
Такие короткие переходы имеют ещё одну особенность - если они расположены на последних байтах страницы ПЗУ, то переход будет выполнен не на адрес внутри страницы с этой инструкцией, а на адрес со следующей страницы. Инструкция по адресу 0x7FF или 0x7FE, выполняющая короткий переход на адрес 0x7BB неожидано передаст управление на 0x8BB.

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

Кроме того, накопилось еще несколько пожеланий к компилятору:

Поддержка двух банков ПЗУ. Программа разбухала, и я опасался, что не уложусь в 4 КиБ.
Необходимость размещения кода по конкретным адресам.

Остановлюсь на этих требованиях подробнее.

Как я уже ранее упоминал, переход между банками ROM не совсем прозрачен для программиста, и требуется самостоятельно переключаться между ними:

__location(0x0:0x00)
entrypoint:
  JMS output
  DB1
  JMS foo
  HLT

output:
  LDM 0xF
  WMP
  BBL 0

__rom_bank(0x1)
foo:
  JMS output
  DB0
  NOP
  BBL 0

В этом листинге можно наблюдать следующее:

entrypoint располагается в первом банке, по адресу 0x000. И отсюда процессор начинает исполнять код после сигнала сброса.
foo будет размещен во втором банке по адресу, который выберет компоновщик.
функция output используется кодом в обоих банках, а значит необходимо её продублировать.
инструкция NOP перед командой возврата BBL необходима для соблюдения нужных таймингов. DBx переключает банк ПЗУ через 2 цикла, так как предполагается использование этой инструкции в связке с командами перехода, которые как раз и исполняются за 2 командных цикла. В отличии от них, BBL занимает только один цикл, так что нам нужно потратить дополнительный цикл для корректного перехода в нужный банк во время возврата из функции.

Мы хотим задавать конкретный адрес размещения кода не только для точки входа, но и для целей оптимизации - создания эффективных таблиц перехода. В ISA присутствует инструкция косвенного переходаJIN rX , которая передаёт управление на адрес, содержащийся в регистровой паре (8-бит). С помощью этой команды я могу переходить на участки кода в зависимости от значения регистра, избежав кучи сравнений. Я активно использовал этот трюк. Вот, к примеру, функция для деления 4-битного числа на 4-битное число. Здесь для каждого значения делителя (1..15) у меня свой блок кода.

__location(00:0x10)
div4x4_1:
  XCH rr12
  LDM 0x0
  XCH rr13
  BBL 0

# INPUT:
#   acc - dividend
#   rr10 - divisor
#   rr11 - 0x0
# OUTPUT:
#   rr12 - quotient
#   rr13 - reminder
__location(00:0x18)
div4x4:
  JIN r5

__location(00:0x20)
div4x4_2:
  RAR
  XCH rr12
  TCC
  XCH rr13
  BBL 0

...

Определившись с необходимостью более изощренной компиляции, я составил такой конвейер превращения исходного кода в образ ROM:

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

Блоки и их связи
На этом этапе у нас есть список блоков и информация о их связях между собой. Теперь можем сформировать более крупные структуры - объединения блоков. Это набор блоков, связанных между собой короткими переходами. Необходимо их рассматривать вместе из-за ограниченний с короткими переходами, описанных выше. Тогда как блоки, ссылающиеся через длинные переходы, могут располагаться в произвольных местах относительно друг друга. К сожалению, мы не можем применить какой-нибудь алгоритм, решающий классическую задачу об упаковке в контейнеры, потому что наши объединения блоков могут пересекать страничную границу ROM, даже несмотря на наличие коротких переходов:

Пример объединения блоков, пересекающий страничную границу
Поэтому я пошёл по более прямолинейному пути - завёл курсор, который указывает на текущую незаполненную позицию в ROM, и который мы двигаем в направлении конца образа. На каждом шаге я пытаюсь вставить какой-либо супер-блок на позицию, заданную курсором, отдавая приоритет наиболее крупным объединениям. Если никакой супер-блок не может быть размещён (скажем, нарушаются внутренние короткие переходы), то двигаем курсор на 1 байт вперёд и опять проверяем оставшиеся блоки.

Блоки, для которых мы указали конкретную позицию, не нуждаются в какой-либо дополнительной обработке компоновщиком. Он их просто вставляет по нужному адресу. Несмотря на то, что основной алгоритм размещения блоков линейно двигает курсор вперёд, я добавил небольшую оптимизацию и стараюсь воткнуть мелкие блоки в свободное пространство между участками кода на фиксированных позициях.

Дополнительное усложнение связано с тем, что мы можем тасовать блоки внутри объединения в поисках варианта, который позволит разместить их по текущему адресу. А одно объединение может содержать весьма большое количество блоков - больше десятка. И количество перестановок у нас будет . Поэтому чтоб не усложнять примитивный алгоритм, я просто-напросто добавил в компоновщик кэш - ключём кэширования является суперблок (какой-нибудь хэш, который зависит от сгенерированного машинного кода), а в качестве значения выступает список смещений внутри страницы ПЗУ и информацией о том, можно ли по данному смещению разместить этот суперблок или нет.

Эмулятор и профилировщик

Не секрет, что отладка и тестирование программ занимают большую часть времени при разработке, а отлаживать код на реальном железе не очень удобно. Я уже писал эмулятор для i4004 и GUI к нему, но потребовался новый функционал:

Хотелось запускать эмулятор не только в браузере, но и иметь программный интерфейс для использования в автоматических тестах.
Адаптация эмулятора под 4040: новые инструкции, расширенный набор регистров, дополнительный банк ПЗУ.
Поддержка снимков ОЗУ. Некоторые функции нашей программы ожидают определенное состояние памяти, поэтому для их тестирования (как автоматического, так и через GUI) нужна возможность указывать начальный снимок RAM.
Профилирование. Без этого мы можем только гадать, что является узким местом нашей программы. При эмуляции мы можем составлять профили на уровне команд, но такая детализация мне была не нужна, и я остановился на уровне функций.
Так как мы добавили компоновщик, то мы потеряли прямое отображение исходного кода в машинный код образа ROM. Раньше мы могли относительно просто транслировать номер строки исходного кода в адрес в ПЗУ. Теперь же требуется поддержка карт кода (source maps) в GUI эмулятора.
Эмулятор уже поддерживал точки останова, но я добавил ещё один полезный инструмент отладки - отображение значения переменных в читаемом виде. Многие переменные являются 16-битными числами. Иногда они хранятся в регистрах, и бывало не так удобно, глядя на список регистров, считывать значение переменной. Или другой пример - переменная может храниться в ОЗУ. Но, имея полную карту памяти перед глазами, не так просто вычленить интересующие ячейки. Особенно, если нужно следить за несколькими переменными.

GUI эмулятора

Первая попытка

Если честно, то к задаче я подошёл расслабленно. 70 часов? Звучит несложно. Поэтому я написал относительно простой код без каких-либо оптимизаций. Тем не менее, стоит упомянуть общие рассуждения, которые пригодились в дальнейшем.

Разрядность переменных

Прежде всего требовалось понять, какой разрядности у нас числа. От этого зависит и используемые алгоритмы примитивных функций, и организация памяти. Мы ограничили искомое количество цифр π () верхней границей в 2043, значит у нас не более 7000, а модуль меньше 14000. Следовательно большинство операций будут над 14-битными числами. Конечно же, мы не будем работать с половиной машинного слова и выровняем нашу разрядность по машинному слову (4 бит). То есть, в основном у нас будут 16-битные переменные.

Напомню о том, как выглядит RAM в системах с Intel 40xx. У нас доступно до 8 банков памяти, каждый банк содержит 16 регистров, состоящих из двух частей: 16 основных 4-битных символа и 4 статусных 4-битных символа. Основное различие между этими двумя категориями заключается в том, как мы получаем к ним доступ.

Пример работы с основными символами:

FIM r0, 0xF0   # loads 0xF0 to r0
SRC r0         # selects register F and main character 0
RDM            # loads the selected main character to the accumulator
XCH rr2        # rr2 = accumulator
FIM r0, 0xF1   # loads 0xF1 to r0
SRC r0         # selects register F and main character 1
RDM            # loads the selected main character to the accumulator
ADD rr2        # accumulator = main character 0 + main character 1

А вот к статусным символам мы можем обращаться через специальные инструкции RDx / WRx, без необходимости предварительно передавать номер символа микросхеме RAM.

FIM r0, 0xF0  # loads 0xF0 to r0
SRC r0        # selects register F and main character 0
RD0           # loads status character 0 to the accumulator
XCH rr2       # rr2 = accumulator
RD1           # loads status character 1 to the accumulator
ADD rr2       # accumulator = status character 0 + status character 1

Как видите, взаимодействие со статусными символами обычно быстрее, чем с основным набором. А так как их всего 4, то нам повезло уложиться в 16 бит и мы сможем разместить переменные в более "быстрой" памяти.

Зная верхние границы для , мы можем также подумать о способе хранения списка простых чисел. Логично получить их только один раз, а не находить их перед каждой новой пачкой цифр π. Примерное количество простых чисел, меньших составляет 2000. Просто записать их все в память не получится, но мы можем хранить их в виде битовой карты: если бит на позиции установлен в 1, это будет значить, что - простое число. Нам доступно 10240 бит в RAM, тогда как нам нужны числа до 14000 (верхняя граница для ), но мы можем схитрить - мы знаем, что чётные числа не являются простыми (за исключением 2), и поэтому можем исключить их из нашей битовой карты. Таким образом, вместо 14000 бит нам требуется только 7000.

Для получения простых чисел используем стандартное решето Эратосфена.

Логарифмы и арифметика чисел с плавающей запятой

Другая тема на подумать - логарифмы. Мне нужно знать значение $log_{2}{10}$ для вычисления . Как вы уже поняли, вычисление логарифма это не такая простая задача для 4040. Конечно же, можно просто использовать предвычисленную константу, но это не укладывается в наши правила. Вместо этого, можем использовать апроксимированное значение: $\log_2{10}\approx\frac{93}{28}$ . Нам должно этого хватить, так как нам не так важно точное значение, лишь бы оно было не слишком маленьким.

Ещё один нужный нам логарифм это $log_{a}{2N}$ . Но, к счастью, нас интересует только целая часть логарифма, поэтому можем использовать простой цикл $v^{max} = \max \{v: a^v < 2N\}$ .

Процессор не поддерживает числа с плавающей точкой из коробки, но по алгоритму должен содержать дробную часть числа. Дабы не реализовывать арифметику чисел с плавающей точкой самостоятельно, я избрал другой путь - будем работать с фиксированной точкой: 1 машинное слово на целую часть и 15 слов на дробную. Наша функция получения цифр π выдаёт по 9 цифр за раз, но еще 6 разрядов необходимо для избежания накопительной погрешности округления.

Модульная арифметика

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

Нам известно мультипликативное свойство модульной арифметики:

$(a \cdot b) \bmod m=(a \bmod m \cdot b \bmod m) \bmod m$

Поэтому мы можем записать такое рекурсивное определение функции модульного умножения $f(x, y) = (a \cdot b) \bmod m$ :

$f(x, y) = \begin{cases} 0 & \text{if $b = 0$}\\ f(2 \cdot a \bmod m, b / 2) & \text{if $b$ is even}\\ (a + f(a, b - 1)) \bmod m & \text{if $b$ is odd}\\ \end{cases}$

Похожая идея используется и для модульного возведения в степень:

$a^b \bmod m = \begin{cases} 1 & \text{if $b = 0$}\\ (a^\frac{b}{2})^2 \bmod m & \text{if $b$ is even}\\ ((a^\frac{b-1}{2})^2 \cdot a) \bmod m & \text{if $b$ is odd}\\ \end{cases}$

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

$A^{\phi(m)} \equiv 1 \pmod m\\A^{\phi(m)-1} \equiv A^{-1} \pmod m$

Так как $m = a^{vmax}$ , то значение функции Эйлера для вычисляется на раз-два:

$\phi(m) = \phi(a^{vmax}) = a^{vmax} - a^{vmax - 1}$

Результаты первого запуска

Насколько мы близки к нашим заветным 70 часам? Увы, очень и очень далеки от них. Я даже не дождался окончания эмуляции, но я использовал полиномиальную аппроксимацию после первых 800 цифр для получения примерной оценки. И увидел ужасающие 14.5 лет. Что ж, видимо намечается весьма долгое и интересное приключение по улучшению производительности...

Оптимизируй это немедленно!

Алгоритм

"Параллельные" вычисления

Прежде чем погружаться в дебри низкоуровневых оптимизаций, я решил поизучать алгоритм более внимательно. Особых чаяний я не питал, ибо Фабрис Беллар очень талантливый программист и математик, и шансы, что он пропустил какой-то очевидный трюк, я оценивал невысоко. Но, тем не менее, он решал немного другую задачу - найти цифру на n-ой позиции в числе π, тогда как я старался найти n цифр числа π. Возможно мы сумеем переиспользовать какие-то вычисления между вызовами основной функции? В памяти мы не можем держать большой кэш - всего порядка 300 записей с 16-битными ключами и значениями. Нужно думать сильнее.

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

$D =\left\{ D + \frac{(f_i \cdot 10^n) \bmod m}{m} \right\}$

В изначальном раскладе мы всё равно будем вычислять для максимального $n_{max}$ , так что это не является дополнительной работой. Да, мы потратим время на ненужные вычисления для $n < n_{max}$ , но зато мы крутим ресурснозатратный внутренний цикл только для $n_{max}$ .

Как обсуждалось выше, хранится как 16-разрядное число, а значит нам нужно 8 байт под него. Если мы хотим проходить внешний цикл по простым числам только один раз и обновлять сразу для всех позиций цифр в π, то нам нужно 1824 байт только под промежуточные значения для . Мы получаем по 9 цифр для заданного , значит для 2048 цифр будем заводить 228 разных значений для - 0, 9, 18, ... А 228 * 8 байт и дают нам требуемые 1824 байт. Это превышает весь доступный объём памяти. Но мы просто разделим всю работу на два прохода - выдадим цифры с позиций [0, 1026], а затем с [1026, 2052]. Тогда потребуется в 2 раза меньше памяти - 912 байт, что, по крайней мере, укладывается в физические ограничения системы.

Код для референсной реализации изменился незначительно:

    for (let i = 0; i < digits.length; i++) {
      d = (f * modularPower(10, startingPosition + i * 9, m)) % m;
      digits[i] = (digits[i] + (d / m)) % 1;
    }

После адаптации версии на ассемблере, я получил впечатляющее ускорение в 85 раз, и программа финиширует за 64 дня вместо нескольких лет. Но нам всё ещё осталось улучшить это время в 20 раз!

Возвращаемся к простым числам

Большая часть памяти теперь отведена под, но у нас же там лежит битовая карта простых чисел. Как же всё это уместить вместе?

Всё просто - вместо обычного решета Эратосфена будем использовать сегментированное. У нас будет начальный сегмент простых чисел меньших $\sqrt{2N}$ (обозначим это значение как длина сегмента). Кроме этого, заведём текущий рабочий сегмент той же длины, и по надобности мы будем его обновлять, используя начальный сегмент.

Произведём оценку требуемой памяти. Длина сегмента у нас 118, а значит мы можем прикинуть, сколько займёт начальный сегмент. Каждое простое число в этом сегменте укладывается в 1 байт, и всего у нас около 30 простых чисел, меньших 118. Следовательно, нам хватит всего 30 байт под начальные данные.

По определению, разница между первым и последним числами в текущем сегменте не превышает длину сегмента. Значит, мы можем организовать хранение теущего сегмента в весьма компактной манере - 16-битная переменная в качестве наименьшего простого числа в сегменте и битовая карта на 59 записей (нам по-прежнему не нужны чётные числа). Таким образом, даже если вместо отдельных битов мы будем использовать целые машинные слова, то нам потребуется всего 32 байта.

Небольшие изменения в формате дробных чисел

Мы не используем BCD формат для (в этом формате, к примеру, число 1234 будет записано в памяти 0x1234 , а не как 0x4D2) потому что он более раздут. А значит десятичная точка для дробных чисел проигрывает шестнадцатеричной. Можно сказать, что в хранилось число, умноженное на $10^{15}$ , а мы меняем этот множитель на $2^{50}$ .

Упрощенная схема для больших простых чисел

Я не слишком верил в то, что только оптимизации кода хватит для достижения нужного результата. Алгоритмические улучшения практически всегда должны идти первым номером, когда не хватает производительности. Поэтому я продолжил исследовать алгоритм. И вскорости подметил, что для простых чисел превышающих $\sqrt{2N}$ , их кратность в всегда будет не больше 1 ( $v^{max} = 1$ ). Следовательно, мы можем упростить кое-какие вычисления.

Значение имеет бинарную природу и чередуется между 0 и 1. И мы даже можем вычислить периоды чередования: начинает иметь значение 1, когда , а возвращается к 0 при . Приведу простую визуализацию такого чередования:

Код стал немного менее читабельным, но зато быстрее.

  do {
    // iterate until next k, that satisfies (2 * k - 1) % a === 0
    for (; (2 * k - 1) % a !== 0; k++) {
      b = (b * k) % a;
      A = (A * (2 * k - 1)) % a;
    }

    // v becomes 1 here
    b = (b * k) % a;
    A = (A * ((2 * k - 1) / a)) % a;
    f = (f + modularInverse(A, a) * b * k) % a;
    k++;

    // iterate until next k, that satisfies k % a === 0
    for (; (k % a) !== 0; k++) {
      b = (b * k) % a;
      A = (A * (2 * k - 1)) % a;
      f = (f + modularInverse(A, a) * b * k) % a;
    }

    // now v becomes 0
    b = (b * (k / a)) % a;
    A = (A * (2 * k - 1)) % a;
    k++;
  } while (k <= N);

В нём мы исключили возведение в степнь и несколько ненужных сравнений, но есть еще простор для улучшений:

Деление тоже оказывается ненужным. Если внимательно присмотреться к циклу, то шагом цикла является текущее простое число . Поэтому мы знаем что внутри цикла значения k / aи (2 * k - 1) / a на каждой итерации цикла увеличиваются на 1 и на 2, соответственно. Целочисленное деление заменяем простым сложением.
Мы можем держать множители для и уже в виде вычетов по модулю . Тогда при модульном умножении нам будет меньше работы.
В начале тела внешнего цикла у нас есть под-цикл для случая, когда . И мы знаем, что перед началом исполнения этого под-цикла $k \bmod a = 1$ , следовательно у нас есть известная последовательность умножений: $b_k \cdot 1 \cdot 2 \cdot 3 ...$ , которая всегда одна и та же на каждой итерации внешнего цикла, поэтому вычисляем этот коэффициент один раз до входа в цикл:
$b_{v=0} = (\frac{a - 1}{2} + 1)! \bmod a$

Делаем референсный код еще менее читабельным.

  let multiplierForZeroV = factorial(((a - 1) / 2) + 1) % a;
  let k = ((a - 1) / 2) + 1;
  let reducedCoefA = 3;

  do {
    // iterate until next k, that satisfies (2 * k - 1) % a === 0
    for (reducedCoefB = 2; reducedCoefB <= (a - 1) / 2; reducedCoefB++) {
      A = (A * reducedCoefA) % a, reducedCoefA += 2;
    }
    b = (b * multiplierForZeroV) % a;

    // v becomes 1 here
    A = (A * multiplierA) % a, multiplierA += 2;
    f = (f + modularInverse(A, a) * b * reducedCoefB) % a;
    k++, reducedCoefB++;

    // iterate until next k, that satisfies k % a === 0
    reducedCoefA = 2;
    for (; reducedCoefB < a; k++, reducedCoefB++) {
      A = (A * reducedCoefA) % a, reducedCoefA += 2;
      b = (b * reducedCoefB) % a;
      f = (f + modularInverse(A, a) * b * reducedCoefB) % a;
    }

    // have a check of the loop boundary now to avoid redundant first loop
    k += (a + 1) / 2;

    // now v becomes 0
    b = (b * multiplierB) % a, multiplierB++;
    A = (A * reducedCoefA) % a, reducedCoefA += 2;

    // try to keep the factor lower than "a"
    if (reducedCoefA > a) {
      // factor is incremented by 2, so it could be bigger than "a" just by 1
      reducedCoefA = 1;
    } else {
      A = (A * reducedCoefA) % a;
      reducedCoefA += 2;
    }
  } while (k <= N);

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

Профилирование

Обычно, если возникают проблемы с производительностью, то необходимо найти узкое горлышко и горячие участки кода. Адаптировал свой профайлер для выдачи трасс, поддерживаемых утилитой flamegraph и получил следующий результат:

Получение обратного числа по модулю

После получения флеймграфа, сразу стало очевидно, какая же операция является узким местом. Судя по всему, "быстрая" реализация, основанная на теореме Эйлера, оказалась не такой уж и быстрой. Набросал вариант с использованием расширенного алгоритма Евклида и получил весьма обнадеживающие цифры. А затем я вспомнил про бинарную версию этого же алгоритма. Как обычно, перед написанием кода на ассемблере, я реализовал алгоритм на языке высокого уровня:

const modularInverse = (a, m) => {
  if (a === 1) {
    return 1;
  }

  let rx = a, ry = m;
  let v = 1, u = 0;

  while (rx % 2 === 0) {
    rx = rx >> 1;
    v = (v % 2 === 0) ? (v >> 1) : (v + m) >> 1;
  }

  while (ry > 1) {
    if (ry % 2 === 0) {
      ry = ry >> 1;
    } else {
      if (ry < rx) {
        [rx, ry] = [ry, rx];
        [v, u] = [u, v];
      }
      ry = (ry - rx) >> 1;
      u = (u - v) % m;
    }

    u = (u % 2 === 0) ? (u >> 1) : ((u + m) >> 1);
  }

  return u < 0 ? (u + m) : u;
}

Перевод данного фрагмента на ассемблер прошёл весьма просто. Единственное, на что можно обратить внимание это то, как производится знаковая операция сдвига (только функция получения обратного оперирует отрицательными числами). На удивление, получилось весьма лаконично:

  LD rr3       # load highest word of number
  RAL          # put sign bit to carry to be able to shift negative numbers
  LD rr3       # again load the highest word of the number
  RAR          # shift it right, but because carry stores sign a bit 
               # it would be propagated back to the highest bit of word
  XCH rr3      # save updated word

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

Модульное возведение в степень

Возводить в степень нам нужно в двух случаях: 1) для вычисления и 2) вычисление $a^{v^{max} - v + u}$ .

Мы знаем что пробегает от 0 до заданного с шагом цикла в 9, следовательно мы можем не вычислять каждый раз, а завести переменную, и на каждой итерации умножать её на . Так мы заменяем возведение в степень тремя умножениями (наши функции работают с 16-битными числами, поэтому мы не можем умножить число на за один вызов функции умножения):

  let powered10 = modularPower(10, startingPosition, m);
    for (let i = 0; i < digits.length; i++) {
      const d = (powered10 * f) % m;
      digits[i] = (digits[i] + (d / m)) % 1;
      powered10 = (powered10 * 1000 * 1000 * 1000) % m;
    }

Второй случай актуален для $v^{max} > 1$ .

Если посмотрим на ограничение сверху: $v^{max} < \lfloor log_{a_0}{2N} \rfloor$ , то нетрудно заметить, что степень, в которую мы возводим, не превосходит 7. Кроме того $a^{v^{max} - v + u} < a^{v^{max}}$ , а значит нам не требуется модульное возведение в степень, достаточно будет обычного.

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

Таким образом мы совсем избавляемся от модульного возведения в степень.

Деление

После обновления алгоритма и отделения ветви кода со случаем $v^{max} = 1$ , у нас осталось не так уж и много операций деления. Тем не менее, если мы можем их оптимизировать, то почему бы и нет?

Деление 4-битного числа на 4-битное

Даже если у нас нет задачи делить 4-битные числа, то всё равно эта операция является базой для длинного деления. Простое вычитание в цикле оказалось не таким уж и медленным и занимало 17 циклов в среднем (для всех комбинаций 4-битного делителя и делимого):

div4x4:
  FIM r6, 0x00
div4x4_loop:
  SUB rr11
  JCN nc, div4x4_return
  CLC
  ISZ rr12, div4x4_loop
div4x4_return:
  ADD rr11
  XCH rr13
  CLC
  BBL 0

В главе про компилятор и компоновщик, я уже упоминал подход с 15 специализированными функциями деления, каждая из которых реализует деление на конкретное число от 1 до 15. Размер такой специализированной функции не может превышать 16 байт, потому что они располагаются по фиксированным адресам 0x10, 0x20, 0x30, ... Для некоторых делителей было невозможно уложиться в данные рамки, но так как у нас был некий общий код, используемый в нескольких блоках и код функций содержал условные переходы, то я смог упаковать код для всех случаев в одну страницу ROM не тратя циклы на сшивание разных участков между собой.

Среднее значение циклов на один вызов операции деления уменьшилось до 11. Казалось бы, 6 циклов не так много, но в относительном выражении это 35% ускорение для этой небольшой функции.

Умножение 4-битных чисел

Другая базовая операция, которая используется повсеместно, в том числе и для длинного деления - это умножение 4-битных чисел в 8-битное произведение. Если наивная версия 4-битного деления уже была относительно быстра до начала оптимизации, то мы не можем сказать того же самого про первую версию операции умножения. Серия сложений одного из множителей в цикле с количеством итераций, заданным другим множителем, показала себя не очень эффективно - около 70 циклов в среднем.

Для выправления ситуации воспользуемся тем же подходом, что и в прошлом разделе. Но придётся еще всё же заиметь таблицу умножения в ОЗУ. Мы можем сократить использование памяти под эту таблицу, написав быстрые функции умножения на 0, 1, 2 и 8. Умножение числа на 2 и 8 можно сделать через сдвиговые операции. Для всех остальных множителей мы обращаемся к памяти.

__location(01:0x00)
mul4x4_0:
  FIM r6, 0x00                                         # low, high = 0
  BBL $BANK_WITH_VARIABLES

# INPUT:
#   acc - first factor (a)
#   rr10 - second factor (b)
#   rr11 - 0x0
# OUTPUT:
#   rr12 - low word of product
#   rr13 - high word of product
__location(01:0x08)
mul4x4:
  JIN r5

__location(01:0x10)
mul4x4_1:
  XCH rr12                                              # low  = a
  LDM 0x0
  XCH rr13                                              # high = 0
  BBL $BANK_WITH_VARIABLES

__location(01:0x20)
mul4x4_2:
  RAL
  XCH rr12
  TCC
  XCH rr13
  BBL $BANK_WITH_VARIABLES

__location(01:0x30)
mul4x4_3:
  XCH rr12
  LDM $BANK_WITH_MUL_TABLE_FOR_FACTOR_3
  DCL
  SRC r6
  RD2
  XCH rr12                                              # low
  RD3
  XCH rr13                                              # high
  BBL $BANK_WITH_VARIABLES

Деление 16-битных чисел

Существует известный алгоритм длинного деления от Кнута [^12], но он весьма обобщенный - вместо этого мы можем написать несколько разных вариаций деления в зависимости от разрядности делителя и делимого. К примеру, вариант функции для деления 12-битного делимого на 8-битный делитель, или 8-битное делимое и 4-битный делитель, и так далее. Всего у нас 10 таких комбинаций (делитель всегда имеет меньшую разрядность чем делимое) и только для 3х из них нам требуется написать реализацию алгоритма D: 16x12, 16x8 и 12x8. Для всех других случаев мы можем разработать более быстрый код.

Модульное умножение

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

Алгоритм умножения

У меня было искушение использовать алгоритм Монтгомери для модульного умножения, причём мы можем держать числа в форме Монтгомери довольно долго - вернуться к нормальной форме нужно будет только перед вычислением . Однако прикинул, что этот алгоритм включает в себя 3 обычных умножения 16-битных чисел (4 разряда, ), и даже со всеми оптимизациями нам потребуется умножений с 4-битными операндами [^13]. Легко посчитать - 30 умножений и каждое требует как минимум 12 циклов, а значит мы потратим 360 циклов только на умножение, а еще не забываем про 150 операций сложения (формула есть в статье выше). Слишком долго.

Ок, а что насчёт обычного длинного умножения с последующим приведением произведения по модулю? Основная операция ведёт к 16 вызовам 4-битного умножения [^14]. А приведение по модулю еще более затратное действие. Тоже не вариант. Алгоритм Карацубы не сильно эффективен для небольших чисел: у нас всё же всего 4 разряда.

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

const modularMultiply = (a, b, m) => {
  let multipliedFactor = a;
  let result = 0;

  for (let shiftedFactor = b; shiftedFactor > 0; shiftedFactor = shiftedFactor >> 1) {
    if (shiftedFactor & 0x1 === 0x1) {
      result = (result + multipliedFactor) % m;
    }

    multipliedFactor = (multipliedFactor + multipliedFactor) % m;
  }

  return result;
};

Ленивое приведение по модулю

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

Однако не всё так просто - если совсем исключить этот этап, то у нас произойдёт целочисленное переполнения (мы же работаем с 16-битными числами). Изначально, приведение по модулю кодировалось всего лишь одной операцией вычитания: x = x > m ? x - m : x ибо оба слагаемых всегда были меньше чем . С новым подходом это уже не так и нам нужно подумать как мы будем эффективно приводить числа по модулю.

Переполнение потенциально может произойти только при операциях с числами больше 0x1000. Поэтому нам не обязательно обеспечивать принадлежность чисел к системе наименьших вычетов, а достаточно просто уменьшать их разрядность корректно. Для этого мы заводим небольшую таблицу поиска (LUT) с числами кратными , подобранными так, чтобы вычитанием числа из таблички уменьшить заданную переменную до 3-х разрядов. Само собой, класс вычета этой переменной не меняется (ведь отнимаем число, кратное модулю), так что всё ок. Математически LUT описывается тривиально, ключем в таблице является старший разряд числа:

$x = \overline{x_3 x_2 x_1 x_0}, x_3 \rightarrow \lfloor \frac{x_3 \cdot 16^{3}}{m} \rfloor \cdot m$

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

$\lfloor \frac{F00_{16}}{m} \rfloor \cdot m$

Еще более ленивое приведение по модулю

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

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

Сложение вместо вычитания

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

# compute [rr0, rr1] = [rr0, rr1] - [rr2, rr3]
LD rr0
SUB rr2    # acc = rr0 - rr2, borrow flag is 0, if there was borrow
CMC        # inverse borrow flag
XCH rr0
LD rr1
SUB rr3    # acc = rr1 - rr3 - borrow, borrow flag is 0, if there was borrow
CMC        # inverse borrow flag
XCH rr3

В свою очередь, сложение не имеет такой особенности:

# compute [rr0, rr1] = [rr0, rr1] + [rr2, rr3]
LD rr0
ADD rr2    # acc = rr0 + rr2, carry flag is 1, if there was carry
XCH rr0
LD rr1
ADD rr3    # acc = rr1 + rr3 + carry, carry flag is 1, if there was carry
XCH rr3

При работе с 16-битными числами мы можем заменить вычитание сложением с (дополнительный код для ) и получить правильный результат, но с меньшим числом инструкций.

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

Нам нужно больше таблиц перехода!

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

  # assume that rr6 has b[0]
  LDM 0x1
  AN6
  JCN z, mulMod_skipMultiplierBit0       # if ((b >> 0) & 0x1)
  JMS mulMod_updateResult                #   res = res + multipliedFactor
mulMod_skipMultiplierBit0:
  JMS mulMod_updateMultipliedFactor      # multipliedFactor = multipliedFactor * 2

  LDM 0x2
  AN6
  JCN z, mulMod_skipMultiplierBit1       # if ((b >> 1) & 0x1)
  JMS mulMod_updateResult                #   res = res + multipliedFactor
mulMod_skipMultiplierBit1:
  JMS mulMod_updateMultipliedFactor      # multipliedFactor = multipliedFactor * 2

  # ... 14 more checks for particular bits

Но если подумать получше, то весь алгоритм заключается в череде сложений - на каждой итерации обновляем multipliedFactor и, в зависимости от текущего бита множителя, увеличиваем result. При этом для каждого 4-битного разряда множителя мы знаем конкретную последовательность операций сложения. Например, если текущий разряд множителя равен 0x5 (0b0101), то нам нужно обновить переменную с промежуточным результатом, дважды удвоить multipliedFactor, снова обновить result и опять учетверить multipliedFactor. В коде это выглядит как-то так:

__location(0x2:0x50)
mulMod_processNibble_5:                   # 0b0101
  JMS mulMod_updateResult
  JMS mulMod_updateMultipliedFactor
  JMS mulMod_updateMultipliedFactor
  JMS mulMod_updateResult
  JMS mulMod_updateMultipliedFactor
  JUN mulMod_updateMultipliedFactor

Опираясь на значение регистра, который содержит значение текущего разряда множителя, мы можем совершить прыжок на нужный фрагмент:

  # rr12 contains current digit for multiplier, rr13 = 0
  JIN r6

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

__location(0x3:0x50)
mulMod_processLastNibble_5:                   # 0b0101
  JMS mulMod_updateResult
  JMS mulMod_updateMultipliedFactor
  JMS mulMod_updateMultipliedFactor
  JUN mulMod_updateResultLast

Выполняем несколько последовательных обновлений multipliedFactor за одну операцию

Иногда (если не возникает переполнение) мы можем вычислять 16 * multipliedFactor и 8 * multipliedFactor более оптимально, с использованием сдвиговых операций.

Если есть риск переполнения, то просто организуем запасной вариант с несколькими вызовами mulMod_updateMultipliedFactor .

Включаем в операцию сложения промежуточного результата еще и удваивание multipliedFactor

Почти всегда за обновлением значения промежуточного результата следует и обновление multipliedFactor :

  JMS mulMod_updateResult
  JMS mulMod_updateMultipliedFactor

И здесь мы можем внедрить микро-оптимизацию, создав функцию mulMod_updateResultAndFactor, которая объединит код из mulMod_updateResultи mulMod_updateMultipliedFactor. Каким же образом это нам поможет? Арифметика простая - раньше у нас было 2 инструкции вызова функции и 2 инструкции возврата, а теперь у нас будет 1 инструкция вызова и 1 инструкция возврата. Каждый раз экономим по 3 машинных цикла. Возможно это звучит не очень впечатляюще, но напомню, что модульное умножение ответственно за 60% времени исполнения программы, а так как сама по себе функция выполняется не так долго, то даже такие единичные сохраненные такты позволяют выиграть минуты и десятки минут в общем табеле.

Более шустрая проверка на потенциальное переполнение переменных

У нас есть контроль за возможным переполнением при удваивании multipliedFactor. Заключается он в проверке на то, что старший разряд этой переменной меньше 0x8. Если это так, то при умножении на 2 мы точно не превысим разрядность:

  XCH rr2                                           # multipliedFactor = multipliedFactor + multipliedFactor
  AN7                                               # check if previous value of multipliedFactor[3] < 4, then new multipliedFactor[3] < 8
  JCN z, mulMod_updateResultAndFactor_return        # if (multipliedFactor[3] < 0x8)

Проверка занимает всего 3 машинных цикла, но, как я писал выше, каждый такт на счету. Поэтому я использовал очередной трюк - заменил AN7/JCN одной инструкцией JIN. Да, будет использована еще одна страница ПЗУ под таблицу перехода, но оно того стоило, тем более небольшой запас ROM у нас был. Если в rr2 у нас содержится старший разряд числа, то используя JIN r1 мы обеспечиваем переход на одну их двух ветвей - когда rr0 < 0x8 и когда rr0 >= 0x8. Да, у нас будут две серии по 8 одинаковых блоков кода, но это уже финальные штрихи, и я смогу пережить такое попрание принципов написания чистого кода.

XCH rr2
JIN r1

__location(0x4:0x00)
mulMod_updateMultipliedFactor_factor0:
  BBL 0

__location(0x4:0x10)
mulMod_updateMultipliedFactor_factor1:
  BBL 0

...

__location(0x4:0x80)
mulMod_updateMultipliedFactor_factor8:
  SRC r1
  RD0
  ADD rr0
  XCH rr0
  RD1
  ADD rr1
  XCH rr1
  RD2
  ADD rr4
  XCH rr4
  RD3
  ADD rr2
  XCH rr2
  CLC
  BBL 0

...

Оптимизации, которые не вошли в конечный код

Было еще несколько рабочих идей, но, по той или иной причине, они не представлены в финальной версии кода:

Ветвь кода для умножения по 12-битному модулю. Это было эффективно, когда произведение всегда было гарантированно меньше модуля. Тогда для небольших модулей мы могли избавиться от кода работы со старшими разрядами.
Быстрое сравнение переменной со значением модуля. Опять же, до внедрения ленивого приведения мы проводили сравнение довольно часто, дабы сразу когда число становится больше мы могли его уменьшить и вернуть в систему наименьших вычетов. Сама операция сравнения относительно затратная - чтоб сравнить два числа, необходимо вычесть одно из другого и проверить флаг заимствования. Ради ускорения я добавлял еще одну таблицу поиска, ключем являлся старший разряд числа, а значением - индикатор того, точно ли переменная больше модуля.
```
SRC rX      # first register in register pair contains the address of LUT, second register - multipliedFactor[3]
RDM         # read isLessThanModulus = LUT[multipliedFactor[3]]
JCN z, ret  # return if (isLessThanModulus)
```
Другой, испробованный мною подход по уменьшению ненужных операций вычитания, заключался в том, чтобы держать промежуточный результат меньше 0. После операции сложения, если число становилось положительным (а проверить это нетрудно), то мы вычитали . До этой оптимизации, псевдокод ассемблерной реализации выглядел примерно так:
```
result = result + multipliedFactor
[tmp, borrow] = result - m

// check if result - m > 0
if (!borrow) {
  result = tmp
}
```
После оптимизации:
```
[result, carry] = result + multipliedFactor

// check if result > 0 now
if (carry) {
  result = result - m
}
```
Заметно отличие в количестве операций. Надо только не забыть скорректировать возравщаемое значение, добавив перед выходом из функции.
Объем работы функции модульного умножения напрямую зависит от позиции старшего единичного бита числа, выступающего в качестве второго операнда функции. Следовательно мы должны стараться сделать так, чтобы второй операнд был меньше чем первый. Мы можем их сравнить и поменять местами, если данное условие не выполняется. Результаты были неоднозначными - сама операция сравнения занимает какое-то время и нужно было очень тонко подбирать в каких случаях мы хотим использовать данный функционал, а в каких сравнение будет лишь тратой времени. В итоге, оптимизация стала напоминать подгонку кода под конкретные входные данные и я убрал её.
Можно было завести еще одну таблицу поиска для более точного приведения переменных по модулю. Даже для "ленивой" версии. Но в ОЗУ оставалось не так много свободного места, да и, навскидку, заметный выигрыш был только для небольших модулей, в остальных случаях код мог работать даже медленее.

Структура памяти

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

Легенда:

Зелёная секция содержит значения , всего 114 регистров из 128 доступных.
Фиолетовая - таблицы поиска для функции модульного умножения.
Красная - таблицы умножения. Каждая ячейка таблицы занимает 8 бит, а в статусных символах регистра памяти хранится 16 бит (4 слова). Так что каждый банк памяти хранит таблицы для двух множителей. Всего у нас 12 множителей (исключаем 0, 1, 2 и 8), а значит у нас заняты статусные символы в 6 банках.
Оранжевая (очень похожа на жёлтую) - таблица со степенями . Может содержать до восьми 16-битных чисел, то есть 32 байта или 2 регистра памяти.
Жёлтая - карта простоты для текущего сегмента простых чисел.
Синяя - начальный сегмент простых чисел.
Белая - свободна для использования под произвольные переменные.

Весьма плотная компоновка, но, к счастью, чуть больше свободы чем в моём прошлом проекте.

Результаты запуска на эмуляторе

После всего этого безобразия с различными хаками и оптимизациями, программа наконец-то уложилась в отведённые рамки - я получил правильную последовательность цифр π за 69 часов, 29 минут и 02 секунды. Отлично, осталось проверить работу на реальном микропроцессоре.

Аппаратная часть

Электрическая принципиальная схема мало чем отличается от схемы платы для i4004. Основные отличия:

Самостоятельно распаял stm32 на плате со всей обвязкой, вместо использования сторонней отладочной платы, совместимой с DIP сокетом.
Питание Intel 4040 процессора берется с USB, тогда как в прошлой версии требовалось внешнее питание на 20 вольт.
Управление системой происходит через USB, а не через UART интерфейс.

В остальном, компонентая база не поменялась - stm32f4 в качестве эмулятора памяти (ROM и RAM), LM339 и CD40109 для конвертации уровней сигналов с 15 вольт для i4040 на 3.3 вольта для stm32, и обратно.

Прошивка

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

Интерфейс памяти у 4040 чуть более продвинутый, чем у 4004: имеется дополнительный пин для выбора банка ПЗУ. Для нас это играет роль, потому что добавляется больше вычислений для определения запрашиваемого адреса. Кроме того, изначально я запускал плату с несколько заниженными значениями тактирующего сигнала, идующего к Intel 4004. 625 килогерц укладывалось в спецификацию процессора, но для того чтобы уложиться в 70 часов, нам нужно брать максимум от процессора, то есть тактировать его частотой в 740 килогерц.

И тогда-то тайминги поплыли - на каждый такт i4040 процессора нужно анализировать входные пины, проводить какие-то вычисления и обновлять сигналы на выходных пинах. При 740 килогерцах, длина такта на i4040 всего 1350нс. При частоте ядра stm32 в 168 мегагерц, у нас всего 200 тактов stm32 на всё про всё.

С оригинальным кодом я перестал укладываться в это узкое окно. Изначально, я использовал прерывания от тактирующего таймера, чтобы знать когда начинается такт на i4040, но с таким подходом были неизбежны задержки между фронтом тактового сигнала и стартом моего кода в обработчике прерывания. Я заменил прерывания на примитивный опрос пина:

  while ((OUT_4040_PHI1_GPIO_Port->IDR & OUT_4040_PHI1_Pin) == OUT_4040_PHI1_Pin);

  OUT_TEST1_GPIO_Port->ODR ^= OUT_TEST1_Pin;
  handleCyclePhi1Falling();
  OUT_TEST1_GPIO_Port->ODR ^= OUT_TEST1_Pin;

Пин TEST1 я использовал для отслеживания сколько времени stm32 тратит на работу между тактами.

Так как тайминги очень жесткие, то я отрубаю все прерывания. После того как stm32 получает команду на старт с дампом ПЗУ, всё взаимодействие с ПК замораживается. Если Intel 4040 использует инструкцию ввода-вывода для взаимодействия с каким-то внешним устройством через порты эмулированной ROM/RAM, то я просто отслеживаю такую активность и записываю информацию в буфер. Как только Intel 4040 переходит в состояние HALT (используется инструкция HLT, которой, кстати, не было в i4004), прошивка на stm32 определяет это и отправляет накопленную статистику и вышеупомянутый буфер на ПК.

Финальные результаты

Само собой, был ряд проблем как в аппаратной, так и в программной частях платы, но в итоге я смог добиться запуска кода на Intel 4040. И тут меня поджидало душераздирающее открытие - я допустил очепятку в коде эмулятора при вычислении времени выполнения программы в секундах. Машинные циклы считались правильно, но для перевода их в секунды я делил на 95000, а должен был на 92500. А это значит, что на тот момент я только думал, что всё почти закончено. Поэтому, скрипя зубами. я вернулся за оптимизацию кода.

Я упомянул этот момент, чтобы показать как важно проводить реальные тесты на настоящем железе, и не довольствоваться оценками, полученными теоретически или даже при симуляции/эмуляции.

В итоге, я добил необходимые несколько процентов скорости выполнения программы, и вернулся к физическому i4040. Не хотелось рисковать и ждать 70 часов без гарантии что код отработает корректно, поэтому сначала я запустил вычисление 255 цифр π, вместо 2048+.

[2023-10-31T12:08:05.527Z] START command has been received, RAM dump size is 8192 bytes
[2023-10-31T13:13:08.988Z] Program has been finished
[2023-10-31T13:13:08.994Z] Instruction cycles processed = 00000000158651E4
[2023-10-31T13:13:08.995Z] Output from i4040 (263 nibbles):
[2023-10-31T13:13:08.995Z]   3 1 4 1 5 9 2 6 5 3 5 8 9 7 9 3 2 3 8 4 6 2 6 4 3 3 8 3 2 7 9 5 0 2 8 8 4 1 9 7 1 6 9 3 9 9 3 7 5 1 0 5 8 2 0 9 7 4 9 4 4 5 9 2 3 0 7 8 1 6 4 0 6 2 8 6 2 0 8 9 9 8 6 2 8 0 3 4 8 2 5 3 4 2 1 1 7 0 6 7
[2023-10-31T13:13:08.995Z]   9 8 2 1 4 8 0 8 6 5 1 3 2 8 2 3 0 6 6 4 7 0 9 3 8 4 4 6 0 9 5 5 0 5 8 2 2 3 1 7 2 5 3 5 9 4 0 8 1 2 8 4 8 1 1 1 7 4 5 0 2 8 4 1 0 2 7 0 1 9 3 8 5 2 1 1 0 5 5 5 9 6 4 4 6 2 2 9 4 8 9 5 4 9 3 0 3 8 1 9
[2023-10-31T13:13:08.995Z]   6 4 4 2 8 8 1 0 9 7 5 6 6 5 9 3 3 4 4 6 1 2 8 4 7 5 6 4 8 2 3 3 7 8 6 7 8 3 1 6 5 2 7 1 2 0 1 9 0 9 1 4 5 6 4 8 5 6 6 9 2 3 F

Как можете увидеть, это заняло всего час времени, в сравнении с 3.5 часами краникового алгоритма на i4004. Само собой, это не является каким-то объективным сравнением: разные процессоры, разные алгоритмы, даже разные тактовые частоты. Более того, мне кажется, что если я перенесу ряд оптимизаций в краниковый алгоритм, то не факт что он будет медленнее. Тем не менее, было любопытно запустить одну и ту же задачу в разных проектах.

Ободрённый тем, что результат совпал с эмулятором, я отправил программу для вычисления 2048+ цифр и стал ждать примерно 70 часов, когда она завершится. Я даже прилепил два небольших радиатора (снятых с dip-8 микросхем) на сам i4040 процессор ради несколько лучшего охлаждения, хотя, если честно, он грелся не так сильно.

Всё прошло успешно и i4040 смог перебить результат ENIAC'a со временем в 69 часов, 28 минут и 31 секунду.

PS: очень внимательный читатель может заметить, что время исполнения программы на железе отличается от времени на эмуляторе, даже если по тактам совпадение очень близкое. Ответ заключается в том, что выходное тактирование с использованием таймеров stm32 не идеально на 100%, и выходная частота на самом деле около 740.1 килогерца, а эмулятор предполагает то, что процессор запущен на частоте ровно в 740 килогерц. Можно сказать, что я даже немного разогнал Intel 4040 (возможно впервые в истории)!

Ссылки

[^1]: Brian J. Shelburne, "The ENIAC's 1949 Determination of π"

[^2]: George W. Reitwiesner, An ENIAC Determination of π and e to more than 2000 Decimal Places

[^3]: Intel Intellec 4/Mod 40 Microcomputer Development System Reference Manual

[^4]: Stanley Rabinowitz and Stan Wagon, "A Spigot Algorithm for the Digits of π"

[^5]: Jeremy Gibbons, "Unbounded Spigot Algorithms for the Digits of Pi "

[^6]: David Bailey, Peter Borwein and Simon Plouffe, "On the Rapid Computation of Various Polylogarithmic Constants"

[^7]: Simon Plouffe, "On the computation of the ݊n'th decimal digit of various transcendental numbers"

[^8]: Fabrice Bellard, "Computation of the n'th digit of pi in any base in O(n^2)"

[^9]: Xavier Gourdon, "Computation of the n-th decimal digit of π with low memory"

[^10]: D. H. Lehmer, "Interesting Series Involving the Central Binomial Coefficient"

[^11]: Intel, "MCS-40 User's Manual For Logic Designers"

[^12]: Donald E. Knuth, "Art of Computer Programming, Volume 2: Seminumerical Algorithms", 4.3. Multiple Precision Arithmetic

[^13]: Hwajeong Seo, Zhe Liu, Yasuyuki Nogami, Jongseok Choi, Howon Kim, "Hybrid Montgomery Reduction"

[^14]: Michael Hutter, Peter Schwabe, "Multiprecision multiplication on AVR revisited"

Комментарии (34)

JordanCpp
14.11.2023 09:32
#26157870
+38
Аплодисменты. Спасибо, неважно насколько это актуально. Сам рассказ очень интересен.

Dark_Purple
14.11.2023 09:32
#26158586
+3
Годнота!

vvbob
14.11.2023 09:32
#26159108
+8
Снимаю виртуальную шляпу!

Интересно, а если заменить эмулированную память на что-то времен создания процессора, насколько и куда поменяются результаты. Что-бы прямо по хардкору все было на той элементной базе.
1. mark_ablov Автор
  14.11.2023 09:32
  #26159588
  +9
  Результат не поменяется - операции с памятью синхронные. Вот 8008 уже имеет READY пин для ожидания когда медленная память сможет выполнить операцию в/в. Но 4004/4040 читают шину данных в следующий такт после передачи адреса. У меня, кстати, есть плата для 4040 с достаточным количеством памяти, но мне было интересно ещё и поработать над железной частью.
  1. vvbob
    14.11.2023 09:32
    #26159704
    Понятно, спасибо!
    
    Всегда когда читал историю процессоров, 4004 казался какой-то малополезной штукой, практически концептом нормального процессора, а оказывается что это не совсем так.
    
    hitomi2500
    14.11.2023 09:32
    #26164570
    +1
    В мемуарах Федерико Фаджина (который проектировал 4004) упоминается, что руководство Intel в те года в целом считало микропроцессоры какой-то малополезной штукой, которая в лучшем случае поможет поднять продажи памяти. И только после того как Фаджин ушёл от них и создал Zilog, зарабатывающий только на процессорах, они спохватились. Возможно Федерико лукавит, ведь борьбу с Intel он проиграл, да и в целом процессоры (хоть и не микро) к тому моменту были уже массовым явлением с понятной областью применения.
    
    vvbob
    14.11.2023 09:32
    #26164676
    Ну, он действительно выглядел очень бледно, даже на фоне тогда уже существующего железа, поэтому вполне возможно что и руководство Intel не заметило тогда перспектив и наверстывать начали когда поняли, что тема очень перспективная. Примеров похожих много было, те-же IBM, неслабо облажались со своей РС, из-за того что на старте считали его чем-то проходным и не особо важным.
  1. sim31r
    14.11.2023 09:32
    #26159988
    Похоже тогда с волновым сопротивлением дорожек не заморачивались, как сейчас это делают.
    
    Dr_Faksov
    14.11.2023 09:32
    #26160172
    +3
    Частоты не те. 740 килогерц, не мегагерц.
    
    sim31r
    14.11.2023 09:32
    #26160688
    Кроме частот это электромагнитные помехи. Возможно требования были ниже, сейчас обязательная сертификация по излучению. Схема работает стабильно, но по шуму может на порядок нормы превышать, пока плату не переделают полностью, иногда переходят с 2 слоёв на 4 или 8 только для этого. Чтобы помех не было, можно с другой стороны зайти сгладить фронт импульса на выходе микросхемы du/dt, но вряд ли с этим заморачивались. Как я понимаю микросхемы старались максимально быстро, насколько выходной каскад позволяет переключить состояние (косвенный признак конденсатор над каждой микросхемой, чтобы питание не просаживалось при переключении). А слабые и медленные транзисторы на выходе косвенно помогали снизить помехи и наводки от переключений. Тут можно исследовать вопрос по старым микросхемам, просто осциллографом посмотреть сигналы на выходе. В современных микроконтроллерах например есть быстрый и медленный режим для переключения ног, 2 МГц и 50 МГц (частота условно как я понимаю, тут речь о сглаживании фронта импульса, будет он пологий в 100 нс или максимально крутой в 5 нс).
    
    mark_ablov Автор
    14.11.2023 09:32
    #26160884
    Вот небольшой пример с оциллографа - SYNC / CM_ROM сигналы генерировались i4004, ph1/ph2 - современным stm32 (ну как современным, 10+ лет уже ему).
    
    sim31r
    14.11.2023 09:32
    #26161068
    Интересно, даже у i4004 переключения жесткие. Только CM-ROM затянут по спаду, возможно на выходе резистор с небольшим сопротивлением к +5В есть, что замедляет переключение.
    
    Если там 5-10 нс переключения плата может шуметь больше, чем современные платы с частотами процессоров в гигагерцы из-за особенностей разводки печатной платы.
    
    mpa4b
    14.11.2023 09:32
    #26161598
    +1
    PMOS логика имеет подтяжки к отрицательному потенциалу и открытые транзисторы тянут к положительному. Поэтому фронт вниз затянутый, вверх -- побыстрее.
    
    Ну и конечно если хочется померять длительность фронтов, картинку бы растянуть по времени раз в 20-30 надо.
    
    mpa4b
    14.11.2023 09:32
    #26161586
    +1
    Шутка что ли такая, супер-медленная PMOS логика, работающая от -15 вольт, с трудом вытягивающая мегагерц, какие нафиг du/dt :)
    
    sim31r
    14.11.2023 09:32
    #26162602
    +1
    Это логические элементы медленные. А переключение состояния на выходе быстрое, SYNC за 10 нс примерно состояние меняет. Этого достаточно чтобы на не оптимизированной плате шуметь эпически в эфир и на соседние дорожки иглами напряжения (спектр до гигагерца). Другое дело что медленная логика на них не отреагирует (на краткие помехи наносекундные), как и норм по помехам тоже не было особых.

VBDUnit
14.11.2023 09:32
#26159494
+17
Хабр торт.

ef_end_y
14.11.2023 09:32
#26160006
+2
Переходы только в тот же банк памяти это жесть, конечно. И сколько времени прошло пока додумались указывать не адрес, а смещение. Даже в 8080 все ещё был адрес. А вообще я не осилил весь текст), но вы слегка затронули мою ностальгию, повспоминал как боролся за уменьшение и тактов и байтов... Это было похоже на решение головоломок

Refridgerator
14.11.2023 09:32
#26160274
+5
Если кто-то всё ещё не понимает, зачем программисту математика - то вот зачем.
1. sim31r
  14.11.2023 09:32
  #26160698
  +1
  Или как вариант, пример того зачем математику уметь программировать )
1. vvbob
  14.11.2023 09:32
  #26160832
  +1
  Что-бы считать Pi на 4004? :)))
  
  Не, респект автору, все круто, но откровенно говоря рядовому программисту очень редко реально оказывается нужен весь этот "матан", даже когда нужна какая-то очень жесткая оптимизация, то обычно оптимизируют алгоритм, без погружения в особо сложную и крутую математику. Но знать все это не помешает, конечно, особенно если приходится что-то писать в таком духе.
  1. mark_ablov Автор
    14.11.2023 09:32
    #26160930
    +9
    Математика нужна, чтобы просто знать о существовании ряда алгоритмов / формул. Не обязательно даже понимать как они работают. Скажем, попалась тебе задача, и ты уже можешь прикинуть как её решать наиболее оптимально. Для многих прикладных задач существуют лучшие практики их решений, но нужно уметь увидеть то, что конкретная проблема в абстрактной форме ложится на известный алгоритм.
    
    PS: само собой, всё это не нужно чтобы перегонять JSON'ы из одного формата в другой, паралельно реализуя CRUD. Но, к примеру, при разработке распределенных систем "матан" уже встречается не так уж и редко.
    
    vvbob
    14.11.2023 09:32
    #26161280
    -2
    Да я и не топлю против математики, она полезна во многих смыслах, но, к сожалению, не так уж часто в обычной практике встречаются задачи, где она нужна, какая-то сложная.
    
    cadmi
    14.11.2023 09:32
    #26161462
    +2
    Просто это именно то, чем различаются собственно "программист" и "кодерок". Но их почему-то постоянно путают, вот и создается иллюзия, что "математика не нужна".
    
    Никакого хейта или обесценивания, просто дежурное напоминание, что это действительно разные профессии.
    
    vvbob
    14.11.2023 09:32
    #26164650
    -1
    Первый раз слышу про профессию "кодерок". Как по мне именно что хейт и обесценивание, вроде если ты в своей работе не решаешь интегралы каждый день и не пользуешься "китайской теоремой", то кодерок ты поганый, а не программист.
    
    sim31r
    14.11.2023 09:32
    #26164752
    Иногда требуется минимальное знание математики, а его у нет у претендентов на вакансию. Для разработки игры искали сотрудника и претенденты не могли найти расстояние между двумя точками на плоскости. Типа оценить в зоне досягаемости ли цель. То есть просто теорему Пифагора применить.

rubinstein
14.11.2023 09:32
#26160338
+8
Читаешь, задержав дыхание! Побольше бы таких статей, а не "как я вошёл в айти".

wmlab
14.11.2023 09:32
#26160374
+8
Отличная работа! Дает представление, чем занимались программисты "Вояджеров".
1. sim31r
  14.11.2023 09:32
  #26164756
  Они еще код меняли чтобы битую ячейку памяти обойти, сожженую высокоэнергетической частицей из дальнего космоса.

Vaitek
14.11.2023 09:32
#26160668
+4
Ничего не понял дальше вступления. Очень круто. Пишите еще!

EvilBeaver
14.11.2023 09:32
#26161318
+2
Безумству храбрых поем мы песню! Ухх... У меня аж привста.. Восторг, чистый восторг!

555www555
14.11.2023 09:32
#26161370
+2
Отличная работа и ее описание. Спасибо, снимаю шляпу. А у нас есть потребность по еще бОльшей оптимизации тока потребления наших самых маленьких в мире диктофонов Edic-mini www.telesys.ru Хотим сделать их еще меньше. Может есть желание поучаствовать?

mordoorg
14.11.2023 09:32
#26161472
+1
Magnum Opus, не иначе. Ещё не дочитал, но в лютом восторге. Спасибо вам за такую классную статью!

Filipp42
14.11.2023 09:32
#26161534
-1
Скажите пожалуйста, а можно ли мыло поднять частоту, скажем, до одного мегагерца? Можно ли было разогнать процессор?
1. mark_ablov Автор
  14.11.2023 09:32
  #26161706
  Моё предположение, что ограничения по частоте шли не столько от самого процессора, как от переферии. У меня был резервный план - попробовать разогнать i4040, если бы не получилось решить задачу только программно. Однако отмечу, что эмулятор памяти на stm32f4, написанный на C (O3, gcc), уже с трудом справлялся с 740кГц. Скорее всего, я сумел бы написать более эффективный код эмулятора, который поддерживает более высокую частоту, но до этого дело не дошло.

Ускоряем программу для 50-летнего процессора на 180000% +245

Введение

Технические ограничения

Теоретическая подготовка

Выбор алгоритма

Исходная формула для π

Вычисление n-ой цифры дроби

Вычисление n-ой цифры числового ряда

Наибольшая кратность простого числа в факторизации делителя

Определение достаточной верхней границы для частичной суммы числового ряда π

Алгоритм

Инструментарий для разработки

Препроцессор

Компилятор и компоновщик

Эмулятор и профилировщик

Первая попытка

Разрядность переменных

Логарифмы и арифметика чисел с плавающей запятой

Модульная арифметика

Результаты первого запуска

Оптимизируй это немедленно!

Алгоритм

"Параллельные" вычисления

Возвращаемся к простым числам

Небольшие изменения в формате дробных чисел

Упрощенная схема для больших простых чисел

Профилирование

Получение обратного числа по модулю

Модульное возведение в степень

Деление

Деление 4-битного числа на 4-битное

Умножение 4-битных чисел

Деление 16-битных чисел

Модульное умножение

Алгоритм умножения

Ленивое приведение по модулю

Еще более ленивое приведение по модулю

Сложение вместо вычитания

Нам нужно больше таблиц перехода!

Выполняем несколько последовательных обновлений multipliedFactor за одну операцию

Включаем в операцию сложения промежуточного результата еще и удваивание multipliedFactor

Более шустрая проверка на потенциальное переполнение переменных

Оптимизации, которые не вошли в конечный код

Структура памяти

Результаты запуска на эмуляторе

Аппаратная часть

Прошивка

Финальные результаты

Ссылки

Комментарии (34)

mark_ablov Автор

mark_ablov Автор

mark_ablov Автор

mark_ablov Автор