多桁円周率の計算(1)

akatuki様が管理人をされているプログ 「高機能電卓の情報」 で、HP の Prime という機種搭載の言語で、円周率を1000桁まで計算してみた、というエントリーがあります。

HP Prime で円周率を計算してみたよ


HP Prime だけでなく、様々なプログラム電卓に移植する試みがゆっくりと進んでいます。

ここで採用した計算方法は、オイラーの arctan(X) を使った展開公式で、X=1 の時の式を使って、以下のように変形したものを使っています。

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

そして、この式を計算するために、Spigot アルゴリズムを使用したプログラムが akatuki様により HP Prime 用に移植されました。
Spigot アルゴリズムは、多倍長変数計算をせずに、計算結果を頭から順次得てゆく方法です。


当該プログで akatuki様が移植されたプログラムから、sentaro様が Casio Basic へ移植されました。それに私も少しお手伝いして出来たのが以下のプログラムです。


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"


fx-5800P用ソースコードのダウンロードは、こちら から。

このプログラムは、メモリ消費量を抑える工夫がされています。

ちなみに、fx-5800P で BASE (結果変数1つあたりの桁数) を 10 として、100桁の円周率を計算させると、10分56秒かかりました。


プログラムを起動すると、桁数 (DIGITS) と 結果変数1つあたりの桁数 (BASE) を聞いてくるので、それぞれ 100 と 10 を入力します。

PICALC_settings 

すると計算が始まり、しばらく(11分弱) 待つと、結果が表示されます。[EXE] キーを押すたびに、10桁の計算結果が表示されます。

小数点以下、30桁 (10桁 x3) が表示されたところ、

PICALC_result_1 

さらに、続く 4桁を表示させると、

PICALC_result_2 

さらに、続き計算結果を表示させると、

PICALC_result_3 

最後に END が表示されたら終わりです。

この画面の1番上の行は、6から始まる9桁です。これは6の前に0があるのですが、数値の最大桁が0の場合は、その0は表示されないので、このようになっています。この行は、0628620899 と読み替えます。

これで、円周率の小数点以下100桁が求まりました。


fx-5800P で条件を色々と変えて遊ぶには時間がかかりすぎるので、fx-9860GII Casio Basic へ移植しました。

fx-9860GII Casio Basic版
ファイル名: PICALC


ClrMat Z
"Digit"?→X
"Base"?→K
Int (X÷K)+1→A
10^K→D
{X+1,1}→Dim Mat Z
X÷log 2→N
"EXE To Go"◢

ClrText
Locate 2+Int log N,1,"/"
Locate 3+Int log N,1,Int N

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
Mat Z[I,1]F+C→W
Int (W÷D)→C
W-CD→Mat Z[I,1]
Next

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

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

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

ClrMat Z
"End"


PC経由で fx-9860GII へ転送できる g1m ファイル のダウンロードは、こちら


Base 10で、100桁の計算を行ったら、fx-9860GII ノーマルクロック(29MHz) で 2分37秒、オーバークロック(280MHz)では 24秒でした。

移植先には、上記の他に、sentaro様により、fx-CG10 / 20 の Casio Basic (fx-9860GIIと同じでOK)、fx-9860 と fx-0860GII の Add-in、TI Nspire CX CAS の BASIC とCネイティブ整数版と実数版、...と色々な移植による速度比較が公開されています。



なお、e-Gadget としては、fx-5800P で、もう少し高速化ができないか?に興味があります。今のところ、10分50秒くらいかかっていたのが、7分50秒くらいまで高速化出来ています。

もう少し検討してから次回紹介しようと思います。




応援クリックをお願いします。励みになるので...

人気ブログランキングへ


FC2ブログランキングへ



 


keywords: fx-5800Pfx-9860GIICasioBasicプログラム関数電卓多桁円周率

リンク集 | ブログ内マップ





関連記事

テーマ : プログラム関数電卓
ジャンル : コンピュータ

コメントの投稿

非公開コメント

お!

ついに新エントリーできましたね!

このネタはハマる人にはたまらないネタだと思います(^^)

