ruby--ruby/lib/mathn.rb

##
# = mathn
#
# mathn is a library for changing the way Ruby does math. If you need
# more precise rounding with multiple division or exponentiation
# operations, then mathn is the right tool. It makes sense to use this
# library if you can use it's late rounding. Mathn does not convert
# Fixnums into Floats as long as you do not convert it yourself.
# Instead of using Float as intermediate value it use Rational as
# value representation.
#
# Example Fixnum with intermediate Float:
#
#   20 / 9 * 3 * 14 / 7 * 3 / 2  # => 18
#
# Example: using mathn Fixnum/Rational:
#
#   require 'mathn'
#   20 / 9 * 3 * 14 / 7 * 3 / 2  # => 20
#
# == Usage
#
# To start using this library, simply:
#
#   require "mathn"
#
# This will change the way division works for Fixnums, specifically
#
#   3 / 2
#
# will return +Rational+ (3/2) instead of the usual +Fixnum+ 1.
#
# == Copyright
#
# Author: Keiju ISHITSUKA(SHL Japan Inc.)
#
#--
# $Release Version: 0.5 $
# $Revision: 1.1.1.1.4.1 $

require "cmath.rb"
require "matrix.rb"
require "prime.rb"

require "mathn/rational"
require "mathn/complex"

unless defined?(Math.exp!)
  Object.instance_eval{remove_const :Math}
  Math = CMath
end

##
# When mathn is required Fixnum's division and exponentiation are enhanced to
# return more precise values in mathematical formulas.
#
#   2/3*3  # => 0
#   require 'mathn'
#   2/3*3  # => 2

class Fixnum
  remove_method :/

  ##
  # +/+ defines the Rational division for Fixnum.
  #
  #   1/3  # => (1/3)

  alias / quo

  alias power! ** unless method_defined? :power!

  ##
  # Exponentiate by +other+

  def ** (other)
    if self < 0 && other.round != other
      Complex(self, 0.0) ** other
    else
      power!(other)
    end
  end

end

##
# When mathn is required Bignum's division and exponentiation are enhanced to
# return more precise values in mathematical formulas.

class Bignum
  remove_method :/

  ##
  # +/+ defines the Rational division for Bignum.
  #
  #   (2**72) / ((2**70) * 3)  # => 4/3

  alias / quo

  alias power! ** unless method_defined? :power!

  ##
  # Exponentiate by +other+

  def ** (other)
    if self < 0 && other.round != other
      Complex(self, 0.0) ** other
    else
      power!(other)
    end
  end

end

##
# When mathn is required Rational changes to simplfy the usage of Rational
# operations.
#
# Normal behaviour:
#
#   Rational.new!(1,3) ** 2 # => Rational(1, 9)
#   (1 / 3) ** 2            # => 0
#
# require 'mathn' behaviour:
#
#   (1 / 3) ** 2            # => 1/9

class Rational
  remove_method :**

  ##
  # Exponentiate by +other+
  #
  #   (1/3) ** 2 # => 1/9

  def ** (other)
    if other.kind_of?(Rational)
      other2 = other
      if self < 0
        return Complex(self, 0.0) ** other
      elsif other == 0
        return Rational(1,1)
      elsif self == 0
        return Rational(0,1)
      elsif self == 1
        return Rational(1,1)
      end

      npd = numerator.prime_division
      dpd = denominator.prime_division
      if other < 0
        other = -other
        npd, dpd = dpd, npd
      end

      for elm in npd
        elm[1] = elm[1] * other
        if !elm[1].kind_of?(Integer) and elm[1].denominator != 1
          return Float(self) ** other2
        end
        elm[1] = elm[1].to_i
      end

      for elm in dpd
        elm[1] = elm[1] * other
        if !elm[1].kind_of?(Integer) and elm[1].denominator != 1
          return Float(self) ** other2
        end
        elm[1] = elm[1].to_i
      end

      num = Integer.from_prime_division(npd)
      den = Integer.from_prime_division(dpd)

      Rational(num,den)

    elsif other.kind_of?(Integer)
      if other > 0
        num = numerator ** other
        den = denominator ** other
      elsif other < 0
        num = denominator ** -other
        den = numerator ** -other
      elsif other == 0
        num = 1
        den = 1
      end
      Rational(num, den)
    elsif other.kind_of?(Float)
      Float(self) ** other
    else
      x , y = other.coerce(self)
      x ** y
    end
  end
end

##
# When mathn is requried the Math module changes as follows:
#
# Standard Math module behaviour:
#   Math.sqrt(4/9)     # => 0.0
#   Math.sqrt(4.0/9.0) # => 0.666666666666667
#   Math.sqrt(- 4/9)   # => Errno::EDOM: Numerical argument out of domain - sqrt
#
# After require 'mathn' this is changed to:
#
#   require 'mathn'
#   Math.sqrt(4/9)      # => 2/3
#   Math.sqrt(4.0/9.0)  # => 0.666666666666667
#   Math.sqrt(- 4/9)    # => Complex(0, 2/3)

module Math
  remove_method(:sqrt)

  ##
  # Computes the square root of +a+.  It makes use of Complex and
  # Rational to have no rounding errors if possible.
  #
  #   Math.sqrt(4/9)      # => 2/3
  #   Math.sqrt(- 4/9)    # => Complex(0, 2/3)
  #   Math.sqrt(4.0/9.0)  # => 0.666666666666667

  def sqrt(a)
    if a.kind_of?(Complex)
      abs = sqrt(a.real*a.real + a.imag*a.imag)
#      if not abs.kind_of?(Rational)
#        return a**Rational(1,2)
#      end
      x = sqrt((a.real + abs)/Rational(2))
      y = sqrt((-a.real + abs)/Rational(2))
#      if !(x.kind_of?(Rational) and y.kind_of?(Rational))
#        return a**Rational(1,2)
#      end
      if a.imag >= 0
        Complex(x, y)
      else
        Complex(x, -y)
      end
    elsif a.respond_to?(:nan?) and a.nan?
      a
    elsif a >= 0
      rsqrt(a)
    else
      Complex(0,rsqrt(-a))
    end
  end

  ##
  # Compute square root of a non negative number. This method is
  # internally used by +Math.sqrt+.

  def rsqrt(a)
    if a.kind_of?(Float)
      sqrt!(a)
    elsif a.kind_of?(Rational)
      rsqrt(a.numerator)/rsqrt(a.denominator)
    else
      src = a
      max = 2 ** 32
      byte_a = [src & 0xffffffff]
      # ruby's bug
      while (src >= max) and (src >>= 32)
        byte_a.unshift src & 0xffffffff
      end

      answer = 0
      main = 0
      side = 0
      for elm in byte_a
        main = (main << 32) + elm
        side <<= 16
        if answer != 0
          if main * 4  < side * side
            applo = main.div(side)
          else
            applo = ((sqrt!(side * side + 4 * main) - side)/2.0).to_i + 1
          end
        else
          applo = sqrt!(main).to_i + 1
        end

        while (x = (side + applo) * applo) > main
          applo -= 1
        end
        main -= x
        answer = (answer << 16) + applo
        side += applo * 2
      end
      if main == 0
        answer
      else
        sqrt!(a)
      end
    end
  end

  class << self
    remove_method(:sqrt)
  end
  module_function :sqrt
  module_function :rsqrt
end

##
# When mathn is required Float is changed to handle Complex numbers.

class Float
  alias power! **

  ##
  # Exponentiate by +other+

  def ** (other)
    if self < 0 && other.round != other
      Complex(self, 0.0) ** other
    else
      power!(other)
    end
  end

end