Avoid Overflow when Calculating π by Evaluating a Series Using 16-bit Arithmetic?

Take a look at related QA:

Its using Wiki: Bailey–Borwein–Plouffe_formula which is more suited for integer arithmetics.

The real challenge however would be:

As you probably want to print the number in dec base …

Also if you need carry in higher level language than asm take a look at this:

You can modify it to handle as many carry bits as you need (if still less than the data type bit-width).

[Edit1] BBP example in C++/VCL

I used this formula (taken from Wiki page linked above):

BBP

converted to fixed point…

//---------------------------------------------------------------------------
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);
    }
//---------------------------------------------------------------------------

The code is using VCL AnsiString which is a self allocating string and mine uint<N> template which is unsigned integer arithmetics of 32*N bitwidth based on mine ALU32. As you can see you only need big integer division addition and multiplication for this (all the other stuff is doable on normal integers).

Here decadic result versus 1000 digit Pi reference:

ref: 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989
BPP: 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259036001133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912983367336244065664308602139494639522473719070217986094370277053921717629317675238467481846766940513200056812714526356082778577134275778960917363717872146844090122495343014654958537105079227968925892354201995611212902196086403441815981362977477130996051870721134999999837297804995105973173281609631859502445945534690830264252230825334468503526193118817101000313783875288658753320838142061717766914730359825349042875546873115956286388235378759375195778185778048187

The computed bigint value is exported to hex string and then converted to decadic base using str_hex2dec from link above. The number of iterations depends on the target bitwidth.

The code is not optimized yet…

Leave a Comment