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