|
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 |
|
-
|