人気ブログランキング |

ブログトップ

電子工作やってみたよ

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)

円周率 DM42で10000桁まで計算出来たのに、そのプログラムの動作がまだ解りません。
プログラム自体はこんなに短くて10000桁でも正しく計算出来るなんてすばらしいのに。
いつまでたっても解らないとそのうち気力が消えてしまいそう。

*最初の方法は、スタックや変数の内容を調べながらプログラムに並べて書いていったのですが、ループで動作するソフトなので追いきれません。 10000桁を出すときは180万回ループします。
試しで桁数10桁の最短にしても100回位まわります。(1ループ50ステップくらいなのでトータルでは最低でも5000ステップくらい。)

*次にやったのは、RPNプログラムを通常の数学の式に変換したものをプログラムに並べて書いてみました。
これで全体としてはかなり解ったのですが、ループで頭に戻って繰り返しの計算になるとまた解らなくなります。
このあたりが連分数の所なのでしょうか。
ある意味、再起呼び出しと同じ動作と言えるのでしょうね。
これになると私の頭では、ごちゃごちゃなってしまいます。

まあ ゆっくり、楽しみながら、諦めないで、あらゆる方向からトライしてみます。

c0335218_21032554.jpg

c0335218_21033187.jpg

c0335218_06272630.jpg

作業部屋の前に新しく作った花壇です。  たたみ一畳ちょっとです。
苗を買ってきて植えたあと一時枯れそうになりましたが、いまは元気になりました。
c0335218_06281295.jpg
不思議です。
花でも野菜の苗でも移植したあとは一度は枯れそうになりますが、そのあとに元気が出てきます。
何故でしょうか?  私の腕が悪いのかな。
c0335218_06281813.jpg




電子工作ランキング

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

by telmic-gunma | 2019-07-06 21:15 | HP電卓 | Comments(4)

DM42で円周率10000桁の計算ができました。
演算時間はUSB電源で10時間 30分前はまだ動いていて10時間目に見たら停止していたので10時間弱ということですね。
これ HP42Sでやると10時間×750倍 => 7500時間 => 312日
電池ももたないし不可能ですね。

演算結果はまた円周率表と比較して正常でした。
(見たのはプログラムの下に示した 前部、中央部、後部についてです。)
ただキャリー処理の部分を停止してあったので直してやり直してみますけれど。
時間のかかるものは、動作の確認するだけでもたいへんですね。


いつもコメントくれるkusakaさんがたまたまそばにいたので「10000桁出来た」と言ったら
「それで」という返事。 
そうなんですよね、電卓や円周率に興味ない人にとっては、「それが どうしたの」という感じなのでしょうね。
わたしにとっては、ゲームセンターの機械よりは、はるかに面白いものですけど。

ここで使う電卓はやはりHP電卓でないとダメですね。
かっこの無いRPN(逆ポーランド法)だからこそ演算のひとつひとつの手順がはっきりとしているのですから。
これ数式通りの電卓やパソコンの数学ソフトでやったら、楽しみが半減すると思います。
こんなこと言っておきながら、わたしまだこのプログラムからどうして円周率が出てくるのか理解できないのです。

まだまだ楽しみは尽きないですね。


c0335218_14441607.jpg




全プログラムリスト 


全プログラムリスト

円周率 1レジスタに10桁表示

Many Digits of Pi by Katie Wasserman - MoHPC
www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/articles.cgi?read=899

オリジナルの計算式

π = 2 + 1/3( 2 + 2/5( 2 + 3/7( 2 + ・・・ ( 2 + k/(2k+1)( 2 + ・・・)))・・・))


計算したのはDM42
  ここまでくると HP42S ではメモリーサイズR880が最大なので動作できない。

01 LBL "AA"
02 SIZE 1005 変数領域の確保 R0~R1004
         (この定義ではR1005をアクセスするとエラーとなる、
         もっと大きく出来るがここまでとする)
03 FIX 00 小数点以下表示しない
04 CLRG 変数領域消去
05 35E3 (35000) ***** 1017=>35000 10000桁計算に概略必要と思われる値
06 STO 01
*

07 LBL 1 mainループの先頭
08 2
09 STO "II"
10 0 CARRYの初期値

