mathn_refinements.rb

# MathN using Ruby 2.0 refinements instead of plain monkey patches
# Patrick OSCITY 2013
#
# Based on original mathn implementation by Keiju ISHITSUKA from Ruby 1.8.7-p72
# http://www.opensource.apple.com/source/ruby/ruby-75/ruby/lib/mathn.rb
# 
# Distributes under the same terms as Ruby.

require "complex"
require "rational"
require "matrix"

module MathN
  refine Integer do

    def gcd2(int)
      a = self.abs
      b = int.abs
      a, b = b, a if a < b

      pd_a = a.prime_division
      pd_b = b.prime_division

      gcd = 1
      for pair in pd_a
        as = pd_b.assoc(pair[0])
        if as
          gcd *= as[0] ** [as[1], pair[1]].min
        end
      end
      return gcd
    end

    def Integer.from_prime_division(pd)
      value = 1
      for prime, index in pd
        value *= prime**index
      end
      value
    end

    def prime_division
      raise ZeroDivisionError if self == 0
      ps = Prime.new
      value = self
      pv = []
      for prime in ps
        count = 0
        while (value1, mod = value.divmod(prime)
               mod) == 0
               value = value1
               count += 1
        end
        if count != 0
          pv.push [prime, count]
        end
        break if prime * prime  >= value
      end
      if value > 1
        pv.push [value, 1]
      end
      return pv
    end
  end

  class Prime
    include Enumerable

    def initialize
      @seed = 1
      @primes = []
      @counts = []
    end

    def succ
      i = -1
      size = @primes.size
      while i < size
        if i == -1
          @seed += 1
          i += 1
        else
          while @seed > @counts[i]
            @counts[i] += @primes[i]
          end
          if @seed != @counts[i]
            i += 1
          else
            i = -1
          end
        end
      end
      @primes.push @seed
      @counts.push @seed + @seed
      return @seed
    end
    alias next succ

    def each
      loop do
        yield succ
      end
    end
  end

  refine Fixnum do
    alias / quo
  end

  refine Bignum do
    alias / quo
  end

  refine Rational do
    Unify = true

    def inspect
      format "%s/%s", numerator.inspect, denominator.inspect
    end

    alias power! **

    def ** (other)
      if other.kind_of?(Rational)
        other2 = other
        if self < 0
          return Complex.new!(self, 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.new!(num, den)
      elsif other.kind_of?(Float)
        Float(self) ** other
      else
        x , y = other.coerce(self)
        x ** y
      end
    end

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

        dem = nil
        x = self.denominator.to_f.to_i
        neard = self.denominator.to_f ** (1.0/other.denominator.to_f)
        loop do
          if (neard**other.denominator == self.denominator)
            dem = neaed
            break
          end
        end
        # nearn = self.numerator.to_f ** (1.0/other.denominator.to_f)
        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.new!(num, den)
      elsif other.kind_of?(Float)
        Float(self) ** other
      else
        x , y = other.coerce(self)
        x ** y
      end
    end
  end

  module Math
    def sqrt(a)
      if a.kind_of?(Complex)
        abs = sqrt(a.real*a.real + a.image*a.image)
        #      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.image >= 0
          Complex(x, y)
        else
          Complex(x, -y)
        end
      elsif a >= 0
        rsqrt(a)
      else
        Complex(0,rsqrt(-a))
      end
    end

    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

    module_function :sqrt
    module_function :rsqrt
  end

  refine Complex do
    Unify = true
  end
end