From: "Masahiro Kanai (CanI)" <cani.m.61st@...>
Date: 2010-01-10T23:36:54+09:00
Subject: [ruby-dev:40019] Re: Almost endless loop of BigMath::atan(x) when 	x.abs >= 1

金井です。
混乱してきたので、整理しながら書いてます。

2010年1月8日14:14 TOYOFUKU Chikanobu <nobu_toyofuku@nifty.com>:
> ということは
>    x = 1 / x if inv = x > 1
> の方も精度指定しないといけないですね。
> 実際 atan(BigDecimal("7"), 100) だと７０数桁目から値が違ってきます。

i).
*[Bug] *[要修正]

atan(BigDecimal("6"), 100)とかもっと顕著ですね。
-
$ echo 'scale=100; a(6)' | bc -l
1.4056476493802697809521934019958079881001980392225250914694378661427625409632993578391302747723584261
$ ./ruby atan6.rb
0.14056476493802697809521930776714836637758737148982182981896716134035368287750380505128119915512264084355727546888185920721779152566674E1


2010年1月8日19:29 TOYOFUKU Chikanobu <nobu_toyofuku@nifty.com>:
>  豊福です。
>
>>  x = 1.div(x, prec+1) if inv = x > 1
>>  x = 2 - 1.div(1+x, prec+1) if inv2 = x > 0.4
>
>  いろいろ間違えてました。式が違うし実行エラーになるし。
>
>  それと pi を使うのにもコストがかかるのですね。100桁？くらいは
> 文字列で準備しておいてもよいんじゃないでしょうか。
>
>
> @@ -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(".314159265358（略）"[0..prec])*10 if prec <= 100
>     n      = prec + BigDecimal.double_fig
>     zero   = BigDecimal("0")

ii).
このパッチをあてたあとでも、atan(BigDecimal('0.7'), 100)は、
精度が保障されている小数第100位が異なる値を返します。

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

速度は、PIの文字列化をしないで、同程度です。
ただ、上の修正は入ると思うので、そうすると、約4割ほど速くなります。

sqrtが重い処理なのは確かなのですが、
BigDecimalで速度のために精度を犠牲にするのは、それこそ本末転倒なので、
もう少し検証が必要かと。

--
Masahiro Kanai (CanI)
http://d.hatena.ne.jp/CanI/