SIMD並列化ライブラリSmartVectorDotNet開発の知見まとめ(3) 初等関数の実装

1. IEEE754浮動小数型の低レベル操作
2. SIMD演算の基礎
3. 初等関数の実装 (本稿)
4. C#とdotnetの最適化

---

関数実装の基本戦略

SmartVectorDotNetでは多くの初等関数のベクトル版を実装しました。
しかし、SIMD命令が提供されているのは四則演算やビットシフト、比較演算など、ほとんどのプログラミング言語において専用の演算子が提供されるような極小数の演算に限定されます。
そこで、計算機的には何らかの近似計算によって本来の計算値と同等の解が得られるようなアルゴリズムを用います。

テイラー展開・マクローリン展開

関数$f(x)$が点$x=a$において無限回微分可能であるとき、$n$次導関数が$f^{(n)}(x)$で表現されるならば、この点付近において

$$f(x) = \sum_{n=0}^\infty\cfrac{f^{(n)}(a)}{n!}(x-a)^n$$

の変換が可能な場合があります。このような変換を$x=a$におけるテイラー展開、これにより得られた級数をテイラー級数といいます。

特に、$a=0$のケース

$$f(x) = = \sum_{n=0}^\infty\cfrac{f^{(n)}(0)}{n!}(x)^n$$

はマクローリン展開/マクローリン級数といいます。
非ゼロ点でのテイラー展開より計算が楽になる場合が多いため非常によく利用されます。

テイラー展開・マクローリン展開が示すのは、これらの変換が可能である関数は加算と乗算のみからなる極めて単純な式で表現できる、ということです。
テイラー級数は一般に次数の若い項ほど解に対し支配的であるため、目的の精度が達成可能であるなら有限個の項を計算した時点で打ち切ることができます。
導関数で表されている各項の係数はあらかじめ計算して定数とすることができる点も計算機的に有利と言えます。

前述のように、テイラー級数がもとの関数と一致するのは展開を実施した点$x=a$の近傍のみであり、この点から離れると大きく精度が低下します。
このため、テイラー展開による計算は定義域を限定する必要があります。

ニュートン法

ニュートン法はテイラー展開とは異なる計算機的近似戦略として知られています。
一部の方程式$f(x)=0$は、1次導関数を用いた漸化式によって解析的に$x$を求めることができます。
具体的には、以下の数式で表現されます。

$$x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}$$

ここで、実装したい関数$x=g(a)$があり、逆関数$a=g^{-1}(x)$が得られている場合を考えます。
これを方程式$0 = g^{-1}(x)-a := f(x)$の解$x$を求める問題と捉えると、$f(x)$が微分可能であるならば前述の漸化式が適用できます。
この性質のため、例えば冪根のように逆関数が平易な形式で表せる関数の実装に対して威力を発揮します。

ニュートン法では$x_1$の選択が解の収束速度に大きく影響します。

各関数の具体的な実装方針

三角関数

正弦・余弦

正弦$\sin x$および余弦$\cos x$のマクローリン級数は以下の形式で表せることが知られています。

$$\begin{aligned} \sin x &= \sum_{n=0}^\infty\cfrac{(-1)^n}{(2n+1)!}x^{2n+1} &= x-\cfrac{x^3}{3!}+\cfrac{x^5}{5!}-\cfrac{x^7}{7!}+\cdots\\ \cos x &= \sum_{n=0}^\infty\cfrac{(-1)^n}{(2n)!}x^{2n} &= 1-\cfrac{x^2}{2!}+\cfrac{x^4}{4!}-\cfrac{x^6}{6!}+\cdots\\ \end{aligned}$$

また、正弦・余弦には$x=0$における対称性および周期$2\pi$での周期性があるため、以下の式により$x$を変換することができます。

$$\begin{aligned} \sin x &= \sin (x + 2n\pi) \\ &= -\sin(-x) \\ &= \sin(\pi-x) \\ \cos x &= \cos (x + 2n\pi) \\ &= \cos(-x) \\ &= -\cos (\pi - x) \\ \end{aligned}$$

この変換により、コアロジックSinCore(x)・CosCore(x)の定義域を$0 \le x \lt \frac{\pi}{2}$に限定することができます。

正接

また、正接$\tan x$は

$$\tan x = \frac{\sin x}{\cos x}$$

ですから、$\sin x$および$\cos x$が実装できればこれを用いて実装することができます。

NOTE1: 三角関数には他にも$\cos x = \sin (\frac{\pi}{2} - x)$などの関係式が存在するため、コアロジックの定義域は更に狭い$0 \le x \lt \frac{\pi}{4}$を採用することもできます。ただし、この戦略はSinCore(x)・CosCore(x)の呼び分けが必要になるためベクトル演算では無駄が多く、SmartVectorDotNetでは採用を見送っています。
NOTE2: NOTE1に示したように$\sin x$と$\cos x$は相互に変換が可能なのでコアロジックの実装はいずれか一方に限ることもできます。しかし、浮動小数では真の円周率を表現できないことから、変換による実装を行った方の関数で$x=0$の場合に厳密解が得られなくなってしまいます。これは精度要求は問題になりませんが、$x=0$のケースは他の$x$よりも利用機会が多い上に用途の特殊性も考えられるため、SmartVectorDotNetでは$\sin x$と$\cos x$それぞれにコアロジックを実装しています。