*
11 LBL 02 レジスター計算ループの先頭
12 1E10 ***** 1レジスタ10桁表示
13 ×     10,000,000,000倍する
14 RCL IND "II"
15 RCL 01
16 ×
17 +
18 ENTER
19 ENTER
20 ENTER
21 RCL 01
22 2
23 ×
24 1
25 +
26 ÷
27 LSTX
28 X<>Y
29 IP (INTのこと)
30 STO IND "II"
31 ×
32 -
33 1
34 STO +"II"
35 R↓
36 RCL "II"
37 1005 ***** 54=> 1005   レジスタR1004まで計算する
38 -
39 X=0?
40 GTO 03
41 R↓
42 GTO 02

*
43 LBL 03
44 2
45 STO 02
46 RCL 01
47 1
48 -
49 STO 01
50 X>0 ?
51 GTO 01
-------------------------------------
52 1004 ***** 53 => 1004  レジスター計算末尾  
53 STO "II"
54 0 最初はキャリー無しから始まる

*
55 LBL 04 ---桁上がり処理ループ先頭------------
56 STO+ IND "II" 上位のレジスタへキャリーを加算 0 or 1
57 RCL IND "II" O or 1 の加算されたレジスタを読み込む
58 1E10      E6 => E10 10,000,000,000で割り算
59 ÷
60 INT 11桁目(キャリー)を取り出す
61 STO 01 キャリーをR1にセーブしておく
62 1E10 ***** キャリーを10,000,000,000倍する
63 ×
64 STO- IND "II" キャリーを含んだ元の数値からキャリーのみ引く
65 RCL 01 キャリーをR1からXレジスターに戻す
66 DSE "II" インデックスポインター1を引いて 0でないならLBL4へ戻る
67 GTO 04
68 RCL 02 円周率の頭の3を読み込んで終わりとする
69 RTN


計算した答え

R02 2,         ミスでキャリー処理を止めていた。
             キャリーくれば3となり正常
R03 11,415,926,535  ここR03の最初の1が繰り上り
             R02に加算される
R04 8,979,323,846
R05 2,643,383,279  


R500 4,263,129,860
R501 8,099,888,687
R502 4,132,604,721
R503 5,695,162,396

|
|
R1000 9,596,881,592
R1001  560,010,165
R1002 5,256,375,678 最後の8が小数点以下 10000桁目
R1003 5,667,227,966 最後10010桁目
R1004 1,279,988,578 最後10020桁目
R1005 ここをアクセスすると サイズエラーとなる



========
修正 2019-7-26
プログラムの表記でHP42Sでは違うところがいくつか在りましたので訂正しました。



今年の田んぼアート まだ田植えしたばかりで苗が小さいですが何とか読めますね 「令和」と川場村の「カワタン」?だったかな
c0335218_07215930.jpg

この花だけは毎年すごく元気ですね。 いつもこうありたいですけれど。
c0335218_07240977.jpg

タマネギ 掘りました。
大きいのから小さいのまでいろいろあります。


c0335218_07254309.jpg


一番大きいやつです。 実物は大きいのですが写真で見るとそれほど大きく感じないですね。
c0335218_07275809.jpg



電子工作ランキング

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





by telmic-gunma | 2019-06-21 04:05 | HP電卓 | Comments(2)

すがわらさんからのアイデアをいただいて1レジスタ 10桁表示をやってみました。
とりあえず答えは2行だけですが、すんなり出来ました。
ただこれはDM42(HP42S)での話でHP15Cで動作するかは試してありません。
全体の動作を理解できていないのに改造ばかりでいいのかな。

改造個所は以下の3か所です。
1E6 を 1E10 に変更しています。
12行、58行、62行

c0335218_14441607.jpg


全プログラムリスト 


全プログラムリスト

円周率 1レジスタに10桁表示(とりあえずの試し実験)

Many Digits of Pi by Katie Wasserman - MoHPC
www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/articles.cgi?read=899

オリジナルの計算式

π = 2 + 1/3( 2 + 2/5( 2 + 3/7( 2 + ・・・ ( 2 + k/(2k+1)( 2 + ・・・)))・・・))


DM42 ( HP42S )

