# BigDecimal(可変長浮動少数点演算用拡張ライブラリ)

English
BigDecimal はオブジェクト指向の強力なスクリプト言語である Ruby に可変長浮動小数点 計算機能を追加するための拡張ライブラリです。 Ruby についての詳しい内容は以下のURLを参照してください。

- http://www.ruby-lang.org/ja/：Ruby公式ページ
- http://ruby.freak.ne.jp/：Rubyに関する情報ページ
- http://kahori.com/ruby/ring/：Rubyに関するページを辿れます

## 目次

- はじめに
- 使用方法とメソッドの一覧
- 無限、非数、ゼロの扱い
- 内部構造
- 2進と10進
- 計算精度について

## はじめに

Ruby には Bignum というクラスがあり、数百桁の整数でも計算することができます。 ただ、任意桁の浮動少数点演算用クラスが無いようです。そこで、 任意桁の浮動少数点演算用拡張ライブラリ BigDecimal を作成しました。 不具合や助言・提案がある場合どしどし、 shigeo@tinyforest.gr.jp までお知らせください。不具合を直す気は大いにあります。ただ、時間などの関係で約束 はできません。また、結果についても保証できるものではありません。 予め、ご了承ください。
このプログラムは、自由に配布・改変して構いません。ただし、著作権は放棄していません。 配布・改変等の権利は Ruby のそれに準じます。詳しくは README を読んでください。

## インストールについて

BigDecimal を含む Ruby の最新版はRuby公式ページからダウンロードできます。 ダウンロードした最新版を解凍したら、通常のインストール手順を実行して下さい。 Ruby が正しくインストールされれば、同時に BigDecimal も利用できるようになるはずです。 ソースファイルは bigdecimal.c,bigdecimal.h の２個のみです。

## 使用方法とメソッドの一覧

「Rubyは既に書ける」という前提で、
`

    require 'bigdecimal'
    a=BigDecimal::new("0.123456789123456789")
    b=BigDecimal::new("123456.78912345678",40)
    c=a+b

`
というような感じで使用します。

### メソッド一覧

以下のようなメソッドが利用可能です。
記述上、BigDecimal オブジェクトを a,b,c,rで、String（文字列）オブジェクトを s、整数を n で表記します。また、「有効桁数」とは BigDecimal が精度を保証する 桁数です。ぴったりではありません、若干の余裕を持って計算されます。また、 例えば３２ビットのシステムでは１０進で４桁毎に計算します。従って、現状では、 内部の「有効桁数」は４の倍数となっています。

