CORDICアルゴリズム〜三角関数等の数値計算法〜

CORDICアルゴリズム

これは、加減算,bitshiftとLUT(Look Up Table)(定数表)のみの計算で　三角関数(sin,cos),逆三角関数(arcsin,arccos),双曲線関数(sinh,cosh),逆双曲線関数(arcsinh,arcosh)などを求めることができる。
数3(複素数)を理解していれば、アイデアは割とすぐに理解できるだろう。

まず前提として、2進数において$2^{-k}(k:整数)$を掛ける操作は数字を小数点に対して右にk個移動することで実現できる。筆算を考えればわかるように、乗算や除算はシフト操作と加減算操作の組み合わせからなるが、加減算に比べて処理が大変である。テーブルの定数表を参照しながら、高速で三角関数の値を計算したいときに有効である。計算方法の導出について見ていくが、複素平面をイメージしながら読むと理解が進むと思う。 CORDIC(Qiita) 参照。
行列を用いた解説 があるが、複素数のほうが本質が理解しやすいように思うので、今回は複素数を用いた説明を試みる。

この表記を用いて、求めるものを表してみる。
$\varepsilon_k\in\{-1,1\},~\tt=\sum_{k=0}^n\ve_k\tt_k$のようにある$\tt_k$を用いて$\tt$が分解できたとする。このとき指数法則より($\prod$は総積記号)

$\begin{align} e^{i\tt}=&\prod_{k=0}^n e^{i\ve_k\tt_k}\\ =&\prod_{k=0}^n (\cos\tt_k+i\ve_k\sin\tt_k)\\ =&\prod_{k=0}^n \cos\tt_k(1+i\ve_k\tan\tt_k) \end{align}$
ここで$\tan\tt_k=2^{-k}$となるよう$0<\tt_k<\pi/2$を選ぶ。そして
$$R_0=\prod_{k=0}^n \cos\tt_k$$とする。すると

となる。この式がアルゴリズムの本質である。$x,y$を実数とすると、
$(1+i\ve_k 2^{-k})(x+iy)=(x-\ve_k 2^{-k}y)+i(y+\ve_k2^{-k}x)$
となるが、よく見ると$ \ve_k 2^{-k}y,\ve_k2^{-k}x$はビットシフトと符号反転だけなので、$1+i\ve_k 2^{-k}$の乗算は加減算とビットシフトでできることがわかる。よって、

とすると分解の結果から、$\cos\tt=Re (Z_{n+1}),~~\sin\tt=Im(Z_{n+1})$のようにして三角関数がn+1回の加減算とビットシフトで計算できる。分解については非回復型除算という手法で求める。
これは2進数の各桁の値が$\{1,-1\}$のどちらかであるとして、除算を筆算で行うやり方とほぼ同じである。$\tt_k=\arctan(2^{-k})$が大まかに指数関数的減少とみなせるので、$\tt_1$から順番に足すか引くかを部分剰余の符号から決定して次の部分剰余にフィードバックする再帰処理なので、詳細は省略する。誤差約$2^{-n}$で分解可能であることの正当性は この記事 で解説している。$\arctan$は$\tan$の逆関数である。

③は特に、複素平面上では点$0,Z_k,Z_{k+1}$という鋭角$\tt_k$の直角三角形が、$\ve_k$の向きに斜辺と辺を共有する形で並んでいる姿を想像する。

10進数化

$\tt_k$を分解基底と称する。分解$ \tt=\sum_{k=0}^n\ve_k\tt_k+\epsilon~~$ である範囲の任意の$\tt$を誤差$\ee$を$\tt_n$未満にできるとき、分解は完全であると称する。
分解の一般論と正当性は 「N進数の一般化について」 を参照してほしい。分解方法を工夫して10進CORDICについて考える。
さて、Minecraftで計算機を制作するにあたり、10進数計算のほうが信号強度の配線上のメリットや変換コスト、伝搬コストが少ないので、CORDICを10進数に改良する必要がある。
しかし安直に

としても、分解を行うことが不可能である(各桁に入る数字が2パターンしか無く、不足する)
今度は冗長10進数で次のように表そうとしてみる