01 LBL AA
02 SIZE 100 変数領域の確保 R0~R99(もっと大きく出来るがここまで)
03 FIX 00 小数点以下表示しない
04 CLRG 変数領域消去
05 70 ***** 1017 =>70 修正 3レジスタしか計算しない
06 STO 01
*


07 LBL 1 mainループの先頭
08 2
09 STO II
10 0 CARRYの初期値

*
11 LBL 02 レジスター計算ループの先頭
12 1E10 ***** 1レジスタ10桁表示
13 ×     10,000,000,000倍する
14 RCL IND II
15 RCL 01
16 ×
17 +
18 ENTER
19 ENTER
20 ENTER
21 RCL 01
22 2
23 ×
24 1
25 +
26 ÷
27 LSTX
28 X<>Y
29 INT
30 STO(i)
31 ×
32 -
33 1
34 STO +II
35 R↓
36 RCL II
37 5 ***** 54=> 5   レジスタ3個しか計算しない
38 -
39 X=0?
40 GTO 3
41 R↓
42 GTO 2

*
43 LBL 3
44 2
45 STO 2
46 RCL 1
47 1
48 -
49 STO 1
50 X>0 ?
51 GTO 1
-------------------------------------
52 4 ***** 53 => 5 修正  
        桁上がり(キャリー)処理
        を後ろのレジスタ(4番目)より順に処理する
53 STO II
54 0 最初はキャリー無しから始まる

*
55 LBL 4 ---桁上がり処理ループ先頭------------
56 STO +(i) 上位のレジスタへキャリーを加算 0 or 1
57 RCL (i) O or 1 の加算されたレジスタを読み込む
58 1E10      E6 => E10 10,000,000,000で割り算
59 ÷
60 INT 7桁目(キャリー)を取り出す
61 STO 1 キャリーをR1にセーブしておく
62 1E10 ***** キャリーを10,000,000,000倍する
63 ×
64 STO -(i) キャリーを含んだ元の数値からキャリーのみ引く
65 RCL 1 キャリーをR1からXレジスターに戻す
66 DSE II インデックスポインター1を引いて 0でないならLBL4へ戻る
67 GTO 4
68 RCL 2 円周率の頭の3を読み込んで終わりとする
69 RTN


答えの 1レジスタに10桁表示。


R02 3,
R03 1,415,926,535
R04 8,979,323,846  

最後小数点以下 20桁目


========




電子工作ランキング

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

by telmic-gunma | 2019-06-19 01:53 | HP電卓 | Comments(0)

まだ 計算式とプログラムの関係がよく解らないのですが、
いろい調べていて演算桁数を広げる方法が見えてきたので挑戦して見ます。
オリジナルのHP15Cでは変数領域はR00~R53を使用し円周率の計算結果をR02~R53まで使って306桁を出力していましたが
DM42の直接アドレス領域一杯のR02~R99までを使って582桁までやってみました。
間接アドレスを使えばさらに広い範囲をアクセスできますが、そこは次回に挑戦して見ます。
修正箇所は以下の3か所だけです。

37行  54 => 100
52行  53 => 99
5行  1017 => 2040

出力された答えは、 先頭部、中央部、終わりの所を円周率表と見比べて確認しました。

ここが同一ならば全体があっていると言ってもいいのですよね。

c0335218_14441607.jpg


全プログラムリスト 

円周率582桁全プログラムリスト

Many Digits of Pi by Katie Wasserman - MoHPC
www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/articles.cgi?read=899

オリジナルの計算式

π = 2 + 1/3( 2 + 2/5( 2 + 3/7( 2 + ・・・ ( 2 + k/(2k+1)( 2 + ・・・)))・・・))


DM42 ( HP42S )

01 LBL AA
02 SIZE 100 変数領域の確保 R0~R99(もっと大きく出来るがここまで)
03 FIX 00 小数点以下表示しない
04 CLRG 変数領域消去
05 2040   ***** 1017 => 2040 修正
06 STO 01
*

07 LBL 1 mainループの先頭
08 2
09 STO II
10 0 CARRYの初期値

