Вычисление числа Пи на ассемблере

[blackmirror]

В википедии есть очень симпатичный алгоритм, N=100 цифр можно вычислить за S=6 шагов, а N=1000 за S=9. Из сложных операций на каждом шаге требуется одно умножение и корень. Еще одно умножение и деление нужно чтобы получить Пи, и один корень для инициализации. То есть требуется одно деление, S+1 умножение и S+1 корень. Что касается деления, то его можно заменить умножением на обратную величину вычисляемую по итеративной формуле типа x=x*(2-x*y) за log(N) шагов. Начинать можно с некоторого приближения x0, на каждой итерации количество верных цифр будет удваиваться, полное умножение NxN потребуется только на последнем шаге, а все предыдущие итерации по времени выполнения примерно равны еще одному умножению. То есть финальное деление примерно эквивалентно 4 умножениям. Что касается корней, то 1/корень можно вычислить используя итеративную формулу x=x*(3-x*x*y). Здесь три умножения, добавив еще 1.5 на предыдущие итерации, и 1 чтобы получить корень, получим что корень эквивалентен 5.5 умножениям, в итоге получается что для 100 цифр требуется примерно 50 умножений, а для 1000 примерно 70. Возможно что для деления и корня есть и более простые варианты. Хорошо бы конечно делать умножения методом Карацубы, но 100 и 1000 десятичных цифр занимают не очень круглое количество байт и не факт, что дробление блоков и лишние проверки не дадут фору какому-либо более примитивному алгоритму.

[Advertiser]

[avivanov76]

А про бизнес-приложения и говорить не приходится

Помню, как-то сгоряча переписал пару функций из бизнес-логики на асме, благо на Дельфи это легко делалось. Потом проследил всю цепочку вызовов и понял, что сэкономил пару микросекунд, в то время как дальше (не в моем коде) шел SQL запрос, выполнявшийся минимум 50 миллисекунд

[perestoronin]

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

Аналогично, могу предположить, и на Бейсиках можно найти такие задачи, которые асму будут не по зубам.

Но у нас то, самая что ни на есть задачка для ассемблеров - расчет числа Пи

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

Ждем помощи математиков.

[b2m]

В википедии есть очень симпатичный алгоритм, N=100 цифр можно вычислить за S=6 шагов, а N=1000 за S=9.

Ну да, а шаги эти включают работу с числами с плавающей точкой, для 100 знаков минимум 45 байт мантиссы, для 1000 - 450. И при этом, вычисление квадратного корня нужно будет в такой длинный ряд разложить (чтобы была нужная точность), что считаться он будет на 8-ми битках вечность. Реализованные в этой теме алгоритмы использовали лишь 16/32 битные числа.

[ivagor]

Потер два предыдущих поста насчет оптимизации деления для z80 и 8085.
В итоге лучший вариант такой: в финальной версии для вектора не трогаем de, а меняем знак самого bc. Фрагменты "jnc DIV0_1\ dad h\ inr l" меняем на adc hl,hl или xchg\ rdel\ xchg. dad d меняем на dad b, а dad b меняем на sbc hl,bc или dsub. Лишнее убираем - древовидность становится не нужна и размер резко сокращается. Также отпадает необходимость в самомодифицирующемся коде (в процедуре деления). Такие версии показывают практически одинаковое быстродействие.

[Titus]

Так а где сам исходник самого быстрого в мире Пи и деления для Z80?)

[blackmirror]

b2m
Для 1000 цифр хватает мантиссы в 420 байт, если длина мантиссы M, а целая часть в нулевом байте, то при перемножении чисел A и B, можно не перемножать байты A[i] и B[j] если i+j>M, в этом случае из-за неучтённых переносов мы теряем только 3-4 младших байта.
Что касается вычисления X=1/Y, то имея приближение X=1/Y+R, где |R|<2^-n, и сделав итерацию x=x(2-xy), получаем X=(1/Y+R)(2-(1/Y+R)*Y)=(1/Y+R)(1-R*Y)=1/Y-R*R*Y, и если Y<1, то ошибка станет менее 2^-2n, иными словами количество верных цифр удваивается. Если Y>1, то можно на предыдущем шаге обрабатывать не половину цифр, а чуть больше. К примеру использовать умножения 45x45, 23x23, 12x12, 7x7 и так далее. С обратным корнем ситуация аналогичная, только умножений там чуть больше.

[ivagor]

Деление 32/16=(16;16) для z80 div320_z80_correct2.zip
Легко адаптируется для 8085 - меняем adc hl,hl на xchg\ rdel\ xchg, sbc hl,bc на dsub.
Или для 580ВМ1 - меняем adc hl,hl на cs\ dad h, sbc hl,bc на dsub b.

- - - Добавлено - - -

При необходимости можно еще добавить в начале изменение знака bc.

UPD: заменил архив - добавил "выпавший" adc a,a. У меня эта процедура в программе в мнемониках 8080 (коды z80 в виде .db), команда выпала при переводе мнемоник 8080->z80.

UPD2: ёпрст, вроде НГ не начался, а пришлось еще раз исправлять

- - - Добавлено - - -

Еще вариант для 8085 - div320_8085.zip
Эту процедуру просто скопировал из программы, ничего не менял, поэтому ошибок сразу нет.

[Titus]

Деление 32/16=(16;16) для z80 div320_z80_correct2.zip

А чего там в какой-то предыдущей версии стек использовался, а тут нет?

[ivagor]

Если речь про использование dad sp (add hl,sp), то оказалось, что без этого можно обойтись и без sp будет даже быстрее (при данной разрядности).
Для z80 можно с использованием sp апгрейдить эту процедуру до получения 32битного остатка. Младшие байты делимого при этом придется тасовать самомодифицирующимся кодом, как в версии для 8080. Зато можно увеличить и разрядность делимого и довести процедуру, например, до 64/32=(32;32)
С другой стороны для расчета пи в пределах первых тысяч знаков это не особо нужно. Если использовать формулу Гаусса, то хватит даже 16 битных процедур (но вроде знаковых, надо уточнить) для расчета 4000 с чем-то знаков (в этой теме про это упоминал).

- - - Добавлено - - -

но вроде знаковых, надо уточнить

Уточнил - достаточно беззнакового деления и умножения.