计算时避免溢出？通过使用16位算术评估系列？

Question

我正在尝试编写一个计算十进制数字的程序？到1000位数或更多。

为了有趣地进行低级编程，最终程序将以汇编形式编写在没有乘法或除法且仅执行16位加法的8位CPU上。为了简化实现，希望只能使用16位无符号整数运算，并使用迭代算法。速度不是主要问题。快速乘法和除法不在此问题的范围内，因此也不要考虑这些问题。

在进行汇编之前，我仍在尝试在台式计算机上的C语言中找到可用的算法。到目前为止，我发现以下系列相当有效并且相对容易实现。

该公式是使用收敛加速技术从Leibniz系列推导而来的，要推导该公式，请参见Carl D.Offner撰写的《计算？中的数字》（https://cs.umb.edu/~offner/files/pi.pdf） ，第19-26页。最终公式显示在第26页。我编写的初始公式有错别字，请刷新页面以查看固定公式。2第54页中解释了最大术语中的常数项。该白皮书还介绍了一种高级迭代算法，但是我在这里没有使用它。

如果一个人使用许多（例如5000）项来评估该系列，则有可能获得千位数的？。容易，并且我发现使用此算法，该系列也易于迭代评估：

算法

1. 首先，重新排列公式以从数组中获取其常数项。

2. 用2填充数组以开始第一次迭代，因此新公式类似于原始公式。

3. 让carry = 0。

4. 从最伟大的任期开始。从数组中获得一项（2），将其乘以PRECISION对进行定点除法2 * i + 1，然后将提示作为新项保存到数组中。然后添加下一项。现在递减i，转到下一学期，重复直到i == 1。最后加上最后一项x_0。

5. 由于使用的PRECISION是16位整数10，因此，将获得2个十进制数字，但只有第一个数字有效。将第二个数字另存为进位。显示第一个数字加进位。

6. x_0 是整数2，不应为后续迭代添加整数，将其清除。

7. 转到第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位数字。如果检测到多余的进位，我们将输出退格键以擦除前一位，执行进位并重新打印。这只是概念验证的丑陋解决方案，与我关于溢出的问题无关，但为了完整性，就在这里。将来会实施更好的方法。

具有重复进位的实现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位CPU，只有16位操作。是否有一个技巧可以解决，我仅使用一个或多个uint16_t解决了我在这里看到的16位整数溢出问题？如果无法避免使用多精度算法，那么在这里实现它的最简单方法是什么？我知道我需要引入一个额外的16位“扩展字”，但是我不确定如何实现它。

并在此先感谢您的耐心配合，以了解此处的详细情况。

Answer 1

有一个技巧：

考虑使用一个数组作为分子，另一个数组作为分母。每个位置代表该数字相乘得到实际数字的次数。

一个例子：

(1 * 2 * 3 * 7 * 7) / (3 * 6 * 8)

将表示为：

num[] = {1, 1, 1, 0, 0, 0, 2};
denom[] = {0, 0, 1, 0, 0, 1, 0, 1};

然后考虑在存储之前将每个数字分解为素数，这样您就可以得到更少的数字。现在您需要另一个数组来存储所有素数：

primes[] = {2, 3, 5, 7};
num[] = {1, 1, 0, 2};
denom[] = {4, 2, 0, 0};

这将允许您存储难以想象的大数字，但您迟早会想要将它们转换回数字，因此您需要首先简化它。方法就是factors[i] += num[i] - denom[i]对数组中的每个字段、系列中的每个分数进行减法。您将希望在每次迭代后进行简化，以便最大限度地减少溢出风险。

factors[] = {-3, -1, 0, 2};

当您需要数字时，只需num *= pow(primes[i], factors[i]);对数组中的每个字段的因子为正或num /= pow(primes, -factors[i]);负数即可。（如果为 0，则不执行任何操作。

num和denom是用于存储分数的临时数组，存储结果的数组是factors。请记住memset每次使用前的临时数组。

这个解释对于任何大分数都是有用的。为了适应您的具体问题，您可能需要使用整数幂函数，并乘以 10^something 将小数部分变成整数部分。这是你的使命，你应该接受它吗:)