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 "complex.rb"
|
---|
14 | require "rational.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 | raise ZeroDivisionError if self == 0
|
---|
47 | ps = Prime.new
|
---|
48 | value = self
|
---|
49 | pv = []
|
---|
50 | for prime in ps
|
---|
51 | count = 0
|
---|
52 | while (value1, mod = value.divmod(prime)
|
---|
53 | mod) == 0
|
---|
54 | value = value1
|
---|
55 | count += 1
|
---|
56 | end
|
---|
57 | if count != 0
|
---|
58 | pv.push [prime, count]
|
---|
59 | end
|
---|
60 | break if prime * prime >= value
|
---|
61 | end
|
---|
62 | if value > 1
|
---|
63 | pv.push [value, 1]
|
---|
64 | end
|
---|
65 | return pv
|
---|
66 | end
|
---|
67 | end
|
---|
68 |
|
---|
69 | class Prime
|
---|
70 | include Enumerable
|
---|
71 |
|
---|
72 | def initialize
|
---|
73 | @seed = 1
|
---|
74 | @primes = []
|
---|
75 | @counts = []
|
---|
76 | end
|
---|
77 |
|
---|
78 | def succ
|
---|
79 | i = -1
|
---|
80 | size = @primes.size
|
---|
81 | while i < size
|
---|
82 | if i == -1
|
---|
83 | @seed += 1
|
---|
84 | i += 1
|
---|
85 | else
|
---|
86 | while @seed > @counts[i]
|
---|
87 | @counts[i] += @primes[i]
|
---|
88 | end
|
---|
89 | if @seed != @counts[i]
|
---|
90 | i += 1
|
---|
91 | else
|
---|
92 | i = -1
|
---|
93 | end
|
---|
94 | end
|
---|
95 | end
|
---|
96 | @primes.push @seed
|
---|
97 | @counts.push @seed + @seed
|
---|
98 | return @seed
|
---|
99 | end
|
---|
100 | alias next succ
|
---|
101 |
|
---|
102 | def each
|
---|
103 | loop do
|
---|
104 | yield succ
|
---|
105 | end
|
---|
106 | end
|
---|
107 | end
|
---|
108 |
|
---|
109 | class Fixnum
|
---|
110 | alias / quo
|
---|
111 | end
|
---|
112 |
|
---|
113 | class Bignum
|
---|
114 | alias / quo
|
---|
115 | end
|
---|
116 |
|
---|
117 | class Rational
|
---|
118 | Unify = true
|
---|
119 |
|
---|
120 | def inspect
|
---|
121 | format "%s/%s", numerator.inspect, denominator.inspect
|
---|
122 | end
|
---|
123 |
|
---|
124 | alias power! **
|
---|
125 |
|
---|
126 | def ** (other)
|
---|
127 | if other.kind_of?(Rational)
|
---|
128 | other2 = other
|
---|
129 | if self < 0
|
---|
130 | return Complex.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.new!(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 |
|
---|
187 | def power2(other)
|
---|
188 | if other.kind_of?(Rational)
|
---|
189 | if self < 0
|
---|
190 | return Complex(self, 0) ** other
|
---|
191 | elsif other == 0
|
---|
192 | return Rational(1,1)
|
---|
193 | elsif self == 0
|
---|
194 | return Rational(0,1)
|
---|
195 | elsif self == 1
|
---|
196 | return Rational(1,1)
|
---|
197 | end
|
---|
198 |
|
---|
199 | dem = nil
|
---|
200 | x = self.denominator.to_f.to_i
|
---|
201 | neard = self.denominator.to_f ** (1.0/other.denominator.to_f)
|
---|
202 | loop do
|
---|
203 | if (neard**other.denominator == self.denominator)
|
---|
204 | dem = neaed
|
---|
205 | break
|
---|
206 | end
|
---|
207 | end
|
---|
208 | nearn = self.numerator.to_f ** (1.0/other.denominator.to_f)
|
---|
209 | Rational(num,den)
|
---|
210 |
|
---|
211 | elsif other.kind_of?(Integer)
|
---|
212 | if other > 0
|
---|
213 | num = numerator ** other
|
---|
214 | den = denominator ** other
|
---|
215 | elsif other < 0
|
---|
216 | num = denominator ** -other
|
---|
217 | den = numerator ** -other
|
---|
218 | elsif other == 0
|
---|
219 | num = 1
|
---|
220 | den = 1
|
---|
221 | end
|
---|
222 | Rational.new!(num, den)
|
---|
223 | elsif other.kind_of?(Float)
|
---|
224 | Float(self) ** other
|
---|
225 | else
|
---|
226 | x , y = other.coerce(self)
|
---|
227 | x ** y
|
---|
228 | end
|
---|
229 | end
|
---|
230 | end
|
---|
231 |
|
---|
232 | module Math
|
---|
233 | def sqrt(a)
|
---|
234 | if a.kind_of?(Complex)
|
---|
235 | abs = sqrt(a.real*a.real + a.image*a.image)
|
---|
236 | # if not abs.kind_of?(Rational)
|
---|
237 | # return a**Rational(1,2)
|
---|
238 | # end
|
---|
239 | x = sqrt((a.real + abs)/Rational(2))
|
---|
240 | y = sqrt((-a.real + abs)/Rational(2))
|
---|
241 | # if !(x.kind_of?(Rational) and y.kind_of?(Rational))
|
---|
242 | # return a**Rational(1,2)
|
---|
243 | # end
|
---|
244 | if a.image >= 0
|
---|
245 | Complex(x, y)
|
---|
246 | else
|
---|
247 | Complex(x, -y)
|
---|
248 | end
|
---|
249 | elsif a >= 0
|
---|
250 | rsqrt(a)
|
---|
251 | else
|
---|
252 | Complex(0,rsqrt(-a))
|
---|
253 | end
|
---|
254 | end
|
---|
255 |
|
---|
256 | def rsqrt(a)
|
---|
257 | if a.kind_of?(Float)
|
---|
258 | sqrt!(a)
|
---|
259 | elsif a.kind_of?(Rational)
|
---|
260 | rsqrt(a.numerator)/rsqrt(a.denominator)
|
---|
261 | else
|
---|
262 | src = a
|
---|
263 | max = 2 ** 32
|
---|
264 | byte_a = [src & 0xffffffff]
|
---|
265 | # ruby's bug
|
---|
266 | while (src >= max) and (src >>= 32)
|
---|
267 | byte_a.unshift src & 0xffffffff
|
---|
268 | end
|
---|
269 |
|
---|
270 | answer = 0
|
---|
271 | main = 0
|
---|
272 | side = 0
|
---|
273 | for elm in byte_a
|
---|
274 | main = (main << 32) + elm
|
---|
275 | side <<= 16
|
---|
276 | if answer != 0
|
---|
277 | if main * 4 < side * side
|
---|
278 | applo = main.div(side)
|
---|
279 | else
|
---|
280 | applo = ((sqrt!(side * side + 4 * main) - side)/2.0).to_i + 1
|
---|
281 | end
|
---|
282 | else
|
---|
283 | applo = sqrt!(main).to_i + 1
|
---|
284 | end
|
---|
285 |
|
---|
286 | while (x = (side + applo) * applo) > main
|
---|
287 | applo -= 1
|
---|
288 | end
|
---|
289 | main -= x
|
---|
290 | answer = (answer << 16) + applo
|
---|
291 | side += applo * 2
|
---|
292 | end
|
---|
293 | if main == 0
|
---|
294 | answer
|
---|
295 | else
|
---|
296 | sqrt!(a)
|
---|
297 | end
|
---|
298 | end
|
---|
299 | end
|
---|
300 |
|
---|
301 | module_function :sqrt
|
---|
302 | module_function :rsqrt
|
---|
303 | end
|
---|
304 |
|
---|
305 | class Complex
|
---|
306 | Unify = true
|
---|
307 | end
|
---|
308 |
|
---|