2010年1月26日火曜日

HP35sで、円周率を計算してみたよ

前のエントリで紹介しました次の記事を参考に、HP35sで円周率を計算してみました。

Many Digits of Pi (HP32sii, HP-16c, HP-12c, HP-12cp, HP-19c, HP-30b)
http://www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/articles.cgi?read=899

この記事のトップにある、HP32sii向けのコードが移植しやすく思われ、試行錯誤を重ねて移植、実行してみた所、およそ01:05時間程度で結果が得られました。
桁数は小数点以下200桁です。以下、オリジナルのコードを移植する際の「差分」について述べておきます。

1) 32siiの行番号は英文字1文字 + 数字2桁ですが、35sは英文字1文字 + 数字3桁と異なります。また、GTOなどのジャンプでは32siiの場合には数字部分を指定しません(出来ないらしい)が、35sでは指定する必要があります。

2) 変数のマッピングが異なります。

2-1) 32siiでは、変数A~Zが配列として、1~26の番号を用いて操作可能ですが、35sではA~Zは-1~-26で操作されます。
2-2) 32siiでは、配列操作のインデクスに、特殊変数「i」を用意していますが、35sでは「I」、「J」を用います。
2-3) 統計機能の変数として「n」(Σ+キーを使って入力した、データの個数)を使いますが、この場所が、32siiでは配列番号28でアクセスされるのに対し、35sでは-27となります。

以上を踏まえて、書き換えたものを載せておきます。

original author ; Katie Wasserman (許可は取っていない ... ; 英語出来ない)


P001 | LBL P
P002 | CLVARS
P003 | CLΣ
P004 | 1 ; Set up index to store variable 'A' (on 35s, it is assigned to index 0)
P005 | STO J
P006 | -27 ; Set up to store the iteration counter in register n
P007 | STO I
P008 | 26 ; calculate number of iterations needed = log-base-2(10) * N
P009 | 8
P010 | *
P011 | 2
P012 | LOG
P013 | /
P014 | IP
P015 | STO(I) ; Store in n
P016 | 2
P017 | STO(J) ; starting value
P018 | *
P019 | 3
P020 | +
P021 | Σ+ ; the SIGMAx register is used to store (iteration_counter x 2 + 1)

B001 | LBL B ; main loop
B002 | 1.026 ; set up for register loop
B003 | STO I
B004 | 0 ; initial carry

C001 | LBL C
C002 | 1E8 ; 100,000,000 this determines the number of digits per register
C003 | *
C004 | n
C005 | RCL*(I)
C006 | + ; multiply and add carry
C007 | ENTER
C008 | ENTER
C009 | Σx
C010 | /
C011 | IP ; divide by (iteration_counter x 2 + 1)
C012 | STO(I)
C013 | Σx
C014 | *
C015 | - ; remainder is the carry
C016 | ISG I
C017 | GTO C001
C018 | 2
C019 | STO+(J)
C020 | Σ- ; decrement iteration_counter by 1; and (2 x iteration_counter + 1) by 2
C021 | x>0?
C022 | GTO B001
C023 | 26 ; set up for register loop to carry overflows up one register,
C024 | STO I ; register value can never be more than 2 x 10^8

E001 | LBL E
E002 | RCL(I)
E003 | 1E8
E004 | -
E005 | x<=0?
E006 | GTO F001 ; no carry
E007 | STO(I)
E008 | 1
E009 | STO- I ; these 3 instructions add the carry forward
E010 | STO+(I)
E011 | STO+ I
E012 | GTO E001 ; this isn't really needed since the overflow can never be more than 1E8, ; but it's nice to have if other algorithms are used
F001 | LBL F
F002 | DSE I
F003 | GTO E001
F004 | RTN

S001 | LBL S
S002 | 1.027
S003 | STO I
S004 | VIEW(I)
S005 | ISG I
S006 | GTO S004
S007 | RTN


実行は[XEQ] P [ENTER]です。実行が終わったら、[XEQ] S [ENTER]で、計算結果が見られます。結果は、配列1~26に、8桁毎に入っていて、それを表示するだけです。

いわゆる「Quick hack」的な書き換えなので、実行効率などは判りません。もっと、良いコードになる可能性はあると思います。
しかし、32sii(40分)より35s(1時間5分)の方が遅いとは ... トホホ。一方、30bでは600桁以上計算するのに、40分掛からないとか。やはり、ARMベースで新しい関数電卓とか、出ないかな ?

4 件のコメント:

hp12c さんのコメント...

金融電卓で検索してたどりきました。
RPN方式のhp12cにはまっています。

akatuki さんのコメント...

hp12c様、今晩は。
blog、拝見しました。12cの使い方について丁寧に解説されて居られます様で、なかなか良いですね。
当方は、金融関係は疎いものですから、余り良く判らんのですが、同じHP電卓の利用者として、頑張りましょう。(何だか、良く判らないような ?)

hiro0903865 さんのコメント...

お久しぶりです。
といっても、以前何をコメントしたのか覚えていないのですが、ははは。
まさかパイ計算を35sで出来るとは考えもしませんでした。
元32siiプログラマーさま、電卓情報更新者さまお疲れさまでした。

35sは主にEQNを利用して簡単な式をプログラム?する程度ですが
暇がありましたら掲載プログラム参考にさせていただきたいと思います。

ところで、私も12cさまと同じく12cノーマルを入手して遊んでおります。
過去発売されていた15cというのがなるほど、確かに使いやすかったのかもと実感しております。

akatuki さんのコメント...

hiro0903865 様、今晩は。遅くなりました。

> といっても、以前何をコメントしたのか覚えていないのですが、ははは。

当blogは余りコメントを戴いていないので、恐らくお名前を変えていらっしゃるのかしらん ?

> 35sは主にEQNを利用して簡単な式をプログラム?する程度ですが
> 暇がありましたら掲載プログラム参考にさせていただきたいと思います。

いやぁね、当方も先日、ようやくこのプログラムの動作について判りかけて来た次第でして。結構トリッキィな感じです。統計メモリをうまく使っていますね。この辺に気を付ければ、他の電卓やポケコンでも動くプログラムに改変できると思います。

> ところで、私も12cさまと同じく12cノーマルを入手して遊んでおります。
> 過去発売されていた15cというのがなるほど、確かに使いやすかったのかもと実感しております。

12cも使って居られますか。当方は12cは使っていないのですが、15Cは名機だと思っています(壊れちゃったけれど)。プログラム領域は余り多くないのですが、行列計算や複素数の計算もOKで、かなり勉強しました。35sに15cの代替を期待したのですが、行列計算がプログラムに依るという所は惜しい、かと。
まあ、いずれ、30bのハードウェアで科学技術計算向けの機種が出ないか、などと期待しまして。