人気ブログランキング |

ブログトップ

電子工作やってみたよ

円周率の計算式 少し解った気がする その2

c0335218_14110076.jpg

c0335218_14110801.jpg

HP電卓のRPN(逆ポーランド記法)で書かれた円周率計算プログラム、動作は10000桁まで正しく行う事を確認してあるのですが、原理と仕組み、そして動作、いまだ理解できません。
このプログラムのオリジナルは
Many Digits of Pi by Katie Wasserman - MoHPC
www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/articles.cgi?read=899
という方が書かれておりHPの各機種(HP32sii,HP-16c,HP-12c,HP-19c,HP-15c)向けにつくられています。
私はこの中のHP-15cのものをDM42(HP42S)用に10000桁まで拡張してみました。

もうかなり考えているのですが(一か月以上かな)脳みそがまったく収束しません。


「e-Gadget-プログラム関数電卓」ブログで、やす様が「多桁円周率の計算(1)」というのを発表されています。
オイラーによる arctan(X) の展開式を使ってカシオの電卓fx-5800P上のカシオBASICで円周率を計算しているのですが、さいわいHP電卓のプログラムも同じ展開式でやっています。
ゆえに このソフトを解読すれば、HP電卓のソフトも解るのではないかと思い挑戦してみました。

ソフトの動作を理解するには自分で動かしてみなくちゃ解りません。
ソフトを動かすにはハードがなくちゃ動きません。
大金をはたいてfx-5800Pを買ってしまいました。
まあ 電卓プログラムはHPのものしか知らなかったので他のプログラムを知るチャンスとも思ったのですけれど。

四苦八苦しながらプログラムを動かしてみたら やす様と同じ結果が画面に出てきました。
3.141592......   バンザイ!

しかし数値の動きを追いかけたり、改造したり、デバッグしたりは慣れた言語、慣れたハードでないとなかなか思い通りに行きません。
そこでこのカシオBASICをN88ベーシックに書き直してやってみました。
そうしたらいつものことですが、やってみてから気付きました。 このN88有効桁数が7桁しかない。
そこですべての変数を倍精度の16桁に変えたのですが、for~nextの所でエラーが出る。
でN88は諦めて、次の手持ちのソフト、F-BASIC97に乗せ換えて倍精度(有効15桁)でやりました。
これはうまく動作して、カシオBASICと同じ 3.141592...がきれいに出てきました。

いま現在の進行状況は ここまでなんですけれど。
いつものごとく「もうちょいだ」と思ってからが長いんですよね。


概略のイメージが解ってきたので、これからまた HP電卓のソフトに挑戦です。

やす様の許可をもらって CasioBasicのソフトを載せます。
実験される方はやす様のブログを見て行ってください。
「e-Gadget-プログラム関数電卓」の「多桁円周率の計算(1)」



fx-5800P CasioBasic
ファイル名: PICALC
0→DimZ
"DIGITS"?→X
"BASE"?→K
Int(X÷K)+1→A
10^(K)→D
X+1→DimZ
X÷log(2)→N
"EXE TO GO"◢

Cls
For Int(N)→F To 1 Step -1
Locate 2,1," "  // スペース2つ
Locate 1,1,F

0→C
For A→I To 1 Step -1
Z[I]F+C→W
Int(W÷D)→C
W-CD→Z[I]
Next

2F+1→G
0→C
For 1→I To A
CD+Z[I]→T
Int(T÷G)→W
T-WG→C
W→Z[I]
Next

Z[1]+2→Z[1]
Next

For 1→I To A
Z[I]◢
Next

0→DimZ
"END"

__________________

以下は F-BASIC97 に書き換えたものです。
変数はすべて倍精度の15桁にしてあります。
動作はOKでした。
90桁の答えが出ます。


'save "PAI015.bas"
dim Z#(20)
for I=0 to 20 : Z#(I)=0 : next
'K#                                      '  計算ループ最大値 固定
for K#=350 to 1 step-1         ' 円周率100dig出すために3.5倍のKから計算する、
' 理由は解らない? 経験則
C#=0 'キャリー最初はゼロ
for I#=11 to 1 step -1 '下位のレジスタ-から上位に向け計算 レジスター10個
W#=Z#(I#)*K#+C# '現在対象レジスタにKをかけ前回のキャリーを加算する
C#=int(W#/(10^10)) '10桁より上をキャリー検出する 1レジスターは10桁
Z#(I#)=W#-(C#*(10^10)) '10桁に抑えるためキャリーを差し引現在値にする
next

G#=(2*K#)+1 '基本式の分母側の計算
C#=0 'キャリー

for I#=1 to 11 ' 上位のレジスタから下位に向けて計算 レジスタ10個分
T#=(C#*(10^10))+Z#(I#) 'キャリーがあったら現在のレジスタの頭に加算
W#=int(T#/G#) '基本式の分母で割り算して整数部分を取り出す。
C#=T#-(W#*G#)
Z#(I#)=W#
next

Z#(1)=Z#(1)+2

?"Z(1_2_3_4_5 )=";Z#(1);" ";Z#(2);" ";Z#(3);" ";Z#(4);" ";Z#(5)
?"Z(6_7_8_9_10)=";Z#(6);" ";Z#(7);" ";Z#(8);" ";Z#(9);" ";Z#(10)
stop

気になったこと
基本公式 2+(K/2K+1)*(1回前の値) で変数の有効桁数が15桁であっても Kの値が3桁であると1レジスタ10桁*3桁で13桁必要となります。
ましてHP電卓のソフトでは10000桁計算するためKの最大値を35000として掛けていますが、これは5桁なので10桁にかけると15桁必要となりHP42Sの規格一杯になってしまいます。
実動作させたのはDM42およびFree42であったのでたまたまうまく行ったのかもしれません。
私の考え違いかも知れませんが。

