lib/rubysl/mathn/mathn.rb in rubysl-mathn-1.0.0 vs lib/rubysl/mathn/mathn.rb in rubysl-mathn-2.0.0
- old
+ new
@@ -1,218 +1,178 @@
+#--
+# $Release Version: 0.5 $
+# $Revision: 1.1.1.1.4.1 $
+
+##
+# = mathn
#
-# 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.)
+# 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.
#
-# --
+# Without mathn:
#
+# 3 / 2 => 1 # Integer
#
+# With mathn:
#
+# 3 / 2 => 3/2 # Rational
+#
+# mathn features late rounding and lacks truncation of intermediate results:
+#
+# Without mathn:
+#
+# 20 / 9 * 3 * 14 / 7 * 3 / 2 # => 18
+#
+# With mathn:
+#
+# 20 / 9 * 3 * 14 / 7 * 3 / 2 # => 20
+#
+#
+# When you require 'mathn', the libraries for Prime, CMath, Matrix and Vector
+# are also loaded.
+#
+# == Copyright
+#
+# Author: Keiju ISHITSUKA (SHL Japan Inc.)
+#--
+# class Numeric follows to make this documentation findable in a reasonable
+# location
-require "complex.rb"
-require "rational.rb"
+class Numeric; end
+
+require "cmath.rb"
require "matrix.rb"
+require "prime.rb"
-class Integer
+unless defined?(Math.exp!)
+ Object.instance_eval{remove_const :Math}
+ Math = CMath # :nodoc:
+end
- def gcd2(int)
- a = self.abs
- b = int.abs
- a, b = b, a if a < b
+##
+# When mathn is required, Fixnum's division and exponentiation are enhanced to
+# return more precise values from mathematical expressions.
+#
+# 2/3*3 # => 0
+# require 'mathn'
+# 2/3*3 # => 2
- pd_a = a.prime_division
- pd_b = b.prime_division
+class Fixnum
+ remove_method :/
- 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
+ ##
+ # +/+ defines the Rational division for Fixnum.
+ #
+ # 1/3 # => (1/3)
- def Integer.from_prime_division(pd)
- value = 1
- for prime, index in pd
- value *= prime**index
- end
- value
- end
+ alias / quo
- def prime_division
- raise ZeroDivisionError if self == 0
- ps = Prime.new
- value = self
- if value < 0
- value = -value
- pv = [[-1, 1]]
- else
- pv = []
- end
- 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
+ alias power! ** unless method_defined? :power!
-class Prime
- include Enumerable
+ ##
+ # Exponentiate by +other+
- 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
+ def ** (other)
+ if self < 0 && other.round != other
+ Complex(self, 0.0) ** other
+ else
+ power!(other)
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
+##
+# When mathn is required Bignum's division and exponentiation are enhanced to
+# return more precise values from mathematical expressions.
class Bignum
- alias / quo
-end
+ remove_method :/
-class Rational
- Unify = true
+ ##
+ # +/+ defines the Rational division for Bignum.
+ #
+ # (2**72) / ((2**70) * 3) # => 4/3
- def inspect
- format "%s/%s", numerator.inspect, denominator.inspect
- end
+ alias / quo
- alias power! **
+ alias power! ** unless method_defined? :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
+ ##
+ # Exponentiate by +other+
- npd = numerator.prime_division
- dpd = denominator.prime_division
- if other < 0
- other = -other
- npd, dpd = dpd, npd
- end
+ def ** (other)
+ if self < 0 && other.round != other
+ Complex(self, 0.0) ** other
+ else
+ power!(other)
+ end
+ 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
+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
+##
+# When mathn is required Rational is changed to simplify the use 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
- num = Integer.from_prime_division(npd)
- den = Integer.from_prime_division(dpd)
+class Rational
+ remove_method :**
- Rational(num,den)
+ ##
+ # Exponentiate by +other+
+ #
+ # (1/3) ** 2 # => 1/9
- 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)
+ def ** (other)
if other.kind_of?(Rational)
+ other2 = other
if self < 0
- return Complex(self, 0) ** other
+ 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
- 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
+ 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
- nearn = self.numerator.to_f ** (1.0/other.denominator.to_f)
+
+ 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
@@ -222,44 +182,75 @@
den = numerator ** -other
elsif other == 0
num = 1
den = 1
end
- Rational.new!(num, den)
+ 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 required, 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.image*a.image)
- # if not abs.kind_of?(Rational)
- # return a**Rational(1,2)
- # end
+ 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.image >= 0
+# 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)
@@ -301,13 +292,36 @@
sqrt!(a)
end
end
end
+ class << self
+ remove_method(:sqrt)
+ end
module_function :sqrt
module_function :rsqrt
end
-class Complex
- Unify = true
+##
+# 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
+module Rubinius
+ def self.mathn_loaded?
+ true
+ end
+end