Pi with dc

pdp11のエミュレータを探してググっていたら、とんでもないものをみつけた。
Pi with dcLink 」Julius Schmidt:German natural science student interested in science。unix に標準的に付属しているスタック型(逆ポーランド記法)電卓 dc で実装した spigot algorithm の円周率計算プログラムである。一桁づつメモリのある限り動作し続ける。遅いけど。
echo '1sq180sr60st2si[3li*1+d1+*3*suli27*12-lq*5lr*+lt5*/d48+Psy10lqlid2*1-***10lulqli5*2-*lr+lylt*-**srsqlult*stli1+silmx]smlmx' | dc
このページからもリンクされているが、原典は「Unbounded Spigot Algorithms for the Digits of PiLink 」JEREMY GIBBONS:University Lecturer in Software Engineering and Continuing Education at the University of Oxford
> piG3 = g(1,180,60,2) where 
>   g(q,r,t,i) = let (u,y)=(3*(3*i+1)*(3*i+2),div(q*(27*i-12)+5*r)(5*t)) 
>                in y : g(10*q*i*(2*i-1),10*u*(q*(5*i-2)+r-y*t),t*u,i+1)
げげ、Haskell じゃぁないか。
root@jesse:~# apt-get install haskell-platform
パッケージリストを読み込んでいます... 完了
依存関係ツリーを作成しています                
状態情報を読み取っています... 完了
提案パッケージ:
  haskell-platform-doc haskell-platform-prof
以下のパッケージが新たにインストールされます:
  haskell-platform