とりあえず、管理人様バージョンのfx-CG10/20版(といっても丸ごと同じですが(^^;)も含めたCasioBasic版です。
http://pm.matrix.jp/Casio_PICALC.zip

おっ !

親分、夜分に恐れ入ります。

円周率計算の記事、かたじけなく。
今度は、肩の力を抜いたネタを準備しているので、お楽しみ(?)に。
(あまり期待しないでネ)

Re: お!

sentaro様

モンテカルロ法は、グラフィックを使ってみる口実でしたが、円周率問題ですね。
で、ひょんなことから、今回は円周率第二弾...って感じになりました。

表題には(1)とあるように、次は高速化を紹介する予定です。


> このネタはハマる人にはたまらないネタだと思います(^^)

はい、ですね。他の方のブログでもかなり濃い内容が多いことも、ハマったら深いことを証明しています。


> とりあえず、管理人様バージョンのfx-CG10/20版(といっても丸ごと同じですが(^^;)も含めたCasioBasic版です。
> http://pm.matrix.jp/Casio_PICALC.zip

ありがとうございます。


Re: おっ !

akatuki様

毎度、お世話になっています。

今回は見事にハマってしまいました。この手の話は、本当に楽しいですね。


> 円周率計算の記事、かたじけなく。

元ネタは、もう1年近い前のことでしたよね。
当時は、多倍長変数を扱うものだと決めてかかっていて、以前 Visual Basic や C で多倍長整数の計算をやったことがあるだけに、fx-5800Pで多倍長整数は、かなり無理だと思っただけでに、移植も性根を入れないと、つまりたっぷり時間が必要だと、勝手に思っていた次第。


まさかakatuki様の取り上げたアルゴリズムが、多倍長を使わない方法だとは知らなかったのです。

今回改めて調べたら、1000なり10000なりをドンと掛け算してから、頭から整数部を取り出すという見事な方法で、多倍長変数を使わない方法だと分かったのです(Spigotアルゴリズム)。

それなら fx-5800Pでもなんとかなりそう、と俄然やる気が出てきました。


> 今度は、肩の力を抜いたネタを準備しているので、お楽しみ(?)に。
> (あまり期待しないでネ)

脱力系ですか?



表示速度

管理人様、こんにちは!

表示のオーバーヘッドですがちょっと調べてみました。

332->N
For 1->I To N
Locate 2,1," "
Locate 1,1,I
Next

表示部分だけを抜き出したものですけど、ループ回数100桁での332回で計測しました。

fx-5800P     14.3秒
fx-9860GII(SH3) 14.3秒
fx-9860GII(SH4A) 11.2秒
fx-CG10      27.7秒

fx-5800Pとfx-9860GII(SH3)が同じ速度というのがちょっと驚きですけど、fx-CG10はやはり遅いです。
内部演算速度が速くなるにつれてこの表示がかなり負担になってきますね。

CasioBasic版が管理人様により完成度が上がってるので再度fx-CG10/20版を更新しました。

http://pm.matrix.jp/Casio_PICALC.zip

追加でPICALC2の計算式最適化バージョンPICALC2Aを加えてます(^^ゞ



akatuki様、こんにちは!

>今度は、肩の力を抜いたネタを準備しているので、お楽しみ(?)に。

楽しみにお待ちしてます!(^^)

効率化→速度向上

sentaro様

関西へ向かう新幹線の中です。

効率化した PICALC2A ですが、確かに効率化の効果がでていて、計算速度が向上しています。
ありがとうございます。

ところで、PICALC2A ですが、試しに 280桁、Base を 6 にして計算させると、以下の箇所でエラーがでています。


Int (T÷B)→W
C+W→Mat Z[F+I,1] <<-- ここでエラー
T-BW→C
Isz I:N-E→N
LpWhile N>0


推測ですが、プログラムの冒頭にある

Int (A÷log2)→F



Int (E(A÷K+A))→F

にすると、エラーが消えました。F が必要数より小さく、多分計算に必要な配列がとれなくなったことによるエラーではないかと思われます。

この部分の計算は、計算項数を求める部分だと思いますが、

PiCalc2 のロジックで、理屈から考えますと、

[計算項数] = ( Int([桁数] / [Base] ) + 1) x ([Base] / log 2 )

となると思うので(100%の自信なし)、F が大きくなり、結果としてより多くの配列表域が確保できて、エラーが消えたのではないかと...


そこで、sentaro様の PICAL2A の効率化部分をそっくりそのままいただいて、PICALC2 と合体させ、上の計算項数の式を使って、PICALC3 にしてみました。

PiCalc3 fx-9860GII 版
----------------------------------------------
Filename:PICALC3

ClrMat Z
Lbl 0
"Digits"?→A
"Base"?→K
10^K→B
Int (K÷log 2)→E◢
Int (E(A÷K+1))→F
Int (A÷K+1)→H
F+H◢
If F+H+1>999:Then
ClrText
Locate 1,4,"Digits exceeds 999"
Goto 0:IfEnd
F→N
{F+H+1,1}→Dim Mat Z
"EXE To Go"◢

ClrText
For 1→I To N
Locate 2,1," "
Locate 1,1,I
Int (B÷5)→Mat Z[I,1]
Next

0→C:1→I
Do
Locate 2,1," "
Locate 1,1,N

0→T
For N-1→J To 1 Step -1
'Locate 2,2," "
'Locate 1,2,J
2J-1→D
TJ+BMat Z[J,1]→W
Int (W÷D)→T
W-DT→Mat Z[J,1]
Next

Int (T÷B)→W
C+W→Mat Z[F+I,1]
T-BW→C
Isz I:N-E→N
LpWhile N>0

For 1→I To H
Mat Z[F+I,1]◢
Next

ClrMat Z
""
Locate 1,7,"Last digit:"
Locate 13,7,KH-1
----------------------------------------------

さて、PiCalc2 で Base 7 で計算すると良いと思っていたのですが、計算結果の照合をしたところ Base 6 でないと途中の桁から結果が違ってしまうことに気がつきました。うっかりしてました。

で、実際の計算結果が何桁までかを、最後に表示するようにして、照合が楽になるようにしました。ここで表示されるのは、頭の 3 を含まず、小数点以下の桁数です。


g1mファイルは、以下の e-Gadgetアーカイブの一番下から...

http://egadget2.web.fc2.com/archives/archives.html#fx9860GII_CasioBasic



ちなみに、PiCalc3 の fx-5800P 版も作って 100桁、Base=6 で計算させると、101桁まで結果が表示されますが、計算は100桁まで正しくて、5分53秒 と PiCalc2 の 6分35秒 よりもさらに高速化しました。(両方とも、内側ループのLocate をコメントアウト/削除しています)

PiCalc2 で削除するだけで1分速くなりますよね(既にご指摘の通り) 。
fx-5800P では、効率化の効果が見た目に大きいですね。


PiCalc3 fx-5800P 版
----------------------------------------------

Cls
Lbl 0
"DIGITS"?→A
"BASE"?→K
10^(K)→B
Int(K÷log(2))→E◢
Int(E(A÷K+1))→F
Int(A÷K+1)→H
F+H◢
F→N
F+H+1→DimZ
"EXE TO GO"◢

Cls
For 1→I To N
Locate 2,1," "
Locate 1,1,I
Int(B÷5)→Z[I]
Next

0→C:1→I
Do
Locate 2,1," "
Locate 1,1,N

0→T
For N-1→J To 1 Step -1
2J-1→D
TJ+BZ[J]→W
Int(W÷D)→T
W-DT→Z[J]
Next

Int(T÷B)→W
C+W→Z[F+I]
T-BW→C
Isz I:N-E→N
LpWhile N>0

For 1→I To H
Z[F+I]◢
Next

0→DimZ
""
Locate 1,4,"LAST DIGIT:"
Locate 12,7,KH-1
----------------------------------------------

sentaro様

fx-9860GII 版で追加した、行列の仕様以上の計算項数になる場合のIf 文の処理で、表示が紛らわしいので、適切な言葉に変えたほうが、良さそうです。

例えば、

Calc items > 999

とか...

精度問題

管理人様、こんにちは!
関西出張おつかれさまです。

>さて、PiCalc2 で Base 7 で計算すると良いと思っていたのですが、計算結果の照合をしたところ Base 6 でないと途中の桁から結果が違ってしまうことに気がつきました。うっかりしてました。

Spigotアルゴリズムでの必要精度のあたりはまだよくわかってないのですけど、CasioBasicも6桁ですがHP Primeも6桁でないと計算できないです。
ただ、TIだと7桁で1000桁でも合うのを確認できました。
ということで、
精度のことがちょっと気になって単純な計算で調べてみたところ、
CASIOの電卓モードでは、
1e14+1-1e14→0
1e13+1-1e13→0
1e12+1-1e12→1
1e11+1-1e11→1
ということで、1e13つまり14桁の整数の計算は上手くいかないことがわかりました。
内部精度が15桁あっても実際にきっちり使えるのは13桁ということになります。
なにやら14桁目と15桁目の最後2桁は四捨五入?されるみたいですね。
ちなみにHP Primeは内部、表示ともに12桁でそのまま12桁精度です。
TI Nspireでは内部14桁、表示12桁ですが14桁まるまる精度があるようです。
今回の円周率計算では計算精度が一番高かったのはTIということに…。


というところで、早速に管理人様の新しいPICALC3試してみました(^^)

280桁の場合ですが、
----------------------
PICALC3

952(905+47) // 確保する配列数(計算項数+計算結果格納配列数)

564856
692346
 34861
 45345
207403
Last digit:281
----------------------

最後の桁ですが正確には、

564856
692346
 34861
 45345
664821

なので、最後の5桁が間違ってます。

調べてみると、280桁の精度を得るには計算項数が931回必要なのですが、PICALC3だと計算項数が905回と少なくなっていたので精度が下がったというわけです。
100桁のときは計算項数が335で必要項数の332を越えてるので精度には問題なかったというわけですね。

ということで、再度精度を決める計算項数のFを算出する計算式は、
管理人様の、
Int (E(A÷K+A))→F
これだと桁数によってはループ回数が少なくなって精度が足りなくなるので、
Int(A÷log 2)+1→F
で正解だろうという結論です。
計算項数は計算時間に直結なので精度が確保できる最小限である必要がありますね。

あと、配列数を求める計算式は、
Int((A+K-1)÷K)+1→H
でよさそうです。

私の改修版PICALC2Aの280桁で配列が足りなくなったのはこの2箇所の部分で確保数が足りなかったということですね。
そこを修正してみるとエラーは出なくなって無事完走しました。

----------------------
PICALC2A補正版

931  // 計算項数
47   // 配列確保数

564856
692346
 34861
 45345
664820
----------------------

fx-5800Pだと配列がZしかないですけど、fx-9860GIIだとそれ以外も使えるので結果格納配列は別に確保した方がより計算桁数増やせそうですが、それ以上に計算桁数を増やすなら最初の配列最小限アルゴリズムに戻すしかなさそうです。



>fx-5800P では、効率化の効果が見た目に大きいですね。

割り算とか配列アクセスとか時間のかかる処理を必要最小限にするだけで演算速度にかなり効果ありますね。
fx-5800Pは表示は比較的速いですが内部がちと遅いので演算量削減は効果的に効いてくると思います(^^)


>例えば、
>Calc items > 999
>とか...

コード節約で、
単純にオーバーした場合はArgument ERRORで済ますというのは?(^^;

No title

管理人様、こんにちは!

PICALCもそろそろ完成系という感じなのでまとめました(^^)
akatuki様のところでもアップしたfx-9860GIIとCG10/20対応のCasioBasic版です。
http://pm.matrix.jp/Casio_PICALC2.zip


fx-9860GIIのCのアドイン版も桁数変更できるようにバージョンアップしてみました。
http://pm.matrix.jp/picalcFX2.zip

No title

管理人様、こんにちは!

先週アップした円周率プログラムは桁数変えていろいろ試してみてたら気が付かなかったバグ?があったりで奥の深さ&詰めの甘さを改めて感じてます(^^;

1変数あたりの桁数を1にするとPICAL1Bでは問題ないもののオリジナル系統のPICAL_Bでは最初の桁からおかしな結果になって上手くいきません。
Spigotアルゴリズム版のPICAL2B、PICAL4Aでは途中の桁でおかしな結果になるようです。
調べてみるとPICAL_Bでは多倍長の掛け算と割り算の間でC=0を省くと上手くいくようになりました。
1桁にしたときは最上位桁の桁上がりが重要な意味を持っていたのですね。

PICAL2B、PICAL4Aでは計算結果を配列に格納していくところで桁上がりが発生して1桁じゃなくなるところがあってそれが原因でした。
これは1桁だけの問題ではないのですが、メモリあたり複数桁だと桁上がりが露見しなかっただけのようです。

----------------
Int (T÷B)→W
C+W→Mat Z[F+I,1] // ←桁上がりが発生する場合がある
T-BW→C
Isz I:N-E→N
LpWhile N>0
----------------

正確な結果を得るには、PICAL1Bと同様の結果配列補正ルーチンを演算終了後に入れると上手くいくようです。

4~6桁で計算する分には何も影響がないのでずっと気がつきませんでした(^^;

ってことで、バグ取れたと思うCasioBasicバージョンを再アップしておきます。
http://pm.matrix.jp/Casio_PICALC2.zip

修正版のfx-9860GIIアドイン版もアップしておきます。
http://pm.matrix.jp/picalcFX2b.zip
最小1桁対応しました。計算結果の保証はないです(^^;
相変わらず計算後にリブートになるエラーが出ることがあるのですが原因不明です…

少ない Base 数について

sentaro様

ご指摘ありがとうございます。

私も Base = 2 や 3 の時には、計算結果が違うことがある点に気付いていました。

具体的な対処方法をご提示頂いて、ありがあとうございます。

いずれ、内容をまとめようと思います。


最新記事
最新コメント
カテゴリ
C# (3)
検索フォーム
Visitors
Online Counter
現在の閲覧者数:
プロフィール

やす (Krtyski)

Author:やす (Krtyski)
since Oct 30, 2013


プログラム関数電卓は、プログラムを作り・使ってナンボ!

実際に触って気づいたこと、自作プログラム、電卓プログラミングについて書いています

おもしろい・役に立つならクリックしてください。励みになります。
にほんブログ村 IT技術ブログ 開発言語へ
にほんブログ村


人気ブログランキングへ


FC2ブログランキングへ


写真: 「4駆で泥んこ遊び@オックスフォード郊外」

リンク
月別アーカイブ
Sitemap

全ての記事を表示する

RSSリンクの表示
最新トラックバック
ブロとも申請フォーム

この人とブロともになる

QRコード
QR