浮動10進少数点四則演算をBigInt演算で記述

多倍長整数演算さえあれば、それらを組み合わせて好きな進数の任意の精度の浮動小数点四則演算を記述することができます。ということで、コード例を。Ruby の BigInt はメモリーが許す限りの精度で整数演算できるので、これを使うと精度も指数の桁も好きなだけ増やせるのですが、ここでは検算を楽におこなうためにIEEE754 の倍精度浮動2進小数点演算結果との比較しやすさを優先し、10進数8桁の小数部をもち、指数部は-1000から+1000の間に制限しています。丸めはゼロ方向への四捨五入だけに対応しています。NaN と無限大を定義してない手抜きな実装です。
小数部の計算は、富豪的に10進数17桁の倍精度プラス1桁の固定小数点でおこない、絶対値が10より小さな10進17桁の固定小数点を整数で表しているとみなして、正規化により小数部の絶対値が1より小さく0.1以上になるように調整します。それから、10進数下8桁の絶対値から同上9桁を繰り上げるかどうかを判定して丸めをおこない、再度正規化を繰り返します。調整がすべて済んだら、小数点以下の上位8桁分を切り出して、新しい浮動小数点数を作ります。このロジックは D. KnuthThe Art of Computer Programming Vol.2 asin:4756145434 記載の正規化アルゴリズムの通りです。MIX 計算機のアセンブラのコードより、何を計算しているのかが、読み取りやすくなっているのは高級言語ならではですけど、普通、この手のコードはアセンブラで記述するものなのでしょう。

2014年4月16日改訂: 除算だけが必ず切り捨てになっていたのを加減乗の丸め方式にあわせる。
2011年12月15日改訂: ruby で 10 ** i で i がゼロ以上のときは Integer になるのに対して、10 ** -1 等の i が負のとき、Float になるのを失念していたので、常に Integer のサブクラスで少数部の計算がおこなわれるように変更しました。

class DecFloat
  BASE = 10
  PREC = 8
  E_MAX = 1000
  E_MIN = -E_MAX

  attr_accessor :e, :f

  def self.[](e, f = nil)
    if f.nil?
      e, f = PREC, e
    end
    new(e, f) + new(E_MIN, 0)
  end

  def initialize(e, f)
    unless Integer === e and Integer === f
      raise ArgumentError, '(Integer, Integer)'
    end
    if e > E_MAX
      raise RangeError, "exponent overflow #{e}"
    elsif e < E_MIN
      raise RangeError, "exponent underflow #{e}"
    end
    if zero?
      e = E_MIN
    end
    @e, @f = e, f
  end

  def to_s() inspect end

  def inspect
    sign = (if @f < 0 then '-' else '' end)
    i, f = @f.abs.divmod(BASE ** (PREC - 1))
    '%s%d.%07de%+03d' % [sign, i, f, self.e - 1]
  end

  def zero?
    @f == 0
  end
  
  def abs
    self.class.new(@e, @f.abs)
  end

  def +(v)
    u = self
    if u.e < v.e
      u, v = v, u
    end
    if v.zero? or u.e - v.e >= PREC + 2
      norm(u.e, u.f * BASE ** PREC)
    elsif u.zero?
      norm(v.e, v.f * BASE ** PREC)
    else
      s = PREC - u.e + v.e
      if s >= 0
        norm(u.e, u.f * BASE ** PREC + v.f * BASE ** s)
      else
        norm(u.e, u.f * BASE ** PREC + v.f / (BASE ** (-s)))
      end
    end
  end

  def -(v)
    self + (-v)
  end

  def -@
    self.class.new(@e, -@f)
  end

  def *(v)
    norm(self.e + v.e, self.f * v.f)
  end

  def /(v)
    norm(PREC - 1 + self.e - v.e, (BASE ** (PREC + 1)) * self.f / v.f)
  end

private

  def norm(e, fg)
    unless Integer === e and Integer === f
      raise ArgumentError, '(Integer, Integer)'
    end
    one = BASE ** (PREC * 2)        # 1.0
    under = BASE ** (PREC * 2 - 1)  # 1.0 / b
    boundary = BASE ** PREC         # fg = f * boundary + g
    sign = fg < 0 ? -1 : +1
    fg = fg.abs
    if fg != 0
      loop do
        while fg >= one
          fg /= BASE
          e += 1
        end
        while fg < under
          fg *= BASE
          e -= 1
        end
        f, g = fg.divmod(boundary)
        if g >= boundary / 2
          f += 1
        end
        fg = f * boundary
        break if fg < one
      end
    end
    self.class.new(e, (fg / boundary) * sign)
  end
end

u = DecFloat[-26, 66256000]   # plank constant 6.6256e-27 erg sec
f = 0.66256000e-26
v = DecFloat[-29, 87654321]
g = 0.87654321e-29
printf "    u = %s,     f = %.7e\n", u.to_s, f
printf "    v = %s,     g = %.7e\n", v.to_s, g
printf "u + v = %s, f + g = %.7e\n", (u + v).to_s, f + g
printf "u - v = %s, f - g = %.7e\n", (u - v).to_s, f - g
printf "v - u = %s, g - f = %.7e\n", (v - u).to_s, g - f
printf "u * v = %s, f * g = %.7e\n", (u * v).to_s, f * g
printf "u / v = %s, f / g = %.7e\n", (u / v).to_s, f / g
printf "v / u = %s, g / f = %.7e\n", (v / u).to_s, g / f

実行結果は次の通りです。

$ ruby decfloat.rb
    u = 6.6256000e-27,     f = 6.6256000e-27
    v = 8.7654321e-30,     g = 8.7654321e-30
u + v = 6.6343654e-27, f + g = 6.6343654e-27
u - v = 6.6168346e-27, f - g = 6.6168346e-27
v - u = -6.6168346e-27, g - f = -6.6168346e-27
u * v = 5.8076247e-56, f * g = 5.8076247e-56
u / v = 7.5587830e+02, f / g = 7.5587831e+02
v / u = 1.3229643e-03, g / f = 1.3229643e-03