arcsinを数値計算する方法を探し求めて彷徨う

2023-07-17
#math

以前数学関数を自前で実装する一環でatanは実装したが、欠けていた逆三角関数asinを実装したかった。

マクローリン展開

まずはなんの工夫もなくマクローリン展開を使ってみる。 \(\arcsin x\) の微分を求めるには、\(y = \arcsin x\) とおくと逆関数だから \(x = \sin y\) となり、全微分 \(dx = (\sin y)^\prime dy = \cos y dy\) から、

$$ \begin{align*} \frac{dy}{dx} &= \frac{1}{\cos y} = \frac{1}{\sqrt{1 - \sin^2 y}} \\ &= \frac{1}{\sqrt{1 - x^2}} \end{align*} $$

以降は愚直に微分すると

$$ \begin{align*} \left( \arcsin x \right) ^{\prime \prime} &= x \cdot \left( 1 - x^2 \right)^{3/2} \\ \left( \arcsin x \right) ^{\prime \prime \prime} &= \left( 2x^2 + 1 \right) \cdot \left( 1 - x^2 \right)^{5/2} \end{align*} $$

（以下略）となる。

微分した関数を使ってマクローリン展開すると一般項は、偶数項が0となり奇数項だけ残り

$$ \begin{align*} f(x) &= \arcsin x \\ &= x + \frac{1}{6}x^3 + \frac{3}{40}x^5 + \frac{5}{112}x^7 + \frac{35}{1152}x^9 + \frac{63}{2816}x^{11} + \cdots \\ &= x + \frac{1}{3} \left( \frac{1}{2} \right) x^3 + \frac{1}{5} \left( \frac{3 \cdot 1}{4 \cdot 2} \right) x^5 + \frac{1}{7} \left( \frac{5 \cdot 3 \cdot 1}{6 \cdot 4 \cdot 2} \right) x^7 + \cdots \\ &= \sum_{n=0}^\infty \frac{1}{4^n (2 n + 1)} \binom{2n}{n} x^{2n + 1} \end{align*} $$

（\(\binom{n}{k}\) は組み合わせ \({}_n C_k = \frac{n!}{(n - k)! k!} \)）。

精度をマシにしたい

マクローリン展開ということもありx=0付近では精度はいいが、遠ざかるにつれ精度が悪くなる。 単に項の数を増やしても精度が上がらなかった。 逆正弦関数なのでx=1ではジャスト\(\pi/2\)を返して欲しい。 じゃあx=1でテイラー展開するか、と思ったが微分した関数に1を与えると0除算で\(\infty\)なので求められない。

ググってみるとFinding the Taylor series of $\arcsin(1-x)$ - Mathematics Stack Exchange が見つかった。 それによると

$$ \arcsin(1 − x) = \frac{\pi}{2} − \sqrt{2}𝑥^{1/2} − \frac{\sqrt{2}}{12} x^{3/2} + O(x^{5/2}) $$

とのこと。 これを使えばx=1付近の値を求められる。

項を増やしたい（sympy使用、失敗編）

誤差が\(x^{5/2}\)のオーダーでは心許ないので、もっと項を増やしたい。 ちょっと手計算で微分を計算してみたが早々に計算間違いをしでかしたので、ググってsympyというPython上で変数を使った微分ができるライブラリがあることを知った：

import numpy as np
import scipy as sc
import sympy as sy
import matplotlib.pyplot as plt
from scipy.misc import derivative

sy.init_printing(order='grevlex')

x = sy.symbols('x', real=True)
t = sy.symbols('t', real=True)  # t = sqrt(x)
y = sy.asin(1 - t ** 2)

for n in range(5):
    deriv = sy.simplify(sy.diff(y, t, n))  # n階微分
    print(f'\n### {n}')
    display(deriv.subs([(t, sy.sqrt(x))]))  # t = sqrt(x)を代入
    display(deriv.subs([(t, 0)]))  # t=0を代入する
  • diffでn階微分、subsでtに0を代入した結果を計算できる

微分を手計算しなくていいのすごい！ でt=0を代入してみるとNaNになってしまった。 分母に\(-(x^2-1)^2+1\)があって0を代入するから0除算で\(\infty\)に発散してしまうからダメか…（と勘違いして別の方法を探した）。