- new
新しい BigDecimal オブジェクトを生成します。
a=BigDecimal::new(s[,n])
s は初期値を文字列で指定します． n は必要な有効桁数（a の最大有効桁数）を整数で指定します。 n が 0 または省略されたときは、n の値は s の有効桁数とみなされます。 s の有効桁数より n が小さいときも n=0 のときと同じです。 a の最大有効桁数は n より若干大い値が採用されます。
- +
加算（c = a + b）
c の精度については「計算精度について」を参照してください。
- -
減算（c = a - b）、または符号反転（c = -a）
c の精度については「計算精度について」を参照してください。
- *
乗算(c = a * b)
cの精度は(aの精度)+(bの精度)程度です。
詳しくは「計算精度について」を参照してください。
- /
除算(c = a / b)
c の精度については「計算精度について」を参照してください。
- assign
以下のように使用します。
c = a.assign(n,f)
f > 0 なら、a を c に、そのまま代入します。 f < 0 なら、-a を c に代入します。 f の絶対値（|f|）は１か２を指定してください。 |f|=2 のときは、c の最大精度が a の実精度より小さいときには 丸められます。|f|=1 のときは切り捨てられます。 n は c の有効桁数です（n 桁以上の精度を持つ c が生成されます）。
- add
以下のように使用します。
c = a.add(b,n)
c = a + b を最大で n 桁まで計算します。 a + b の精度が n より大きいときは丸められます。
- sub
以下のように使用します。
c = a.sub(b,n)
c = a - b を最大で n 桁まで計算します。 a - b の精度が n より大きいときは丸められます。
- mult
以下のように使用します。
c = a.mult(b,n)
c = a * b を最大で n 桁まで計算します。 a * b の精度が n より大きいときは丸められます。
- div
以下のように使用します。
c,r = a.div(b,n)
c=a/b の計算をします。 r には剰余が代入されます。a/bは 必要ならn 桁まで計算されます。divmod メソッド と異なり、c は整数とは限りません。 また、 c は丸められることはありません。 a = c*b + r の関係は成立します。
- %
r = a%b
a/b の余りを計算します。以下の計算と同じものです。
r = a-((a/b).floor)*b
- fix
a の小数点以下の切り捨て。
c = a.fix
- frac
a の整数部分の切り捨て。
c = a.frac
- floor
a 以下の最大整数を表す値（BigDecimal 値）を返します。
c = a.floor
以下のように引数を与えて、小数点以下 n+1 位の数字を操作することもできます (少数点以下を、最大 n 桁にします)。
c = BigDecimal("1.23456")
d = c.floor(4) # d = 1.2345 になります。
n が負のときは小数点以上 n 桁目を操作します。
c = BigDecimal("15.23456")
d = c.floor(-1) # d = 10.0 になります。
- ceil
a 以上の整数のうち、最も小さい整数を計算し、その値（BigDecimal 値）を返します。
c = a.ceil
以下のように引数を与えて、小数点以下 n+1 位の数字を操作することもできます (少数点以下を、最大 n 桁にします)。
c = BigDecimal::new("1.23456")
d = c.ceil(4) # d = 1.2346 になります。
n が負のときは小数点以上 n 桁目をを操作します。
c = BigDecimal::new("15.23456")
d = c.ceil(-1) # d = 20.0 になります。
- round
小数点以下第一位の数を四捨五入して整数（BigDecimal 値）にします。
c = a.round
以下のように引数を与えて、小数点以下 n+1 位の数字を操作することもできます (少数点以下を、最大 n 桁にします)。
n が正の時は、小数点以下 n+1 位の数字を四捨五入します。 c = BigDecimal::new("1.23456")
d = c.round(4) # d = 1.235 になります。
n が負のときは小数点以上 n 桁目をを操作します。
c = BigDecimal::new("15.23456")
d = c.round(-1) # d = 20.0 になります。
- truncate
小数点以下の数を切り捨てて整数（BigDecimal 値）にします。
c = a.truncate
以下のように引数を与えて、小数点以下 n+1 位の数字を操作することもできます (少数点以下を、最大 n 桁にします)。
n が正の時は、小数点以下 n+1 位の数字を切り捨てます。 c = BigDecimal::new("1.23456")
d = c.truncate(4) # d = 1.2345 になります。
n が負のときは小数点以上 n 桁目をを操作します。
c = BigDecimal::new("15.23456")
d = c.truncate(-1) # d = 10.0 になります。
- divmod
商と剰余の配列を返します。
c,r = a.divmod(b) # a = c*b + r
divmodメソッドは a = c * b + r となる a / b の浮動小数点型の商 c と剰余 r を 計算します。ここで c は整数（少数部分のない実数）になります。
c = (a/b).floor
r = a - c*b
で計算されます。
- remainder
r=a.remainder(b)
a/b の剰余 r を計算します。
c = (a/b).fix
r = a - c*b
で計算されます。
- abs
ａの絶対値
c = a.abs
- to_i
少数点以下を切り捨てて整数に変換します。
i = a.to_i
i は値に応じて Fixnum か Bignum になります。 a が Infinity や NaN のとき、i は nil になります。
- to_f
dup と全く同じです。 同じ値の BigDecimal オブジェクトを生成します。
- to_s
文字列に変換します("0.xxxxxEn"の形になります）。
s = a.to_s
- to_s2
文字列に変換します。仮数部分を n 桁毎に空白で区切ります。
s = a.to_s2(n)
- exponent
指数部を整数値で返します。 n = a.exponent
は a の値が 0.xxxxxxx*10**n を意味します。
- E
自然対数の底e(=2.718281828....)を計算します（正直にテイラー展開で）。
e = BigDecimal::E(n)
nは必要な有効桁数を整数で指定します。
- PI
円周率(=3.14159265358979....)を計算します（J．Machinの公式を用います）。
e = BigDecimal::PI(n)
n は必要な有効桁数を整数で指定します。
- BASE
内部で使用される基数の値です。整数が 32 ビットの処理系では10000です。
b = BigDecimal::BASE
- mode
BigDecimalの実行結果を制御します。以下の使用方法が定義されています。
f = BigDecimal::mode(BigDecimal::EXCEPTION_NaN,flag)
f = BigDecimal::mode(BigDecimal::EXCEPTION_INFINITY,flag)
f = BigDecimal::mode(BigDecimal::EXCEPTION_UNDERFLOW,flag)
f = BigDecimal::mode(BigDecimal::EXCEPTION_OVERFLOW,flag)
f = BigDecimal::mode(BigDecimal::EXCEPTION_ZERODIVIDE,flag)
f = BigDecimal::mode(BigDecimal::EXCEPTION_ALL,flag)
EXCEPTION_NaN は結果が NaN になったときの指定です。 EXCEPTION_INFINITY は結果が無限大(±Infinity) になったときの指定です。 EXCEPTION_UNDERFLOW は指数部がアンダーフローするときの指定です。 EXCEPTION_OVERFLOW は指数部がオーバーフローするときの指定です。 EXCEPTION_ZERODIVIDE はゼロによる割り算を実行したときの指定です。 EXCEPTION_ALL は、可能な全てに対して一括して設定するときに 使用します。 flag が true のときは、指定した状態になったときに例外を発行 するようになります。 flag が false（デフォルト）なら、例外は発行されません。計算結果は 以下のようになります。
EXCEPTION_NaN のとき、非数(NaN)
EXCEPTION_INFINITY のとき、無限(+ or -Infinity)
EXCEPTION_UNDERFLOW のとき、ゼロ
EXCEPTION_OVERFLOW のとき、+Infinity か -Infinity
EXCEPTION_ZERODIVIDE のとき、+Infinity か -Infinity
EXCEPTION_INFINITY、EXCEPTION_OVERFLOW、EXCEPTION_ZERODIVIDE は今のところ同じです。
戻り値は、設定後の値です。「値」の意味は、例えば BigDecimal::EXCEPTION_NaNと「値」の & が ゼロ以外ならば EXCEPTION_NaNが設定されているという意味です。 flag が nil、または、true と false 以外なら現在の設定値が返ります。
- limit([n])
生成されるBigDecimalオブジェクトの最大桁数をn桁に制限します。戻り値は 設定する前の値です。設定値のデフォルト値は０で、桁数無制限という意味です。 nを指定しない場合は、現状の最大桁数が返ります。
mf = BigDecimal::limit(n)
- sign
値の属性を返します。 n = a.sign
としたとき n の値は a が以下のときを意味します。
() の中の数字は、実際の値です(「内部構造」を参照)。
n = BigDecimal::SIGN_NaN(0) : a は NaN
n = BigDecimal::SIGN_POSITIVE_ZERO(1) : a は +0
n = BigDecimal::SIGN_NEGATIVE_ZERO(-1) : a は -0
n = BigDecimal::SIGN_POSITIVE_FINITE(2) : a は正の値
n = BigDecimal::SIGN_NEGATIVE_FINITE(-2) : a は負の値
n = BigDecimal::SIGN_POSITIVE_INFINITE(3) : a は+Infinity
n = BigDecimal::SIGN_NEGATIVE_INFINITE(-3) : a は-Infinity
- nan?
a.nan? は a がNaNのとき真を返します。
- infinite?
a.infinite? は a が+∞または-∞のとき真を返します。
- finite?
a.finite? は a が∞または NaN でないとき真を返します。
- to_parts
BigDecimal 値を 0.xxxxxxx*10**n と表現したときに、符号（NaNのときは 0、それ以外は+1か-1になります）、 仮数部分の文字列（"xxxxxxx"）と、基数（10）、更に指数 n を配列で 返します。
a=BigDecimal::new("3.14159265",10)
f,x,y,z = a.to_parts
とすると、f=+1、x="314159265"、y=10、z=1になります。
従って、
s = "0."+x
b = f*(s.to_f)*(y**z)
で Float に変換することができます。
- inspect
デバッグ出力に使用されます。
p a=BigDecimal::new("3.14",10)
とすると、[0x112344:'0.314E1',4(12)]のように出力されます。 最初の16進数はオブジェクトのアドレス、次の '0.314E1' は値、 次の4は現在の有効桁数(表示より若干大きいことがあります)、 最後はオブジェクトが取り得る最大桁数になります。
- dup
同じ値の BigDecimal オブジェクトを生成します。
- sqrt
aの有効桁 n 桁の平方根（n の平方根ではありません）。 これまた、正直にニュートン法で計算します。
c = a.sqrt(n)
- sincos
a の有効桁 n 桁の sin と cos を同時に（テイラー展開で）計算して、 sin と cos の配列を返します。 n は必要な有効桁数です（ n の sin や cos を計算するわけではありません）。
sin,cos = a.sincos(n)
- exp
自然対数の底e(=2.718281828....)の a 乗を計算します。
c = a.exp(n)
n は必要な有効桁数です。
- power
a の n 乗を計算します。ｎは整数。
c = a.power(n)
結果として c の有効桁は a の n 倍以上になるので注意。
- zero?
a が 0 なら true になります。
c = a.zero?
- nonzero?
a が 0 なら false、0 以外なら a そのものが返ります。
c = a.nonzero?
- <=>
a==b なら 0、a > b なら 1、a < b なら -1 になります。
c = a <=> b

後は、読んで字の如くです。

- ==
- ===
「==」と同じですが case 文で使用されます。
- !=
- <
- <=
- >
- >=

### coerceについて

BigDecimal オブジェクトが算術演算子の左にあるときは、BigDecimal オブジェクトが 右にあるオブジェクトを(必要なら) BigDecimal に変換してから計算します。 従って、BigDecimal オブジェクト以外でも数値を意味するものなら右に置けば 演算は可能です。
文字列で数値を与える場合は注意が必要です。数値に変換できない文字があると、 単に変換を止めるだけでエラーにはなりません。"10XX"なら１０、"XXXX"は０ と扱われます。
`

       a = BigDecimal.E(20)
       c = a * "0.123456789123456789123456789" # 文字を BigDecimal に変換してから計算

` 無限大や非数を表す文字として、"Infinity"、"+Infinity"、"-Infinity"、"NaN" も使用できます(大文字・小文字を区別します)。ただし、mode メソッドで false を 指定した場合は例外が発生します。
また、BigDecimalクラスは coerce（Ruby本参照）をサポートしています。 従って、BigDecimal オブジェクトが右にある場合も大抵は大丈夫です。 ただ、現在の Ruby インタプリタの仕様上、文字列が左にあると計算できません。
`

      a = BigDecimal.E(20)
      c = "0.123456789123456789123456789" * a # エラー

` 必要性があるとは思いませんが、どうしてもと言う人は String オブジェクトを継承した新たなクラスを作成してから、 そのクラスで coerce をサポートしてください。

## 無限、非数、ゼロの扱い

「無限」とは表現できないくらい大きな数です。特別に扱うために +Infinity（正の無限大）や -Infinity（負の無限大）という ように表記されます。 無限は 1.0/0.0 のようにゼロで割るような計算をしたときに生成されます。
「非数」は 0.0/0.0 や Infinity-Infinity 等の結果が定義できない 計算をしたときに生成されます。非数は NaN（Not a Number）と表記されます。 NaN を含む計算は全て NaN になります。また NaN は自分も含めて、どんな数 とも一致しません。
ゼロは +0.0 と -0.0 が存在します。ただし、+0.0==-0.0 は true です。
Infinity、NaN、 +0.0 と -0.0 等を含んだ計算結果は組み合わせに より複雑です。興味のある人は、以下のプログラムを実行して結果を 確認してください（結果について、疑問や間違いを発見された方は お知らせ願います）。

    require "bigdecimal"

    aa  = %w(1 -1 +0.0 -0.0 +Infinity -Infinity NaN)
    ba  = %w(1 -1 +0.0 -0.0 +Infinity -Infinity NaN)
    opa = %w(+ - * / <=> > >=  < == != <=)

    for a in aa
      for b in ba
        for op in opa
          x = BigDecimal::new(a)
          y = BigDecimal::new(b)
          eval("ans= x #{op} y;print a,' ',op,' ',b,' ==> ',ans.to_s,\"\n\"")
        end
      end
    end

## 内部構造

BigDecimal内部で浮動小数点は構造体(Real)で表現されます。 そのうち仮数部は unsigned long の配列(以下の構造体要素frac)で管理されます。 概念的には、以下のようになります。
<浮動小数点数> = 0.xxxxxxxxx*BASE**n
ここで、xは仮数部を表す数字、BASEは基数（１０進なら１０）、nは指数部を表す 整数値です。BASEが大きいほど、大きな数値が表現できます。つまり、配列のサイズを 少なくできます。BASEは大きいほど都合がよいわけですが、デバッグのやりやすさなどを 考慮して、10000になっています（BASEはVpInit()関数で自動的に計算します）。 これは、32ビット整数の場合です。64ビット整数の場合はもっと大きな値になります。 残念ながら、64ビット整数でのテストはまだやっていません（もし、やられた方がいれば 結果を教えていただければありがたいです）。 BASEが10000のときは、以下の仮数部の配列(frac)の各要素には最大で４桁の 数字が格納されます。
浮動小数点構造体(Real)は以下のようになっています。
`

      typedef struct {
         unsigned long MaxPrec; // 最大精度(frac[]の配列サイズ)
         unsigned long Prec;    // 精度(frac[]の使用サイズ)
         short    sign;         // 以下のように符号等の状態を定義します。
                                //  ==0 : NaN
                                //    1 : +0
                                //   -1 : -0
                                //    2 : 正の値
                                //   -2 : 負の値
                                //    3 : +Infinity
                                //   -3 : -Infinity
         unsigned short flag;   // 各種の制御フラッグ
         int      exponent;     // 指数部の値(仮数部*BASE**exponent)
         unsigned long frac[1]; // 仮数部の配列(可変)
      } Real;

` 例えば 1234.56784321 という数字は(BASE=10000なら)

        0.1234 5678 4321*(10000)**1

ですから frac[0]=1234、frac[1]=5678、frac[2]=4321、 Prec=3、sign=2、exponent=1 となります。MaxPrecは Prec より大きければいくつでもかまいません。flag の 使用方法は実装に依存して内部で使用されます。

## 2進と10進

BigDecimal は <浮動小数点数> = 0.xxxxxxxxx*10**n という10進形式で数値を保持します。 しかし、計算機の浮動小数点数の内部表現は、言うまでもなく <浮動小数点数> = 0.bbbbbbbb*2**n という 2進形式が普通です(x は 0 から 9 まで、b は 0 か 1 の数字)。 BigDecimal がなぜ10進の内部表現形式を採用したのかを以下に説明します。

#### 10進のメリット

デバッグのしやすさ
  まず、プログラム作成が楽です。frac[0]=1234、frac[1]=5678、frac[2]=4321、 exponent=1、sign=2 なら数値が 1234.56784321 であるのは見れば直ぐに分かります。
10進表記された数値なら確実に内部表現に変換できる
  例えば、以下のようなプログラムは全く誤差無しで 計算することができます。以下の例は、一行に一つの数値 が書いてあるファイル file の合計数値を求めるものです。

         file = File::open(....,"r")
         s = BigDecimal::new("0")
         while line = file.gets
            s = s + line
         end

  この例を2進数でやると誤差が入り込む可能性があります。 例えば 0.1 を2進で表現すると 0.1 = b1*2**(-1)+b1*2**(-2)+b3*2**(-3)+b4*2**(-4).... と無限に続いてしまいます(b1=0,b2=0,b3=0,b4=1...)。ここで bn(n=1,2,3,...) は 2進を表現する 0 か 1 の数字列です。従って、どこかで打ち切る必要があります。 ここで変換誤差が入ります。もちろん、これを再度10進表記にして印刷するような 場合は適切な丸め操作（四捨五入）によって再び "0.1" と表示されます。しかし、 内部では正確な 0.1 ではありません。
有効桁数は有限である（つまり自動決定できる）
  0.1 を表現するための領域はたった一つの配列要素（ frac[0]=1 ）で済みます。 配列要素の数は10進数値から自動的に決定できます。これは、可変長浮動小数点演算では 大事なことです。逆に 0.1 を2進表現したときには2進の有効桁をいくつにするのか 0.1 を 見ただけでは決定できません。

### 10進のデメリット

実は今までのメリットは、そのままデメリットにもなります。 そもそも、10進を2進、2進を10進に変換するような操作は変換誤差 を伴う場合を回避することはできません。 既に計算機内部に取り込まれた2進数値を BigDecimal の内部表現に 変換するときには誤差が避けられない場合があります。

### 最初は何か？

自分で計算するときにわざわざ2進数を使う人は極めてまれです。 計算機にデータを入力するときもほとんどの場合、 10進数で入力します。その結果、double 等の計算機内部 表現は最初から誤差が入っている場合があります。 BigDecimal はユーザ入力を誤差無しで取り込むことができます。 デバッグがしやすいのと、データ読みこみ時に誤差が入らない というのが実際のメリットです。

## 計算精度について

c = a op b という計算(op は + - * /)をしたときの動作は 以下のようになります。
１．乗算と除算は(a の有効桁数)+(a の有効桁数)分の最大桁数（実際は、余裕を持って、 もう少し大きくなります）を持つ変数 c を新たに生成します。 加減算の場合は、誤差が出ないだけの精度を持つ c を生成します。例えば c = 0.1+0.1*10**(-100) のような場合、c の精度は１００桁以上の精度を 持つようになります。
２．次に c = a op b の計算を実行します。
このように、加減算と乗算での c は必ず「誤差が出ない」だけの精度を 持って生成されます。除算は(a の有効桁数)+(a の有効桁数)分の最大桁数 を持つ c が生成されますが、c = 1.0/3.0 のような計算で明らかなように、 c の最大精度を超えるところで計算が打ち切られる場合があります。
いずれにせよ、c の最大精度は a や b より大きくなりますので c が必要とする メモリー領域は大きくなることに注意して下さい。
注意：「+,-,*,/」では結果の精度（有効桁数）を自分で指定できません。 精度をコントロールしたい場合は、以下の add,sub 等のメソッド を使用します。

### 自分で精度をコントロールしたい場合

自分で精度(有効桁数)をコントロールしたい場合は assign、add、sub、mult、div 等のメソッド が使用できます。 以下の円周率を計算するプログラム例のように、 求める桁数は自分で指定することができます。
`

    #!/usr/local/bin/ruby

    #
    # pai.rb
    #  USAGE: ruby pai.rb n
    #   where n is the number of digits required.
    #  EX.: ruby pai.rb 1000
    #

    require "bigdecimal"
    #
    # Calculates 3.1415.... using J. Machin's formula.
    #
    def pai(sig) # sig: Number of significant figures
      exp    = -sig
      pi     = BigDecimal::new("0")
      two    = BigDecimal::new("2")
      m25    = BigDecimal::new("-0.04")
      m57121 = BigDecimal::new("-57121")

      u = BigDecimal::new("1")
      k = BigDecimal::new("1")
      w = BigDecimal::new("1")
      t = BigDecimal::new("-80")
      while (u.exponent >= exp)
        t   = t*m25
        u,r = t.div(k,sig)
        pi  = pi + u
        k   = k+two
      end

      u = BigDecimal::new("1")
      k = BigDecimal::new("1")
      w = BigDecimal::new("1")
      t = BigDecimal::new("956")
      while (u.exponent >= exp )
        t,r = t.div(m57121,sig)
        u,r = t.div(k,sig)
        pi  = pi + u
        k   = k+two
      end
      pi
    end

    if $0 == __FILE__
      print "PAI("+ARGV[0]+"):\n"
      p pai(ARGV[0].to_i)
    end

`

小林 茂雄 (E-Mail:<shigeo@tinyforest.gr.jp>)