エピ１７ 円周率を計算しましょう！ ＢＡＳＩＣでね
＊＊＊
エピ１７ 円周率を計算しましょう！ ＢＡＳＩＣでね
＊＊＊
ＳＰ５０３０なら "PRINT π" で、
3.1415927
と表示されます。ですが、この値は 1624h 番地にある C2 A1 DA 0F C9 という定数なのです。
それではなくて、ですね。ちゃんと円周率を計算してみたいのです。ＢＡＳＩＣでね。本日はお日柄もよくて。円周率の日（2026.3.14）ですし。
＊＊＊
円周率の計算は、マチンの公式かその類似の式で求めるのがポピュラーらしいです。ＥＮＩＡＣもマチンの公式で２０３７桁を計算しています（１９４９年）。
でもこれは昔のことで、現在の円周率の記録は別の方式で計算したものですって。
とはいえ。手軽に計算できることから数千桁くらいの小さな計算（！）をするなら、今でもマチンの公式を使うことがあるそうです。大きな計算ってば、５兆桁とか５０兆桁、...。あ、はい。小さな計算がしたいので、それでいいです。
＊
マチンの公式とは。ググれば沢山の解説が見つかりますけど。簡単に説明します。
三角関数の定義から、
tan(45度) = tan(π/4) = 1
ですから、tan() の逆関数を atn() とすると、....難しいことはおいておいて、
π/4 = atn(1)
つまり、
π = 4*atn(1)
です。
atn() はテイラ級数（あるいはグレゴリ級数）として展開して計算できます。atn() をテイラ展開すると、
atn(x) = x^1/1 - x^3/3 + x^5/5 - x^7/7 ...
です。x=1 を代入して沢山の項の合計を求めれば円周率になります（４倍してね）。ただし、そのままでは収束が遅くて効率が悪いのです。分母は奇数つまり 2*n+1 で、分子は ±1 ですから、n=10 のときは -1/21 = -0.0476 なのです。3.14 の 4 すら怪しいのです。
で、...そこを工夫をしたのがマチンの公式です。
atn(1) = 4*atn(1/5) - atn(1/239)
これで x=1 の代りに x=1/5 や x=1/239 を代入することになって、... -x^3/3, +x^5/5, -x^7/7, ... なのですから、その効果は分かりますね。x=1/5 を代入して n=10 なら、-(1/5)^21/21 = -0.0000000000000000998 です。とても小さいネ。
これで円周率を実用的に計算することができるのです。
＊
ＳＰ５０３０の関数 ATN() を使って 4*ATN(1) を求めたりすると、3.1416231 となって悲惨です。マチンの公式を使えば、3.1415927 となります。
＊
ですが、足したり引いたりするのがちょっとな、...と感じて、ここではオイレルの公式を使います。テイラ級数もね。加算だけで計算できるように変形します。
これは感性、...音楽性とかのアレ、...の違いなので異論は認めます。普通に書いてマチンの方が速いはず。
＊＊＊
はい。プログラムを書きました。
10 X=0:S=0:T=10:A=0:B=0:D=25:E=3125:Y=0:Z=0:W=50
11 DIM S(W),A(W),B(W)
12 A(0)=2:A(1)=8:B(0)=0:B(1)=3:B(2)=0:B(3)=3:B(4)=3:B(5)=6
20 FOR N=0 TO W
31 Y=Y+1:Z=Z+9:D=D+50:E=E+6250
41 S=0:A=0:B=0
42 FOR X=W TO N STEP -1
43 S(X)=S+S(X)+A(X)+B(X):S=INT(S(X)/T):S(X)=S(X)-S*T
44 A(X)=A+A(X)*Y:A=INT(A(X)/T):A(X)=A(X)-A*T
45 B(X)=B+B(X)*Z:B=INT(B(X)/T):B(X)=B(X)-B*T
46 NEXT
51 A=0:B=0
52 FOR X=N TO W
53 S=A*T+A(X):A(X)=INT(S/D):A=S-A(X)*D
54 S=B*T+B(X):B(X)=INT(S/E):B=S-B(X)*E
55 NEXT
60 NEXT
70 PRINT S(0);".";:FOR X=1 TO W:PRINT STR$(S(X));:NEXT:PRINT
＊
でも、ＢＡＳＩＣだとアルゴリズムが分かり辛いですね。スクリプトで普通に書くと、...
w = 50 -6; s = 0; a = 280000*(10**w); b = 30336*(10**w)
for n in range(1,w):
s = s + a+b
a = a * (2*n)/(2*n+1)*1/50
b = b * (2*n)/(2*n+1)*9/6250
print s
> 31415926535897932384626433832795028841971693993730
となります。変数 a と b が級数の項で、値を更新しながら変数 s に加算しているのです。スクリプトは多倍長の演算を実装しているので簡潔に書けますね。
変数 a の更新は (2*n)/(2*n+1)*1/50 を掛けていて、n が大きくなると、
2n 1
―― ・ ―― ⇒ 0.01999...
2n+1 50
となりますから、１回の更新で１桁以上落ちます。変数 b は 0.00143999... ですから、もっと速いです。
それで w 回ループすれば w 桁は求まるでしょう、...というザックリした処理です。
なお、正しい円周率は、
3.14159 26535 89793 23846 26433 83279 50288 41971 69399 37510
58209 74944 59230 78164 06286 20899 86280 34825 34211 70679
.....
ですから、最後の２桁は違います。w=52 とすれば 50 桁までを正しく計算できますけど。
計算時間は、w=1000 ならほぼ瞬時で、w=10000 だと少し待つかな、数秒くらい？
＊＊＊
ＢＡＳＩＣのプログラムも、LOAD して RUN すれば、
3.1415926535897932384626433832795028841971693993749
と表示します。最後の２桁は違います。ステップ数は 44,997,784 です（LOAD と RUN と BYE を含む）。
計算時間は、...。後で評価しますね。
＊
ＢＡＳＩＣのプログラムについても少し触れます。
10 X=0:S=0:T=10:A=0:B=0:D=25:E=3125:Y=0:Z=0:W=50
文番号 10 で変数 T の使用を宣言して数値 10 を設定していますが、その後に変更していません。
定数なら変数にする必要はないのですが、あえて処理速度のために変数にしています（...、あ。マジックナンバをソースコードに埋め込むな、とするコーディング規約があることは知ってますよ。ソレとは違う話です）。
数値 10 のままだと、数式の処理の際に 31h, 30h と調べて、これを５バイトの浮動小数点数に変換して数値ワークに転送することになります。数値変数なら変数領域を検索して５バイトの数値（の番地）を得るだけなので、こちらの方が速いのです。
数値 10 は、文番号 43, 44, 45 と 53, 54 で計８回登場して、更に FOR ループの中にあるので、変数にしました。
その効果は（W=10 の場合で）
数値の場合： 2,562,246 ステップ
変数の場合： 2,331,510 ステップ
ステップ数での見積もりですが、2331510/2562246 = 0.9099 で約１割削減です。
＊
余談になりますが。
「地底最大の作戦」（Ｉ／Ｏ、１９８０年７月号）では、ＢＡＳＩＣのスピード性として、
＊ マルチ・ステートメントの乱用
＊ 数字の定数化
＊ REM 文の排除
＊ 条件文の順序に対する考慮
の４点を挙げています。
数値 10 を変数 T にするとは「数字の定数化」のことに他なりません。他の３点もスピードを追求するなら、その通りですね。
マルチ・ステートメントにすれば別々の文にするより速いですし、文の数が減ると GOTO, GOSUB の処理も速くなります。
プログラムリストを見ると「FORX=0TO799:」や「IFY=109THEN450」などと、スペースを削除していますし。当然に "LET" はなくて。"NEXT" のループ変数名は省略です。変数名は１文字で、時々２文字です。...、これらは敢て言うまでもない常識的なこと、...なのでしょう。
でも、あと１つ追加したら良い点があると思います。それは、...
＊
10 X=0:S=0:T=10:A=0:B=0:D=25:E=3125:Y=0:Z=0:W=10
文番号 10 では、変数 X や S の使用を宣言して 0 を設定しています。これらの変数を宣言することは、無意味な、あるいは宗教的な信仰だけではなくて、実は処理速度にも関係することなのです。
試しに、
10 T=10:D=25:E=3125:W=10
として、X や S などの変数の宣言を削除すると、プログラムは正しく動作しますが、
宣言あり： 2,331,510 ステップ
宣言なし： 2,518,363 ステップ
...、なんと、2518363/2331510 = 1.0801 と８％も増えてしまいます。
ＳＰ５０３０は変数を参照するとき、変数領域の先頭から順番に検索するのです。変数領域の変数の順番は単に変数の宣言順です。だから使用頻度の高い変数を先に宣言すると、変数の参照が効率的になるのです。
ループ変数やループ内で参照する変数は先に宣言すると良いでしょう。
効果はプログラム次第ですが、８％なら小さくないよ（消費税相当だよ！）
＊
ＳＰ５０３０の数値変数は仮数が４バイトありますから、１桁（０〜９）だけに使うのはもったいないのです。
変数 T を T=100 として、
12 A(0)=2:A(1)=80:B(0)=0:B(1)=30:B(2)=33:B(3)=60
20 FOR N=0 TO W STEP 0.5
とすれば２桁（０〜９９）毎に処理できますし、変数 T を T=10000 として、
12 A(0)=2:A(1)=8000:B(0)=0:B(1)=3033:B(2)=6000
20 FOR N=0 TO W STEP 0.25
とすれば４桁（０〜９９９９）毎です。
20桁）T= 10, W=20： 7,934,095 ステップ
20桁）T= 100, W=10： 4,249,935 ステップ
20桁）T=10000, W= 5： 2,400,440 ステップ
圧倒的にイイですね！ 20 桁の計算が 2400440/7934095 = 0.3025 と７割オフです。営業終了時間が近づいた頃のお刺身みたいです。
それからループの終了判定を追加します。文番号 53 の末尾で、変数 S を変数 R に加算します。文番号 56 で判定して R=0 なら終了です。文番号 10 での変数 R の宣言は W よりは先ですかね。
51 A=0:B=0:R=0
53 S=A*T+A(X):A(X)=INT(S/D):A=S-A(X)*D:R=R+S
56 IF R=0 THEN 70
これを試すと、
20桁）T= 10, W=20： 7,934,095 --> 6,868,540 ステップ
20桁）T= 100, W=10： 4,249,935 --> 3,768,198 ステップ
20桁）T=10000, W= 5： 2,400,440 --> 1,777,719 ステップ
となって、1777719/2400440 = 0.7405 と更に２５％オフ！ 賞味期限が間際の...。
＊
アルゴリズムの見直しによる計算量の改善は、コーディング技法による効率化とは異次元です。
ループ変数が５バイトの浮動小数点数なんてこともあるけどさ。ループを回さない、ループ回数を減らすことができるのはアルゴリズムです。
「地底最大の作戦」には、スピード性の追求ではなくて、...それも大切だけど、...ゲームとしてのアルゴリズムというか、デザイン？、アイデア？、...何と言うべきか分からないけど、そんな感じの何かがあると思う。
＊＊＊
＊＊
＊
エピ５で志した円周率を表示するプログラム、...これが動く雰囲気にたっぷりと浸れました。
LOAD したプログラムが RUN で動き始めて、FOR ループが回って、HL が命令文を渡って中間コードで分岐して、LET で数値変数を参照して、算術演算子が加算や乗算をして結果が数値ワークにスタックされて、領域管理が忙しなく動き、...そうして、級数が収束して円周率が求まっていく、...。ああ、楽しいねー
＊
＊＊
＊＊＊
間違いの指摘とか疑問とか、ご意見・ご感想とかありましたら、どうぞ感想欄に！
＊＊＊
2026.3.15 加筆と推敲
2026.3.19 微推敲
2026.3.20 推敲
2026.4.3 推敲
2026.5.3 推測
2026.5.4 微推敲
