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

525

CORDICアルゴリズム

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

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

記号

絶対値 $r$ ,偏角 $θ$ の複素数を $r e^{i θ} (= r \cos θ + i r \sin θ)$ と表記する。

指数法則(三角関数の加法定理・複素数の乗除算)

$r_{1} e^{i θ_{1}} r_{2} e^{i θ_{2}} = r_{1} r_{2} e^{i (θ_{1} + θ_{2})}$
$(r e^{i θ})^{- 1} = r^{- 1} e^{- i θ}$

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

$\begin{aligned} e^{i θ} = & \prod_{k = 0}^{n} e^{i ε_{k} θ_{k}} \\ = & \prod_{k = 0}^{n} (\cos θ_{k} + i ε_{k} \sin θ_{k}) \\ = & \prod_{k = 0}^{n} \cos θ_{k} (1 + i ε_{k} \tan θ_{k}) \end{aligned}$
ここで $\tan θ_{k} = 2^{- k}$ となるよう $0 < θ_{k} < π / 2$ を選ぶ。そして
$R_{0} = \prod_{k = 0}^{n} \cos θ_{k}$ とする。すると

計算手順の本質

$e^{i θ} = \cos θ + i \sin θ = R_{0} \prod_{k = 0}^{n} (1 + i ε_{k} 2^{- k})$

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

漸化式

$Z_{0} = R_{0} + i 0, Z_{k + 1} = (1 + i ε_{k} 2^{- k}) Z_{k}$
$\cos θ = R e (Z_{n + 1}), \sin θ = I m (Z_{n + 1})$

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

CORDICまとめ

　①事前に計算して表に格納しておく:
$k = 0, . . ., n, θ_{k} = \arctan 2^{- k}, R_{0} = \prod_{k = 0}^{n} \cos θ_{k}$
②入力 $θ$ を非回復型除算によって次のように分解する:
$θ = \sum_{k = 0}^{n} ε_{k} θ_{k} + ϵ (ε_{k} \in {1, - 1}, 誤差 ϵ, | ϵ | < 2^{- n})$
③加減算とビットシフトの繰り返しで複素数の乗算を行い、逐次 $Z_{k}$ を計算する
$Z_{0} = R_{0} + i 0, Z_{k + 1} = (1 + i ε_{k} 2^{- k}) Z_{k}$
$\cos θ = R e (Z_{n + 1}), \sin θ = I m (Z_{n + 1})$

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

10進数化

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

$k = 0, \dots, n, θ_{k} = \arctan 10^{- k}, R_{0} = \prod_{k = 0}^{n} \cos θ_{k}$
$θ = \sum_{k = 0}^{n} ε_{k} θ_{k} + ϵ (ε_{k} \in {1, - 1}, 誤差 ϵ)$

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

$k = 0, \dots, n, θ_{k} = \arctan 10^{- k}, R_{0} = \prod_{k = 0}^{n} \cos θ_{k}$
$θ = \sum_{k = 0}^{n} ε_{k} θ_{k} + ϵ (ε_{k} \in {5, 4, \dots, - 4, - 5}, 誤差 ϵ)$

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

$n = 4 m, k = 4 s + t, t \in {1, 2, 3, 4}$
$(c_{1}, c_{2}, c_{3}, c_{4}) = (4, 2, 2, 1)$
$θ_{0} = \arctan (1)$
$θ_{k} = \arctan (\frac{c_{t}}{10^{s + 1}})$
正当性,誤差評価については「N進数展開の一般化と分解について」で既に説明している。

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

実装の工夫

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

線形近似

$e^{i θ} = (\cos ϵ + i \sin ϵ) (\cos Θ + i \sin Θ) = (\cos Θ - ϵ \sin Θ) + i (\sin Θ + ϵ \cos Θ)$

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

定数削減

$10^{- N / 2} < \frac{c_{t}}{10^{s + 1}} < 10^{- N / 3}$ のとき
$θ_{k} = \arctan (\frac{c_{t}}{10^{s + 1}}) = \frac{c_{t}}{10^{s + 1}} - \frac{1}{3} {(\frac{c_{t}}{10^{s + 1}})}^{3} = \frac{c_{t}}{10^{s + 1}}$

となる。今回はラジアンではなくデグリーでの実装を試みたので、実際には
$θ_{k} d e g = \frac{180 c_{t}}{10^{s + 1} π}$ をsについてシフトさせて流用している。

双曲線CORDIC

実は複素数を用いると双曲線CORDICもセットで理解できる。
$j$ を分解型複素数 $j^{2} = 1$ とする。分解型複素数の場合も同様の計算法則が成立する。
$e^{j θ} = \cosh θ + j \sinh θ$
$θ_{k} = a r c t a n h 2^{- m_{k}}$
$θ = \sum_{k = 0}^{n} ε_{k} θ_{k} + ϵ (ε_{k} \in {1, - 1}, 誤差 ϵ)$
$e^{ｊ θ} = \cosh θ + j \sinh θ = R_{0} \prod_{k = 0}^{n} (1 + j ε_{k} 2^{- m_{k}})$
$R_{0} = \prod_{k = 0}^{n} \cosh θ_{k}$
$Z_{0} = R_{0} + i 0, Z_{k + 1} = (1 + j ε_{k} 2^{- m_{k}}) Z_{k}$
$\cos θ = R e (Z_{n + 1}), \sin θ = I m (Z_{n + 1})$
$(1 + j ε_{k} \tanh θ_{k}) (x + j y) = (x + ε_{k} 2^{- m_{k}} c_{k} y) + j (y + ε_{k} 2^{- m_{k}} c_{k} x)$
ここで注意しなければならないのは、三角関数の場合と違い、
$θ_{k} = a r c t a n h 2^{- k}$
を分解基底に取ると
$\sum_{k = 2}^{n} \arctan (2^{- k}) > \arctan (2^{- 1})$
$\sum_{k = 2}^{n} a r c t a n h (2^{- k}) < a r c t a n h (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な算法)