pdp11のエミュレータを探してググっていたら、とんでもないものをみつけた。
「
Pi with dc
」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 Pi 」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.pdf
// 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/pi
// 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 |
１をｑレジスタに積む |
g(1,180,60,2) |
| 180sr |
１８０をｒレジスタに積む |
g(1,180,60,2) |
| 60st |
６０をｔレジスタに積む |
g(1,180,60,2) |
| 2si |
２をｉレジスタに積む |
g(1,180,60,2) |
| [ |
ここから |
|
| 3li* |
3 とｉレジスタを掛けてメインスタックに積む |
3*(3*i+1)*(3*i+2) |
| 1+ |
メインスタックに１を足す |
3*(3*i+1)*(3*i+2) |
| d |
メインスタックを複製して積む |
3*(3*i+1)*(3*i+1+1) |
| 1+ |
メインスタックに１足す |
3*(3*i+1)*(3*i+1+1) |
| * |
メインスタック（top と next）を掛ける |
3*(3*i+1)*(3*i+2) |
| 3* |
メインスタックに３を掛ける |
3*(3*i+1)*(3*i+2) |
| su |
メインスタックをｕレジスタに積む |
let (u,y)= |
| li27* |
ｉレジスタと２７を掛けてメインスタックに積む |
div(q*(27*i-12)+5*r)(5*t) |
| 12- |
メインスタックから１２引く |
div(q*(27*i-12)+5*r)(5*t) |
| lq* |
メインスタックにｑレジスタを掛ける |
div(q*(27*i-12)+5*r)(5*t) |
| 5lr* |
５とｒレジスタを掛けメインスタックに積む |
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+ |
３０ｈ（４８ｄ）の下駄を履かせて |
asciiコード変換 |
| P |
プリント |
y : |
| sy |
メインスタックをｙレジスタに積む |
let (u,y)= |
| 10lqlid2*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*-** |
１０とｕとｑとｉと５をメインスタックに積んで２を引いて掛けてｒを足してｙとｔを掛けて引いて掛けて掛ける |
g(10*q*i*(2*i-1),10*u*(q*(5*i-2)+r-y*t),t*u,i+1) |
| srsq |
メインスタックのtopをｒレジスタに、nextをｑレジスタに積む |
g(10*q*i*(2*i-1),10*u*(q*(5*i-2)+r-y*t),t*u,i+1)
|
| lult*st |
ｕレジスタとｔレジスタを掛けてｔレジスタに積む |
t*u |
| li1+si |
ｉレジスタと１を足してｉレジスタに積む |
i+1 |
| lmx |
ｍレジスタのマクロを無条件で実行 |
再帰マクロ呼び出し |
| ]sm |
ここまでをマクロとしてｍレジスタに積む |
マクロ定義 |
| lmx |
ｍレジスタのマクロを無条件で実行 |
初回マクロ呼び出し |
以上の文字列をechoしてdcに食わせるという寸法だ。
It's a great implementation.
pi
RSS2.0
Comments