Я пытаюсь написать программу, которая вычисляет десятичные цифры от π до 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 * i + 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-разрядное «слово расширения», но я не уверен, как его реализовать.
И заранее спасибо за ваше терпение, чтобы понять длинный контекст здесь.