mathn.rb

#
#   mathn.rb - 
#   	$Release Version: 0.5 $
#   	$Revision: 1.1.1.1.4.1 $
#   	$Date: 1998/01/16 12:36:05 $
#   	by Keiju ISHITSUKA(SHL Japan Inc.)
#
# --
#
#   
#

require "complex.rb"
require "rational.rb"
require "matrix.rb"

class Integer

  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

class Fixnum
  alias / quo
end

class Bignum
  alias / quo
end

class Rational
  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

class Complex
  Unify = true
end