Almost endless loop of BigMath::atan(x) when x.abs >= 1

e$B6b0fe(B e$B?N90$H?=$7$^$9!#e(B

e$B0J2<$N$h$&$K$9$k$H!"7W;;NL$NA}Bg$K$h$j!"e(B
e$B$[$\L58B%k!<%W$KF~$j$^$9!#e(B

$ ./ruby -v
ruby 1.9.2dev (2009-09-23 trunk 25052) [i686-linux]
$ ./ruby -e ’
require “bigdecimal”
require “bigdecimal/math”
N = 20
p BigMath::atan(BigDecimal(“1.0”), N)

e$B$3$l$O!"e(BBigMath::atan(x)e$B$Ne(Bx.abse$B$,e(B1e$B0J>e$N;~$K!"H/@8$9$kLdBj$N$h$&$G$9!#e(B
e$B$3$l$K$h$j!"e(Bmake test-alle$B$bESCf$GL58B%k!<%W$KF~$j$^$9!#e(B

x.abs >= 1e$B$N$H$-!"e(B
atan(x) = pi/2 - atan(1/x)
e$B$rMQ$$$F!"0J2<$N$h$&$K=q$-D>$7$^$7$?!#e(B

e$B$^$?!">e$N<0$G$b!"e(Bxe$B$,e(B1e$B$KHs>o$K6a$$?tCM$N>l9g$O!"e(B
e$B6a;wCM$r5a$a$k$?$a$KG|Bg$J7W;;NL$,I,MW$J$?$a!"e(B
x.abs.round(prec) == 1e$B$N$H$-!"e(B
atan(x) = atan((-1+sqrt(1+x**2, prec))/x)
e$B$H$7$^$7$?!#e(B

e$B$h$m$7$/$*4j$$$$$?$7$^$9!#e(B

— /ext/bigdecimal/lib/bigdecimal/math.rb (revision 25052)
+++ /ext/bigdecimal/lib/bigdecimal/math.rb (working copy)
@@ -122,11 +122,18 @@
raise ArgumentError, “Zero or negative precision for atan” if prec
<= 0
return BigDecimal(“NaN”) if x.infinite? || x.nan?
n = prec + BigDecimal.double_fig

  • y = x
  • if 1 > x.abs
  •    x1 = x
    
  • elsif 1 == x.abs.round(prec)
  •    x1 = (-1 + sqrt(1 + x**2, prec))/x
    
  • else
  •    x1 = 1 / x
    
  • end
  • y = x1
    d = y
  • t = x
  • t = x1
    r = BigDecimal(“3”)
  • x2 = x.mult(x,n)
  • x2 = x1.mult(x1,n)
    while d.nonzero? && ((m = n - (y.exponent - d.exponent).abs) > 0)
    m = BigDecimal.double_fig if m < BigDecimal.double_fig
    t = -t.mult(x2,n)
    @@ -134,7 +141,13 @@
    y += d
    r += 2
    end
  • y
  • if 1 > x.abs
  •   y
    
  • elsif 1 == x.abs.round(prec)
  •   y*2
    
  • else
  •   PI(prec)*y.sign/4 - y
    
  • end
    end

Computes the value of e (the base of natural logarithms) raised to

the

e$BK-J!$G$9!#CY$$H?1~$G$9$,!#e(B

Date: Wed, 23 Sep 2009 19:12:28 +0900
Subject: [ruby-dev:39367] Almost endless loop of BigMath::atan(x) when x.abs >= 1

atan(x) = atan((-1+sqrt(1+x**2, prec))/x)
e$B$H$7$^$7$?!#e(B

sqrt e$B;H$&$H=E$/$J$k$N$Ge(B
atan(pi/4 - x) == (1 - atan(x))/(1 + atan(x))
e$B$r;H$&$N$O$I$&$G$7$g$&!#e(B
atan(BigDecimal(“0.7”), 100) e$B$G<B83$7$?$H$3$mG$/$i$$B.$/e(B
e$B$J$k$h$&$G$9!#e(B

— /ext/bigdecimal/lib/bigdecimal/math.rb (revision 25184)
+++ /ext/bigdecimal/lib/bigdecimal/math.rb (working copy)
@@ -8,3 +8,3 @@

cos (x, prec)

-# atan(x, prec) Note: |x|<1, x=0.9999 may not converge.
+# atan(x, prec)

exp (x, prec)

@@ -128,3 +128,3 @@
x = 1 / x if inv = x > 1

  • x = (-1 + sqrt(1 + x**2, prec))/x if dbl = x > 0.5
  • x = (1-x)/(1+x) if inv2 = x > 0.4
    n = prec + BigDecimal.double_fig
    @@ -142,3 +142,3 @@
    end
  • y *= 2 if dbl
  • y = pi / 4 - y if inv2
    y = pi / 2 - y if inv

e$B!V1_<~N($N8x<0$H7W;;K!!We(B
http://www.kurims.kyoto-u.ac.jp/~ooura/pi04.pdf
e$B$Ne(B (8)e$B<0!Je(Bx - x^3/3 + x^5/5 - …
e$B$r%*%$%i!<2CB.$7$?<0!K$re(B
e$B;H$C$?<BAu$b;n$7$F$_$?$+$C$?$N$G$9$,8m:[email protected]$b$jJ}$,$o$+$ie(B
e$B$J$$$N$G$"$-$i$a$^$7$?!#e(B

May Haskell be with you.

e$BK-J!$G$9!#e(B

atan(pi/4 - x) == (1 - atan(x))/(1 + atan(x))

tan(pi/4 - x) == (1 - tan(x))/(1 + tan(x))

e$B$G$7$?!#e(B

e$B$D$$$G$KJs9p$7$F$*$/$He(B revision 25184 e$B$Ge(B
atan(BigDecimal(“0.7”), 100) e$B$r7W;;$9$k$H>[email protected]<#3#0?t7eL\e(B
e$B$+$iCM$,0c$C$F$$$^$9!#e(B
e$B!Je(Becho ‘scale=100; a(0.7)’ | bc -l e$B$N7k2L$HHf$Y$^$7$?!Ke(B

May Haskell be with you.

e$BK-J!$G$9!#e(B

e$B$D$$$G$KJs9p$7$F$*$/$He(B revision 25184 e$B$Ge(B
atan(BigDecimal(“0.7”), 100) e$B$r7W;;$9$k$H>[email protected]<#3#0?t7eL\e(B
e$B$+$iCM$,0c$C$F$$$^$9!#e(B

e$B!!Lse(B2e$BG$NB.EY8~>e$O:F8=$7$^$7$?$,!"@:EY$Oe(Bpatche$B$r$"$F$?8e$NJ}$,Mn$A$F$$$k5$$,$7$^$9!#e(B

e$B$9$$$^$;$s!#>[email protected]<#3#0?t7eL$+$iCM$,0c$C$F$k$N$O;d$Ne(B
e$BJ}$G$7$?!#e(Borz

e$B$I$&$b=|;;[email protected]:EY$r;XDj$7$F$J$$$N$,0l0x$C$][email protected]$7$^$7$?!#e(B

  • x = (1-x).div(1+x, prec) if inv2 = x > 0.4

e$B$H$$$&$3$H$Oe(B
x = 1 / x if inv = x > 1
e$B$NJ}[email protected]:EY;XDj$7$J$$$H$$$1$J$$$G$9$M!#e(B
e$B<B:]e(B atan(BigDecimal(“7”), 100)
[email protected]$H#7#0?t7eL$+$iCM$,0c$C$F$-$^$9!#e(B

e$B$=$l$G$b!"[email protected]:EY$O2<2s$C$F$^$9$M!&!&!&e(B

e$B<}B+$9$k$^$G$N9?tJ,8m:9$,N/$C$F$/$k$N$G$=$NJ,@:EY$re(B e$BA}$d$5$J$$$H$$$1$J$$5$$,$7$^$9!#9?t$r=q$+$;$F$$k$H:GBge(B
e$B$G$b#2#0#09`$r1[$($J$$$h$&$J$N$Ge(B
x = (1-x).div(1+x, prec+3) if inv2 = x > 0.4
e$B$G$7$g$&$+!#e(B
e$BB>$N?t3X4X?t$K$D$$$F$bF1MM$N$3$H$,5/$-$F$$$J$$$+D4$Y$F$
$^$9!#e(B

May Haskell be with you.

e$BK-J!$G$9!#e(B

e$B<}B+$9$k$^$G$N9`?tJ,8m:9$,N/$C$F$/$k$N$G$=$NJ,@:EY$re(B
e$BA}$d$5$J$$$H$$$1$J$$5$$,$7$^$9!#e(B

e$B9M$(D>$7$?$i#29`L\0J9_==J,>.$5$/$J$k$N$G$=$NI,MW$O$J$5$=$&e(B
e$B$G$9$M!#[email protected](B pi/2 e$B$de(B pi/4
e$B$+$i$N8:;;$,$"[email protected]:EY$re(B
e$BA}$d$9$Y$-5$$,$7$^$9!#e(B

pi = PI(prec+1)

x = 1.div(x, prec+1) if inv = x > 1

x = (1-x).div(1+x, prec+1) if inv2 = x > 0.4
e$B$^$?$Oe(B
x = 2 - 1.div(1+x, prec+1) if inv2 = x > 0.4
e$B!J8e<T$NJ}$,$[$s$N$o$:$+B.$$!)!Ke(B

May Haskell be with you.

e$BK-J!$G$9!#e(B

x = 1.div(x, prec+1) if inv = x > 1
x = 2 - 1.div(1+x, prec+1) if inv2 = x > 0.4

e$B$$$m$$$m4V0c$($F$^$7$?!#<0$,0c$&$7<B9T%(%i!<$K$J$k$7!#e(B

e$B$=$l$He(B pi
e$B$r;H$&$N$K$b%3%9%H$,$+$+$k$N$G$9$M!#e(B100e$B7e!)$/$i$$$Oe(B
e$BJ8;zNs$G=`Hw$7$F$*$$$F$b$h$$$s$8$c$J$$$G$7$g$&$+!#e(B

@@ -122,14 +122,12 @@
raise ArgumentError, “Zero or negative precision for atan” if prec
<= 0
return BigDecimal(“NaN”) if x.nan?

  • pi = PI(prec)
    x = -x if neg = x < 0
  • pi = PI(prec+1) unless x <= 0.4
    return pi.div(neg ? -2 : 2, prec) if x.infinite?
    return pi / (neg ? -4 : 4) if x.round(prec) == 1
  • x = 1 / x if inv = x > 1
  • x = (-1 + sqrt(1 + x**2, prec))/x if dbl = x > 0.5
  • x = BigDecimal(“1”).div(x, prec+1) if inv = x > 1
  • x = BigDecimal(“2”).div(1+x, prec+1) - 1 if inv2 = x > 0.4
    n = prec + BigDecimal.double_fig
  • y = x
  • d = y
  • t = x
  • t = d = y = x
    r = BigDecimal(“3”)
    x2 = x.mult(x,n)
    @@ -141,5 +139,5 @@
    r += 2
    end
  • y *= 2 if dbl
  • y = pi / 4 - y if inv2
    y = pi / 2 - y if inv
    y = -y if neg
    @@ -216,4 +214,5 @@
    def PI(prec)
    raise ArgumentError, “Zero or negative argument for PI” if prec <=
    0
  • return BigDecimal(".314159265358e$B!JN,!Ke(B"[0…prec])*10 if prec
    <= 100
    n = prec + BigDecimal.double_fig
    zero = BigDecimal(“0”)

May Haskell be with you.

e$B6b0f$G$9!#e(B

2010e$BG/e(B1e$B7ne(B7e$BF|e(B19:46 TOYOFUKU Chikanobu
[email protected]:

e$B$D$$$G$KJs9p$7$F$*$/$He(B revision 25184 e$B$Ge(B
atan(BigDecimal(“0.7”), 100) e$B$r7W;;$9$k$H>[email protected]<#3#0?t7eL\e(B
e$B$+$iCM$,0c$C$F$$$^$9!#e(B
e$B!Je(Becho ‘scale=100; a(0.7)’ | bc -l e$B$N7k2L$HHf$Y$^$7$?!Ke(B

e$BLse(B2e$BG$NB.EY8~>e$O:F8=$7$^$7$?$,!"@:EY$Oe(Bpatche$B$r$"$F$?8e$NJ}$,Mn$A$F$$$k5$$,$7$^$9!#e(B

$ ./ruby -v
ruby 1.9.2dev (2010-01-07 trunk 26248) [i386-linux]

$ cat atan.rb
#!/usr/local/bin/ruby
require “bigdecimal”
require “bigdecimal/math.rb”
include BigMath

print atan(BigDecimal(“0.7”), 100)

$cat atan_patched.rb
#!/usr/local/bin/ruby
require “bigdecimal”
require “bigdecimal/math_patched.rb”
include BigMath

print atan(BigDecimal(“0.7”), 100)

$ echo ‘scale=100; a(0.7)’ | bc -l
.6107259643892086165437588764902360938185030661288276158428677300002315242905175259244566500045283295
$ ./ruby atan.rb
0.61072596438920861654375887649023609381850306612882761584286773000023152429051752592445665000452832950428011235007529724968872659238685714285714285714285714285714286E0
$ ./ruby atan_patched.rb
0.6107259643892086165437588764902360935332681667999685554401831662418422625599085443470502268704352706330580953308279960143351279436517E0

e$B$I$&$b=|;;[email protected]:EY$r;XDj$7$F$J$$$N$,0l0x$C$][email protected]$7$^$7$?!#e(B

— /ext/bigdecimal/lib/bigdecimal/math.rb (revision 25184)
+++ /ext/bigdecimal/lib/bigdecimal/math.rb (working copy)
@@ -8,3 +8,3 @@

cos (x, prec)

-# atan(x, prec) Note: |x|<1, x=0.9999 may not converge.
+# atan(x, prec)

exp (x, prec)

@@ -128,3 +128,3 @@
x = 1 / x if inv = x > 1

  • x = (-1 + sqrt(1 + x**2, prec))/x if dbl = x > 0.5
  • x = (1-x).div(1+x, prec) if inv2 = x > 0.4
    n = prec + BigDecimal.double_fig
    @@ -142,3 +142,3 @@
    end
  • y *= 2 if dbl
  • y = pi / 4 - y if inv2
    y = pi / 2 - y if inv

$ echo ‘scale=100; a(0.7)’ | bc -l
.6107259643892086165437588764902360938185030661288276158428677300002315242905175259244566500045283295
$ ./ruby atan_repatched.rb
0.6107259643892086165437588764902360938185030661288276158428677300002315242905175259244566500045283294757566224171893928143351279436517E0

e$B$=$l$G$b!"[email protected]:EY$O2<2s$C$F$^$9$M!&!&!&e(B
e$BB.EY$O!">/!9Mn$A$Fe(B1.5e$BG\DxEY$G$9!#e(B

e$BK-J!$G$9!#e(B

e$B$3$N%Q%C%A$r$"$F$?$"$H$G$b!"e(Batan(BigDecimal(‘0.7’), 100)e$B$O!"e(B
[email protected]:EY$,J]>c$5$l$F$$$k>.?tBhe(B100e$B0L$,0[$J$kCM$rJV$7$^$9!#e(B

e$B8m:9$,e(B 10^-100 e$B$h$j2<2s$C$F$$$k$N$G;[email protected]$H;W$$$^$9!#e(B
e$BNc$($Pe(B 0.50000000001 e$B$He(B 0.49999999999
e$B$O>[email protected]<Bh#17e$Ge(B
e$B0c$$$^$9$,8m:9$Oe(B 0.00000000002 e$B$7$+$"$j$^$;$s!#e(B
sqrte$BHG$NJ}$,8m:9$,>.$5$$$N$Oe(B sqrt e$B$r5a$a$k$H$-$K7W;;$Ne(B
e$BET9g>eMW5a$5$l$?M-8z7e?t0J>[email protected]:EY$G7W;;$9$k$?$a$=$l$,e(B
atan [email protected]$H;W$$$^$9!#e(B

e$B$H$3$m$GA0$Ke(B
|> e$B<}B+$9$k$^$G$N9?tJ,8m:9$,N/$C$F$/$k$N$G$=$NJ,@:EY$re(B |> e$BA}$d$5$J$$$H$$$1$J$$5$$,$7$^$9!#e(B | | e$B9M$(D>$7$?$i#29L\0J9_==J,>.$5$/$J$k$N$G$=$NI,MW$O$J$5$=$&e(B

e$B$H=q$-$^$7$?$,!V$=$NI,MW$O$J$$!W$N$Oe(B atan [email protected]$1$G$7$?!#e(B
e$B!Je(Batan e$B$Oe(B x -> 1/x e$B$J$I$7$F0z?t$r#10J2<$KM^$($F$$$k!Ke(B
sin,cos,exp e$B$O0z?t$,Bg$-$/[email protected]:EY$,$+$J$j0-$/$J$j$^$9!#e(B
e$B$=$l$r9MN8$7$?9)IW$r$7$J$$$H$$$1$J$$$h$&$G$9!#e(B

May Haskell be with you.

e$B6b0f$G$9!#e(B
e$B:.Mp$7$F$-$?$N$G!"@0M}$7$J$,$i=q$$$F$^$9!#e(B

2010e$BG/e(B1e$B7ne(B8e$BF|e(B14:14 TOYOFUKU Chikanobu
[email protected]:

e$B$H$$$&$3$H$Oe(B
x = 1 / x if inv = x > 1
e$B$NJ}[email protected]:EY;XDj$7$J$$$H$$$1$J$$$G$9$M!#e(B
e$B<B:]e(B atan(BigDecimal(“7”), 100) [email protected]$H#7#0?t7eL$+$iCM$,0c$C$F$-$^$9!#e(B

i).
*[Bug] *[[email protected](B]

atan(BigDecimal(“6”), 100)e$B$H$+$b$C$H82Cx$G$9$M!#e(B

$ echo ‘scale=100; a(6)’ | bc -l
1.4056476493802697809521934019958079881001980392225250914694378661427625409632993578391302747723584261
$ ./ruby atan6.rb
0.14056476493802697809521930776714836637758737148982182981896716134035368287750380505128119915512264084355727546888185920721779152566674E1

2010e$BG/e(B1e$B7ne(B8e$BF|e(B19:29 TOYOFUKU Chikanobu
[email protected]:

  • x = BigDecimal(“1”).div(x, prec+1) if inv = x > 1
    end
  • y *= 2 if dbl
  • y = pi / 4 - y if inv2
    y = pi / 2 - y if inv
    y = -y if neg
    @@ -216,4 +214,5 @@
    def PI(prec)
    raise ArgumentError, “Zero or negative argument for PI” if prec <= 0
  • return BigDecimal(".314159265358e$B!JN,!Ke(B"[0…prec])*10 if prec <= 100
    n = prec + BigDecimal.double_fig
    zero = BigDecimal(“0”)

ii).
e$B$3$N%Q%C%A$r$"$F$?$"$H$G$b!"e(Batan(BigDecimal(‘0.7’), 100)e$B$O!"e(B
[email protected]:EY$,J]>c$5$l$F$$$k>.?tBhe(B100e$B0L$,0[$J$kCM$rJV$7$^$9!#e(B

$ echo ‘scale=100; a(0.7)’ | bc -l
.6107259643892086165437588764902360938185030661288276158428677300002315242905175259244566500045283295
$ ./ruby atan07.rb
0.6107259643892086165437588764902360938185030661288276158428677300002315242905175259244566500045283294757566224171893928143351279436517E0

e$BB.EY$O!“e(BPIe$B$NJ8;zNs2=$r$7$J$$$G!“F1DxEY$G$9!#e(B
[email protected]!”>[email protected]$OF~$k$H;W$&$N$G!”$=$&$9$k$H!"Lse(B4e$B3d$[$IB.$/$J$j$^$9!#e(B

sqrte$B$,=E$$=hM}$J$N$O3N$+$J$N$G$9$,!“e(B
[email protected]:EY$r5>@7$K$9$k$N$O!”$=$l$3$=K\KvE>E]$J$N$G!"e(B
e$B$b$&>/$78!>Z$,I,MW$+$H!#e(B

e$BK-J!$G$9!#e(B

BigDecimal BigMath e$B$Ne(B sin(BigDecimal(“1E50”), 1) e$B$,e(B PI e$B$Ke(B
[email protected]:EY$NITB-$GL58B%k!<%W$K$J$j$^$9!#e(B

[email protected]:EY$rA}$d$5$J$$$H$$$1$J$$$3$H$O$$$1$J$$$N$G$9$,L58B%k!<%We(B
e$B!JK\Ev$OL58B%k!<%W$G$O$J$/$b$N$9$4$/;~4V$N$+$+$k%k!<%W$G$7$?!K$Ke(B
e$B$J$kD>@$N860x$Oe(B
BigDecimal(“1E50”) % BigDecimal(“3”)
e$B$,e(B 0.1E45 e$B$K$J$k$3$H$G$7$?!#e(B
e$BB?J,e(B y = q*x + r e$B$G>&e(B q e$B$re(B x
[email protected]:EY$G5a$a$F$$$k$+$ie(B
[email protected]$H;W$$$^$9!#$G$bIaDL$Oe(B r < x e$B$r4|BT$7$^$9$h$M!#e(B

May Haskell be with you.

e$BK-J!$G$9!#e(B

BigDecimal BigMath e$B$Ne(B sin(BigDecimal(“1E50”), 1) e$B$,e(B PI
e$B$Ke(B
[email protected]:EY$NITB-$GL58B%k!<%W$K$J$j$^$9!#e(B
e$B$5$5$d$+$J:GE,2=$bF~$l$^$7$?!#e(B
e$B%F%9%H%W%m%0%i%`$bIU$1$F$*$-$^$9!#e(B

— /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

======== e$B%F%9%He(B ========

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.

小林です。

----- Original Message -----
From: “TOYOFUKU Chikanobu” [email protected]
Subject: [ruby-dev:40085] Re: BigMath の sin(BigDecimal(“1E50”), 1)
が無限ループ

多分 y = q*x + r で商 q を x と同程度の精度で求めているから
だと思います。でも普通は r < x を期待しますよね。
ですね。
パッチを添付しましたので、どなたかコミットしていただけるとありがたいです。

それでは

e$BK-J!$G$9!#e(B

e$B8m:9$NI>2A$,IT==J,$G$7$?!#e(B
e$B&2e(B(-1)^nx^e$B!Je(B2n+1)/e$B!Je(B2n+1)! e$B$N7W;;$G8m:9$Oe(B
e$B&2e(Bx^e$B!Je(B2
n+1)/e$B!Je(B2*n+1)! = sinh(x)
e$BG$K$J$k$H;W$$$^$9!#e(B

sinh(2e$B&Pe(B) e$B$be(B sinh(e$B&Pe(B) e$B$bCM$,e(B 10
e$B$rD6$([email protected]:EY$rJ]$D$?$a$Ke(B
x e$B$re(B e$B&Pe(B/2 e$BDxEY0J2<$K$7$^$7$?!#e(Bcos e$B$bF1MM$G$9!#e(B

— /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.

This forum is not affiliated to the Ruby language, Ruby on Rails framework, nor any Ruby applications discussed here.

| Privacy Policy | Terms of Service | Remote Ruby Jobs