はじめに
Ｒ言語を使ってJPL(ジェット推進研究所)の天体暦DE405/DE430のパラメータ[ftp://ssd.jpl.nasa.gov/pub/eph/planets/ascii/de405/ascp2000.405]
を使用して計算しようとしてますが、18桁の精度が必要です。
精度の扱いがよく分からないため、調べた結果のメモになります。
mpfrの仕様
Ｒのコマンドプロンプトから「mpfr」と入力します。
mpfr
function (x, precBits, base = 10, rnd.mode = c("N", "D", "U", "Z", "A"))
baseの意味がよく分からない。
rnd.modeは「丸めモード」か？
10進18桁を保障するには
19桁目で丸めるとして、
log 2^64 = 19.2
2進64bitあればよいということになります。
確認してみます。
Const("pi", 64)
1 'mpfr' number of precision 64 bits
[1] 3.14159265358979323851
3.14159 26535 89793 23851
238まで正しいので1桁＋小数点以下18桁の計19桁の精度があることが分かります。
DE405/DE430の数値を読み込む
実際には'0.245153650000000000D+07'といった倍精度指数表記なので、一工夫が必要です。
mpfr('0.113078768223014772D+08',64)
mpfr("0.113078768223014772D+08", 64) でエラー:
str2mpfr1_list(x, *): x[1] cannot be made into MPFR
エラーになります。
'D'を'E'に置き換えます。置換にはsub()/gsub()/chartr()がありますが、chartr()を使ってみます。
mpfr(chartr('D','E','0.113078768223014772D+08'),64)
1 'mpfr' number of precision 64 bits
[1] 11307876.8223014772002
大丈夫そうです。
scanを使ってファイルから読み込みます。全て文字列として読み込むためにwhat=''を指定します。
val <- mpfr(chartr('D','E',scan('/proj/R/ascp2000.405',what='')),64)
Read 234038 items
val[1:2]
2 'mpfr' numbers of precision 64 bits
[1] 1 1018
val[3:8]
6 'mpfr' numbers of precision 64 bits
[1] 2451536.5 2451568.5 -33800878.7742210924989
[4] 11307876.8223014772002 397860.220581820933006 -25563.1650982068168005
読み込むことが出来ました。実際の値は以下です。
1 1018
0.245153650000000000D+07 0.245156850000000000D+07 -0.338008787742210925D+08
0.113078768223014772D+08 0.397860220581820933D+06 -0.255631650982068168D+05
10進18桁で出力する。
formatMpfr()を使う。
formatMpfr( val[6], digits=18 )
[1] "11307876.8223014772"
指数表示にする。
formatMpfr( val[6], digits=18, scientific=TRUE )
[1] "1.13078768223014772e7"formatMpfr( val[3:8], digits=18, scientific=TRUE )
[1] "1.13078768223014772e7"formatMpfr( val[3:8], digits=18, scientific=TRUE )
[1] "2.45153650000000000e6" "2.45156850000000000e6" "-3.38008787742210925e7"
[4] "1.13078768223014772e7" "3.97860220581820933e5" "-2.55631650982068168e4"
これでなんとかなったか。
実際の値は以下（再掲）
0.245153650000000000D+07 0.245156850000000000D+07 -0.338008787742210925D+08
0.113078768223014772D+08 0.397860220581820933D+06 -0.255631650982068168D+05