зеркало из https://github.com/github/ruby.git
315 строки
6.7 KiB
Ruby
315 строки
6.7 KiB
Ruby
#
|
|
# 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
|
|
|
|
remove_method :inspect
|
|
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
|
|
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
|
|
|