[ruby-dev:4978] Re: compare to NaN (Re: Re: call for finite/isnan testers)

えぐち＠エスアンドイー です。

>>> On Wed, 10 Feb 1999 00:44:57 +0900, えぐち said:

えぐち> えぐち＠エスアンドイー です。
えぐち> 
えぐち> >>> On Tue, 9 Feb 1999 21:50:25 +0900, [e-mail] (石塚圭樹 ) said:
えぐち> 
えぐち> keiju> けいじゅ＠日本ラショナルソフトウェアです.
...
えぐち> keiju> (^^;;; これって, Cで実装しはじめているって意味ですかね?
えぐち> 
えぐち> いやぁ、いきなり C って事はさすがにないです。 ;-)
えぐち> IEEE 754 の巧妙なデータパッキングを見た影響でしょう。
えぐち> 
えぐち> Float#<=> は実は、Integer#<=>そのものだと言うあたりの
えぐち> 影響は受けました。

仮数部固定長(１０進３２桁)で Bigfloat を実装しました。
Ruby のコードです。  添付します。

BASE はともかく、PREC を可変出来るように設計したつもりです。

あと、無限／非数／正規化への配慮もしました。

このレベルのコードをいぢるといろいろわかるもので

 def inf?(n)
     n == n + 1.0
 end

は誤りで、ちょっと n が大きくなると
すぐ無限と早とちりします。^^;;;

 def inf?(n)
     n == n * 1.1
 end

が良い様です。(1.01 でも 1.000000000001 でもいいけど、、)

また、
	a = 1.0/0.0)
	a.to_i 

するとインタプリタがかえって来なくなるような
雰囲気ですが、まだ追求していません。

	えぐち

#
# bigfloat.rb:
#
require "e2mmap"

def Bigfloat(num, exp = 0)
  Bigfloat.new(num, exp)
end

# dummy definition
class Rational<Numeric;end
class Complex<Numeric;end

def inf?(n)
  n == n * 1.1
end
def nan?(n)
  n != n
end

class Bigfloat < Numeric
  extend Exception2MessageMapper
  def_e2message(TypeError, "wrong argument type %s (expected %s)")

  BASE = 10 # XXX base of @exp
  PREC = 32 # XXX
  def Bigfloat.count_prec_of_float
    i, a, b = 0, 1, 1.0 / BASE
    begin
      o, a, b, i = a, a + b, b / BASE, i + 1
    end until(a == o)
    i
  end
  FLOAT_PREC = Bigfloat.count_prec_of_float

  attr :sign # plus(and zero):+1, minus:-1, +Inf:+Inf, -Inf:-Inf, NaN:NaN
  attr :prov # PREC'digits prov.
  attr :exp  # base BASE of exp.

  def normalize!
    if @prov == 0 then
      @sign, @prov, @exp = +1, 0, 0
      return
    end
    @sign, @prov = +1, @sign * @prov
    @sign, @prov = -@sign, -@prov if @prov < 0
    while @prov < BASE ** PREC
      @prov *= BASE
      @exp -= 1
    end
    while @prov >= BASE ** (PREC + 1)
      @prov /= BASE
      @exp += 1
    end
  end

  def initialize(num, exp = 0)
    @sign, @prov, @exp = +1, 0, 0
    case num
    when Integer
      @sign, @prov, @exp = +1, num, exp
    when Rational
      Bigfloat.new(num.numerator) / Bigfloat.new(num.denominator)
      return
    when Float
      if inf?(num) or nan?(num) then
	@sign, @prov, @exp = num, 0, 0
	return
      end
      @sign = +1;
      @sign, num = -1, -num if num < 0
      exp10 = Integer(Math.log10(num))
      @exp = exp10 - FLOAT_PREC
      @prov = Integer(num / (BASE ** @exp))
    when String
      case num 
      when /^\s*([+]?Inf)$/ 
	@sign = +1.0 / 0.0
	return
      when /^\s*(-Inf)$/ 
	@sign = -1.0 / 0.0
	return
      when /^\s*NaN$/ 
	@sign = 1.0 % 0.0
	return
      when /^\s*([-+]?[0-9]*)\.?([0-9]*)([eE]([+-]?[0-9]+))?$/
	@prov = $1.to_i * (10 ** $2.size) 
	@sign, @prov = -1, -@prov if @prov < 0 
	@prov += $2.to_i
	@exp = -$2.size
	@exp += $4.to_i if $4
      end
    when Bigfloat
       @sign, @prov, @exp = num.instance_eval("[@sign, @prov, @exp ]")
    else
      Bigfloat.fail(TypeError, 
		    num.type,
		    "Inetger, Float, String, Bigfloat or Rational")
    end

    normalize!
  end

  def inspect
    "Bigfloat(#{to_a})"
  end

  def to_s
    return @sign.inspect if inf?(@sign) or nan?(@sign)
    s = ""
    s += "-" if @sign < 0
    n, d = @prov.divmod(10 ** PREC)
    s += n.inspect
    s += "."
    s += format("%0*d", PREC, d.inspect)