これで終わりです。




追記
以下はレジスタを分割せずに一つのままの単純な形で計算したものです。

変数レジスターを分割せずに一つとして扱ったときのF-BASIC97でのプログラムです。
1変数の有効桁数が大きければ、この単純プログラムで円周率の桁数が増やせるのですね。

'save "PAI001.bas"
BF#=1 'このBF#の初期値は1でも0でも結果は変わらない。
for K#=55 to 1 step -1 '1変数の有効桁15とし3.5倍+余裕でK=55から出発
ANS#=2+(K#/((2*K#)+1))*BF#
BF#=ANS#
?"K=";K#;" ANS#=";ANS#
next
stop

これで終わりです。



下は上記プログラムの数学的な式、
c0335218_13040846.jpg

下の写真は1 レジスタのプログラムを実行した結果です。
Kの数値が減少するにしたがって円周率の値が3.141に近づいていくのがわかります。
不思議なのは、Kが大きいときはなかなか近づかなくて最後のkが2から1になるときに、ピタと正解になることですね。
これ Kを1から増加させる方向で計算させるプログラム出来るのでしょうか。
そうすれば、答えが限りなく円周率に近づいていく様子が見れると思うのですが。
c0335218_14111779.jpg



我が家の軒先の鉢に咲いた「あさがお」
去年からそのままにしてありました。
昔からやはり夏は「あさがお」ですね

c0335218_14112809.jpg

c0335218_14113726.jpg

c0335218_14114877.jpg



電子工作ランキング

にほんブログ村 その他趣味ブログ 電子工作へ
にほんブログ村

Commented by すがわら at 2019-08-18 10:15 x
古い人間なのでFORループを浮動小数で制御するのは怖いです。ただK値が35000となるとこれは単精度整数では収まらないので2重ループにするといった面倒な処理が必要になりますが...
それにしてもFX-5800ってBASIC風なんですね。FX-520Pから続いていたキーボタンだけで組んでいく簡易プログラム風も味が有って好きなんですけどね。
Commented by telmic-gunma at 2019-08-18 11:51
> すがわらさん
すごいね。この暑さ 生きていますか
そちら 海抜低くて町の中だから、暑さもすごいでしょう。
こちら 海抜600m位で山の中だけど猛烈な暑さです。
クーラーの効いた部屋にコモっています。

「FORループを浮動小数で制御するのは怖いです」
いままで 何にも気にしないときは平気でしたけれど、それ聞いた後では気になっちゃいますね。

「キーボタンだけで組んでいく簡易プログラム風も味が有って好き」
これ プログラム電卓の本質かも
今回 fx-5800Pのプログラムを組んでいて何か違和感を感じたのは、コマンドをメニューから選ぶということからきているのかも。
こういこと感じるのも、古い人間ということ?
HP電卓もコマンドが増えたせいか、HP32S以後そうなっちゃったものね。
むかしの表示が7セグメントだけの物は、コマンドがキートップやパネルに印刷されており
プログラムが楽だったものね。
今までモヤモヤしていたものが、少し晴れてきたような気がします。
ありがとうございます。
7セグメントの手持ちの機種というと
HP19C, HP97, HP15C かな  
もう一度みなおして、可愛がってあげようかな。
Commented by Kusaka at 2019-08-18 12:47 x
こちらは海抜1000mで室温28℃ぐらいなので扇風機とビールで頑張っています。 大金をはたいてとあったので、FX-5800Pを検索したら7000円ぐらいですね。HP電卓と比べたらめちゃくちゃ安いじゃないですか。 検索結果にCasioのFX-CG50っていうものも出てきて、Pythonでプログラミング出来るとか。こっちも面白そうな。まあ、円周率への興味が終わったらCasio電卓は不要と思いますので社食の昼飯代¥310で買い取りしますよ。
Commented by すがわら at 2019-08-18 15:51 x
まず訂正... FX-520PでなくてFX-502Pでした。私が最初にいじったプログラム電卓です。
当地は平地というほど標高は低くないですけど(約150mくらい)、気温は上がりますね。日が照れば当たり前のように36℃に達します。東京都内における最高温度記録も近所のアメダスポイントで記録してますし(40.8℃)。結構窓全開でお過ごしのお宅を見かけるのですが、いったい何かと戦っているのか、はたまた何かの修行なのか謎です。
Commented by telmic-gunma at 2019-08-18 18:28
> Kusakaさんすがわらさん
今週一週間山梨にいます いま着いたところ
ノートパソコンの外部電源忘れた
もうすぐ 電池切れると思う
今週はブログの見るのも書くのも出来ないと思います
アポロ13号の気分だね
さようなら
Commented by あけび at 2019-08-23 10:33 x
オシロイバナはもうちょいで開花ですね~~
朝顔といえば夏!! 花びらを見るだけで猛暑でも何気に涼しく感じられるから不思議です!!
我が家に朝顔はありませんが これを見て 朝顔でグリーンカーテンを作りたくなりました。
ゴーヤーで作っているお宅がありますね。
Commented by telmic-gunma at 2019-08-25 08:46
> あけびさん
ありがとうございます。
おしろい花 出張から帰ってきてみたらちょっと咲いていました。
日当たりの良いところのものは、たくさん咲いています。
次のブログに載せますね。
朝顔でグリーンカーテンなんて綺麗でしょうね。
来年 挑戦してみようかな。
名前
URL
画像認証
削除用パスワード

※このブログはコメント承認制を適用しています。ブログの持ち主が承認するまでコメントは表示されません。

by telmic-gunma | 2019-08-17 16:16 | HP電卓 | Comments(7)