Я пытаюсь написать программу, которая рассчитывает десятичные цифры от π до 1000 цифр и более.
Чтобы практиковать низкоуровневое программирование для развлечения, конечная программа будет написана на ассемблере, на 8-битном процессоре, который не имеет умножения или деления и выполняет только 16-битные сложения. Чтобы упростить реализацию, желательно иметь возможность использовать только 16-разрядные целочисленные операции без знака и использовать итерационный алгоритм. Скорость не является серьезной проблемой. А быстрое умножение и деление выходят за рамки этого вопроса, поэтому не стоит также рассматривать эти проблемы.
Перед тем, как реализовать его в сборке, я все еще пытаюсь выяснить пригодный для использования алгоритм на C на моем настольном компьютере. До сих пор я обнаружил, что следующие серии достаточно эффективны и относительно просты в реализации.
Формула получена из ряда Лейбница с использованием метода ускорения сходимости. Чтобы получить ее, см. "Вычисление цифр в π" Карла Д. Оффнера (https://cs.umb.edu/~offner/files/pi.pdf), стр. 19-26. Окончательная формула показана на странице 26. Первоначальная формула, которую я написал, содержала некоторые опечатки, обновите страницу, чтобы увидеть фиксированную формулу. Постоянный член 2 для наибольшего члена объясняется на странице 54. В статье также описан продвинутый итерационный алгоритм, но я здесь не использовал его.
Если оценить ряд, используя много (например, 5000) терминов, можно легко получить тысячи цифр числа π, и я обнаружил, что этот ряд также легко оценить итеративно, используя этот алгоритм:
Алгоритм
- Сначала измените формулу, чтобы получить ее постоянные члены из массива.
-
Заполните массив 2, чтобы начать первую итерацию, поэтому новая формула напоминает исходную.
-
Пусть
carry = 0. -
Начните с самого большого срока. Получите одно слагаемое (2) из массива, умножьте его на
PRECISIONчтобы выполнить деление с фиксированной точкой на2 * я + 1, и сохраните напоминание как новый термин в массиве. Затем добавьте следующий термин. Теперь уменьшите значениеi, переходите к следующему члену, повторяйте, покаi == 1. Наконец добавьте последний членx_0. -
Поскольку используется 16-битное целое число,
PRECISIONравно10, следовательно, получаются 2 десятичных знака, но действительна только первая цифра. Сохраните вторую цифру как перенос. Показать первую цифру плюс перенос. -
x_0- это целое число 2, его нельзя добавлять для последовательных итераций, очистите его. -
Перейдите к шагу 4, чтобы вычислить следующую десятичную цифру, пока у нас не появятся все нужные цифры.
Реализация 1
Перевод этого алгоритма в C:
#include <stdio.h>
#include <stdint.h>
#define N 2160
#define PRECISION 10
uint16_t terms[N + 1] = {0};
int main(void)
{
/* initialize the initial terms */
for (size_t i = 0; i < N + 1; i++) {
terms[i] = 2;
}
uint16_t carry = 0;
for (size_t j = 0; j < N / 4; j++) {
uint16_t numerator = 0;
uint16_t denominator;
uint16_t digit;
for (size_t i = N; i > 0; i--) {
numerator += terms[i] * PRECISION;
denominator = 2 * i + 1;
terms[i] = numerator % denominator;
numerator /= denominator;
numerator *= i;
}
numerator += terms[0] * PRECISION;
digit = numerator / PRECISION + carry;
carry = numerator % PRECISION;
printf("%01u", digit);
/* constant term 2, only needed for the first iteration. */
terms[0] = 0;
}
putchar('\n');
}
Код может вычислять от π до 31 десятичного знака, пока не произойдет ошибка.
31415926535897932384626433832794
10 <-- wrong
Иногда digit + carry больше 9, поэтому требуется дополнительный перенос. Если нам очень не повезло, может быть даже двойной перенос, тройной перенос и т.д. Мы используем кольцевой буфер для хранения последних 4 цифр. Если обнаружен дополнительный перенос, мы выводим клавишу возврата, чтобы стереть предыдущую цифру, выполнить перенос и перепечатать их. Это просто уродливое решение Proof-of-Concept, которое не имеет отношения к моему вопросу о переполнении, но для полноты вот оно. Что-то лучшее будет реализовано в будущем.
Реализация 2 с повторным переносом
#include <stdio.h>
#include <stdint.h>
#define N 2160
#define PRECISION 10
#define BUF_SIZE 4
uint16_t terms[N + 1] = {0};
int main(void)
{
/* initialize the initial terms */
for (size_t i = 0; i < N + 1; i++) {
terms[i] = 2;
}
uint16_t carry = 0;
uint16_t digit[BUF_SIZE];
int8_t idx = 0;
for (size_t j = 0; j < N / 4; j++) {
uint16_t numerator = 0;
uint16_t denominator;
for (size_t i = N; i > 0; i--) {
numerator += terms[i] * PRECISION;
denominator = 2 * i + 1;
terms[i] = numerator % denominator;
numerator /= denominator;
numerator *= i;
}
numerator += terms[0] * PRECISION;
digit[idx] = numerator / PRECISION + carry;
/* over 9, needs at least one carry op. */
if (digit[idx] > 9) {
for (int i = 1; i <= 4; i++) {
if (i > 3) {
/* allow up to 3 consecutive carry ops */
fprintf(stderr, "ERROR: too many carry ops!\n");
return 1;
}
/* erase a digit */
putchar('\b');
/* carry */
digit[idx] -= 10;
idx--;
if (idx < 0) {
idx = BUF_SIZE - 1;
}
digit[idx]++;
if (digit[idx] < 10) {
/* done! reprint the digits */
for (int j = 0; j <= i; j++) {
printf("%01u", digit[idx]);
idx++;
if (idx > BUF_SIZE - 1) {
idx = 0;
}
}
break;
}
}
}
else {
printf("%01u", digit[idx]);
}
carry = numerator % PRECISION;
terms[0] = 0;
/* put an element to the ring buffer */
idx++;
if (idx > BUF_SIZE - 1) {
idx = 0;
}
}
putchar('\n');
}
Отлично, теперь программа может правильно вычислить 534 цифры π, пока не сделает ошибку.
3141592653589793238462643383279502884
1971693993751058209749445923078164062
8620899862803482534211706798214808651
3282306647093844609550582231725359408
1284811174502841027019385211055596446
2294895493038196442881097566593344612
8475648233786783165271201909145648566
9234603486104543266482133936072602491
4127372458700660631558817488152092096
2829254091715364367892590360011330530
5488204665213841469519415116094330572
7036575959195309218611738193261179310
5118548074462379962749567351885752724
8912279381830119491298336733624406566
43086021394946395
22421 <-- wrong
16-разрядное целочисленное переполнение
Оказывается, при вычислении самых больших слагаемых в начале, слагаемое ошибки становится довольно большим, поскольку делители в начале находятся в диапазоне ~ 4000. При оценке ряда numerator фактически начинает сразу переполняться при умножении.
Целочисленное переполнение незначительно при вычислении первых 500 цифр, но начинает ухудшаться и ухудшаться, пока не даст неверный результат.
Изменение uint16_t numerator = 0 на uint32_t numerator = 0 может решить эту проблему и вычислить цифры от π до 1000+.
Однако, как я упоминал ранее, моей целевой платформой является 8-битный процессор, и он выполняет только 16-битные операции. Есть ли хитрость для решения проблемы переполнения 16-битного целого числа, которую я здесь вижу, используя только один или несколько uint16_t? Если невозможно избежать арифметики с множественной точностью, какой самый простой способ реализовать это здесь? Я как-то знаю, что мне нужно ввести дополнительное 16-разрядное "слово расширения", но я не уверен, как это реализовать.
И заранее спасибо за ваше терпение, чтобы понять длинный контекст здесь.
