@@ -200,52 +200,55 @@ def self.log(x, prec)
200
200
return BigDecimal ::Internal . infinity_computation_result if x . infinite?
201
201
return BigDecimal ( 0 ) if x == 1
202
202
203
- if x > 10 || x < 0.1
204
- log10 = log ( BigDecimal ( 10 ) , prec )
205
- exponent = x . exponent
206
- x = x * BigDecimal ( "1e#{ -x . exponent } " )
207
- if x < 0.3
208
- x *= 10
209
- exponent -= 1
203
+ BigDecimal . save_limit do
204
+ BigDecimal . limit ( 0 )
205
+ if x > 10 || x < 0.1
206
+ log10 = log ( BigDecimal ( 10 ) , prec )
207
+ exponent = x . exponent
208
+ x *= BigDecimal ( "1e#{ -x . exponent } " )
209
+ if x < 0.3
210
+ x *= 10
211
+ exponent -= 1
212
+ end
213
+ return log10 * exponent + log ( x , prec )
210
214
end
211
- return log10 * exponent + log ( x , prec )
212
- end
213
215
214
- x_minus_one_exponent = ( x - 1 ) . exponent
215
- prec += BigDecimal . double_fig
216
+ x_minus_one_exponent = ( x - 1 ) . exponent
217
+ prec += BigDecimal . double_fig
216
218
217
- # log(x) = log(sqrt(sqrt(sqrt(sqrt(x))))) * 2**sqrt_steps
218
- sqrt_steps = [ 2 * Integer . sqrt ( prec ) + 3 * x_minus_one_exponent , 0 ] . max
219
+ # log(x) = log(sqrt(sqrt(sqrt(sqrt(x))))) * 2**sqrt_steps
220
+ sqrt_steps = [ 2 * Integer . sqrt ( prec ) + 3 * x_minus_one_exponent , 0 ] . max
219
221
220
- # Reduce sqrt_step until sqrt gets fast
221
- # https://github.com/ruby/bigdecimal/pull/323
222
- # https://github.com/ruby/bigdecimal/pull/343
223
- sqrt_steps /= 10
222
+ # Reduce sqrt_step until sqrt gets fast
223
+ # https://github.com/ruby/bigdecimal/pull/323
224
+ # https://github.com/ruby/bigdecimal/pull/343
225
+ sqrt_steps /= 10
224
226
225
- lg2 = 0.3010299956639812
226
- prec2 = prec + [ -x_minus_one_exponent , 0 ] . max + ( sqrt_steps * lg2 ) . ceil
227
+ lg2 = 0.3010299956639812
228
+ prec2 = prec + [ -x_minus_one_exponent , 0 ] . max + ( sqrt_steps * lg2 ) . ceil
227
229
228
- sqrt_steps . times do
229
- x = x . sqrt ( prec2 )
230
+ sqrt_steps . times do
231
+ x = x . sqrt ( prec2 )
230
232
231
- # Workaround for https://github.com/ruby/bigdecimal/issues/354
232
- x = x . mult ( 1 , prec2 + BigDecimal . double_fig )
233
- end
233
+ # Workaround for https://github.com/ruby/bigdecimal/issues/354
234
+ x = x . mult ( 1 , prec2 + BigDecimal . double_fig )
235
+ end
234
236
235
- # Taylor series for log(x) around 1
236
- # log(x) = -log((1 + X) / (1 - X)) where X = (x - 1) / (x + 1)
237
- # log(x) = 2 * (X + X**3 / 3 + X**5 / 5 + X**7 / 7 + ...)
238
- x = ( x - 1 ) . div ( x + 1 , prec2 )
239
- y = x
240
- x2 = x . mult ( x , prec )
241
- 1 . step do |i |
242
- n = prec + x . exponent - y . exponent + x2 . exponent
243
- break if n <= 0 || x . zero?
244
- x = x . mult ( x2 . round ( n - x2 . exponent ) , n )
245
- y = y . add ( x . div ( 2 * i + 1 , n ) , prec )
246
- end
237
+ # Taylor series for log(x) around 1
238
+ # log(x) = -log((1 + X) / (1 - X)) where X = (x - 1) / (x + 1)
239
+ # log(x) = 2 * (X + X**3 / 3 + X**5 / 5 + X**7 / 7 + ...)
240
+ x = ( x - 1 ) . div ( x + 1 , prec2 )
241
+ y = x
242
+ x2 = x . mult ( x , prec )
243
+ 1 . step do |i |
244
+ n = prec + x . exponent - y . exponent + x2 . exponent
245
+ break if n <= 0 || x . zero?
246
+ x = x . mult ( x2 . round ( n - x2 . exponent ) , n )
247
+ y = y . add ( x . div ( 2 * i + 1 , n ) , prec )
248
+ end
247
249
248
- y . mult ( 2 ** ( sqrt_steps + 1 ) , prec )
250
+ y . mult ( 2 ** ( sqrt_steps + 1 ) , prec )
251
+ end
249
252
end
250
253
251
254
# call-seq:
@@ -266,30 +269,33 @@ def self.exp(x, prec)
266
269
return BigDecimal ( 1 ) if x . zero?
267
270
return BigDecimal ( 1 ) . div ( exp ( -x , prec ) , prec ) if x < 0
268
271
269
- # exp(x * 10**cnt) = exp(x)**(10**cnt)
270
- cnt = x > 1 ? x . exponent : 0
271
- prec2 = prec + BigDecimal . double_fig + cnt
272
- x *= BigDecimal ( "1e-#{ cnt } " )
273
- xn = BigDecimal ( 1 )
274
- y = BigDecimal ( 1 )
275
-
276
- # Taylor series for exp(x) around 0
277
- 1 . step do |i |
278
- n = prec2 + xn . exponent
279
- break if n <= 0 || xn . zero?
280
- x = x . mult ( 1 , n )
281
- xn = xn . mult ( x , n ) . div ( i , n )
282
- y = y . add ( xn , prec2 )
283
- end
272
+ BigDecimal . save_limit do
273
+ BigDecimal . limit ( 0 )
274
+ # exp(x * 10**cnt) = exp(x)**(10**cnt)
275
+ cnt = x > 1 ? x . exponent : 0
276
+ prec2 = prec + BigDecimal . double_fig + cnt
277
+ x *= BigDecimal ( "1e-#{ cnt } " )
278
+ xn = BigDecimal ( 1 )
279
+ y = BigDecimal ( 1 )
280
+
281
+ # Taylor series for exp(x) around 0
282
+ 1 . step do |i |
283
+ n = prec2 + xn . exponent
284
+ break if n <= 0 || xn . zero?
285
+ x = x . mult ( 1 , n )
286
+ xn = xn . mult ( x , n ) . div ( i , n )
287
+ y = y . add ( xn , prec2 )
288
+ end
284
289
285
- # calculate exp(x * 10**cnt) from exp(x)
286
- # exp(x * 10**k) = exp(x * 10**(k - 1)) ** 10
287
- cnt . times do
288
- y2 = y . mult ( y , prec2 )
289
- y5 = y2 . mult ( y2 , prec2 ) . mult ( y , prec2 )
290
- y = y5 . mult ( y5 , prec2 )
291
- end
290
+ # calculate exp(x * 10**cnt) from exp(x)
291
+ # exp(x * 10**k) = exp(x * 10**(k - 1)) ** 10
292
+ cnt . times do
293
+ y2 = y . mult ( y , prec2 )
294
+ y5 = y2 . mult ( y2 , prec2 ) . mult ( y , prec2 )
295
+ y = y5 . mult ( y5 , prec2 )
296
+ end
292
297
293
- y . mult ( 1 , prec )
298
+ y . mult ( 1 , prec )
299
+ end
294
300
end
295
301
end
0 commit comments