added tag v1_9_0_4

git-svn-id: svn+ssh://ci.ruby-lang.org/ruby/tags/v1_9_0_4@18845 b2dd03c8-39d4-4d8f-98ff-823fe69b080e

1 files changed, 264 insertions, 0 deletions

diff --git a/trunk/lib/mathn.rb b/trunk/lib/mathn.rb
new file mode 100644
index 0000000000..ce1acc927f
--- /dev/null
+++ b/trunk/lib/mathn.rb
@@ -0,0 +1,264 @@
+#
+#   mathn.rb - 
+#   	$Release Version: 0.5 $
+#   	$Revision: 1.1.1.1.4.1 $
+#   	by Keiju ISHITSUKA(SHL Japan Inc.)
+#
+# --
+#
+#   
+#
+
+require "complex.rb"
+require "rational.rb"
+require "matrix.rb"
+
+class Integer
+
+  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
+  # These are included as class variables to cache them for later uses.  If memory
+  #   usage is a problem, they can be put in Prime#initialize as instance variables.
+
+  # There must be no primes between @@primes[-1] and @@next_to_check.
+  @@primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101]
+  # @@next_to_check % 6 must be 1.  
+  @@next_to_check = 103            # @@primes[-1] - @@primes[-1] % 6 + 7
+  @@ulticheck_index = 3            # @@primes.index(@@primes.reverse.find {|n|
+                                   #   n < Math.sqrt(@@next_to_check) })
+  @@ulticheck_next_squared = 121   # @@primes[@@ulticheck_index + 1] ** 2
+
+  class << self
+    # Return the prime cache.
+    def cache
+      return @@primes
+    end
+    alias primes cache
+    alias primes_so_far cache
+  end
+  
+  def initialize
+    @index = -1
+  end
+  
+  # Return primes given by this instance so far.
+  def primes
+    return @@primes[0, @index + 1]
+  end
+  alias primes_so_far primes
+  
+  def succ
+    @index += 1
+    while @index >= @@primes.length
+      # Only check for prime factors up to the square root of the potential primes,
+      #   but without the performance hit of an actual square root calculation.
+      if @@next_to_check + 4 > @@ulticheck_next_squared
+        @@ulticheck_index += 1
+        @@ulticheck_next_squared = @@primes.at(@@ulticheck_index + 1) ** 2
+      end
+      # Only check numbers congruent to one and five, modulo six. All others
+      #   are divisible by two or three.  This also allows us to skip checking against
+      #   two and three.
+      @@primes.push @@next_to_check if @@primes[2..@@ulticheck_index].find {|prime| @@next_to_check % prime == 0 }.nil?
+      @@next_to_check += 4
+      @@primes.push @@next_to_check if @@primes[2..@@ulticheck_index].find {|prime| @@next_to_check % prime == 0 }.nil?
+      @@next_to_check += 2 
+    end
+    return @@primes[@index]
+  end
+  alias next succ
+
+  def each
+    return to_enum(:each) unless block_given?
+    loop do
+      yield succ
+    end
+  end
+end
+
+class Fixnum
+  remove_method :/
+  alias / quo
+end
+
+class Bignum
+  remove_method :/
+  alias / quo
+end
+
+class Rational
+  Unify = true
+
+  alias power! **
+
+  def ** (other)
+    if other.kind_of?(Rational)
+      other2 = other
+      if self < 0
+	return Complex.__send__(: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(num, den)
+    elsif other.kind_of?(Float)
+      Float(self) ** other
+    else
+      x , y = other.coerce(self)
+      x ** y
+    end
+  end
+end
+
+module Math
+  remove_method(:sqrt)
+  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