|
ソフトウェア実験室 |
|
逆三角関数を求める |
| 最終更新：2006/01/30 新規 |
三角関数ができれば、その逆も欲しくなる。ATan(Arctan) が基本なので、ATan が欲しいところであるが・・・・
●原理
○ATanの級数展開
長い円周率 で既に紹介しているが、ATan の級数展開は、|x| ≦ 1 でのみ定義されているので、一般値では実用にならない。
○ASinの級数展開
そこで、ASin を求め、ATan に変換することになる。
上図より、A = ATan(X) = ASin(X/Y)、ところで、
Y2 = 1 + X2
なので、ASin(X/Y) = ASin(X/Sqrt(1 + X2))となる。つまり、
ATan(X) = ASin(X/Sqrt(1 + X2))
ASinは、
ASin(X) = ∑ (2 * k - 1)!!/((2 * k)!! * (2 * k + 1))*X(2*k+1) [k = 0 to ∞]
で表される。!! は2重階乗であるが、これに怖気づかずに、素直に定義に基づいて紐解けば、簡単である。2重階乗を通常の階乗に変換する方法*もあるが、実計算上、余計にややこしくなる。
→* 2重階乗を通常の階乗に変換し、整理すると、 (2 * k)!/(2(2*k) * (k!)2)
2重階乗は、偶数の場合は、偶数のみの乗積、奇数の場合は、奇数のみの乗積となる。k = 3 と 4 のときの階乗部分の係数の状態は、
k = 3 のとき、 (1*3*5)/(2*4*6)
k = 4 のとき、 (1*3*5*7)/(2*4*6*8)
となる。k が、1増加するごとに、分子は、2 * k -1、分母は、2 * k の積が追加されて行くことになる。従って、各項の係数Ckは、初期値を 1 として、k = 1 から開始し、
Ck = Ck * (2 * k - 1)/(2 * k)
と、累算していけば良いことがわかる。これで、毎回階乗の計算はしなくて良い。また、Ckは常に 1 以下の安定した数値になるので、都合が良い。総合の係数は、
Ck / (2 * k + 1)
となる。
○級数の収束を予測する。
この級数の最悪の場合は、X=1 ときで、係数で収束が決まる。これを実測してみると、
| k | 係数のオーダ |
| 1 | -1 |
| 100 | -4 |
| 500 | -5 |
| 1,000 | -6 |
| 5,000 | -7 |
係数のオーダ(指数)は表のようになる。減少はするも、5,000項まで計算しても、7桁程度の精度しか出ない。X(2*k+1) を早く 1 より小さくなるような工夫が必要となる。
○工夫(余角を使う)
A = π / 2 - B
であるから、
ASin(X) = π / 2 - ASin(Y)
ところで、
X2 + Y2 = 1
だから、
ASin(Y) = ASin(Sqrt(1 - X2))
となる。Sqrt(1 - X2) は、X が、1 に近いほど、小さくなるので、収束が期待できる。今、X = 0.99 とすれば、 Y = 0.141067 となる。k = 100 のとき、X(2*k+1) は、-171乗であり、100項で170桁の精度がでる。この場合、高精度の開平が必要となる。 余角を採用する判定は、
X > Sqrt(1 - X2)
を解けばよい。これを解くと、
X > 1 / Sqrt(2) = 0.7
つまり、絶対値が0.7以上の場合は、余角を採用する。
●実験結果
UltraMathにて、約500桁の精度で、ASin(0.5) を計算してみた。上半分が、ラジアン、下半分が度で表示している。30度に限りなく近い値が得られた。