逆三角関数

逆正接

逆正接$\arctan x$のマクローリン級数は以下の形式で表せることが知られています。

$$\arctan x = \sum_{n=0}^{\infty}\cfrac{(-1)^n}{2n+1}x^{2n+1} = x - \cfrac{x^3}{3} + \cfrac{x^5}{5} - \cfrac{x^7}{7} + \cdots$$

更に、逆正接には次の性質があります。

$$\begin{aligned} &\arctan(-x) = -\arctan x && (\mathrm{i})\\ &\arctan\cfrac{1}{x} = \mathrm{sign}\,x\cdot\cfrac{\pi}{2} - \arctan x && (\mathrm{ii})\\ &\arctan x + \arctan y = \arctan\cfrac{x + y}{1 - xy} \\ &\Rightarrow \arctan x + \arctan (-1) = \arctan x - \cfrac{\pi}{4} = \arctan\cfrac{1 - x}{1 + x} \\ &\Rightarrow \arctan x = \arctan\cfrac{1 - x}{1 + x} + \cfrac{\pi}{4} && (\mathrm{iii}) \end{aligned}$$

この性質について、コアロジックAtanCore(x)の定義域に関して次のことが言えます。

• $|x|$と$|-x|$の大小が入れ替わる点は$x=0$であるから、式(i)を用いることでAtanCore(x)の定義域を$x \ge 0$に限定できます。
• $|x|$と$\left|\frac{1}{x}\right|$の大小が入れ替わる点は$x=\pm 1$であるから、式(ii)を用いることでAtanCore(x)の定義域を$-1 \le x \le 1$に限定できます。
• $|x|$と$|\frac{1-x}{1+x}|$の大小が入れ替わる点は$x=-1\pm\sqrt{2}$であるから、式(iii)を用いることでAtanCore(x)の定義域を$-1-\sqrt{2} \le x \le -1+\sqrt{2}$に限定できます。

以上の3つを組み合わせれば、AtanCore(x)の定義域を$0 \le x \le \sqrt{2} - 1$に限定することができます。

逆正弦・逆余弦

また、逆正弦及び逆余弦は逆正接を用いて次の式で求めることができます。

$$\begin{aligned} \arcsin x &= 2\arctan\cfrac{x}{1 + \sqrt{1 - x^2}} \\ \arccos x &= 2\arctan\cfrac{\sqrt{1 - x^2}}{1 + x} \end{aligned}$$

指数関数

ネイピア数を底とする指数関数

ネイピア数を底とする指数関数$e^x$のマクローリン級数は以下の形式で表せることが知られています。

$$e^x = \sum_{n=0}^{\infty}\cfrac{1}{n!}x^n = 1 + x + \cfrac{x^2}{2} + \cfrac{x^3}{6} + \cfrac{x^4}{24} + \cdots$$

正の実数$a,\ b$に対して次の関係が成り立つことが知られています。

$$b^x = a^{x\log_a b}$$

ここで、特に$a=2,\ b=e$の場合を考えます。

$$e^x=2^{x\log_2 e}$$

更に、$n:=\mathrm{round}(x\log_2 e),\ y:= x\log_2 e - n$と置き換えると、更に次のような変換が可能です。

$$e^x = 2^n \cdot 2^p = 2^n \cdot e^{y\log_e 2}$$

この式について、以下のことが言えます。

• $\log_2 e$および$\log_e 2$は定数です。
• $n\in\mathbb{Z}$であることから、$2^n$はビットシフトで置換できます。
• $y$について、$-0.5 \le y \le +0.5$が成り立ちます。

故にコアロジックExpCore(x)の定義域は$-0.5 \le x \le +0.5$に限定できます。

対数関数

対数関数のマクローリン展開は、$\ln x$よりも$\ln (1+x)$について考える方が平易な形式を得ることが知られています。

$$\ln(1 + x) = \sum_{n=0}^{\infty}\cfrac{(-1)^n}{n+1}x^{n+1} = x - \cfrac{x^2}{2} + \cfrac{x^3}{3} - \cfrac{x^4}{4} + \cdots \ (\Leftarrow -1 \lt x \le 1)$$

式中に示したとおり、このマクローリン級数は$-1 \lt x \le 1$の定義域でしか収束しないことに注意してください。

さて、自然対数について、正の実数$a,\ b$に対し次の関係が成り立つことが知られています。

$$\ln(a \cdot b) = \log a + \log b$$

