lib/mathn.rb

   1 #
   2 #   mathn.rb -
   3 #       $Release Version: 0.5 $
   4 #       $Revision: 1.1.1.1.4.1 $
   5 #       by Keiju ISHITSUKA(SHL Japan Inc.)
   6 #
   7 # --
   8 #
   9 #
  10 #
  11
  12 require "complex.rb"
  13 require "rational.rb"
  14 require "matrix.rb"
  15
  16 class Integer
  17
  18   def Integer.from_prime_division(pd)
  19     value = 1
  20     for prime, index in pd
  21       value *= prime**index
  22     end
  23     value
  24   end
  25
  26   def prime_division
  27     raise ZeroDivisionError if self == 0
  28     ps = Prime.new
  29     value = self
  30     pv = []
  31     for prime in ps
  32       count = 0
  33       while (value1, mod = value.divmod(prime)
  34              mod) == 0
  35         value = value1
  36         count += 1
  37       end
  38       if count != 0
  39         pv.push [prime, count]
  40       end
  41       break if prime * prime  >= value
  42     end
  43     if value > 1
  44       pv.push [value, 1]
  45     end
  46     return pv
  47   end
  48 end
  49
  50 class Prime
  51   include Enumerable
  52   # These are included as class variables to cache them for later uses.  If memory
  53   #   usage is a problem, they can be put in Prime#initialize as instance variables.
  54
  55   # There must be no primes between @@primes[-1] and @@next_to_check.
  56   @@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]
  57   # @@next_to_check % 6 must be 1.
  58   @@next_to_check = 103            # @@primes[-1] - @@primes[-1] % 6 + 7
  59   @@ulticheck_index = 3            # @@primes.index(@@primes.reverse.find {|n|
  60                                    #   n < Math.sqrt(@@next_to_check) })
  61   @@ulticheck_next_squared = 121   # @@primes[@@ulticheck_index + 1] ** 2
  62
  63   class << self
  64     # Return the prime cache.
  65     def cache
  66       return @@primes
  67     end
  68     alias primes cache
  69     alias primes_so_far cache
  70   end
  71
  72   def initialize
  73     @index = -1
  74   end
  75
  76   # Return primes given by this instance so far.
  77   def primes
  78     return @@primes[0, @index + 1]
  79   end
  80   alias primes_so_far primes
  81
  82   def succ
  83     @index += 1
  84     while @index >= @@primes.length
  85       # Only check for prime factors up to the square root of the potential primes,
  86       #   but without the performance hit of an actual square root calculation.
  87       if @@next_to_check + 4 > @@ulticheck_next_squared
  88         @@ulticheck_index += 1
  89         @@ulticheck_next_squared = @@primes.at(@@ulticheck_index + 1) ** 2
  90       end
  91       # Only check numbers congruent to one and five, modulo six. All others
  92       #   are divisible by two or three.  This also allows us to skip checking against
  93       #   two and three.
  94       @@primes.push @@next_to_check if @@primes[2..@@ulticheck_index].find {|prime| @@next_to_check % prime == 0 }.nil?
  95       @@next_to_check += 4
  96       @@primes.push @@next_to_check if @@primes[2..@@ulticheck_index].find {|prime| @@next_to_check % prime == 0 }.nil?
  97       @@next_to_check += 2
  98     end
  99     return @@primes[@index]
 100   end
 101   alias next succ
 102
 103   def each
 104     return to_enum(:each) unless block_given?
 105     loop do
 106       yield succ
 107     end
 108   end
 109 end
 110
 111 class Fixnum
 112   remove_method :/
 113   alias / quo
 114 end
 115
 116 class Bignum
 117   remove_method :/
 118   alias / quo
 119 end
 120
 121 class Rational
 122   Unify = true
 123
 124   alias power! **
 125
 126   def ** (other)
 127     if other.kind_of?(Rational)
 128       other2 = other
 129       if self < 0
 130         return Complex.__send__(:new!, self, 0) ** other
 131       elsif other == 0
 132         return Rational(1,1)
 133       elsif self == 0
 134         return Rational(0,1)
 135       elsif self == 1
 136         return Rational(1,1)
 137       end
 138
 139       npd = numerator.prime_division
 140       dpd = denominator.prime_division
 141       if other < 0
 142         other = -other
 143         npd, dpd = dpd, npd
 144       end
 145
 146       for elm in npd
 147         elm[1] = elm[1] * other
 148         if !elm[1].kind_of?(Integer) and elm[1].denominator != 1
 149          return Float(self) ** other2
 150         end
 151         elm[1] = elm[1].to_i
 152       end
 153
 154       for elm in dpd
 155         elm[1] = elm[1] * other
 156         if !elm[1].kind_of?(Integer) and elm[1].denominator != 1
 157          return Float(self) ** other2
 158         end
 159         elm[1] = elm[1].to_i
 160       end
 161
 162       num = Integer.from_prime_division(npd)
 163       den = Integer.from_prime_division(dpd)
 164
 165       Rational(num,den)
 166
 167     elsif other.kind_of?(Integer)
 168       if other > 0
 169         num = numerator ** other
 170         den = denominator ** other
 171       elsif other < 0
 172         num = denominator ** -other
 173         den = numerator ** -other
 174       elsif other == 0
 175         num = 1
 176         den = 1
 177       end
 178       Rational(num, den)
 179     elsif other.kind_of?(Float)
 180       Float(self) ** other
 181     else
 182       x , y = other.coerce(self)
 183       x ** y
 184     end
 185   end
 186 end
 187
 188 module Math
 189   remove_method(:sqrt)
 190   def sqrt(a)
 191     if a.kind_of?(Complex)
 192       abs = sqrt(a.real*a.real + a.image*a.image)
 193 #      if not abs.kind_of?(Rational)
 194 #       return a**Rational(1,2)
 195 #      end
 196       x = sqrt((a.real + abs)/Rational(2))
 197       y = sqrt((-a.real + abs)/Rational(2))
 198 #      if !(x.kind_of?(Rational) and y.kind_of?(Rational))
 199 #       return a**Rational(1,2)
 200 #      end
 201       if a.image >= 0
 202         Complex(x, y)
 203       else
 204         Complex(x, -y)
 205       end
 206     elsif a >= 0
 207       rsqrt(a)
 208     else
 209       Complex(0,rsqrt(-a))
 210     end
 211   end
 212
 213   def rsqrt(a)
 214     if a.kind_of?(Float)
 215       sqrt!(a)
 216     elsif a.kind_of?(Rational)
 217       rsqrt(a.numerator)/rsqrt(a.denominator)
 218     else
 219       src = a
 220       max = 2 ** 32
 221       byte_a = [src & 0xffffffff]
 222       # ruby's bug
 223       while (src >= max) and (src >>= 32)
 224         byte_a.unshift src & 0xffffffff
 225       end
 226
 227       answer = 0
 228       main = 0
 229       side = 0
 230       for elm in byte_a
 231         main = (main << 32) + elm
 232         side <<= 16
 233         if answer != 0
 234           if main * 4  < side * side
 235             applo = main.div(side)
 236           else
 237             applo = ((sqrt!(side * side + 4 * main) - side)/2.0).to_i + 1
 238           end
 239         else
 240           applo = sqrt!(main).to_i + 1
 241         end
 242
 243         while (x = (side + applo) * applo) > main
 244           applo -= 1
 245         end
 246         main -= x
 247         answer = (answer << 16) + applo
 248         side += applo * 2
 249       end
 250       if main == 0
 251         answer
 252       else
 253         sqrt!(a)
 254       end
 255     end
 256   end
 257
 258   module_function :sqrt
 259   module_function :rsqrt
 260 end
 261
 262 class Complex
 263   Unify = true
 264 end