[ruby-core:25600] [Bug #2101] Better numerical approximation for bigdecimal/math.rb

Bug #2101: Better numerical approximation for bigdecimal/math.rb
http://redmine.ruby-lang.org/issues/show/2101

Author: Makoto Yamashita
Status: Open, Priority: Normal
ruby -v: ruby 1.8.7 (2008-08-11 patchlevel 72) [universal-darwin10.0]

The functions in the BigMath module tend to be slow and/or give erroneous result when you supply "large" numbers as their arguments, e.g.
$ ruby -rbigdecimal -rbigdecimal/math -e 'include BigMath;p BigMath.exp(BigDecimal("-50"),10).to_s'
"-0.956441768151696446193043446622932686846145E-4"
(the correct answer would be e**(-50) ~ 1.928749847964E-22; anyway this should remain always positive)

$ ruby -rbigdecimal -rbigdecimal/math -e 'include BigMath;p BigMath.sin(BigDecimal("100"),10).to_s'
"0.285414041009447759758871460612315212089477E18"
(the correct answer would be sin(100) = -0.50636564...)

In addition, it takes too long to calculate things like BigMath.log(BigDecimal("1E-10"),10).

There are two reasons for this:
1) approximation by the Taylor series f(x) = f(0) + f'(0)x + f''(0)x**2/2 + ... is error-prone when |x| >> |f(x)|, since the first a few terms in the series have too large absolute values and the errors coming from them make the final estimate of f(x) very inaccurate (this is visible for exp(x) when x ~ -50, and sin(x), cos(x) for |x| ~ 100.
2) the sum we compute in the current code becomes inadequate when it becomes similar to the harmonic series 1 + 1/2 + 1/3 + ...  The convergence is going to be really slow and again error-prone.  This is the case for log(x) when x ~ 1E-10, atan(x) for x ~ 1 - 1E-5.

The proposed change
1) uses exp(x) = 1/exp(-x) when x is negative (since y -> 1/y does not affect the relative magnitude of the error) and computes sin(x + 2k pi), cos(x + 2k pi) with k satisfying |x + 2k pi| <= 2pi.
2) computes log(x*(10**a)) (and return log(x) = log(x*(10**a)) - a*log(10)) when 10**(-a-1) < x < 10**(-a) < 1.  For the arch tangent, it uses Euler's formula atan(x) = sum 2**(2n)(n!)**2/(2n+1)! x**(2n+1)/(1+x**2)**(n+1).  This has the byproduct of (rightfully!) allowing |x| >= 1 as the argument.


----------------------------------------
http://redmine.ruby-lang.org