lib/mathn.rb
DEFINITIONS
This source file includes following functions.
1 #
2 # mathn.rb -
3 # $Release Version: 0.5 $
4 # $Revision: 1.1.1.1.4.1 $
5 # $Date: 1998/01/16 12:36:05 $
6 # by Keiju ISHITSUKA(SHL Japan Inc.)
7 #
8 # --
9 #
10 #
11 #
12
13 require "rational.rb"
14 require "complex.rb"
15 require "matrix.rb"
16
17 class Integer
18
19 def gcd2(int)
20 a = self.abs
21 b = int.abs
22 a, b = b, a if a < b
23
24 pd_a = a.prime_division
25 pd_b = b.prime_division
26
27 gcd = 1
28 for pair in pd_a
29 as = pd_b.assoc(pair[0])
30 if as
31 gcd *= as[0] ** [as[1], pair[1]].min
32 end
33 end
34 return gcd
35 end
36
37 def Integer.from_prime_division(pd)
38 value = 1
39 for prime, index in pd
40 value *= prime**index
41 end
42 value
43 end
44
45 def prime_division
46 ps = Prime.new
47 value = self
48 pv = []
49 for prime in ps
50 count = 0
51 while (value1, mod = value.divmod(prime)
52 mod) == 0
53 value = value1
54 count += 1
55 end
56 if count != 0
57 pv.push [prime, count]
58 end
59 break if prime * prime >= value
60 end
61 if value > 1
62 pv.push [value, 1]
63 end
64 return pv
65 end
66 end
67
68 class Prime
69 include Enumerable
70
71 def initialize
72 @seed = 1
73 @primes = []
74 @counts = []
75 end
76
77 def succ
78 i = -1
79 size = @primes.size
80 while i < size
81 if i == -1
82 @seed += 1
83 i += 1
84 else
85 while @seed > @counts[i]
86 @counts[i] += @primes[i]
87 end
88 if @seed != @counts[i]
89 i += 1
90 else
91 i = -1
92 end
93 end
94 end
95 @primes.push @seed
96 @counts.push @seed + @seed
97 return @seed
98 end
99 alias next succ
100
101 def each
102 loop do
103 yield succ
104 end
105 end
106 end
107
108 class Fixnum
109 alias divmod! divmod
110 alias / rdiv
111 def divmod(other)
112 a = self.div(other)
113 b = self % other
114 return a,b
115 end
116 end
117
118 class Bignum
119 alias divmod! divmod
120 alias / rdiv
121 end
122
123 class Rational
124 Unify = true
125
126 def inspect
127 format "%s/%s", @numerator.inspect, @denominator.inspect
128 end
129
130 alias power! **
131
132 def ** (other)
133 if other.kind_of?(Rational)
134 if self < 0
135 return Complex(self, 0) ** other
136 elsif other == 0
137 return Rational(1,1)
138 elsif self == 0
139 return Rational(0,1)
140 elsif self == 1
141 return Rational(1,1)
142 end
143
144 npd = @numerator.prime_division
145 dpd = @denominator.prime_division
146 if other < 0
147 other = -other
148 npd, dpd = dpd, npd
149 end
150
151 for elm in npd
152 elm[1] = elm[1] * other
153 if !elm[1].kind_of?(Integer) and elm[1].denominator != 1
154 return Float(self) ** other
155 end
156 elm[1] = elm[1].to_i
157 end
158
159 for elm in dpd
160 elm[1] = elm[1] * other
161 if !elm[1].kind_of?(Integer) and elm[1].denominator != 1
162 return Float(self) ** other
163 end
164 elm[1] = elm[1].to_i
165 end
166
167 num = Integer.from_prime_division(npd)
168 den = Integer.from_prime_division(dpd)
169
170 Rational(num,den)
171
172 elsif other.kind_of?(Integer)
173 if other > 0
174 num = @numerator ** other
175 den = @denominator ** other
176 elsif other < 0
177 num = @denominator ** -other
178 den = @numerator ** -other
179 elsif other == 0
180 num = 1
181 den = 1
182 end
183 Rational.new!(num, den)
184 elsif other.kind_of?(Float)
185 Float(self) ** other
186 else
187 x , y = other.coerce(self)
188 x ** y
189 end
190 end
191
192 def power2(other)
193 if other.kind_of?(Rational)
194 if self < 0
195 return Complex(self, 0) ** other
196 elsif other == 0
197 return Rational(1,1)
198 elsif self == 0
199 return Rational(0,1)
200 elsif self == 1
201 return Rational(1,1)
202 end
203
204 dem = nil
205 x = self.denominator.to_f.to_i
206 neard = self.denominator.to_f ** (1.0/other.denominator.to_f)
207 loop do
208 if (neard**other.denominator == self.denominator)
209 dem = neaed
210 break
211 end
212 end
213 nearn = self.numerator.to_f ** (1.0/other.denominator.to_f)
214 Rational(num,den)
215
216 elsif other.kind_of?(Integer)
217 if other > 0
218 num = @numerator ** other
219 den = @denominator ** other
220 elsif other < 0
221 num = @denominator ** -other
222 den = @numerator ** -other
223 elsif other == 0
224 num = 1
225 den = 1
226 end
227 Rational.new!(num, den)
228 elsif other.kind_of?(Float)
229 Float(self) ** other
230 else
231 x , y = other.coerce(self)
232 x ** y
233 end
234 end
235 end
236
237 module Math
238 def sqrt(a)
239 if a.kind_of?(Complex)
240 abs = sqrt(a.real*a.real + a.image*a.image)
241 # if not abs.kind_of?(Rational)
242 # return a**Rational(1,2)
243 # end
244 x = sqrt((a.real + abs)/Rational(2))
245 y = sqrt((-a.real + abs)/Rational(2))
246 # if !(x.kind_of?(Rational) and y.kind_of?(Rational))
247 # return a**Rational(1,2)
248 # end
249 if a.image >= 0
250 Complex(x, y)
251 else
252 Complex(x, -y)
253 end
254 elsif a >= 0
255 rsqrt(a)
256 else
257 Complex(0,rsqrt(-a))
258 end
259 end
260
261 def rsqrt(a)
262 if a.kind_of?(Float)
263 sqrt!(a)
264 elsif a.kind_of?(Rational)
265 rsqrt(a.numerator)/rsqrt(a.denominator)
266 else
267 src = a
268 max = 2 ** 32
269 byte_a = [src & 0xffffffff]
270 # ruby's bug
271 while (src >= max) and (src >>= 32)
272 byte_a.unshift src & 0xffffffff
273 end
274
275 answer = 0
276 main = 0
277 side = 0
278 for elm in byte_a
279 main = (main << 32) + elm
280 side <<= 16
281 if answer != 0
282 if main * 4 < side * side
283 applo = main.div(side)
284 else
285 applo = ((sqrt!(side * side + 4 * main) - side)/2.0).to_i + 1
286 end
287 else
288 applo = sqrt!(main).to_i + 1
289 end
290
291 while (x = (side + applo) * applo) > main
292 applo -= 1
293 end
294 main -= x
295 answer = (answer << 16) + applo
296 side += applo * 2
297 end
298 if main == 0
299 answer
300 else
301 sqrt!(a)
302 end
303 end
304 end
305
306 module_function :sqrt
307 module_function :rsqrt
308 end
309
310 class Complex
311 Unify = true
312 end
313