From: TOYOFUKU Chikanobu <nobu_toyofuku@...>
Date: 2010-01-15T19:51:02+09:00
Subject: [ruby-dev:40104] Re: BigMath の sin(BigDecimal("1E50"), 1) が無限ループ

  豊福です。

  誤差の評価が不十分でした。
Σ(-1)^n*x^（2*n+1)/（2*n+1)! の計算で誤差は
Σx^（2*n+1)/（2*n+1)! = sinh(x) 倍になると思います。

sinh(2π) も sinh(π) も値が 10 を超えるので精度を保つために
x を π/2 程度以下にしました。cos も同様です。

--- /ext/bigdecimal/lib/bigdecimal/math.rb      (revision 25184)
+++ /ext/bigdecimal/lib/bigdecimal/math.rb      (working copy)
@@ -54,5 +54,7 @@
     two  = BigDecimal("2")
     x = -x if neg = x < 0
-    if x > (twopi = two * BigMath.PI(prec))
+    if x >= 6.3
+      pi = PI(x.exponent+prec+1)
+      twopi = two * pi
       if x > 30
         x %= twopi
@@ -61,7 +63,17 @@
       end
     end
+    if x >= 1.6
+      pi ||= PI(x.exponent+prec+1)
+      while x >= 1.6
+        x -= pi
+        neg = !neg
+      end
+      if x < 0
+        x = -x
+        neg = !neg
+      end
+    end
     x1   = x
-    x2   = x.mult(x,n)
-    sign = 1
+    x2   = -x.mult(x,n)
     y    = x
     d    = y
@@ -70,9 +82,8 @@
     while d.nonzero? && ((m = n - (y.exponent - d.exponent).abs) > 0)
       m = BigDecimal.double_fig if m < BigDecimal.double_fig
-      sign = -sign
       x1  = x2.mult(x1,n)
       i  += two
       z  *= (i-one) * i
-      d   = sign * x1.div(z,m)
+      d   = x1.div(z,m)
       y  += d
     end
@@ -90,5 +101,7 @@
     two  = BigDecimal("2")
     x = -x if x < 0
-    if x > (twopi = two * BigMath.PI(prec))
+    if x >= 6.3
+      pi = PI(x.exponent+prec+1)
+      twopi = two * pi
       if x > 30
         x %= twopi
@@ -97,7 +110,15 @@
       end
     end
+    neg = false
+    if x >= 1.6
+      pi ||= PI(x.exponent+prec+1)
+      while x >= 1.6
+        x -= pi
+        neg = !neg
+      end
+      x = -x if x < 0
+    end
     x1 = one
-    x2 = x.mult(x,n)
-    sign = 1
+    x2 = -x.mult(x,n)
     y = one
     d = y
@@ -106,12 +127,11 @@
     while d.nonzero? && ((m = n - (y.exponent - d.exponent).abs) > 0)
       m = BigDecimal.double_fig if m < BigDecimal.double_fig
-      sign = -sign
       x1  = x2.mult(x1,n)
       i  += two
       z  *= (i-one) * i
-      d   = sign * x1.div(z,m)
+      d   = x1.div(z,m)
       y  += d
     end
-    y
+    neg ? -y : y
   end

---
May Haskell be with you.