*
11 LBL 02 レジスター計算ループの先頭
12 1E6
13 ×
14 RCL IND II
15 RCL 01
16 ×
17 +
18 ENTER
19 ENTER
20 ENTER
21 RCL 01
22 2
23 ×
24 1
25 +
26 ÷
27 LSTX
28 X<>Y
29 INT
30 STO(i)
31 ×
32 -
33 1
34 STO +II
35 R↓
36 RCL II
37 100 ***** 54=> 100 ここを修正
38 -
39 X=0?
40 GTO 3
41 R↓
42 GTO 2

*
43 LBL 3
44 2
45 STO 2
46 RCL 1
47 1
48 -
49 STO 1
50 X>0 ?
51 GTO 1
-------------------------------------
52 53 ***** 53 => 99 修正  
        桁上がり(キャリー)処理
        を後ろのレジスタ(99番目)より順に処理する
53 STO II
54 0 最初はキャリー無しから始まる

*
55 LBL 4 ---桁上がり処理ループ先頭------------
56 STO +(i) 上位のレジスタへキャリーを加算 0 or 1
57 RCL (i) O or 1 の加算されたレジスタを読み込む
58 E6 1000000で割り算
59 ÷
60 INT 7桁目(キャリー)を取り出す
61 STO 1 キャリーをR1にセーブしておく
62 E6 キャリーを1000000倍する
63 ×
64 STO -(i) キャリーを含んだ元の数値からキャリーのみ引く
65 RCL 1 キャリーをR1からXレジスターに戻す
66 DSE II インデックスポインター1を引いて 0でないならLBL4へ戻る
67 GTO 4
68 RCL 2 円周率の頭の3を読み込んで終わりとする
69 RTN


答えは R3 から R99 までひとつのメモリーに6桁分づつ入っています。
この計算でキャリーのあるところは R32 R23 R19 R17 R11 R3 です。
円周率582桁の答えは以下の通りです。(書くまでもないですね。)
+
R02 3
R03 141592
R04 653589
R05 793238
R06 462643
R07 383279
R08 502884
R09 197169
R10 399375
R11 105820
R12 974944
R13 592307
R14 816406
|
|
R50 393607
R51 260249
R52 141273
R53 724586
|
|
R70 057270
R71 365759
R72 591953
|
|
R97 217176
R98 293176
R99 752384 最後小数点以下 582桁目

========



ジャガイモの花
それなりにきれいですね。
c0335218_14594850.jpg
これもジャガイモ

c0335218_14595700.jpg

日陰のミョウガ   元気です。  まだ移植して数週間ですから。
c0335218_15015816.jpg

日向のミョウガ  こっちのが元気かな?  何年もこの場所だからね。
c0335218_15021878.jpg

これは何でしょうか?  地面がバックだと引き立たないですね。
c0335218_15081407.jpg

三角の形の白のアジサイ、まだ途中なのでこれからですね。
c0335218_15094637.jpg


電子工作ランキング

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

by telmic-gunma | 2019-06-15 14:45 | HP電卓 | Comments(2)

最近、毎日がジタバタしてシヌーという感じ。
仕事では客先で新型商品が出るというのでそれに対応した試験機のハード、ソフトを請け負ってアップアップ
個人的には村の役を引き受けたので、人間関係を含めてガチャガチャ
わたしの脳みそは、ボケかかり。  どうしたらいいんだ。

そうそう ブログでは 「人の悪口と愚痴は書くな」と言われたけどつい出ちゃう。  
やばいやばい。

c0335218_10495804.jpg
パイを計算するのに級数展開というのがあるそうです。
今回 比較的わかり易そうな下図二つ目のライプニッツ級数の式でやってみました。
左辺のパイの下にある4を右辺に移して、左辺は パイ= の形にしました。
しかし こうするとパイがでる、なんてどうして考え付いたのでしょうね。 すごい頭ですね。


c0335218_19244143.jpg



全プログラムリスト 
DM42(HP42Sと同じ) 動作はUSBで電源供給なので電池の3倍速で動いています。

パイを求めるライプニッツ級数
使用変数
00  現時点でのパイ
01 現在の項の分母
02 演算回数 2000 を減算して=> 0
03 1回前のパイ
04 1回前の分母

