А вот это что:
Имхо, это добавление переноса от умножения предыдущего разряда на 10 и деления на 2i-1.Код:x = 10*a[i-1] + q * i;
Вид для печати
Есть формула, по которой считается очередная 16-ичная цифра числа Pi, т.е. фактически очередные 4 бита двоичного представления. И она не требует вычисления предыдущих цифр. Проблема только в том, что чтобы перевести бинарное число в десятичное, нужно много раз умножать на 10 всё число целиком. Т.е. от наличия дополнительных цифр в двоичном числе зависят цифры в десятичном представлении этого числа. А после перевода первых N цифр в десятичный вид уже нельзя добавить новые двоичные данные.
Получается, сначала нужно рассчитать все цифры 16-ичного (двоичного) числа, и только потом заняться переводом его в десятичное. Для 8-мибитного компьютера перевод двоичного числа в десятичный вид не менее трудоёмок (в плане времени вычисления), чем рассчёт каких либо других формул.
- - - Добавлено - - -
a[i] это цифра в другом счислении. Понятно, что умножая это число на 10 нужно учитывать ещё и перенос из другого разряда. Очередной множитель цифры этой системы сисления как раз равен i/(2*i-1), а q, как ты правильно заметил - это перенос. Просто умножение на i в одной строке, а деление - в другой.
Забавно, получается что:
просто задает число пи равное 2.222222...., а все остальное ниже - просто перевод в десятичную систему? А в какой системе исчисления пи равно 2.222...?Код:for(j=0;j<LEN;j++) a[j]=2;
Тут надо читать про spigot-алгоритм вычисления числа pi:
http://zx-pk.ru/attachment.php?attac...8&d=1447061792
А, тогда не совсем позиционная система исчисления (что предполагает фиксированное основание), а все-таки обычный ряд. Просто тут его вычисление удачно совместили с переводом в десятичную систему. А в массиве находятся не цифры, а обновляемые члены ряда.
Update: кстати, добавление нового члена ряда прикольно сделано, с отложенным переносом. Если новая полученная цифра не девятка, то переноса в предыдущую цифру не предполагается (неочевидно, как по мне - требует доказательства, но ладно), значит предыдущую цифру можно вывести. А если девятка - то копим число девяток подряд, вдруг произойдет перенос и его распространение.
Очередная и, надеюсь, последняя версия spigota - Вложение 54829
100 знаков - 19897238 тактов - 11.19 сек (это с расчетом таблицы!)
535 знаков - 551337754 тактов - 5 мин 10 сек
Обращаю внимание, что запускать нужно G4200
Можно еще чуть-чуть ускорить, если размещать таблицу для умножения с 0000h, но я не стал заниматься крохоборством.
Ну, молодец, молодец, ускорил в 5 раз :)
Теперь напиши это на псевдоассемблере и ещё ассемблер, генерирующий прогу для любого из компов, что есть у меня в эмуляторе :)
[/sarcasm detection off]
А тебя не смущает, что по ссылке, которую приводил litwr, программа на си считает 800 знаков за 1,458,354,526 тактов? Если попробовать сравнить, то 535 знаков она посчитала бы примерно за 1,458,354,526/2,2=662888420 тактов. Т.е. сейчас программа на асме 8080 наконец-то стала считать быстрее программы на си для z80. Как то не очень круто, если учесть, что в той проге используются и 32 битные операции. Алгоритм, который используется в сишном варианте по ссылке очень похож на spigot, но, если я правильно понимаю, переводит не по одной цифре, а по 4.
[/sarcasm detection on]
А это и есть тот же spigot только переводит в систему с основанием 10000. Ну и массив 32-битный, чтобы не переполнялись члены. Кстати, возможно, там ошибка - если будет 9999 и перенос из следующего разряда. Вероятность небольшая, но пи - число иррациональное и трансцендентное :), в нем все комбинации группы цифр рано или поздно встретятся.
LLVM с его IR (intermediate representation, по факту это псевдоассемблер), возможно "спасет отца русской демократии" :)
Если совсем грубо, то вероятность появления "цифры" означающей 9999 - порядка 1/10000, да еще и перенос должен быть из следующего разряда, то есть вероятность около 1/100К, поэтому ошибка от 100К знаков может вылезти. Это если она там есть (вдруг ряд сходится так быстро что перенос не возникает), пока это только мое подозрение.
http://bio.fizteh.ru/student/spravoc...rpfelw5bsu.zip
999999,пятаяседьмая строка
Если верить архиву с миллиардом знаков числа пи, который скачал на рутракере, та программа правильно считает максимум 16256 знаков (за 1.7 секунды). Для 8биток выше крыши
- - - Добавлено - - -
Если заменить long на int64_t, то правильно считает до 54932 знака включительно (за 25 секунд).
Смущает. Там есть одна хитрость, которая позволит ускорить наш spigot ( уже наш :) ) ровно в два раза.
- - - Добавлено - - -
Может просто тупо тот алгоритм на асме переписать? Хотя, основная масса времени это всё равно умножение/деление, не думаю, что в сишном рантайме они какие-то неоптимальные.
- - - Добавлено - - -
Если только использовать не чисто 32-битные умножение/деление, а смешанные 32/16 битные. В том алгоритме явно указано, что нужно умножать uint32 на uint32, ну и делить также. А на самом деле нужно uint16*uint16 -> uint32, а затем uint32/uint16 -> uint16.
У тебя есть готовые арифметические процедуры для 32 разрядных? Если да, то я за "тупое переписывание". Пусть даже будет медленнее, чем pirk20, зато больше цифр можно считать.
У меня были процедурки для 32битных, но они, насколько помню, в архиве на другом компе, который сдох (но винт целый).
Была старинная математическая библиотека для 8080, но не могу вспомнить как назывался архив или хотя бы кто автор.
Возможно у vinxru есть готовые процедуры для 32-разрядных для его си?
Конечно можно и заново написать.
- - - Добавлено - - -
В Гуртовцеве Гудыменко есть готовые процедуры для 16*16=32 и 32/16
Вот это можно как-то ускорить (не разворачивая цикл)?
Код:MUL32: MOV B,H ; DEHL = DE*HL
MOV C,L
LXI H,0
MVI A,17
MUL321: DCR A
RZ
XCHG
DAD H
XCHG
JC MUL322
DAD H
JNC MUL321
INX D
JMP MUL321
MUL322: DAD H
JNC $+4
INX D
DAD B
JNC MUL321
INX D
JMP MUL321
Эх, для меня это все очень сложно, качать мегабайтные архивы, сверять :)
Есть путь попроще:
Сначала оно прикольно упало, когда я NUMS поставил 280000 :), стек коротковат оказался, пришлось static вкрутить.Код:#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdarg.h>
#include <windows.h>
#define NUMS 280000
int main()
{
static long r[NUMS + 1];
long i, k;
long b, d;
long c;
long x;
c = 0;
for (i = 0; i < NUMS; i++)
r[i] = 2000;
for (k = NUMS; k > 0; k -= 14)
{
d = 0;
i = k;
while (1)
{
d += r[i] * 10000;
b = i * 2 - 1;
r[i] = d % b;
d /= b;
i--;
if (i == 0) break;
d *= i;
}
x = c + d / 10000;
if (x >= 10000)
{
printf("\r\nERROR (%d, %d)", x, k);
break;
}
printf("%.4d", (int)(c + d / 10000));
c = d % 10000;
}
return 0;
}
А потом оно секунд за 10 привело к ошибке:
ERROR (10000, 119266)
Ожидаем четыре цифры, а оказывается их пять, и пятая-то, вот досада, должна быть перенесена в уже выведенное, косяк-с, мда.
Итого - ошибка есть, возникает на 45-тысячной с копейками цифре, если кому надо меньше цифр (ограничен там восьмибитками или еще чего) - то и так сгодится
Вот из вышеупомянутой книги, как ни странно - работает :)
Код:DIV32: MOV A,B ; DE = HLDE/BC, HL = HLDE%BC
CMA
MOV B,A
MOV A,C
CMA
MOV C,A
INX B
XRA A
DIV321: DAD H
RAR
XCHG
DAD H
XCHG
JNC $+4
INX H
PUSH H
DAD B
JNC DIV322
RAL
DIV323: INX D
INX SP
INX SP
ADI 10h
JNC DIV321
RET
DIV322: RAL
JC DIV323
POP H
ADI 10h
JNC DIV321
RET
Подобные вещи и в pirk катят, но убрал, чтобы работало и при разрешенных прерываниях (там где они есть).
- - - Добавлено - - -
Это я к тому, что желательно этот фрагмент потом заменить
- - - Добавлено - - -
Пардон, здесь проблем не будет (у меня несколько другая ситуация была).
Я тоже не согласен с 45000 :biggrin:
Взял Пи вот тут
и сравнил с вычисленным по 4-цифровому шпигату (по ссылке что давал litwr и в моем листинге), оно навернулось на 16278 цифре - начиная с нее не совпадает с эталонными цифрами. То есть, найденная мной ошибка - там не единственная.
Update:
там банальное целочисленное переполнение происходит, то есть - не совсем ошибка, просто неприменимость в данных условиях. Поменял типы на 64-битные - доехало до 50К, потом опять ERROR
Update2:
Правильная последовательность с 54920-ой цифры:
61527146690058147000026330
Вычисленная по 4-цифровому шпигату (с int64):
61527146690058146000026330
То есть - перенос 10000 не отработал, ошибка в наличии.
Нехватает однако 16*16=32 и 32/16=16. Наткнулся на итерацию, в которой после деления 17 бит получается, дальше можно не проверять. Нужны нормальные 32-битные процедуры, но тут регистров уже не хватит. Да и на первый взгляд, медленновато пока получается (пусть даже и неверный результат). Точно не засекал, но даже на глаз - медленнее. Так что первый шпигат - для восьмибиток пока лучше.
Я так понимаю, что нужна процедура переваривающая 32 битные: делимое, делитель и частное. Остаток м.б. 16 битным
Cделал код для 6502. На Коммодоре+4 800 знаков отсчитало за 6 мин 32 сек. Там по cсылке есть указания на алгоритм, по которому миллиарды знаков отсчитал некто Чудновский, с 32-разрядными базовыми константами. Но кода не искал, может где и есть в сети.
Для затравки процедура 32 битного деления для 8080 - Вложение 54857
Все 32битное - делитель, делимое, частное и остаток. Очень медленно. Единственный плюс - легко масштабировать хоть до 256 бит.
С делением я разобрался, доделал до 32-битного частного и 16-битного остатка. Сделал так: если надо, деление будет выполнено дважды (для старших двух байт и для младших двух байт). Умножение оставил как есть, но сделал цикл сложением для старшего байта 24-битных чисел. Вобщем, пока считает 100 знаков за 12,5 сек.
Ну и для 536 и 800 знаков результат очень хороший, особенно если учесть, как реализованы умножение и деление.
Но но, папрашу без фамильярнасти!
Фактически у тебя получилось вручную откомпилировать для 8080 с несколько большей эффективностью, чем тот си для z80. И это еще без всяких развертываний циклов и таблиц. Интересно, будет ли эффективно дальнейшее увеличение основания? Хотя для 8080 вряд ли, с учетом роста сложности реализации умножения и деления.
Развернул циклы умножения и деления, теперь 100 знаков 10.644 сек.
Видимо нет. Если считать по 5 знаков за раз, то делить/умножать нужно будет уже на 100000, а оно в 16 бит не полезет. Эффективные процедуры для "честных" 32-битных умножения/деления я вряд-ли сделаю. Да и неизвестно, что там будет после умножения на 100000, может произведение и в 32 бита не влезет.
Тут похоже никого кроме как с РК86 и нет. Кроме бкашников, возможно никто и не знает. Но прикрепляю программку и снимок экрана.
800 знаков за 5м 41с
100 знаков за 5.8c
Использовалось умножение с 2К таблицей и деление с раскрученными циклами.
На БК должно хорошо деление получится, а с умножением будет плохо - медленно там с байтовыми таблицами.
Вложение 54870
Вложение 54869
Результат 100 знаков за 10 сек. на Радио-РК впечатляет. Это на 2 МГц? Популярнейший Коммодор 64 (как и АГАТ) с этой программкой не должен справиться быстрее чем за сек 14.
Проверил программку на типовом ПК c AMD 3.4 МГц и N=210000, никаких ошибок переполнения d не было, но расхождение на 32375-й цифре. Пи брал с программы Reduce (https://en.wikipedia.org/wiki/Reduce...lgebra_system)).
Насчет программы с арктангенсами. Она работает, но последние цифры выдает неправильно. Т.е. например, чтобы правильно посчитать 1000 цифр, нужно задать расчет 1004 цифр. И, похоже, с определенной цифры она начинает совсем неправильно считать. Но, повторюсь, 1000 цифр посчитать ею можно. Причем, если уменьшить константы
short B=10; /* Working base */
short LB=1; /* Log10(base) */
short MaxDiv=57; /* about sqrt(2^15/B) */
то все long можно заменить на short
- - - Добавлено - - -
C short считает правильно до 4138 цифр
- - - Добавлено - - -
В исходном виде (с long) считает заметно лучше spigotа с основанием 10000, который с 64-битной арифметикой считал правильно только до 54392 цифр. Арккотангенсная программа правильно посчитала и 110000, наверняка и больше правильно посчитает, но это уже долговато.
- - - Добавлено - - -
Ну и 200000 знаков правильно посчитала почти за 5 минут. Дальше не буду пробовать, это уже и для 16биток многовато
8080 совсем неплох для арифметики, z80 тут ничего почти не улучшает, только частота повыше.
Прогнал вашу программку (с минимальным подгоном под систему) на Amstrad CPC6128 c z80 на 4 мгц (реально 3.2).
Вложение 54874Вложение 54875
Результаты: 100 знаков - 6.14 с, 800 - 6м 4.91c, 1000 - 9м 26.82c
Интересно будет сравнить со Спектрумом...
Ha Commodore+4 1000 знаков за 9м 1.04c (программку чуть улучшил - Вложение 54876).
Получается, что на "тяжёлой" арифметике 6502 только раза в полтора побыстрее на одинаковой частоте.
Под z80 можно заметно оптимизировать, да и для 8080 предел не достигнут.
Для каждой точности есть свой вариант уточненной формулы Мэчина.
Самая простая формула обладает невысокой точностью, но при этом является более быстрой.
https://en.wikipedia.org/wiki/Approximations_of_π
Думаю, что упущен момент с расчетом неких констант, выраженный в этой программе на бейсике:
- - - Добавлено - - -Код:; This BASIC program was used to calculate the equates for PIF.ASM.
; DEFDBL A-Z
; INPUT "Digits required"; A
; MPVSize = 32 * ((A / (LOG(2) / LOG(10)) + 255) \ 256)
; PRINT "NumDigits = "; A + 64
; PRINT "MPVSize = "; MPVSize
; k = A / .69897
; k = INT(15 + (k + k * .1)) 'floor 1
; PRINT "Last1 = "; k
; k = INT(15 + (A / 2.37)) 'floor 2
; PRINT "Last2 = "; k
;
numDigits=164
mpvSize=74
last1=172
last2=57
Реализация на Си далека от идеала. Смотрите реализацию на асме 286 и возможно ли её перенести на 8080 без потери эффективности ?
https://sites.google.com/site/richge...-pi-calculator