人気ブログランキング |

ブログトップ

電子工作やってみたよ


c0335218_13030539.jpg
円周率の計算について 「e-Gadget-プログラム関数電卓」ブログで、やす様が「多桁円周率の計算(1)」というのを発表されています。
これは カシオのfx-5800上のカシオBASICで書かれていますが、実験しやすいように F-BASIC97 に書き換えて動作を調べてみました。

一言で言ってしまえば、 
オイラーによる arctan(X) の展開式をそのまま使って電卓が計算できる桁数(今回は10桁)を使って多倍長計算をさせて結果を出しています。
結果は以下の物です。


注) ミスあり ”Kを3500から計算”とあるのは、 ”Kを35000から計算”が正しいです。

c0335218_20204731.jpg



c0335218_20205809.jpg



プログラムから元の式や構造、動作、など調べるのは大変。

「RPN(逆ポーランド記法)で書かれていてもボトムアップでコツコツ調べて行けばいつかは全体が見えてくるはずだ」と思っていたのですが、今回かなり難しい(以前 石取りゲームを解読した時もそうだったような)
次にやす様のCASIO BASICでやったけれど解らず、F-BASIC97に乗せ換えて何度も数値の流れを見てやっと解ってきました。
やはり 物事はトップダウンでやるべきなのでしょうね。
全体が解ってないのにすみっこからチマチマやっても進みません。
ゲームとして楽しむためならそれでもいいのかな。
でも 最初は全体なんて解っていないのだから、どうどうめぐりですね。

プログラムは短いのに、繰り返しの回数が多くとにかく全体の計算量はべら棒に多いのです。
こういうのは、動作を調べるのは難しいですね。
かえって、プログラムは長くても一回しか通過しないやつならば、流れがすぐ解る気がします。


以下は 実際のプログラムの繰り返し数値計算しているところ調べたものの一部です。
1レジスタで計算するのを3桁にしてトータル10桁までの計算の流れを調べました。
ここは多倍長計算の繰り返しされるところこの数値の流れをA4の紙で20枚以上書いておぼろげながらわかってきました。

c0335218_21043692.jpg
これはまだ BASICでの計算手順が見えただけです。
本式のHP電卓のRPNでの動作解析はまだこれから、また挑戦です。

解っている人、あるいは解る人にとっては、何と言う事無いでしょうが、私のように解らない人にとって、この説明でわかるでしょうか。
たぶん半年後に私がこのブログを見ても解らなくなっている気がしてなりません。

他人に(明日の自分でもある)プログラムを解ってもらうには、日本語で普通の文章で「ここはこれこれのことをしていますよ。」 とていねいに書いてあるのが良いのかな。
あと「図」ですね。
「図」と「説明文」 まずこれさえあれば、プログラムは後でいいのでしょうね。

なに言ってるんだか。 
昔のわたしは、説明など何も書かず使い捨てソフトばかり書きなぐっていたような気がしますが。


わが村の田んぼアート 
今が一番 緑が濃いのかな。 
これからだんだん黄色くなるのでしょうね。
c0335218_16233875.jpg


電子工作ランキング

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

by telmic-gunma | 2019-08-31 04:38 | HP電卓 | Comments(0)

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



電子工作ランキング

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

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


c0335218_13030539.jpg
円周率の計算する以下の式、最後は延々と繋がっており、どこからどのようにして計算を開始すればいいのか、見当がつきませんでした。
そんな時ある本(参照1)を見ていたら式の最終としたいところの(2+ ・・・)を1とおいてやるように書いてありました。
早速 その手順でやると円周率らしき値が出てきました。
Kの値の大きいところから計算を開始することにより有効桁がどんどん増えていきます。
今回はHP42Sを使ってこの計算式を普通に計算するだけなので得られる答えは12桁しかありませんが、Kの値により有効桁が増えていくのがわかりました。


解っている人にとっては何ともない事でしょうが、解らない人間にとってはこんな事でも進めなくなってしまいます。
c0335218_13040846.jpg
今回 ノートに書いたものをカメラで撮ったのでちょっと見にくい感じです。
後でスキャナーで撮り直してみます。



スタートする  円周率の答え  誤差(%)  有効桁数
Kの値


10 3.14084209564 0.023896 2桁
20 3.14159211321 0.000017201 6桁
30 3.14159265315 0.000000014 9桁
35 3.14159265358 0.00000000001 10桁
40 3.14159265359 0.00000000000 11桁
50 3.14159265359 0.00000000000 11桁

10桁の答えを出すのにKは35必要でした。
故に10000桁の答えを出すとき35000を設定したのは比例関係が成り立つならばちょうど良かったということになりますね。
理屈は解らなかったのですけれどね。

これで終わりです。

c0335218_13044686.jpg
参照1   (数学問題へのコンピュータアプローチ 昭和51年5月発行)

いま我が家の周りに咲いている花
c0335218_15030304.jpg

いま我が家の周りに咲いている花

c0335218_15025808.jpg

いま我が家の周りに咲いている花
c0335218_15031417.jpg

いま我が家の周りに咲いている花
この花だけは、春から秋まで、いつでも、どんな時でも綺麗に精一杯咲いています。 不思議です。
見習いたい。
c0335218_15025358.jpg






電子工作ランキング

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

by telmic-gunma | 2019-08-06 15:17 | HP電卓 | Comments(13)