#    return s if @prov == 0 or (@exp + PREC) == 0 
    s += "E"
    s += (@exp + PREC).inspect
  end

  def abs
    sign, prov, exp = instance_eval("[@sign, @prov, @exp]")
    Bigfloat.new(sign, prov, exp)
  end

  def to_i
    Integer(@sign * @prov * 10 ** @exp)
  end

  def to_f
    @sign * @prov * 10.0 ** @exp
  end

  def <=>(other)
    case other
    when Bigfloat
      #break
    when Integer, Rational, Float, String
      other = Bigfloat(other)
    when Complex
      return Complex(self) <=> other
    when Array
      Bigfloat.fail(TypeError, 
		    other.type,
		    "Inetger, Float, String, Bigfloat or Rational")
      return -1;
    else
      x, y = other.coerce(self)
      return x <=> y
    end
    result = @sign <=> other.sign 
    return result if result != 0
    result = @exp  <=> other.exp
    return result if result != 0
    return @prov <=> other.prov
  end

  def +(other)
    case other
    when Bigfloat
      #break
    when Integer, Rational, Float
      other = Bigfloat(other)
    when Complex
      return Complex(self) + other
    else
      x, y = other.coerce(self)
      return x + y
    end

    left  = @sign * @prov
    right = other.sign * other.prov
    if @exp > other.exp then
      right *= BASE ** (@exp - other.exp)
      exp = @exp
    else
      left *= BASE ** (other.exp - @exp)
      exp = other.exp
    end
    Bigfloat.new(left + right, exp)
  end

  def -(other)
    case other
    when Bigfloat
      #break
    when Integer, Rational, Float
      other = Bigfloat(other)
    when Complex
      return Complex(self) - other
    else
      x, y = other.coerce(self)
      return x - y
    end

    left  = @sign * @prov
    right = other.sign * other.prov
    if @exp > other.exp then
      right *= BASE ** (@exp - other.exp)
      exp = @exp
    else
      left *= BASE ** (other.exp - @exp)
      exp = other.exp
    end
    Bigfloat.new(left - right, exp)
  end

  def *(other)
    case other
    when Bigfloat
      #break
    when Integer, Rational, Float
      other = Bigfloat(other)
    when Complex
      return Complex(self) * other
    else
      x, y = other.coerce(self)
      return x * y
    end

    sign = @sign * other.sign
    prov = @prov * other.prov
    exp  = @exp  + other.exp
    Bigfloat.new(sign * prov, exp)
  end

  def /(other)
    case other
    when Bigfloat
      #break
    when Integer, Rational, Float
      other = Bigfloat(other)
    when Complex
      return Complex(self) / other
    else
      x, y = other.coerce(self)
      return x / y
    end

    sign = @sign * other.sign
    prov = @prov * BASE ** PREC / other.prov
    exp  = @exp  - PREC - other.exp 
    Bigfloat.new(sign * prov, exp)
  end

  def %(other)
    Bigfloat.new ( self - Integer ( self / other ) * other )
  end

  def coerce(other)
    case other
    when Integer, Rational, Float
      return Bigfloat(other), self
    when Complex
      return other, Complex(self)
    else
      super
    end
  end
end