こんにちわ。
>円周率を求めるソフトを作成しようと思っているのですが、
>莫大な桁数の少数を保管できる方法がわかりません。
GENKIさんが仰るように、dim または sdim を使用するのが良いかと思います。
具体的な方法は、円周率を求める方法によって変わってくるのではないでしょうか。
ただ、アークタンジェント法にしろ、算術幾何平均法にしろ、
無限に続く円周率を求めるのはＨＳＰには向かないかと思います。
しかし、もしそのアルゴリズムを理解なさっていて、
円周率の近似値を求めるのが目的でないのであれば、近いことはできそうです。
一番簡単な方法を例に挙げます。
HSPは、標準では double型が限界で、円周率の近似値を sin関数で求めると、
mes " PI = 3.14159265\n"
repeat 5
x = 3.
repeat ( cnt )
x += sin( x )
loop
mes str( cnt + 1 ) + "回目 : " + x
loop
stop
３回目で限界に達します。
これは sin(x) にも double型 にも原因があるためです。
sin関数は初等関数なので、実質上は無限回連続微分が可能なはずです。
そこで、sin関数をTaylor展開によって、置き換えてみます。
sin(x) = ... + { ( -1 )~( n - 1 ) } * { ( x~( 2n - 1 ) )/( ( 2n - 1 )! ) } + ...
（ チルダ(~) は累乗を意味し、{ } は大括弧を意味します）
そうすると、
[A] = ( -1 )~( n - 1 )
[B] = ( x~( 2n - 1 ) )
[C] = ( ( 2n - 1 )! )
に分解できるので、それぞれを計算します。
n (int型・正の整数)を求める項数、t (double型)を合計とすると、
repeat ( n )
_cnt = cnt + 1
lp = ( 2*_cnt - 1 )
// [A]
sign = ( _cnt\2 )*2 - 1
// [B]
vp = 1.0
repeat ( lp )
vp *= x
loop
// [C]
vw = 0
repeat ( lp )
vw += lp
lp --
loop
// sin(x)
t += ( vp / vw * sign )
loop
で求まるはずです。
これを見ると、実際に無限大の実数が必要なのは、
vp *= x で“実数と整数の乗算”
t += ( vp / vw * sign ) で“実数と整数の除算”と“実数と実数の加減算”
の４点だけで済むことになります。
勿論、この方法でも double型の限界があります。
なので、double型もsdim型(str)に置き換えてみます。
無限長の型は、sdim命令で実現します。
これを「XLong型」とします。
以下は、
vp *= x で使うための、XLong × Int のサンプルモジュールです。
#define ctype XLongOut( %1, %2 ) ( %1 + "." + %2 )
#module XLONG_
// ■■■■■■■■■■■■■■■■■■/* by Lisa*/■
// ■ 無限長 XLong型 のサンプル
// ■■■■■■■■■■■■■■■■■■■■■■■■■
#deffunc XLongDeclare var p1, var p2, str p3, str p4
sdim p1, strlen( p3 )
sdim p2, strlen( p4 )
p1 = p3
p2 = p4
return
#deffunc XLongMultipleXI var p1, var p2, int p3
cry = 0 // 繰上げ用変数
repeat 2
if ( cnt ) {
// 計算済みの小数部を、元の小数部変数に戻す
sdim p2, strlen( val )
p2 = val
// 受け取った整数部変数をコピー
len = strlen( p1 )
sdim val, len
val = p1
} else {
// 受け取った小数部変数をコピー
len = strlen( p2 )
sdim val, len
val = p2
}
// 乗算のループ
repeat ( len )
// 後尾から先頭へインデックスを移動する
idx = len - cnt - 1
// 読み出し、数値に置き換え、乗算をして、繰り上がりを足す
num = ( peek( val, idx ) - '0' ) * p3 + cry
// 一桁目を文字コードに戻して、書き換える
poke val, idx, ( ( num \ 10 ) + '0' )
// 二桁目を繰り上げ変数へ保存する
cry = int( num / 10 )
loop
loop
// 整数部での繰り上がりがあるか
if ( cry ) :inc = str( cry ) :else :inc = ""
// 計算済みの整数部を、元の整数部変数に戻す
sdim p1, ( strlen( val ) + strlen( inc ) )
p1 = inc + val
return
#global
/*// XLong型を使用してみる //*/
// 321.823546879012345678900001 を XLong型に代入
XLongDeclare i, d, "821", "323546879012345678900001"
mes
mes "元になる値\n" + XLongOut( i, d )
// XLong変数に 2 を掛ける
XLongMultipleXI i, d, 2
mes
mes "計算した値\n" + XLongOut( i, d )
stop
しかし、ご存知かと思いますが、円周率の世界では
このように筆算アルゴリズムで処理する方法はタブーとされています。
（処理が桁数の二乗になりますので）
これ以上詳しくは、調べてみないと分かりませんが、
このサンプルが何かのヒントになればと思い、長々と書かせて貰いました。
長文・乱文、本当に申し訳ありません。
莫大な少数の扱い
