diff options
Diffstat (limited to 'ruby_1_8_6/lib/mathn.rb')
| -rw-r--r-- | ruby_1_8_6/lib/mathn.rb | 308 |
1 files changed, 308 insertions, 0 deletions
diff --git a/ruby_1_8_6/lib/mathn.rb b/ruby_1_8_6/lib/mathn.rb new file mode 100644 index 0000000000..a5a121c6c6 --- /dev/null +++ b/ruby_1_8_6/lib/mathn.rb @@ -0,0 +1,308 @@ +# +# 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 + |