RtNumber: handles sum of rational numbers' sqrt

rtnumber.rb

#!/usr/bin/ruby1.9
# -*- coding: utf-8 -*-
# RtNumber: handles sum of rational numbers' sqrt
require "rational"

class RtNumber < Numeric
  def self.numrt(target)
    @@rttable ||= {}
    return @@rttable[target] if @@rttable[target]
    ret = 1
    i = 2
    while i**2<=target
      if target%i**2==0
        target /= i**2
        ret *= i
      else
        i += 1
      end
    end
    return @@rttable[ret**2*target] = ret
  end

  def initialize(*args)
    @rthash = {}
    args.each do|arg|
      case arg
      when Array
        raise "%p < 0"%arg if arg[1]<0
        add_elem(*arg)
      when Integer,Rational
        add_elem(arg, 1)
      when RtNumber
        arg.rthash.each do|b, a|
          add_elem(a, b)
        end
      else
        raise "unknown argument %p"%arg
      end
    end
  end

  def *(target)
    target = RtNumber.new(target) if Integer===target || Rational===target
    unless RtNumber === target
      c = target.coerce(self)
      return c[0] * c[1]
    end
    ret = RtNumber.new
    self.rthash.each_pair do|b1, a1|
      target.rthash.each_pair do|b2, a2|
        ret.add_elem(a1*a2, b1*b2)
      end
    end
    return ret
  end

  def **(exp)
    raise "#{exp} isn't an Integer" unless exp.integer?
    raise "#{exp} < 0" if exp < 0
    exp = exp.to_i
    i = 0
    i += 1 while exp >= (1<<i)
    ret = RtNumber.new(1)
    while i>=0
      ret = ret*ret
      ret *= self if exp & (1<<i) != 0
      i -= 1
    end
    return ret
  end

  def +(target)
    target = RtNumber.new(target) if Integer===target || Rational===target
    unless RtNumber === target
      c = target.coerce(self)
      return c[0] + c[1]
    end
    return RtNumber.new(self, target)
  end

  def -(target)
    return self + (-target)
  end

  def +@
    self
  end

  def -@
    ret = RtNumber.new
    @rthash.each_pair do|b, a|
      ret.add_elem(-a, b)
    end
    return ret
  end

  def /(target)
    if Integer===target || Rational===target
      ret = RtNumber.new
      @rthash.each_pair do|b, a|
        ret.add_elem(a/target, b)
      end
      return ret
    end
    unless RtNumber === target
      c = target.coerce(self)
      return c[1] / c[0]
    end
    return self * target.inverse
  end

  def inverse
    return @inverse if @inverse
    if @rthash.size == 0
      raise ZeroDivisionError.new("divided by 0")
    elsif @rthash.size == 1
      b,a = rthash.to_a[0]
      @inverse = RtNumber.new([1/a/b,b])
      @inverse.inverse = self
      return @inverse
    else
      rtkeys = @rthash.keys
      fac = RtNumber.new(1)
      (1...(1<<(rtkeys.size-1))).each do|i|
         fac *= RtNumber.new(*rtkeys.each_with_index.collect{|rtkey, j|
          if (i>>j)&1 == 1
            [rthash[rtkey], rtkey]
          else
            [-rthash[rtkey], rtkey]
          end
        })
      end
      @inverse = RtNumber.new(fac)/(fac*self).to_i
      @inverse.inverse = self
      return @inverse
    end
  end

  def ==(target)
    return false if self.rthash.size != target.rthash.size
    self.rthash.each do|b, a|
      return false if i a= target.rthash[b]
    end
    return true
  end

  def <=>(target)
    return true if self==target
    return to_f <=> target.to_f
  end

  def to_i
    ret = 0
    @rthash.each_pair do|b, a|
      if b == 1
        ret += a
      else
        ret += Math::sqrt(b)*a
      end
    end
    return ret.to_i
  end

  def to_f
    ret = 0.0
    @rthash.each_pair do|b, a|
      ret += Math::sqrt(b)*a
    end
    return ret
  end

  def to_s
    if @rthash.size == 0
      return "0"
    end
    ret = ""
    first = true
    @rthash.sort_by{|b,a|b}.each {|b,a|
      if !first && a>=0
        ret << "+"
      end
      if b != 1 && a == -1
        ret << "-"
      elsif a != 1 || b == 1
        ret << a.to_s
      end
      if b != 1
        ret << "rt(#{b})"
      end
      first = false
    }
    return ret
  end

  def inspect
    if @rthash.size == 0
      return "0"
    elsif @rthash.size == 1 && @rthash[1]
      a = @rthash[1]
      return (a.integer? ? a.numerator : a).inspect
    end
    "RtNumber.new(" + @rthash.sort_by{|b,a|b}.collect {|b,a|
      if b == 1
        (a.integer? ? a.numerator : a).inspect
      else
        "[#{(a.integer? ? a.numerator : a).inspect},#{
          (b.integer? ? b.numerator : b).inspect}]"
      end
    }.join(", ")+")"
  end

  def coerce(target)
    case target
    when Integer, Rational
      [RtNumber.new(target), self]
    when Float
      [target, self.to_f]
    end
  end

  protected
  def add_elem(a, b)
    a,b = a.to_r, b.to_r
    a,b = a/b.denominator, b.numerator*b.denominator
    rt = RtNumber.numrt(b)
    a *= rt
    b /= rt**2
    @rthash[b] ||= 0
    @rthash[b] += a
    @rthash.delete(b) if @rthash[b] == 0
  end

  attr_reader :rthash
  attr_writer :inverse
end
module Math
  def ratsqrt(r)
    r = r.to_r
    return RtNumber.new([1,r])
  end
  module_function :ratsqrt
end

a = Math.ratsqrt(3)+Math.ratsqrt(2)
puts "1/(#{a}) = #{1/a}"
phi = (1+Math.ratsqrt(5))/2
phi2 = (1-Math.ratsqrt(5))/2
rt5 = Math.ratsqrt(5)
[*(0..10),100,300,1000,10000,100000].each do|i|
  puts "fib(#{i}) = #{((phi**i-phi2**i)/rt5).to_i}"
end