Many Digits of Pi by Katie Wasserman - MoHPC
www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/articles.cgi?read=899
のプログラムをHP Prime向けに焼き直してみました。このプログラムでは小数点以下1000桁程度までの円周率を計算します。
------
EXPORT ARRYS := 126;
EXPORT DIGIT := 100000000;
// C5(B) -> C4(A)
mv_A_B()
BEGIN
LOCAL i;
FOR i FROM 1 TO ARRYS DO
C4(i) := C5(i);
END;
END;
// div C4(A)/f -> C5(B)
div_A_f(f)
BEGIN
LOCAL i;
LOCAL cy := 0;
FOR i FROM 1 TO ARRYS DO
C5(i) := floor((cy*DIGIT+C4(i)) / f);
cy := (cy*DIGIT+C4(i)) - C5(i)*f;
END;
END;
// mul C4(A)/f -> C5(B)
mul_A_f(f)
BEGIN
LOCAL i, w;
LOCAL cy := 0;
FOR i FROM ARRYS DOWNTO 1 DO
w := C4(i)*f+cy;
C5(i) := w-floor(w/DIGIT)*DIGIT;
cy := floor(w/DIGIT);
END;
END;
EXPORT PI_CALC()
BEGIN
LOCAL n, f, w;
// array initialize
C4:=MAKELIST(0.0,X,1,ARRYS,1);
C5:=MAKELIST(0.0,X,1,ARRYS,1);
// loop countre
n := log(10)/log(2)*1000;
FOR f FROM floor(n) DOWNTO 1 DO
// A * f -> B
mul_A_f(f);
// B -> A
mv_A_B();
// A / f -> B
div_A_f(f*2+1);
// B -> A
mv_A_B();
// A + 2 -> A
C4(1) := C4(1)+2;
END;
END;
--------実行前に、CAS setting の「Exact」フラグを外しておくのが吉です (Exactフラグが付いていると、時間が掛かってしまいます)。
変更 on 2015/03/18
最新のエミュレータで実行した所、何と6分半も掛かってしまいました。また、実行後、画面表示が「凍てつく」状態になり、一旦エミュレータを終了してからでないと、2変数統計表示が出て来ません。
計算結果は2変数統計機能のC4配列変数に収蔵されます。C4(1)には3, 以下、小数点以下を8桁ごとに分割して配列変数に収容しています。
変数DIGIT は、配列要素1つに割り当てる桁数を規定するもので、この数値の場合、配列要素1つにつき、8桁の計算を行います。残念ながら、これ以上にすると計算精度が損なわれる様です。
仕方がないので、1000桁の計算のため、125個 (=1000桁/8桁) の配列要素を用意しています。
AmazonのJulyさん、「日本語クィックガイドなし(正規輸入)」というカテゴリになりました。未だ、「日本語クィックガイドあり」にはモノがありませんが、準備中との事なので、首を長くして待っております。
追記 on 2014/05/27
よくコードを検討したら、若干の変更で、配列変数のコピー部分は不要でしたネ。後で修正版をupしようかしらん ?
