
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 