テイラー展開を試す

x=0だけで展開してるし、x=1では傾きが無限大で展開できないので、もう少し違う点をいくつか持っておいて展開すればいいんじゃないか？ と試してみたが、結局1付近では精度が悪かった。

二分探索を試す

級数展開が使えないならどうしたものか…と見ていたらxが正の場合には\(x \le \arcsin x \le \frac{\pi}{2} x\)のようなので、二分探索はどうかと考えた。 asinは単調増加なので中点のsinの値を計算してその値より大きいか小さいかで探索していく：

double my_asin_bisect(double x) {
  if (x < -1 || x > 1)
    return NAN;
  if (x < 0)
    return -asin(-x);

  // Binary search
  double min = x, max = M_PI_2 * x;
  for (int i = 0; i < 16; ++i) {
    double m = (min + max) * 0.5;
    double y = sin(m);
    if (x >= y)
      min = m;
    else
      max = m;
  }
  return (min + max) * 0.5;
}

計算したところ、16回ループで誤差（標準のライブラリで求めた値との差）は1e-6くらいだった。 sin(m)も独自実装ではマクローリン展開なので相当重そうだけど、最悪計算できるだけマシだからこれでもいいか…と一瞬思った。

Wolfram|Alphaで係数を調べる

その後ググってるとWolfram|Alphaでもっと項の係数を調べられた： taylor series arcsin(1-x) - Wolfram|Alpha 一般項はわからないけど誤差が\(x^{11}\)なら十分か、と思った。

ていうか係数は∞にならずに計算できるんだっけ？と気づいた。 どうもsympyでうまく約分できてないだけで、１階微分\(\frac{-2}{\sqrt{2 - x^2}}\)から計算していけばちゃんと計算できた。

組み合わせる

二分探索で計算するよりも、x=0近辺では\(\arcsin x\)、x=1近辺では\(\arcsin (1-x)\)のマクローリン展開を使うのが精度がよかった。 項をどこまで計算するかによるけどx=0.55あたりで切り替えるのがよさそうで、誤差は最大で1.0e-10のオーダーだった。

muslのソースを見てみる

答え合わせ的にオープンソースのライブラリ、muslのソースを見てみた：musl/src/math/asin.c at master · cloudius-systems/musl コメントによれば、

  • \(R(x^2)\)を\(x^3\)以降の有理近似（誤差2^-58.75以下）
  • [0, 0.5] では \(x + x \cdot x^2 \cdot R(x^2)\)
  • [0.5, 1] では \(\frac{\pi}{2} - 2 \arcsin {\sqrt{(1 - x) / 2}}\)

と書かれている。[0.5, 1.0]の場合の計算式

$$ \arcsin x = \frac{\pi}{2} - 2 \arcsin {\sqrt{\frac{1 - x}{2}}} $$

という公式 がどうやったら導けるのか結局わからなかった…。 （これはsinの半角の公式 \(\sin^2 \frac{\theta}{2} = \frac{1 - \cos \theta}{2}\) を使うようで、

$$ \begin{align*} \sin^2 \frac{\theta}{2} &= \frac{1 - \cos \theta}{2} \\ \therefore \theta &= 2 \arcsin \sqrt{\frac{1 - \cos \theta}{2}} \end{align*} $$

\(\cos \theta = x\) とおき、 \(\arccos x = \pi / 2 - \arcsin x\) を使って

$$ \begin{align*} \theta = \arccos x &= 2 \arcsin \sqrt{\frac{1 - x}{2}} \\ \frac{\pi}{2} - \arcsin x &= 2 \arcsin \sqrt{\frac{1 - x}{2}} \\ \arcsin x &= \frac{\pi}{2} - 2 \arcsin \sqrt{\frac{1 - x}{2}} \end{align*} $$

xが0.5より大きければ平方根の値は0.5以下になるのでその範囲の計算で求められる。）

muslはMITライセンスとのことなので、使用する分には問題なさそう。

arccosは？

arcsinを使って、 \(\frac{\pi}{2} - \arcsin x\) で求められる。

リンク

関連記事

新着記事