今度は誤差$|\epsilon|<10^{-n}$程度で分解することができる。しかし$\sin(\ve_k\tt_k)=^?\ve_k\sin\tt_k$が成立しなくなるので、
$$e^{i\tt}=\cos\tt+i\sin\tt=^?R_0\prod_{k=0}^n (1+i\ve_k10^{-k})$$
が使えなくなる。$|1+i\ve_k\tan\tt_k|$が$\ve_k$に依らないためには、$\ve_k\in\{1,-1\}$にせざるを得ないのである。
そこで条件を少し緩める。$c_k$を事前に決めた自然数として、$\theta_k=\arctan(c_k10^{-k})$とする。$x,y$を実数とすると、
$(1+i\ve_k \tan\tt_k)(x+iy)=(x-\ve_k 10^{-k}c_ky)+i(y+\ve_k10^{-k}c_kx)$
は$x,y$のビットシフト及び、実部虚部ともに$c_k+1$回の加減算で実現できる。
次のように$c$を構成すると完全に分解できる。

LUTに格納しておくのは$ \tt_k~~(t=1,2,4),~~\tt_0,$
$$R_0=\prod_{k=0}^n \cos\tt_k$$である。

実装の工夫

さらに線形近似を用いて必要なテーブル数を減らす。
$$ \tt=\sum_{k=0}^n\ve_k\tt_k+\epsilon=\Theta+\epsilon$$
として$\Theta$部分をCORDICにより計算し、$\epsilon$に対して線形近似を行う。
有効桁数をNとして$|\ee|<10^{-N/2}$のとき次のように近似できる:
$\cos\ee=1-\frac12\ee^2=1$
$\sin\ee=\ee.$
これにより

となる。$\ee\sin\Theta,\ee\cos\Theta$は乗算が必要である。CORDICはレジスタのコピーや多数のビットシフトが必要なので、高速なバレルシフタが無いとシフトの律速になってしまう。乗算1回で済む場合はそちらを使うほうが有利である。

となる。今回はラジアンではなくデグリーでの実装を試みたので、実際には
$\tt_k \rm deg=\dfrac {180c_t}{10^{s+1}\pi}$をsについてシフトさせて流用している。

双曲線CORDIC

実は複素数を用いると双曲線CORDICもセットで理解できる。
$j$を分解型複素数$j^2=1$とする。分解型複素数の場合も同様の計算法則が成立する。
$e^{j\theta}=\cosh\tt+j\sinh\tt$
$$\tt_k=\rm arctanh~ 2^{-m_k}$$
$$\tt=\sum_{k=0}^n\ve_k\tt_k+\epsilon~~(\ve_k\in\{1,-1\},~~ 誤差\epsilon)$$
$$e^{ｊ\tt}=\cosh\tt+j\sinh\tt=R_0\prod_{k=0}^n (1+j\ve_k2^{-m_k})$$
$$R_0=\prod_{k=0}^n \cosh\tt_k$$
$$Z_0={R_0}+i0,~~Z_{k+1}=(1+j\ve_k2^{-m_k})Z_k$$
$\cos\tt=Re (Z_{n+1}),~~\sin\tt=Im(Z_{n+1})$
$(1+j\ve_k \tanh\tt_k)(x+jy)=(x+\ve_k 2^{-m_k}c_ky)+j(y+\ve_k2^{-m_k}c_kx)$
ここで注意しなければならないのは、三角関数の場合と違い、
$\tt_k=\rm arctanh ~2^{-k}$
を分解基底に取ると
$$\sum_{k=2}^n\arctan(2^{-k})>\arctan(2^{-1})$$
$$\sum_{k=2}^n\rm arctanh(2^{-k})< arctanh(2^{-1})$$
なので分解が完全でない点である。 命題:基底の完全性 を参照。そこで、途中でkをダブらせて2回同じ分解を行うなどする必要がある。その意味で$m_k$を用意した。
たとえば
$$(m_k)=[1,2,3,4,4,5,6,7,8,9,10,11,12,13,13,14,15,16,17,18,19,20,21,22,23]$$
などとすれば完全である(greedyな算法)