この性質を踏まえ、$y, n :\Leftrightarrow x=y\cdot 2^n, n \in \mathbb{Z}, 1 \le |y| \lt 2$となるような$y,\ n$を置くと次の変換が成り立ちます。

$$\begin{aligned} \ln x =& \ln y\cdot 2^n \\ =& \ln y + \ln 2^n \\ =& \ln y + n\ln 2 &(\mathrm{i}) \\ =& \left(\ln\cfrac{y}{2} + \ln 2\right) + n\ln 2 &(\mathrm{ii}) \end{aligned}$$

この式について以下のことが言えます。

• $\ln 2$は定数です。
• IEEE754 binary32/64の正規化数においては、$y$および$n$はビット演算のみで求めることができます。
• $y\le\sqrt{2}$であるならば$(\mathrm{i})$の式を、そうでなければ$(\mathrm{ii})$の式を利用するようにすれば、より少ない項でも精度の良い計算結果を得ることができます。

双曲線関数

双曲線関数は次のように定義されています。

$$\begin{aligned} \sinh x &= \cfrac{e^x - e^{-x}}{2} \\ \cosh x &= \cfrac{e^x + e^{-x}}{2} \\ \tanh x &= \cfrac{e^x - e^{-x}}{e^x + e^{-x}} = \cfrac{e^{2x} - 1}{e^{2x} + 1} \\ \end{aligned}$$

このことから、指数関数の実装が完了していれば双曲線関数は容易に実装できます。

逆双曲線関数

逆双曲線関数は次のように定義されています。

$$\begin{aligned} \mathrm{arsinh}\,x &= \ln(x + \sqrt{x^2+1}) \\ \mathrm{arcosh}\,x &= \ln(x + \sqrt{x^2-1}) \\ \mathrm{artanh}\,x &= \cfrac{1}{2}\ln\cfrac{1+x}{1-x} \\ \end{aligned}$$

このうち、$\mathrm{artanh}\,x$は定義域が$-1 \lt x \lt 1$に限定されているため対数関数の実装をそのまま用いて実装しても問題なく動作します。

一方、$\mathrm{arcsinh}\,x$および$\mathrm{arccosh}\,x$は$x$の2乗項を含んでいるため、
浮動小数演算においては本来の定義域よりも大幅に狭い範囲でしか適切に解を得ることができません。
そこで、次のように変形することで本来の定義域による計算が可能なようにします。

• $\mathrm{arccosh}\,x$に対して、

  $$\begin{aligned} \mathrm{arcosh}\,x &= \ln(x + \sqrt{x^2-1}) \\ &= \ln\left\{x\left(1 + \frac{\sqrt{(x - 1)}\sqrt{(x + 1)}}{x}\right)\right\} \\ &= \ln x + \ln\left(1 + \frac{\sqrt{(x - 1)}\sqrt{(x + 1)}}{x}\right) \end{aligned}$$

• $\mathrm{arcsinh}\,x$に対して、$|x|$が十分に大きいならば

  $$\begin{aligned} \mathrm{arsinh}\,x &= \ln(x + \sqrt{x^2+1}) \\ &\sim \mathrm{Sign}(x) \ln(2|x|) \sim \mathrm{Sign}(x) (2 + \ln|x|) \end{aligned}$$

立方根

立方根はこれまでの初等関数とは異なり、ニュートン法での実装を検討します。
$x$の立方根$y=\sqrt[3]{x}$を求めるに当たり、方程式$f(x)=x^3-y=0$を考えます。

$$\begin{aligned} &f(x) = x^3 - y = 0 \\ &\begin{aligned} \Rightarrow x_{n+1} &= x_n - \cfrac{{x_n}^3 - y}{3{x_n}^2} \\ &= \cfrac{2}{3}x_n + \cfrac{y}{3{x_n}^2} \end{aligned} \end{aligned}$$

この漸化式は非常に平易であるため容易に実装できます。

$x_1$の選び方

前述のとおり、ニュートン法では初期値$x_1$の選定が収束速度に大きく影響します。
SmartVectorDotNetでは次の戦略に基づいて初期値の決定を行っています。

計算の対象がIEEE754 binary32/64であることを踏まえると、$a, n :\Leftrightarrow x = a\cdot 2^n, n \in \mathbb{Z}, 1 \le a \lt 2$であるような$a,\ n$はビット演算のみで求めることができます。
ここで、$a\cdot2^{\left\lfloor\frac{n}{3}\right\rfloor}$を考えると、この値は$x$の立方根$\sqrt[3]{x} = \sqrt[3]{a\cdot 2^n} = \sqrt[3]{a}\cdot 2^{\frac{n}{3}}$に対し、少なくとも$x$よりは近い値であると期待できます。
この事実に基づき、$a\cdot2^{\left\lfloor\frac{n}{3}\right\rfloor}$を初期値$x_1$とします。