アップグレード: 0 個、新規インストール: 1 個、削除: 0 個、保留: 3 個。
7,340 B 中 0 B のアーカイブを取得する必要があります。
この操作後に追加で 7,168 B のディスク容量が消費されます。
以前に未選択のパッケージ haskell-platform を選択しています。
(データベースを読み込んでいます ... 現在 187646 個のファイルとディレクトリがインストールされています。)
.../haskell-platform_2013.2.0.0.debian12_all.deb を展開する準備をしています ...
haskell-platform (2013.2.0.0.debian12) を展開しています...
haskell-platform (2013.2.0.0.debian12) を設定しています ...
root@jesse:~# ^dログアウト
hiroaki@jesse:~$ cd pi
hiroaki@jesse:~/pi$ cat piG3.hs
>piG3 = g(1,180,60,2) where
>   g(q,r,t,i) = let (u,y)=(3*(3*i+1)*(3*i+2),div(q*(27*i-12)+5*r)(5*t))
>                in y : g(10*q*i*(2*i-1),10*u*(q*(5*i-2)+r-y*t),t*u,i+1)
main = print piG3
hiroaki@jesse:~/pi$ ghc --make piG3
[1 of 1] Compiling Main             ( piG3.hs, piG3.o )
Linking piG3 ...
hiroaki@jesse:~/pi$ ./piG3
[3,1,4,1,5,9,2,6,5,3,5,8,9,7,9,3,2,3,8,4,6,2,6,4,3,3,8,3,2,7,9,5,0,2,8,8,4,1,9,7,1,6,9,3,9,9,3,7,5,1,0,5,8,2,0,9,7,4,9,4,4,5,9,2,3,0,7,8,1,6,4,0,6,2,8,6,2,0,8,9,9,8,6,2,8,0,3,4,8,2,5,3,^c
げげ、リスト表示だ。こうなったらC言語で実装してやる。
hiroaki@jesse:~/pi$ cat piG3.c
// Pi by spigot algorithm.
//
// [1] JEREMY GIBBONS, 
// "Unbounded Spigot Algorithms for the Digits of Pi"(2005/07/29)
// http://www.cs.ox.ac.uk/people/jeremy.gibbons/publications/spigot.pdfLink 
// JEREMY GIBBONS is a University Lecturer in Software Engineering and Continuing Education at the University of Oxford.
// > piG3 = g(1,180,60,2) where 
// > g(q,r,t,i) = let (u,y)=(3*(3*i+1)*(3*i+2),div(q*(27*i-12)+5*r)(5*t)) 
// > in y : g(10*q*i*(2*i-1),10*u*(q*(5*i-2)+r-y*t),t*u,i+1)
//
// [2] Julius Schmidt, "AIJU'S HOMEPAGE/Code/Pi"
// http://aiju.de/code/piLink 
// Pi with dc
// echo '1sq180sr60st2si[3li*1+d1+*3*suli27*12-lq*5lr*+lt5*/d48+Psy10lqlid2*1-***10lulqli5*2-*lr+lylt*-**srsqlult*stli1+silmx]smlmx' | dc
// It's a great implementation.
// 
// Implementation of C with GMP by Hiroaki Nitobe. 2016/02/07
//

#include        <stdio.h>
#include        <signal.h>
#include        <gmp.h>
#include        "time_tools.h"

static int detect_sigint = 0;

void signal_handler(int signo){
        if(signo==SIGINT){
        fprintf(stderr, "\nReceived SIGINT\n");
        detect_sigint = 1;
}
}

void g(void){
        mpz_t        q, r, t, y, temp;
        unsigned long int        i, u, d;

        mpz_init(q);
        mpz_init(r);
        mpz_init(t);
        mpz_init(y);
        mpz_init(temp);

        // q=1
        mpz_set_ui(q, 1UL);
        // r=180
        mpz_set_ui(r, 180UL);
        // t=60
        mpz_set_ui(t, 60UL);
        // i=2
        i = 2UL;
        // d=0
        d = 0UL;

        while(1) {
        // u=3*(3*i+1)*(3*i+2)
                u=3UL*(3UL*i+1UL)*(3UL*i+2UL);
        // y=(q*(27*i-12)+5*r)/(5*t)
                mpz_mul_ui(y, q, 27UL*i-12UL);
                mpz_mul_ui(temp, r, 5UL); 
                mpz_add(y, y, temp);
                mpz_mul_ui(temp, t, 5UL);
                mpz_fdiv_q(y, y, temp);
                mpz_out_str(stdout, 10, y);
                                if(d%100UL==0UL){
                                        if(d==0){
                                                printf(".\n");
                                        }else{
                                                printf("\n");
                                        }
                                        if(d%10000UL==0UL){
                                                fprintf(stderr, "%lu, ", d);
                                                my_timestamp(stderr);
                                                if(detect_sigint==1)break;
                                        }
                                }

        // next r = 10*u*(q*(5*i-2)+r-y*t)
                mpz_mul_ui(temp, q, 5UL*i-2UL);
                mpz_add(r, r, temp);
                mpz_mul(temp, y, t);
                mpz_sub(r, r, temp);
                mpz_mul_ui(r, r, 10UL*u);
        // next q = 10*q*i*(2*i-1)
                mpz_mul_ui(q, q, 10UL*i*(2UL*i-1));
        // next t = t*u
                mpz_mul_ui(t, t, u);
        // next i = i+1
                i = i + 1;
                d = d + 1;
        }
        mpz_clear(q);
        mpz_clear(r);
        mpz_clear(t);
        mpz_clear(y);
        mpz_clear(temp);
        fprintf(stderr, "Normal end.\n");
}

int main(int argc, char *argv[]){
        if(signal(SIGINT, signal_handler) == SIG_ERR){
                fprintf(stderr, "Can't catch SIGINT\n");
        }

        fprintf(stderr, "with gcc-%d.%d.%d, gmp-%s\n",
                        __GNUC__, __GNUC_MINOR__, __GNUC_PATCHLEVEL__, gmp_version);
                                my_date(stderr);
                                my_set_start(stderr);

        g();
                                my_date(stderr);
        return 0;
}
hiroaki@jesse:~/pi$ 
hiroaki@jesse:~/pi$ cat time_tools.c
#include        <stdio.h>
#include        <stdlib.h>
#include        <time.h>
#include        <sys/time.h>
#include        <sys/types.h>
#include        <unistd.h>

void my_date(FILE *fp);
double my_get_time(void);
void my_set_start(FILE *fp);
void my_timestamp(FILE *fp);

static double        start, lap;

void my_date(FILE *fp) {
        time_t        timer;
        time(&timer);
        fprintf(fp, "%s", ctime(&timer));
}

double my_get_time(void) {
        struct timeval        tv;

        gettimeofday(&tv, NULL);

        return (tv.tv_sec + tv.tv_usec * 1e-6);
}

void my_set_start(FILE *fp) {
        start = lap = my_get_time();
        fprintf(fp, "split: %f, lap: %f\n", start - start, lap - lap);
}

void my_timestamp(FILE *fp) {
        double        now;

        now = my_get_time();

        fprintf(fp, "split: %f, lap: %f\n", now - start, now - lap);

        lap = now;
}
hiroaki@jesse:~/pi$ 
hiroaki@jesse:~/pi$ cat time_tools.h
extern void my_date(FILE *fp);
extern double my_get_time(void);
extern void my_set_start(FILE *fp);
extern void my_timestamp(FILE *fp);
hiroaki@jesse:~/pi$ 
hiroaki@jesse:~/pi$ gcc -l gmp -o piG3 piG3.c time_tools.c
hiroaki@jesse:~/pi$ ./piG3 > piG.txt
with gcc-4.9.2, gmp-6.0.0
Thu Feb 11 17:50:18 2016
split: 0.000000, lap: 0.000000
0, split: 0.000048, lap: 0.000048
10000, split: 0.260359, lap: 0.260311
20000, split: 1.232651, lap: 0.972292
30000, split: 3.073861, lap: 1.841210
40000, split: 6.088222, lap: 3.014361
^C
Received SIGINT
50000, split: 10.240260, lap: 4.152038
Normal end.
Thu Feb 11 17:50:28 2016
hiroaki@jesse:~/pi$ 
んん。動いてるみたいだ。
さて、最初の dc 版をみてみよう。
 
echo '1sq180sr60st2si[3li*1+d1+*3*suli27*12-lq*5lr*+lt5*/d48+Psy10lqlid2*1-***10lulqli5*2-*lr+lylt*-**srsqlult*stli1+silmx]smlmx' | dc

1sq 1をqレジスタに積む g(1,180,60,2)
180sr 180をrレジスタに積む g(1,180,60,2)
60st 60をtレジスタに積む g(1,180,60,2)
2si 2をiレジスタに積む g(1,180,60,2)
[ ここから  
  3li* 3 とiレジスタを掛けてメインスタックに積む 3*(3*i+1)*(3*i+2)
  1+ メインスタックに1を足す 3*(3*i+1)*(3*i+2)
  d メインスタックを複製して積む 3*(3*i+1)*(3*i+1+1)
  1+ メインスタックに1足す 3*(3*i+1)*(3*i+1+1)
  * メインスタック(top と next)を掛ける 3*(3*i+1)*(3*i+2)
  3* メインスタックに3を掛ける 3*(3*i+1)*(3*i+2)
  su メインスタックをuレジスタに積む let (u,y)=
  li27* iレジスタと27を掛けてメインスタックに積む div(q*(27*i-12)+5*r)(5*t)
  12- メインスタックから12引く div(q*(27*i-12)+5*r)(5*t)
  lq* メインスタックにqレジスタを掛ける div(q*(27*i-12)+5*r)(5*t)
  5lr* 5とrレジスタを掛けメインスタックに積む div(q*(27*i-12)+5*r)(5*t)
  +   div(q*(27*i-12)+5*r)(5*t)
  lt5*   div(q*(27*i-12)+5*r)(5*t)
  /   div(q*(27*i-12)+5*r)(5*t)
  d 印刷用にメインスタックを複製
  48+ 30h(48d)の下駄を履かせて asciiコード変換
  P プリント y :
  sy メインスタックをyレジスタに積む let (u,y)=
  10lqlid2*1-*** 10とqとiをメインスタックに積んでiを複製して2を掛けて1を引いて掛けて掛けて掛ける g(10*q*i*(2*i-1),10*u*(q*(5*i-2)+r-y*t),t*u,i+1)
  10lulqli5*2-*lr+lylt*-** 10とuとqとiと5をメインスタックに積んで2を引いて掛けてrを足してyとtを掛けて引いて掛けて掛ける g(10*q*i*(2*i-1),10*u*(q*(5*i-2)+r-y*t),t*u,i+1)
  srsq メインスタックのtopをrレジスタに、nextをqレジスタに積む g(10*q*i*(2*i-1),10*u*(q*(5*i-2)+r-y*t),t*u,i+1)
  lult*st uレジスタとtレジスタを掛けてtレジスタに積む t*u
  li1+si iレジスタと1を足してiレジスタに積む i+1
  lmx mレジスタのマクロを無条件で実行 再帰マクロ呼び出し
]sm ここまでをマクロとしてmレジスタに積む マクロ定義
lmx mレジスタのマクロを無条件で実行 初回マクロ呼び出し
以上の文字列をechoしてdcに食わせるという寸法だ。
It's a great implementation.

 

— posted by nitobe at 05:01 pm   commentComment [0] 

T: Y: ALL: Online:
ThemeSwitch
  • Basic
Created in 0.0168 sec.
prev
2016.2
next
  1 2 3 4 5 6
7 8 9 10 11 12 13
14 15 16 17 18 19 20
21 22 23 24 25 26 27
28 29          
 
strawberry-linux geigercounter Ver.2
Sibasaki, Cyofu City, Tokyo, JAPAN
blogBar