Syn BASICのサンプルプログラム(4)

このページをスマホなどでご覧になる場合は、画面を横長にする方が読みやすくなります。
目次へ  前のページへ (1) (2) (3) (4) (5) 次のページへ
2022年11月24日 公開。以前はSyn BASIC(オンラインBASICインタプリタ)のページのコンテンツであったSyn BASICで作ったプログラムの例を、独立させて当ページにした。

前のページから引き続き、Syn BASICで作ったプログラムの例を紹介します。

2-6.円周率を求めるプログラム(1)

リスト11は、式(1)の公式を使って円周率を求めるプログラムです。

π=0141+x2dx ・・・ (1)

参考:式(1)が成立する理由については、後のコラムで説明します。

式(1)の定積分をシンプソン法で近似計算します。シンプソン法では、積分区間を等分割して近似計算を行いますが、リスト10のプログラムでは、0~1の区間を100分割しています。

リスト11、シンプソン法の定積分で円周率を求めるプログラムCOPY
10 'シンプソン法で円周率を求める
20 DIV=100 :'積分区間の分割数
30 SUM=0   :'区間積分値の合計
40 FOR I=0 TO DIV
50 X=I/DIV  :'区間のX
60 IF I=0 OR I=DIV THEN K=1 ELSE IF I MOD 2=1 THEN K=4 ELSE K=2  :'Kは重み(1か2か4)
70 SUM=SUM+K*(4/(1+X^2))  :'関数4/(1+X^2)を積分する
80 NEXT
90 PRINT SUM/(3*DIV) :'円周率を表示
100 END

リスト11のプログラムは、画面10のRUNボタン(青い三角形のボタン)をクリックすると、実際に実行できます。

円周率の正しい値は3.141592653589793238…ですが、プログラムを実行すると、表示される小数点以下11桁まで正確に円周率が求まっている事が分かります。

10 'シンプソン法で円周率を求める
20 DIV=100 :'積分区間の分割数
30 SUM=0 :'区間積分値の合計
40 FOR I=0 TO DIV
50 X=I/DIV :'区間のX
60 IF I=0 OR I=DIV THEN K=1 ELSE IF I MOD 2=1 THEN K=4 ELSE K=2 :'Kは重み(1か2か4)
70 SUM=SUM+K*(4/(1+X^2)) :'関数4/(1+X^2)を積分する
80 NEXT
90 PRINT SUM/(3*DIV) :'円周率を表示
100 END
CLS:LIST
画面10、シンプソン法の定積分で円周率を求めるプログラム
【コラム】式(1)の証明

式(1)の両辺を14倍して得られる式(2)を証明します。

π4=0111+x2dx ・・・ (2)

まず、式(3)の不定積分を計算してみましょう。

11+x2dx ・・・ (3)

式(4)の様にθをおいて、置換積分をします。

x=tanθ ・・・ (4)

式(4)の両辺をθで微分すると、式(5)を得ます。

dxdθ=1cos2θ
∴ dx=1cos2θdθ ・・・ (5)

式(3)に式(4)と式(5)を代入すると、式(6)を得ます。

11+x2dx =11+tan2θ1cos2θdθ =11+sin2θcos2θ1cos2θdθ =1cos2θ+sin2θdθ =dθ =θ+C ・・・ (6)

ただし、Cは積分定数

さらに、式(4)より式(7)が成立します。

θ=arctanx ・・・ (7)

式(7)を式(6)に代入すると、式(8)を得ます。

11+x2dx =arctanx+C ・・・ (8)

式(8)を利用すれば、式(2)は次の様に証明できます。

0111+x2dx =arctanx01 =arctan 1−arctan 0 =π4 ・・・ (9)

2-7.円周率を求めるプログラム(2)

リスト12は、円周率.jpというサイトに載っていた、C言語で記述された、Spigotアルゴリズムで円周率を求めるプログラムを、BASICに移植したものです。円周率を2400桁計算します。

Spigotアルゴリズムについては、下記のサイトに分かりやすい解説が載っています。

リスト12、Spigotアルゴリズムによって円周率を求めるプログラムCOPY
100 'SPIGOT.BAS
110 'Spigotアルゴリズムによる円周率の計算
120 BASE=10000               :'基底
130 N=8400                   :'計算項数
140 DIM NUMERATOR(8400)      :'分子
150 FOR I=0 to N-1
160   NUMERATOR(I)=BASE/5
170 NEXT
180 OUT=0                    :'出力値
190 FOR N=8400 TO 1 STEP -14
200   TEMP=0                 :'一時変数/繰り上がり
210   FOR I=N-1 TO 1 STEP -1
220     DENOM=2*I-1          :'分母
230     TEMP=TEMP*I+NUMERATOR(I)*BASE
240     NUMERATOR(I)=TEMP MOD DENOM
250     TEMP=INT(TEMP/DENOM)
260   NEXT
270   S$=RIGHT$("000"+STR$(OUT+INT(TEMP/BASE)),4)
280   IF N=8400 THEN S$=LEFT$(S$,1)+"."+RIGHT$(S$,3)
290   PRINT S$;
300   OUT=TEMP MOD BASE
310 NEXT
320 PRINT

リスト12のプログラムは、画面11でRUNボタン(青い三角形のボタン)をクリックすれば、実際に実行する事ができます。

100 'SPIGOT.BAS
110 'Spigotアルゴリズムによる円周率の計算
120 BASE=10000 :'基底
130 N=8400 :'計算項数
140 DIM NUMERATOR(8400) :'分子
150 FOR I=0 to N-1
160 NUMERATOR(I)=BASE/5
170 NEXT
180 OUT=0 :'出力値
190 FOR N=8400 TO 1 STEP -14
200 TEMP=0 :'一時変数/繰り上がり
210 FOR I=N-1 TO 1 STEP -1
220 DENOM=2*I-1 :'分母
230 TEMP=TEMP*I+NUMERATOR(I)*BASE
240 NUMERATOR(I)=TEMP MOD DENOM
250 TEMP=INT(TEMP/DENOM)
260 NEXT
270 S$=RIGHT$("000"+STR$(OUT+INT(TEMP/BASE)),4)
280 IF N=8400 THEN S$=LEFT$(S$,1)+"."+RIGHT$(S$,3)
290 PRINT S$;
300 OUT=TEMP MOD BASE
310 NEXT
320 PRINT
画面11、Spigotアルゴリズムによって円周率を求めるプログラム

次のページでは、ゲームのプログラムを紹介します。

目次へ  前のページへ (1) (2) (3) (4) (5) 次のページへ

関連ページ

Arduino 電子工作
このサイトの記事が本になりました。
書名:Arduino 電子工作
ISBN:978-4-7775-1941-5
工学社の書籍の内容の紹介ページ
本のカバーの写真か書名をクリックすると、Amazonの書籍購入ページに移動します。