c - 在通过使用 16 位算术评估一个系列来计算 π 时避免溢出？

Question

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

为了练习低级编程的乐趣，最终的程序将用汇编编写，在没有乘法或除法的 8 位 CPU 上，只执行 16 位加法。为了简化实现，最好只能使用 16 位无符号整数运算，并使用迭代算法。速度不是主要问题。而快速乘除法超出了这个问题的范围，所以不要考虑这些问题。

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

该公式是使用收敛加速技术从莱布尼茨级数推导出来的，要推导出它，请参阅计算 π 中的数字，Carl D. Offner ( 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 位操作。有没有什么技巧可以解决我在这里看到的 只使用一个或多个 uint16_t 的 16 位整数溢出问题？如果无法避免多精度算术，那么在这里实现它的最简单方法是什么？我知道我需要引入一个额外的 16 位“扩展词”，但我不确定如何实现它。

并提前感谢您耐心理解这里的长篇大论。

Answer 1

看看相关的QA:

Baking-Pi Challenge - Understanding & Improving

它使用更适合整数算术的 Wiki: Bailey–Borwein–Plouffe_formula 。

然而，真正的挑战是:

How do I convert a very long binary number to decimal? 。

因为您可能想以 dec base 打印数字...

另外，如果您需要使用比 asm 更高级别的语言，请查看以下内容:

Cant make value propagate through carry

您可以修改它以根据需要处理尽可能多的进位位(如果仍然小于数据类型位宽)。

[Edit1] C++/VCL 中的 BBP 示例

我使用了这个公式(取自上面链接的维基页面):



转换为定点...

//---------------------------------------------------------------------------
AnsiString str_hex2dec(const AnsiString &hex)
    {
    char c;
    AnsiString dec="",s;
    int i,j,l,ll,cy,val;
    int  i0,i1,i2,i3,sig;
    sig=+1; l=hex.Length();
    if (l) { c=hex[l]; if (c=='h') l--; if (c=='H') l--; }
    i0=0; i1=l; i2=0; i3=l;
    for (i=1;i<=l;i++)      // scan for parts of number
        {
        char c=hex[i];
        if (c=='-') sig=-sig;
        if ((c=='.')||(c==',')) i1=i-1;
        if ((c>='0')&&(c<='9')) { if (!i0) i0=i; if ((!i2)&&(i>i1)) i2=i; }
        if ((c>='A')&&(c<='F')) { if (!i0) i0=i; if ((!i2)&&(i>i1)) i2=i; }
        if ((c>='a')&&(c<='f')) { if (!i0) i0=i; if ((!i2)&&(i>i1)) i2=i; }
        }

    l=0; s=""; if (i0) for (i=i0;i<=i1;i++)
        {
        c=hex[i];
             if ((c>='0')&&(c<='9')) c-='0';
        else if ((c>='A')&&(c<='F')) c-='A'-10;
        else if ((c>='a')&&(c<='f')) c-='A'-10;
        for (cy=c,j=1;j<=l;j++)
            {
            val=(s[j]<<4)+cy;
            s[j]=val%10;
            cy  =val/10;
            }
        while (cy>0)
            {
            l++;
            s+=char(cy%10);
            cy/=10;
            }
        }
    if (s!="")
        {
        for (j=1;j<=l;j++) { c=s[j]; if (c<10) c+='0'; else c+='A'-10; s[j]=c; }
        for (i=l,j=1;j<i;j++,i--) { c=s[i]; s[i]=s[j]; s[j]=c; }
        dec+=s;
        }
    if (dec=="") dec="0";
    if (sig<0) dec="-"+dec;

    if (i2)
        {
        dec+='.';
        s=hex.SubString(i2,i3-i2+1);
        l=s.Length();
        for (i=1;i<=l;i++)
            {
            c=s[i];
                 if ((c>='0')&&(c<='9')) c-='0';
            else if ((c>='A')&&(c<='F')) c-='A'-10;
            else if ((c>='a')&&(c<='f')) c-='A'-10;
            s[i]=c;
            }
        ll=((l*1234)>>10);  // num of decimals to compute
        for (cy=0,i=1;i<=ll;i++)
            {
            for (cy=0,j=l;j>=1;j--)
                {
                val=s[j];
                val*=10;
                val+=cy;
                s[j]=val&15;
                cy=val>>4;
                }
            dec+=char(cy+'0');
            for (;;)
                {
                if (!l) break;;
                if (s[l]) break;
                l--;
                }
            if (!l) break;;
            }
        }

    return dec;
    }
//---------------------------------------------------------------------------
AnsiString pi_BBP() // https://en.wikipedia.org/wiki/Bailey–Borwein–Plouffe_formula
    {
    const int N=100;        // 32*N bit uint arithmetics
    int sh;
    AnsiString s;
    uint<N> pi,a,b,k,k2,k3,k4;

    for (pi=0,sh=(N<<5)-8,k=0;sh>=0;k++,sh-=4)
        {
        k2=k*k;
        k3=k2*k;
        k4=k3*k;
        a =k2* 120;
        a+=k * 151;
        a+=     47;
        b =k4* 512;
        b+=k3*1024;
        b+=k2* 712;
        b+=k * 194;
        b+=     15;
        a<<=sh;
        pi+=a/b;
        }
    pi<<=4;
    s=pi.strhex();
    s=s.Insert(".",2);
    return str_hex2dec(s);
    }
//---------------------------------------------------------------------------

该代码使用 VCL AnsiString，这是一个自分配字符串和我的 uint<N> 模板，它是基于我的 ALU32 的 32*N 位宽的无符号整数算法。正如你所看到的，你只需要大整数除法加法和乘法(所有其他东西都可以在普通整数上进行)。

这里十进制结果与 1000 位 Pi 引用:

ref: 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989
BPP: 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259036001133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912983367336244065664308602139494639522473719070217986094370277053921717629317675238467481846766940513200056812714526356082778577134275778960917363717872146844090122495343014654958537105079227968925892354201995611212902196086403441815981362977477130996051870721134999999837297804995105973173281609631859502445945534690830264252230825334468503526193118817101000313783875288658753320838142061717766914730359825349042875546873115956286388235378759375195778185778048187

计算出的 bigint 值导出为十六进制字符串，然后使用上面链接中的 str_hex2dec 转换为十进制基数。迭代次数取决于目标位宽。

代码还没有优化...