#
# An extension to the standard library's Rational
#
# Includes conversions from String and Float, inversion, reduction methods
# (trim and approximate), and a to_s with a base conversion parameter.
#
# A bunch of this code is taken from the Python project:
#   http://cvs.sourceforge.net/viewcvs.py/python/python/nondist/sandbox
#   /rational/Rational.py?rev=1.3&view=markup
#
# Author: Dave Burt <dave at burt.id.au>
# Created: 5 May 2005
# Last modified: 8 May 2005
#
# Methods:
#
# Rational(n, d = 1)           # replaced - now accepts Floats and Strings
#                              # like Rational(1.3e15) or Rational("4/5")
#
# Rational::to_s(n = 10)       # replaced - now accepts a radix
# Rational::to_f               # replaced - doesn't return NaN any more
# Rational::hash               # replaced - r.hash==r.to_i.hash if r==r.to_i
# Rational::trim(max_d)        # simplify approximately, set max denominator
# Rational::approximate(err = 1e-12)  # simplify approximately, within +/-err
# Rational::inv                # invert
#
# Float::to_r   # converts to Rational, exactly
# String::to_r  # converts to Rational if possible, or returns Rational(0)
#

# Rational() has been improved to handle Floats and Strings, so you can do
# +Rational(1.3e15)+ and +Rational("2/3")+
alias std_Rational Rational
def Rational(a, b = 1)
  if b == 1
    case a
    when Rational
      a
    when String
      case a
      when /\//
        n, d = a.split('/', 2)
        Rational(n) / Rational(d)
      when /e/
        mant, exp = a.split('e', 2)
        Rational(mant) * (10 ** Integer(exp))
      when /\./
        i, f = a.split('.', 2)
        Rational(Integer(i)) + Rational(Integer(f), 10 ** f.length)
      else
        Rational(Integer(a))
      end
     when Float
       a.to_r
     when Integer
       Rational.new!(a)
    end
  elsif a.kind_of?(Integer) and b.kind_of?(Integer)
    Rational.reduce(a, b)
  else
    Rational(a) / Rational(b)
  end
rescue ArgumentError
  raise ArgumentError.new("invalid value for Rational: #{num.inspect}")
end


class Rational

  # Cast to Float. This improved version won't return NaN for Rationals with
  # large numerators or denominators.
  def to_f
    r = if @denominator > (1 << 1022)  # Where did the extra bit go?
      self.trim(1 << 1022)           # (Python handles 1 << 1023)
    else
      self
    end
    r.numerator.to_f / r.denominator.to_f
  end
  
  # This version returns the same hash as the numerator itself if
  # the denominator is 1.
  def hash
    @numerator.hash ^ @denominator.hash ^ 1.hash
  end
  
  # Return the closest rational number such that the denominator is at most
  # +max_d+
  def trim(max_d)
    n, d = @numerator, @denominator
    if max_d == 1
      return Rational(n/d, 1)
    end
    last_n, last_d = 0, 1
    current_n, current_d = 1, 0
    begin
      div, mod = n.divmod(d)
      n, d = d, mod
      before_last_n, before_last_d = last_n, last_d
      next_n = last_n + current_n * div
      next_d = last_d + current_d * div
      last_n, last_d = current_n, current_d
      current_n, current_d = next_n, next_d
    end until mod == 0 or current_d >= max_d
    if current_d == max_d
      return Rational(current_n, current_d)
    end
    i = (max_d - before_last_d) / last_d
    alternative_n = before_last_n + i*last_n
    alternative_d = before_last_d + i*last_d
    alternative = Rational(alternative_n, alternative_d)
    last = Rational(last_n, last_d)
    if (alternative - self).abs < (last - self).abs
      alternative
    else
      last
    end
  end
  
  # Return the simplest rational number within +err+
  def approximate(err = Rational(1, 1e12))
    r = self
    n, d = @numerator, @denominator
    last_n, last_d = 0, 1
    current_n, current_d = 1, 0
    begin
      div, mod = n.divmod(d)
      n, d = d, mod
      next_n = last_n + current_n * div
      next_d = last_d + current_d * div
      last_n, last_d = current_n, current_d
      current_n, current_d = next_n, next_d
      app = Rational(current_n, current_d)
    end until mod == 0 or (app - r).abs < err
    p "mod==0" if mod == 0
    app
  end
  
  # Return the inverse
  def inv
    Rational(@denominator, @numerator)
  end
  
  # Represent the fraction as a string, in the given base.
  def to_s(base = 10)
    if @denominator == 1
      @numerator.to_s(base)
    else
      @numerator.to_s(base) + "/" + @denominator.to_s(base)
    end
  end
end

class Float
  #  Return Rational(top, bot), where top and bot are a pair of co-prime
  #  integers such that x = top/bot.
  #  
  #  The conversion is done exactly, without rounding.
  #  bot > 0 guaranteed.
  #  Some form of binary fp is assumed.
  #  Pass NaNs or infinities at your own risk.
  #  
  #  >>> rational(10.0)
  #  rational(10L, 1L)
  #  >>> rational(0.0)
  #  rational(0L)
  #  >>> rational(-.25)
  #  rational(-1L, 4L)
  def to_r
    x = self
    
    if x == 0.0
      return Rational(0, 1)
    end
    signbit = false
    if x < 0
      x = -x
      signbit = true
    end
    f, e = Math.frexp(x)
    raise unless 0.5 <= f and f < 1.0
    # x = f * 2**e exactly
    
    # Suck up chunk bits at a time; 28 is enough so that we suck
    # up all bits in 2 iterations for all known binary double-
    # precision formats, and small enough to fit in an int.
    chunk = 28
    top = 0
    # invariant: x = (top + f) * 2**e exactly
    while f > 0.0
      f = Math.ldexp(f, chunk)
      digit = f.to_i
      raise unless digit >> chunk == 0
      top = (top << chunk) | digit
      f = f - digit
      raise unless 0.0 <= f and f < 1.0
      e = e - chunk
    end
    raise if top == 0
    
    # now x = top * 2**e exactly; fold in 2**e
    r = Rational(top, 1)
    if e > 0
      r *= 2**e
    else
      r /= 2**-e
    end
    if signbit
      -r
    else
      r
    end
  end
end # class Float

class String
  # Convert to Rational, returning Rational(0) if it's can't be converted
  def to_r
    Rational(self) rescue Rational(0)
  end
end



if $0 == __FILE__
  
  require 'rational_ext_test'

end