マチンの公式は
グレゴリー級数は
したがって
これを 
とおくと
漸化式は
で表され、
まで計算されたら、それに 
をかけ、
で割ると 
の値を得る。
も同様である。
もう、４０年も前のことになる。大学１年生のとき、数学の教員免許を取るためにはどうしても、コンピュータの講義を受けなければならなかった。１年を前期と後期に分け、前期だけの講義で２単位を取得する。私の通う大学にはコンピュータがなかったため、近くの専門学校のコンピュータをお借りしての授業となる。がしかし、借りられるのは１週間で、そのときに初めてコンピュータを目の当たりにすることとなる。そのため、それまでは、座学でコンピュータ言語の文法を学ぶ。コンピュータ言語の習得は実践が伴ってこそなのに、それはもう大変であった。言語はFORTRAN、その講義を担当してくださったのは猪野教授で、「変数は１文字とは限らず…」から始まって「A=A+1は式としてはおかしいだろう？この意味は…」最後に教授は「今わからなくても良いからな、そんなもんだと思っておけば良い。」何とも実の入らない講義を受けた。やがて、実習、大きな部屋に「鉄腕アトム」によく出てきたコンピュータ本体、動くとパネルがピカピカ光る。傍らにはカードリーダと呼ばれる機器、命令は１枚のシートに穴を開けたものを数枚重ね読み取らせる。シートに穴を開ける人をキーパンチャーと言うんだそうです。だから、こちらはコーディング用紙にプログラムを書き、キーパンチャーにお願いしてカードを作ってもらう。話はずれたが、もう一方の傍らには人の背丈もあるプリンター、１３５個(？)のドラムがガチャン、ガチャンと回転しながら、アルファベットを用紙に印字する。
猪野教授から、「以前この授業を受けた学生が、円周率の計算をさせ、コンピュータが２，３日止まったようになった」という話を聞いた。自分も挑戦してみたくなり、単位修得後もその授業を受けたが夢果たせずに卒業。
マチンの公式やグレゴリー級数を知っていたとしても、小数をどこまでも計算させるにはどうしたらよいか。5÷7＝0.714285714285…、有効数字が16桁のコンピュータでどうやって16桁以上を計算させるというのか。悩んだ。でも、人間が計算させるように計算させれば良いということにようやく気づいた。10倍して商を求めて貯める、余りを求め10倍して商を求めて貯める。これを繰り返せば良い。貯める場所は配列。初心者の私には、漸化式と無限小数のプログラムが１つにならないでいた。
１９８２年、日立の「ベーシックマスターJr」というパソコンを購入し、BASICを学んだ。FORTRANとよく似ているBASICは覚えやすかった。まもなく富士通の「FM-8」に乗り換え、その頃から円周率のコンピュータによる計算プログラムを作成していったが、何度となく失敗。ところが、１９８３年１２月３１日、NHK紅白歌合戦を見ていた最中に、ふと頭をよぎった。早速コンピュータに向かい、プログラムを手直しし、「RUN」、１００桁まで、時間は忘れてしまった。そんなにかかっていないと思う。正しい円周率の値を得た。その後、２０００桁に挑戦。夜にコンピュータを走らせ、朝方２０００桁までの値がプリンターで打ち出された。１３時間かかったと思う。
私のノートパソコンでは１００万桁を計算するのに６時間以上もかかっているが、「東京大学 金田研究室」の「スーパーπ」(フリーソフト)では１５秒で１０４万桁の計算を終える。何故半端な4万が付くんだ？速さも追求するためプログラムに工夫が施されているのでしょうね。
私が使った式は、マチンの公式で
マチンの公式は
グレゴリー級数は
したがって
これを とおくと
漸化式は
で表され、 まで計算されたら、それに をかけ、 で割ると の値を得る。 も同様である。
当時のプログラム(F-BASIC)がまだ残っている。幼稚なプログラムだな～。
10 '***************************************
15 '*
20 '* File name : PI
30 '*
40 '* 1985/02/14
50 '*
60 '* Programed by F-NISHINO
70 '*
80 '***************************************
100 DEFDBL A-Z
110 DEFINT C,F,I,J
120 INPUT "何桁まで求めますか。";FIG
130 TIME$="00:00:00"
140 N=10000000000#
150 COUNT=INT(FIG/10)+1
160 DIM AN(COUNT+1),SN(COUNT+1),TN(COUNT+1),PI(COUNT+1)
170 'SN1
180 CNST=5
190 GOSUB 1000
200 AN(1)=2000000000#
210 SN(1)=AN(1)
220 'SN
230 FOR I=2 TO ITEM
240 FOR J=1 TO COUNT+1
250 AN(J)=AN(J)*(2*I-3)
260 NEXT J
270 DVS=25*(2*I-1)
280 GOSUB 2000
290 FOR J=1 TO COUNT+1
300 IF I MOD 2=0 THEN SN(J)=SN(J)-AN(J) ELSE SN(J)=SN(J)+AN(J)
310 NEXT J
320 NEXT I
330 'TN1
340 CNST=239
350 GOSUB 1000
360 BN=1
370 FOR I=1 TO COUNT+1
380 AN(I)=INT(N*BN/239)
390 BN=N*BN-239*AN(I)
400 TN(I)=-AN(I)
410 NEXT I
420 'TN
430 FOR I=2 TO ITEM
440 FOR J=1 TO COUNT+1
450 AN(J)=AN(J)*(2*I-3)
460 NEXT J
470 DVS=239*239
480 GOSUB 2000
490 DVS=2*I-1
500 GOSUB 2000
510 FOR J=1 TO COUNT+1
520 IF I MOD 2=0 THEN TN(J)=TN(J)+AN(J) ELSE TN(J)=TN(J)-AN(J)
530 NEXT J
540 NEXT I
550 FOR I=1 TO COUNT+1
560 PI(I)=16*SN(I)+4*TN(I)
570 NEXT I
580 FOR I=COUNT+1 TO 1 STEP -1
590 QUO=INT(PI(I)/N)
600 PI(I-1)=PI(I-1)+QUO
610 PI(I)=PI(I)-QUO*N
620 IF PI(I)<0 THEN PI(I)=PI(I)+N:PI(I-1)=PI(I-1)-1
630 NEXT I
640 T$=TIME$
650 'PRINT OUT
660 '
670 LPRINT RIGHT$(STR$(PI(0)),1)+".";
680 FOR I=1 TO COUNT-1
690 PI$=STR$(PI(I))
700 LPRINT RIGHT$("000000000"+RIGHT$(PI$,LEN(PI$)-1),10);
710 IF I MOD 5=0 THEN LPRINT
720 LPRINT " ";
730 NEXT I
740 LPRINT
750 LPRINT TAB(20);"円周率 ";FIG;" 桁まで 計算時間 ";T$
760 END
1000 'Item caluculator
1010 ITEM=1
1020 WHILE (FIG+10)*LOG(10)>=LOG(2*ITEM-1)+(2*ITEM-1)*LOG(CNST)
1030 ITEM=ITEM+1
1040 WEND
1050 RETURN
2000 ' N term caluculator
2010 FOR J=1 TO COUNT
2020 QUO=INT(AN(J)/DVS)
2030 AN(J+1)=AN(J+1)+N*(AN(J)-DVS*QUO)
2040 AN(J)=QUO
2050 NEXT J
2060 AN(COUNT+1)=INT(AN(COUNT+1)/N)