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

Masahiro_Kanai_SCan · September 23, 2009, 12:13pm

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

Masahiro_Kanai_SCan · January 7, 2010, 10:45am

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:9$N8+@Q$b$jJ}$,$o$+$ie(B
e$B$J$$$N$G$"$-$i$a$^$7$?!#e(B

May Haskell be with you.

Masahiro_Kanai_SCan · January 7, 2010, 11:46am

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>.?tE@0J2<#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.

Masahiro_Kanai_SCan · January 8, 2010, 6:15am

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>.?tE@0J2<#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!#>.?tE@0J2<#3#0?t7eL$+$iCM$,0c$C$F$k$N$O;d$Ne(B
e$BJ}$G$7$?!#e(Borz

e$B$I$&$b=|;;$N@:EY$r;XDj$7$F$J$$$N$,0l0x$C$]$$$N$G=$@5$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}$b@:EY;XDj$7$J$$$H$$$1$J$$$G$9$M!#e(B
e$B<B:]e(B atan(BigDecimal(“7”), 100)
e$B$@$H#7#0?t7eL$+$iCM$,0c$C$F$-$^$9!#e(B

e$B$=$l$G$b!"$b$H$N@: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.

Masahiro_Kanai_SCan · January 8, 2010, 7:58am

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!#$?$@e(B pi/2 e$B$de(B pi/4
e$B$+$i$N8:;;$,$"$k$N$G0l7e$O@: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.

Masahiro_Kanai_SCan · January 8, 2010, 11:30am

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.

Masahiro_Kanai_SCan · January 7, 2010, 7:22pm

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>.?tE@0J2<#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=|;;$N@:EY$r;XDj$7$F$J$$$N$,0l0x$C$]$$$N$G=$@5$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!“$b$H$N@:EY$O2<2s$C$F$^$9$M!&!&!&e(B
e$BB.EY$O!”>/!9Mn$A$Fe(B1.5e$BG\DxEY$G$9!#e(B

Masahiro_Kanai_SCan · January 10, 2010, 4:30pm

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
e$B@: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;EMMDL$j$@$H;W$$$^$9!#e(B
e$BNc$($Pe(B 0.50000000001 e$B$He(B 0.49999999999
e$B$O>.?tE@0J2<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>e$N@:EY$G7W;;$9$k$?$a$=$l$,e(B
atan e$B$N7k2L$K$bH?1G$7$F$$$k$N$@$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 e$B$@$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$-$/$J$k$H@: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.

Masahiro_Kanai_SCan · January 10, 2010, 3:37pm

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}$b@:EY;XDj$7$J$$$H$$$1$J$$$G$9$M!#e(B
e$B<B:]e(B atan(BigDecimal(“7”), 100) e$B$@$H#7#0?t7eL$+$iCM$,0c$C$F$-$^$9!#e(B

i).
*[Bug] *[e$BMW=$@5e(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
e$B@: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
e$B$?$@!”>e$N=$@5$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
BigDecimale$B$GB.EY$N$?$a$K@: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

Masahiro_Kanai_SCan · January 14, 2010, 10:26am

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
e$BMW5a$9$k@:EY$NITB-$GL58B%k!<%W$K$J$j$^$9!#e(B

e$B@: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
e$B$HF1DxEY$N@:EY$G5a$a$F$$$k$+$ie(B
e$B$@$H;W$$$^$9!#$G$bIaDL$Oe(B r < x e$B$r4|BT$7$^$9$h$M!#e(B

May Haskell be with you.

Masahiro_Kanai_SCan · January 14, 2010, 9:54am

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
e$BMW5a$9$k@: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.

Masahiro_Kanai_SCan · January 15, 2010, 1:20pm

å°æž—ã§ã™ã€‚

----- 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 ã‚’æœŸå¾…ã—ã¾ã™ã‚ˆãã€‚
ã§ã™ãã€‚
ãƒ‘ãƒƒãƒã‚’æ·»ä»˜ã—ã¾ã—ãŸã®ã§ã€ã©ãªãŸã‹ã‚³ãƒŸãƒƒãƒˆã—ã¦ã„ãŸã ã‘ã‚‹ã¨ã‚ã‚ŠãŒãŸã„ã§ã™ã€‚

ãã‚Œã§ã¯

Masahiro_Kanai_SCan · January 15, 2010, 11:56am

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(B2n+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$($k$N$G@: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.

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

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

cos (x, prec)

exp (x, prec)

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:9$N8+@Q$b$jJ}$,$o$+$ie(B e$B$J$$$N$G$"$-$i$a$^$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>.?tE@0J2<#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$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

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

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

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

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

cos (x, prec)

exp (x, prec)

e$B$H=q$-$^$7$?$,!V$=$NI,MW$O$J$$!W$N$Oe(B atan e$B$@$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$-$/$J$k$H@: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

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

sin test

PI/2

cos test

PI

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:9$N8+@Q$b$jJ}$,$o$+$ie(B
e$B$J$$$N$G$"$-$i$a$^$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>.?tE@0J2<#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$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

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

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

e$B$H=q$-$^$7$?$,!V$=$NI,MW$O$J$$!W$N$Oe(B atan e$B$@$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$-$/$J$k$H@: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