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

遅うなりました。HP Primeで円周率を計算してみました。

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フラグが付いていると、時間が掛かってしまいます)。

PCのシミュレータでも20秒くらいは掛かります。
変更 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しようかしらん ?