HP電卓で遊ぶ ⇒ π（パイ）の計算

円周率（パイ　π）はいろいろな事で人々を楽しませていますね。
もちろん楽しむ為に存在しているのではなく学問や土木、工業などあらゆる方面で役に立っているわけですけど。
円周率の楽しみ方のなかで一番人気があるのは、いかに多くの数字を記憶するかということでしょうか。
日本人で10万桁も覚えたひとがいるそうですね。
私も中学生のとき友達と何桁まで覚えられるか競争したことがあります。
勉強へのやる気はなかなか起きないのに、こういう事にはやたら一生懸命になれるんですね。
そしてこのように遊びながら本気でやったことはいつまでも覚えていられるようです。
いま紙に書きながらテストしたら26桁出来ました。（3. 14159 26535 89793 23846 26433)
なのに　昨日言われたことはすぐ忘れてしまうのですよ。　これはなんなのでしょうか。
今までの人生でこの円周率覚えていたことで役立ったことないですね。
だからと言って　ゆとり世代の「円周率は３」というのはどうだかなーと思いますが。

もう一つ円周率の楽しみは、いかにして計算するかですね。
ネットで調べたら10種類くらいありそうです。
その中で　一番電卓向きと思われる　モンテカルロ法をHP電卓を使ってやってみます。
HP32S を使いました。　　最近はこればかりになっちゃいました。



モンテカルロ法をやるには　まず乱数の発生のばらつきが均一でなければ正しい測定はできないのでまず乱数を計ってみます。

全プログラムリスト 

A01  LBL A        **************************
A02  123.4567　　乱数発生器の初期値　何でもよい
A03  SEED　　　　乱数初期値セット
A04  CR VARS　　すべての変数をゼロにする

B01  LBL B        **************************
B02  RANDOM　　乱数の発生　0　~　1まで
B03  10　　　　　I 乱数を10倍して　０～１０までとする
B04  X                I
B05  1　　　　　 乱数に１を足して 0,1,2,  を　1,2,3,とする
B06  +
B07  STO i　　　乱数データが配列のポインターとしてセットされる（少数部分に何かあっても整数部分のみ自動的にポインターになる）
B08  1
B09  STO +(i)　　発生した乱数に対応する変数を一つ加算する　０～1の時はAが+1  1～2の時はBが+1   2～3の時はCが+1 
B10  GTO B


出力されたデータ　（30分間測定）
変数名　　データー範囲　　　出力回数
　A　　　　　0～1　　　　　3920
　B　　　　　1～2　　　　　3913
　C　　　　　2～3　　　　　4075
　D　　　　　3～4　　　　　3937
　E　　　　　4～5　　　　　3937
　F　　　　　5～6　　　　　3965
　G　　　　　6～7　　　　　3925
　H　　　　　7～8　　　　　4007
　I　　　　　8～9　　　　　 4017
　J　　　　　9～１0　　　　 3983

最大値　　4075
最小値　　3913
差　　　　162

生データの合計　　39677　
データの平均値　   39677 / 10 = 3967.7
差　/  平均値　　  162 / 3867.7 ＝　0.040829  
 乱数の発生頻度のバラツキ 　4.08%

ばらつき　4%   こんなものでしょうか。

最初測定１分間でデータの数150個でやりましたが　ばらつきが１６％にもなったので　30分間やり直しました。

これで終わりです。

乱数の発生のばらつきはほどほどに働いているので　次　モンテカルロ法でやってみます。
これは　x,y　二つの値を乱数でだして座標としそれが円の内部にあるか外部なのかで円周率を求めるものです。　

全プログラムリスト 


C01  LBL Ｃ        **************************
C02  123.4567　　乱数発生器の初期値　何でもよい
C03  SEED　　　　乱数初期値セット
C04  CR VARS　　すべての変数をゼロにする

D01  LBL D        **************************
D02  1
D03  STO +C　　全体の回数　カウントアップ
D04  RANDOM 　乱数作成　Ｘ
D05  X^2　　　　2乗してaにしまう
D06  STO A
D07  RANDOM 　乱数作成　Ｙ
D08  X^2　　　  2乗する
D09  RCL A 
D10  +　　　　　乱数のＸとＹを足す
D11  1　　　　　乱数を2乗して足したものが１を超えたか
D12  X ＜ Y ?
D13  GTO D　　 超えていたらデータ無視して次の乱数の計算に入る。
D14  1
D15  STO +D　 １以内ならばＤをカウントアップする
D16  GTO D　　 次の乱数の計算に戻る


測定は30分間やりました。
変数名　　データー
　A　　　　Ｘのランダム数
　B　　　　未使用
　C　　　　全体のポイントの回数　　　　30分ごのデータ　　23999
　D　　　　１以下のポイントの回数　　　30分ごのデータ　　18861

円周率　　（　18861  X   4   )    /   23999    =   3.14363098    　　すごい　　これ　嘘じゃないの？

誤差の計算　　（測定値　ー　　真の円周率　）/  真の円周率
誤差の計算　　（3.14363098  -  3.14159265) / 3.14159265  =  0.00064882

これで終わりです。

：
計算した円周率　3.14363098 　　
誤差　　　　　　0.064%

乱数のばらつきが４％あるのに最終誤差がどうしてこんなに減るのだろう。
モンテカルロ法　こんなので精度は出るはずないと思っていたのだけれど。
測定誤差か　計算ミスか　　これが正しければ　モンテカルロ法ってすごいんだ。　

HP32Sの乱数発生は品質いいみたい。　これソフトで同等の精度のものを作ろうとしたら、結構大変じゃないかな。
配列データの間接指定 ( i ) はうまく使うとコンパクトで解りやすいですね。これ　IF文　でやったら、ぐちゃぐちゃになりそう。
HP電卓でプログラムを組むのは　もうゲームですね。　役立つかどうかは二の次。

この次は　10種ほどあるという　円周率の計算　別のものでやって見ます。


追記
再現性があるか確認するために次の日（2017年11月１日）に乱数の種データを変えてモンテカルロ法（同じソフト）でやってみました。

乱数の種データ　　2.718281828
電卓のRUN時間　　30分　（HP32Sにて）
結果　　　　　　　C　＝　全体のポイントの回数　23734
　 D　＝１以下のポイントの回数　18638

円周率　　（　18638 X 4 ) / 23734 = 3.141147721

誤差の計算　　（測定値　ー　　真の円周率　）/ 真の円周率
誤差の計算　　（ 3.141147721 - 3.14159265) / 3.14159265 = 0.000141627

計算した円周率　3.141147721 　　
誤差　　　　　　0.01416%

再現性はいいですね。　しかし気分はデータ良すぎて納得できないな。


これ　次回ブログ作るときの参考の枠です。　関係ないけどここに置きます。

全プログラムリスト 

これで終わりです。