LBL "AA" ソフトスタート
LCRG   変数オールクリア
1
STO 01  最初の分母 1を設定
4     最初のパイの値を設定     
STO 00
LBL "BB"
2000   途中経過表示のカウンタ(演算2000回ごとに表示)
STO 02
4      項の分子の読み込み
ENTER
2      項の分母の階差の読み込み
RCL + 01 今までの分母に階差を加算する
STO 01  現在の分母をしまう
÷     分子(4)÷分母
STO - 00  過去のパイの演算結果から今の値を引く
RCL 00     現時点のパイの値を一時退避
STO 03
RCL 01  現時点の項の分母を一時退避
STO 04
LBL "CC"
4      項の分子の読み込み
ENT
2      項の分母の階差の読み込み
RCL + 01 今までの分母に階差を加算する
STO 01  現在の分母をしまう
÷     分子(4)÷分母
STO + 00 過去のパイの演算結果から今の値を引く
RCL 03      1回前のパイ
RCL 04  1回前の項の分母
PSE    2000回ごとに引き算の時の値を表示する
RCL 00  現時点のパイ
RCL 01  現時点の項の分母
PSE    2000回ごとに足し算の時の値を表示する
RCL 02  演算回数を読み込む
1
-     演算回数を
X キ 0   演算回数がゼロかチェック
GTO "CC"  ゼロでないなら"CC"へ行く
GTO "BB"  ゼロならば"BB"へ行き2000をセットし直し
END    ソフト 終わり

結果
R/S で演算を停止させたときのレジスタ X,Y,Z,Tを書きます。

T 3.14159254855
Z 19,039,999 元のミスの値 1903999
Y 3.14159275863
X 19,040,001 元のミスの値 1904001

TとYの値は加算する項と減算する項の違いで正解を挟んで上下の値。
計算ごとに段々近づくのですね。
Xはこの時の項の分母となります。
すなわち 演算回数はXの値の約半分の950万回(元ミス95万回)となります。
これだけ計算しても有効桁数は 3.141592 の 7桁しかありません。
とても手計算で出来るやり方ではないですね。

=====================================
追記 2019-04-30
計算を続けてやってみました。

T 3.14159258965
Z 31,279,999
Y 3.14159271753
X 31,280,001

精度は上がらず上下の位置がずれました。
ゆえに 演算桁数を上げずに計算項数を増やしても精度は上がらないと思います。


=====================================
=====================================
追記 2019-04-30
すがわらさんからコメントいただいたのでここに返事を書いてみます。

すがわら氏
横入りすいません。VS入れたのでテストとしてC#で計算させました。FORTRANと違って桁数が少ないので大して精度が出ないのですね。ライプニッツ式で1000万項計算させて、 3.1415926035897882384621752928が得られました。100万項だと、 3.1415921535892932379938932832でした。 無限級数の和は積分に置換できるので、積分からの逆変換じゃないですかね。π/4だからarctanの計算が怪しいかな。

高井
すがわらさんの ライプニッツ式で1000万項計算させて、
3.1415926035897882384621752928
これの9桁目の0は5の移し間違えでないですか、
このままだと3.1415926の計8桁しか有効でありませんね。
ここに5を置くと 3.1415926535897 までの計14桁まで正しくなりますね。

私思うのですがこのライプニッツ式ではいくら演算の項を増やしても精度は上がらないきがします。
右辺2項目の1/3(あるいは4/3)の所で0.333333333333で最後切り捨てられるので演算ソフトの有効桁数以上の精度は出ないでしょう。
HP42Sの場合 12桁(内部演算15桁)でやってますから良くてそこどまりでしょう。

外国のサイトでHP15CやHP32SiiのRPN法でパイを級数を解いているところがありました。
そこ調べてどんなことやってるのか調べて載せてみますね。
ひとの作ったソフトの解析するのを楽しいと感じているうちはまだボケていないのかな。


参考サイト 2個
https://akatuki-724.blogspot.com/2014/05/hp-prime.html

www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/articles.cgi?read=899

========


もう 疲れたので とりあえずここまで。 あとで追加します。
****************************************************************


電子工作ランキング

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

by telmic-gunma | 2019-04-29 19:28 | HP電卓 | Comments(6)