[ruby-dev:40083] BigMath の sin(BigDecimal("1E50"), 1) が無限ループ

  豊福です。

  BigDecimal BigMath の sin(BigDecimal("1E50"), 1) が PI に
要求する精度の不足で無限ループになります。
ささやかな最適化も入れました。
テストプログラムも付けておきます。

--- /ext/bigdecimal/lib/bigdecimal/math.rb      (revision 25184)
+++ /ext/bigdecimal/lib/bigdecimal/math.rb      (working copy)
@@ -54,5 +54,6 @@
     two  = BigDecimal("2")
     x = -x if neg = x < 0
-    if x > (twopi = two * BigMath.PI(prec))
+    if x >= 6.3
+      twopi = two * PI(x.exponent+prec+1)
       if x > 30
         x %= twopi
@@ -62,6 +63,5 @@
     end
     x1   = x
-    x2   = x.mult(x,n)
-    sign = 1
+    x2   = -x.mult(x,n)
     y    = x
     d    = y
@@ -70,9 +70,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 +89,6 @@
     two  = BigDecimal("2")
     x = -x if x < 0
-    if x > (twopi = two * BigMath.PI(prec))
+    if x >= 6.3
+      twopi = two * PI(x.exponent+prec+1)
       if x > 30
         x %= twopi
@@ -98,6 +98,5 @@
     end
     x1 = one
-    x2 = x.mult(x,n)
-    sign = 1
+    x2 = -x.mult(x,n)
     y = one
     d = y
@@ -106,9 +105,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

======== テスト ========

require 'bigdecimal'
require 'bigdecimal/math'
include BigMath

# sin test
[
  ["1", ".8414709848078965066525023216302989996225630607983710656727517099919104043912396689486397435430526958"],
  ["-1", "-.8414709848078965066525023216302989996225630607983710656727517099919104043912396689486397435430526958"],
  ["1.5707963", ".9999999999999996410167575823529650844152981034590074856559180469605393157919430532545927996479289687"], # PI/2
  ["6.2", "-.08308940281749657800057928909836718528109967293845540803530361987239105120353659681491938223261595192"],
  ["10", "-.5440211108893698134047476618513772816836430129162238915741840126167572096404934257070756738949832161"],

  ["100", "-.5063656411097587936565576104597854320650327212906573234433924735943579134194766964992366645129273922"],
].each {|sx,sy|
  m = sy.gsub(/^-?0*\.?0*|E.*|\./,'').size
  x = BigDecimal(sx)
  y = BigDecimal(sy)
  n = y == 0 ? 1: y.exponent
  yeps = BigDecimal("1E#{n-m}")
  (1..m).each {|i|
    if (sin(x, i) - y).abs >= BigDecimal("1E#{n-i}") + yeps
      print "sin fail i: #{i}, x: #{x}\n"
      print sin(x, i), "\n"
      print y, "\n"
      break
    end
  }
}

# cos test
[
  ["1", ".5403023058681397174009366074429766037323104206179222276700972553811003947744717645179518560871830893"],
  ["-1", ".5403023058681397174009366074429766037323104206179222276700972553811003947744717645179518560871830893"],
  ["3.1415926", "-.9999999999999985640670303294121180755978657881108015358856694424326400998221703599107746795702773875"], # PI
  ["6.2", ".9965420970232174751394026238692639525278846240987172016105014589865551509745000082886442339885138615"],
  ["10", "-.8390715290764524522588639478240648345199301651331685468359537310487925868662707684009337127604221389"],
  ["100", ".8623188722876839341019385139508425355100840085355108292801621126927210880509266241030951056842772850"],
].each {|sx,sy|
  m = sy.gsub(/^-?0*\.?0*|E.*|\./,'').size
  x = BigDecimal(sx)
  y = BigDecimal(sy)
  n = y == 0 ? 1: y.exponent
  yeps = BigDecimal("1E#{n-m}")
  (1..m).each {|i|
    if (cos(x, i) - y).abs >= BigDecimal("1E#{n-i}") + yeps
      print "cos fail i: #{i}, x: #{x}\n"
      print cos(x, i), "\n"
      print y, "\n"
      break
    end
  }
}
---
May Haskell be with you.