|
1 |
##
|
|
2 |
# = mathn
|
1 |
3 |
#
|
2 |
|
# mathn.rb -
|
3 |
|
# $Release Version: 0.5 $
|
4 |
|
# $Revision: 1.1.1.1.4.1 $
|
5 |
|
# by Keiju ISHITSUKA(SHL Japan Inc.)
|
|
4 |
# mathn is a library for changing the way Ruby does math.
|
6 |
5 |
#
|
7 |
|
# --
|
|
6 |
# == Usage
|
|
7 |
#
|
|
8 |
# To start using this library, simply:
|
|
9 |
#
|
|
10 |
# require "mathn"
|
|
11 |
#
|
|
12 |
# This will change the way division works for Fixnums, specifically
|
|
13 |
#
|
|
14 |
# 3 / 2
|
8 |
15 |
#
|
|
16 |
# will return (3/2) instead of the usual 1.
|
9 |
17 |
#
|
|
18 |
# == Copyright
|
10 |
19 |
#
|
|
20 |
# Author: Keiju ISHITSUKA(SHL Japan Inc.)
|
|
21 |
#
|
|
22 |
# --
|
|
23 |
# $Release Version: 0.5 $
|
|
24 |
# $Revision: 1.1.1.1.4.1 $
|
11 |
25 |
|
12 |
26 |
require "cmath.rb"
|
13 |
27 |
require "matrix.rb"
|
... | ... | |
27 |
41 |
|
28 |
42 |
alias power! ** unless method_defined? :power!
|
29 |
43 |
|
|
44 |
##
|
|
45 |
# exponentiate by +other+
|
30 |
46 |
def ** (other)
|
31 |
47 |
if self < 0 && other.round != other
|
32 |
48 |
Complex(self, 0.0) ** other
|
... | ... | |
43 |
59 |
|
44 |
60 |
alias power! ** unless method_defined? :power!
|
45 |
61 |
|
|
62 |
##
|
|
63 |
# exponentiate by +other+
|
46 |
64 |
def ** (other)
|
47 |
65 |
if self < 0 && other.round != other
|
48 |
66 |
Complex(self, 0.0) ** other
|
... | ... | |
55 |
73 |
|
56 |
74 |
class Rational
|
57 |
75 |
remove_method :**
|
|
76 |
|
|
77 |
##
|
|
78 |
# exponentiate by +other+
|
58 |
79 |
def ** (other)
|
59 |
80 |
if other.kind_of?(Rational)
|
60 |
81 |
other2 = other
|
61 |
82 |
if self < 0
|
62 |
|
return Complex(self, 0.0) ** other
|
|
83 |
return Complex(self, 0.0) ** other
|
63 |
84 |
elsif other == 0
|
64 |
|
return Rational(1,1)
|
|
85 |
return Rational(1,1)
|
65 |
86 |
elsif self == 0
|
66 |
|
return Rational(0,1)
|
|
87 |
return Rational(0,1)
|
67 |
88 |
elsif self == 1
|
68 |
|
return Rational(1,1)
|
|
89 |
return Rational(1,1)
|
69 |
90 |
end
|
70 |
91 |
|
71 |
92 |
npd = numerator.prime_division
|
72 |
93 |
dpd = denominator.prime_division
|
73 |
94 |
if other < 0
|
74 |
|
other = -other
|
75 |
|
npd, dpd = dpd, npd
|
|
95 |
other = -other
|
|
96 |
npd, dpd = dpd, npd
|
76 |
97 |
end
|
77 |
98 |
|
78 |
99 |
for elm in npd
|
79 |
|
elm[1] = elm[1] * other
|
80 |
|
if !elm[1].kind_of?(Integer) and elm[1].denominator != 1
|
81 |
|
return Float(self) ** other2
|
82 |
|
end
|
83 |
|
elm[1] = elm[1].to_i
|
|
100 |
elm[1] = elm[1] * other
|
|
101 |
if !elm[1].kind_of?(Integer) and elm[1].denominator != 1
|
|
102 |
return Float(self) ** other2
|
|
103 |
end
|
|
104 |
elm[1] = elm[1].to_i
|
84 |
105 |
end
|
85 |
106 |
|
86 |
107 |
for elm in dpd
|
87 |
|
elm[1] = elm[1] * other
|
88 |
|
if !elm[1].kind_of?(Integer) and elm[1].denominator != 1
|
89 |
|
return Float(self) ** other2
|
90 |
|
end
|
91 |
|
elm[1] = elm[1].to_i
|
|
108 |
elm[1] = elm[1] * other
|
|
109 |
if !elm[1].kind_of?(Integer) and elm[1].denominator != 1
|
|
110 |
return Float(self) ** other2
|
|
111 |
end
|
|
112 |
elm[1] = elm[1].to_i
|
92 |
113 |
end
|
93 |
114 |
|
94 |
115 |
num = Integer.from_prime_division(npd)
|
... | ... | |
98 |
119 |
|
99 |
120 |
elsif other.kind_of?(Integer)
|
100 |
121 |
if other > 0
|
101 |
|
num = numerator ** other
|
102 |
|
den = denominator ** other
|
|
122 |
num = numerator ** other
|
|
123 |
den = denominator ** other
|
103 |
124 |
elsif other < 0
|
104 |
|
num = denominator ** -other
|
105 |
|
den = numerator ** -other
|
|
125 |
num = denominator ** -other
|
|
126 |
den = numerator ** -other
|
106 |
127 |
elsif other == 0
|
107 |
|
num = 1
|
108 |
|
den = 1
|
|
128 |
num = 1
|
|
129 |
den = 1
|
109 |
130 |
end
|
110 |
131 |
Rational(num, den)
|
111 |
132 |
elsif other.kind_of?(Float)
|
... | ... | |
119 |
140 |
|
120 |
141 |
module Math
|
121 |
142 |
remove_method(:sqrt)
|
|
143 |
|
|
144 |
##
|
|
145 |
# compute the square root of +a+
|
122 |
146 |
def sqrt(a)
|
123 |
147 |
if a.kind_of?(Complex)
|
124 |
148 |
abs = sqrt(a.real*a.real + a.imag*a.imag)
|
125 |
149 |
# if not abs.kind_of?(Rational)
|
126 |
|
# return a**Rational(1,2)
|
|
150 |
# return a**Rational(1,2)
|
127 |
151 |
# end
|
128 |
152 |
x = sqrt((a.real + abs)/Rational(2))
|
129 |
153 |
y = sqrt((-a.real + abs)/Rational(2))
|
130 |
154 |
# if !(x.kind_of?(Rational) and y.kind_of?(Rational))
|
131 |
|
# return a**Rational(1,2)
|
|
155 |
# return a**Rational(1,2)
|
132 |
156 |
# end
|
133 |
157 |
if a.imag >= 0
|
134 |
|
Complex(x, y)
|
|
158 |
Complex(x, y)
|
135 |
159 |
else
|
136 |
|
Complex(x, -y)
|
|
160 |
Complex(x, -y)
|
137 |
161 |
end
|
138 |
162 |
elsif a.respond_to?(:nan?) and a.nan?
|
139 |
163 |
a
|
... | ... | |
144 |
168 |
end
|
145 |
169 |
end
|
146 |
170 |
|
147 |
|
def rsqrt(a)
|
|
171 |
def rsqrt(a) # :nodoc:
|
148 |
172 |
if a.kind_of?(Float)
|
149 |
173 |
sqrt!(a)
|
150 |
174 |
elsif a.kind_of?(Rational)
|
... | ... | |
155 |
179 |
byte_a = [src & 0xffffffff]
|
156 |
180 |
# ruby's bug
|
157 |
181 |
while (src >= max) and (src >>= 32)
|
158 |
|
byte_a.unshift src & 0xffffffff
|
|
182 |
byte_a.unshift src & 0xffffffff
|
159 |
183 |
end
|
160 |
184 |
|
161 |
185 |
answer = 0
|
162 |
186 |
main = 0
|
163 |
187 |
side = 0
|
164 |
188 |
for elm in byte_a
|
165 |
|
main = (main << 32) + elm
|
166 |
|
side <<= 16
|
167 |
|
if answer != 0
|
168 |
|
if main * 4 < side * side
|
169 |
|
applo = main.div(side)
|
170 |
|
else
|
171 |
|
applo = ((sqrt!(side * side + 4 * main) - side)/2.0).to_i + 1
|
172 |
|
end
|
173 |
|
else
|
174 |
|
applo = sqrt!(main).to_i + 1
|
175 |
|
end
|
176 |
|
|
177 |
|
while (x = (side + applo) * applo) > main
|
178 |
|
applo -= 1
|
179 |
|
end
|
180 |
|
main -= x
|
181 |
|
answer = (answer << 16) + applo
|
182 |
|
side += applo * 2
|
|
189 |
main = (main << 32) + elm
|
|
190 |
side <<= 16
|
|
191 |
if answer != 0
|
|
192 |
if main * 4 < side * side
|
|
193 |
applo = main.div(side)
|
|
194 |
else
|
|
195 |
applo = ((sqrt!(side * side + 4 * main) - side)/2.0).to_i + 1
|
|
196 |
end
|
|
197 |
else
|
|
198 |
applo = sqrt!(main).to_i + 1
|
|
199 |
end
|
|
200 |
|
|
201 |
while (x = (side + applo) * applo) > main
|
|
202 |
applo -= 1
|
|
203 |
end
|
|
204 |
main -= x
|
|
205 |
answer = (answer << 16) + applo
|
|
206 |
side += applo * 2
|
183 |
207 |
end
|
184 |
208 |
if main == 0
|
185 |
|
answer
|
|
209 |
answer
|
186 |
210 |
else
|
187 |
|
sqrt!(a)
|
|
211 |
sqrt!(a)
|
188 |
212 |
end
|
189 |
213 |
end
|
190 |
214 |
end
|
... | ... | |
199 |
223 |
class Float
|
200 |
224 |
alias power! **
|
201 |
225 |
|
|
226 |
##
|
|
227 |
# exponentiate by +other+
|
202 |
228 |
def ** (other)
|
203 |
229 |
if self < 0 && other.round != other
|
204 |
230 |
Complex(self, 0.0) ** other
|
205 |
|
-
|