Javaによる任意精度演算

javaによる任意精度演算のプログラムです。
fortranでやろうとしましたがサポートされてないようです。

私はjavaに精通していません。とりあえず何のプログラムでもいいので任意精度演算を使おうとした際、javaが一番手軽だっただけです。

ここでは、ベルヌーイ数を計算することを考えます。
擬コードとして、wikipediaより、

for m in 0, 1, …, n do
A[m] ← 1 / (m + 1)
for j in m, m − 1, …, 1 do
A[j − 1] ← j × (A[j − 1] - A[j])
end for
end for
return A[0] # which is Bn

でベルヌーイ数が得られます。
倍精度演算(10進数で16桁)では桁落ちのため\(B_{14}\)で6桁しか一致しません。

計算過程で100桁まで用いて(割り算で100桁目で四捨五入)、\(B_{50}\)を計算してみましょう。

import java.math.BigDecimal;

public class bernoulli {

final static int divideN=100;

public static void main(String[] args) {
int m,j;
BigDecimal b = new BigDecimal("0");
BigDecimal bm = new BigDecimal("0");
BigDecimal a[]= new BigDecimal[51];

a[0]=new BigDecimal("0");

for(m=0;m<=50;m++){
bm = BigDecimal.valueOf(m);
a[m]=(new BigDecimal("1")).divide(bm.add(new BigDecimal("1")), divideN, BigDecimal.ROUND_HALF_UP);
for(j=m;j>=1;j--){
a[j-1]=(BigDecimal.valueOf(j)).multiply(a[j-1].subtract(a[j]));
}
b=a[0];
String text=String.format("%.36e",b);
System.out.println(text);
}
}
}

コンパイルと実行は

javac bernoulli.java
java bernoulli

です。

出力は36桁としています。上記プログラムによる結果は、
7.500866746076964366855720075757575758e+24
となり、wolframによる値Table[bernoulli b(i), {i,0,50}] with 36 digits
7.50086674607696436685572007575757576e+24
です。全桁一致しました。
[adsense1]

以下、javaの任意精度演算に関するメモです。

メモ

変数の定義方法

整数型
int i,j;
倍精度実数型
double x;
任意精度型
BigDecimal a,b,c;

配列の定義方法

BigDecimal r[]= new BigDecimal[51];

↑で配列r[0]~r[50]が確保される。

配列の初期化

for文でまわしましょう。

for(i=0;i<=N;i++){
r[i]=new BigDecimal("0");
}

変数へ値の代入

aに任意精度型で\(1.2000000000\ldots\)を代入する。

a=new BigDecimal("1.2");

足し算

\(
a=b+1.2
\)

a=b.add(new BigDecimal("1.2"))

(整数型の型変換)
\(
a=b+i
\)
\(i\)は整数型

a=b.add(BigDecimal.valueOf(i));

引き算

\(
a=b-1.2
\)

a=b.subtract(new BigDecimal("1.2"));

掛け算

\(
a=b\times c
\)

a=b.multiply(c);

割り算

\(
a=b/c
\)
丸めの桁を100桁にする時

a=b.divide(c, 100, BigDecimal.ROUND_HALF_UP);

足し算と割り算

\(
a=b/(c+9.92)
\)
丸めの桁を100桁にする時

a=b.divide(c.add(new BigDecimal("9.92")), 100, BigDecimal.ROUND_HALF_UP);

累乗

\(
a=b^c
\)

a=b.pow(c);

ループ

\(
\displaystyle s=\sum_{i=0}^{N-1}a_i
\)

b=new BigDecimal("0");
for(i=0;i<=N-1;i++){
b=b.add(r[i]);
}

端末への指数表記,有効桁数36桁で出力

String text=String.format("%.36e",a);
System.out.println(text);

[adsense2]

参考文献

【Java】任意精度整数(多倍長整数)で演算する -ヽ｜∵｜ゝ(Fantom) の開発blog？
Java BigDecimalの四則演算とフォーマット処理 -ITのおもちゃ箱
 intとBigDecimalの型変換 -Javaプログラミング講座

（更新停止）

Javaによる任意精度演算