電子ブロック工房 IC・トランジスタで出来たコンピューターを設計・製作するためのブログ

C++で。

#include "stdafx.h"

#include <stdlib.h>

#define FIG 37

void prt(unsigned char* v1){
    unsigned char i;
    for (i=0;i<FIG-1;i++) {
        if (v1[i]<10) printf("0%d",v1[i]);
        else printf("%d",v1[i]);
    }
    printf("\n");
}

void copy(unsigned char* res, unsigned char* v1){
    unsigned char i;
    unsigned short t,c=0;
    v1[FIG-1]=10*((v1[FIG-1]+5)/10);
    for (i=FIG-1;i<FIG;i--) {
        t=v1[i]+c;
        c=t/100;
        if (i) t-=c*100;
        res[i]=t;
    }
}

void clear(unsigned char* res){
    unsigned char i;
    for (i=0;i<FIG;i++) res[i]=0;    
}

void add(unsigned char* res, unsigned char* v1, unsigned char* v2){
    unsigned char cache[FIG];
    unsigned char i;
    unsigned short t,c=0;
    for (i=FIG-1;i<FIG;i--) {
        t=v1[i]+v2[i]+c;
        c=t/100;
        if (i) t-=c*100;
        cache[i]=t;
    }
    copy(res,(unsigned char*)&cache);
}

void mul(unsigned char* res, unsigned char* v1, unsigned char v2){
    unsigned char cache[FIG];
    unsigned char i;
    unsigned short t,c=0;
    for (i=FIG-1;i<FIG;i--) {
        t=v1[i]*v2+c;
        c=t/100;
        if (i) t-=c*100;
        cache[i]=t;
    }
    copy(res,(unsigned char*)&cache);
}

void mul(unsigned char* res, unsigned char* v1, unsigned char* v2){
    unsigned char cache[FIG*3];
    unsigned char* r=&cache[0];
    unsigned char* t=&cache[FIG*2];
    unsigned char i;
    clear(r);
    clear(r+FIG);
    for (i=FIG-1;i<FIG;i--) {
        mul(t,v1,v2[i]);
        add(r+i,r+i,t);
    }
    copy(res,r);
}
void div(unsigned char* res, unsigned char* v1, unsigned char v2){
    unsigned char cache[FIG];
    unsigned char* r=&cache[0];
    unsigned char i;
    unsigned short c=0;
    for (i=0;i<FIG;i++) {
        c=c*100+v1[i];
        r[i]=c/v2;
        c=c-r[i]*v2;
    }
    copy(res,r);
}

int _tmain(int argc, _TCHAR* argv[])
{
    unsigned char* pi=new unsigned char[FIG];
    unsigned char* y=new unsigned char[FIG];
    unsigned char* t=new unsigned char[FIG];
    clear(pi);
    prt(pi);
    for (int n=0;n<127;n++) {
        clear(y);
        y[0]=3;
        div(y,y,2*n+1);
        for (int x=1;x<=n;x++) {
            clear(t);
            t[0]=x+n;
            div(t,t,x);
            div(t,t,16);
            mul(y,y,t);
        }
        add(pi,pi,y);
        prt(pi);
    }
    return 0;
}

copy関数内で丸めの処理をし、下２桁は表示しないようにしてある。このアルゴリズムだとこの辺が限界らしく、桁数をさらに上げると誤差が出るようだ。PICのpi calculatorにも入れてみよう。