我正在尝试编写一个程序，将 π 的十进制数字计算为 1000 位或更多位。
为了有趣地练习低级编程，最终程序将用汇编编写，在没有乘法或除法的 8 位 CPU 上，只执行 16 位加法。为了简化实现，希望能够仅使用 16 位无符号整数运算，并使用迭代算法。速度不是主要问题。而快速乘法和除法超出了这个问题的范围，所以也不要考虑这些问题。
在汇编中实现它之前，我仍然试图在我的台式计算机上找出一个可用的 C 算法。到目前为止，我发现以下系列相当有效且相对容易实现。
该公式是使用收敛加速技术从莱布尼茨系列推导出来的，要推导它，请参阅 Carl D. Offner 的 Computing the Digits in π ( 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 位数字。如果检测到一个额外的进位，我们输出一个退格来擦除前一个数字，执行一个进位，然后重新打印它们。这只是概念证明的一个丑陋的解决方案，这与我关于溢出的问题无关，但为了完整起见，这就是它。将来会实施更好的东西。
重复进位的实现 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 位操作。是否有一个技巧可以解决我在这里看到的 16 位整数溢出问题，只使用一个或多个 uint16_t？如果无法避免多精度算术，那么在这里实现它的最简单方法是什么？我知道我需要引入一个额外的 16 位“扩展字”，但我不确定如何实现它。
并提前感谢您耐心地理解这里的长篇大论。
