CORDICアルゴリズムという、乗算なしにsin/cos/sinh/coshを求めるアルゴリズムがあります。乗算器すらないような非常に小規模なプロセッサでsinやcosが必要になった場合などに有用です。

sin/cosを求める場合

戦略

加法定理を思い出します。

$\sin(\alpha + \beta) = \sin\alpha\cos\beta + \cos\alpha\sin\beta$
$\cos(\alpha + \beta) = \cos\alpha\cos\beta - \sin\alpha\sin\beta$

$\cos\beta$ を因数としてくくりだすと、以下のようになります。

$\sin(\alpha + \beta) = \cos\beta(\sin\alpha + \cos\alpha\tan\beta)$
$\cos(\alpha + \beta) = \cos\beta(\cos\alpha - \sin\alpha\tan\beta)$

したがって、「 $\tan\beta$ と $\cos\beta$ が分かっていれば、 $\sin\alpha$ と $\cos\alpha$ から $\sin(\alpha+\beta)$ と $\cos(\alpha+\beta)$ が求められる」ということになります。

$\tan\beta$ が二のべき乗になっていれば、 $\tan\beta$ を掛けている部分は単純なシフトで済みます。したがって、 $\tan\beta$ が二のべき乗になるような $\beta$ をうまく選び、これを足したり引いたりすることで求めたい角度に近づけていくという戦略を取ります。

$\cos\beta$ については、全体にかかっているので、最後にまとめて計算します。実は、この「最後にまとめて計算」するときにかけるべき値は、求めたい角度に近づけていくときの経路によらず事前に計算した値とすることができます。これは、 $\cos\beta$ が $\beta$ の符号によらず同じ値になることを利用しています。求めたい角度に近づけていく際、「 $\beta$ を足す、または足さない」のように求めたい角度に近づけていくのではなく、「 $\beta$ を足す、または $-\beta$ を足す」のように求めたい角度に近づけていきます。このようにすることで、「最後にまとめて計算」するときにかけるべき値が、どう近づけていったかによらずに決まります。

準備

CORDICアルゴリズムでは、以下の値を事前に計算して持っておきます。

$\beta_k = \arctan(1/2^k) (k=0,1,\dots N)$
$M = \prod_{k=0}^{N}\cos(1/2^k)$

計算

CORDICアルゴリズムでは、sin/cosを求めたい角度 $\theta$ を以下のように分解します。

$\theta = \pm\beta_0\pm\beta_1\pm\dots\pm\beta_N+\varepsilon$

ここで、 $\pm$ については、誤差 $\varepsilon$ が十分小さくなるよう、うまく符号を選びます。とはいっても、上位桁から貪欲に選んでいく、という方法が実用上十分です。選んだ符号を $s_k$ とすると、以下のように書けます。

$\theta = s_0\beta_0+s_1\beta_1+\dots+s_N\beta_N+\varepsilon$

すると、以下のように漸化式を順に計算していけば、 $\sin \theta$ や $\cos \theta$ が求まります。

$x_0 = 0.0$
$y_0 = 0.0$
$x_{k+1} = x_k + s_k y_k2^{-k}$
$y_{k+1} = y_k - s_k x_k2^{-k}$
$\sin\theta \simeq Mx_N$ 、 $\cos\theta \simeq My_N$

コード

まずは説明のために分かりやすいC++コードを紹介します。

std::pair<double, double> cordic_sincos( double a ) {
        double angle = 0.0;
        double x = 1.0;
        double y = 0.0;
        double M = 1.0;

        for( int i = 0; i < 53; ++i ) {
                double t = std::ldexp(1.0, -i);
                double b = std::atan(t);
                if( angle > a ) {
                        std::tie(x, y) = std::tuple(x + y * t, y - x * t);
                        angle += -b;
                        M *= std::cos(-b);
                } else {
                        std::tie(x, y) = std::tuple(x - y * t, y + x * t);
                        angle += b;
                        M *= std::cos(b);
                }
                M *= std::cos(b);
        }
        return { y * M, x * M };
}

このコードでは毎回Mを計算していますが、実際には分岐の方向によらず同じ値になるため、事前計算可能です。また、xを初期化する際にMで初期化してしまえば、最後の乗算も不要になります。

それらの最適化を行い、また整数演算だけで計算できるようにしたコードが以下になります。

std::pair<double, double> int64_cordic_sincos( double a ) {
        double angle = 0.0;
        std::int64_t X = 0x4dba76d421af2d34; // Πcos(2^-k)
        std::int64_t Y = 0;
        for( int i = 0; i < 53; ++i ) {
                double t = std::ldexp(1.0, -i);
                double b = std::atan(t);
                if( angle > a ) {
                        std::tie(X, Y) = std::tuple(X + (Y>>i), Y - (X>>i));
                        angle -= b;
                } else {
                        std::tie(X, Y) = std::tuple(X - (Y>>i), Y + (X>>i));
                        angle += b;
                }
        }
        return { Y/0x1.p+63, X/0x1.p+63 };
}

このように、確かに加減算とシフトのみで計算できていることが分かります。 CORDIC法は、固定小数点数で求めたい場合には非常に有力な手法です。

sinh/coshを求める場合

事前に計算する値が $\cos$ や $\arctan$ ではなく $\cosh$ や $\mathrm{arctanh}$ になること、加法定理の符号が一部異なることを除いて、ほぼ同じ……ではありません。

$\beta_k = \arctan(1/2^k) (k=0,1,\dots N)$ はこれらを組み合わせて適当に足し引きすることで入力 $\theta$ を精度良く表すことができるのですが、 $\beta_k = \mathrm{arctanh}(1/2^k) (k=1,2,\dots N)$ はそうではありません。具体的には、例えば0.0付近の数を精度良く表すことができません。 $\mathrm{arctanh}(1/2) - \sum_{k=2}^{\infty}\mathrm{arctanh}(1/2^k) > 0.043$ であり、これより絶対値の小さい値を精度良く表すことができません。

この問題に対しては、 $\beta_k$ の複数使用を許した分解を考えればよいです。例えば、全部二回使って、 $\theta = \pm\beta_1\pm\beta_2\pm\beta_2\pm\beta_3\pm\beta_3\pm\dots\pm\beta_N\pm\beta_N+\varepsilon$ のようにします。参考文献[2]によれば、重複利用するのは三回に一回程度でよいようです。これは、 $\mathrm{arctanh}(1/2^n) - \sum_{k=n+1}^{\infty}\mathrm{arctanh}(1/2^k) \lt \mathrm{arctanh}(1/2^{3n+1})$ が成り立つためです。具体的なコードは以下のようになります。

std::pair<double, double> cordic_sinhcosh( double a ) {
        double angle = 0.0;
        double x = 1.0;
        double y = 0.0;
        double M = 1.0;

        for( int i = 1; i < 53; ++i ) {
                double t = std::ldexp(1.0, -i);
                double b = std::atanh(t);
                if( angle > a ) {
                        std::tie(x, y) = std::tuple(x - y * t, y - x * t);
                        angle += -b;
                        M *= std::cosh(-b);
                } else {
                        std::tie(x, y) = std::tuple(x + y * t, y + x * t);
                        angle += b;
                        M *= std::cosh(b);
                }
                if( i != 1 && i % 3 == 1 ) {
                        if( angle > a ) {
                                std::tie(x, y) = std::tuple(x - y * t, y - x * t);
                                angle += -b;
                                M *= std::cosh(-b);
                        } else {
                                std::tie(x, y) = std::tuple(x + y * t, y + x * t);
                                angle += b;
                                M *= std::cosh(b);
                        }
                }
        }
        return { y * M, x * M };
}

CORDICを逆に使う場合

CORDICを逆に使うと、 $\arctan$ や $\mathrm{arctanh}$ を求めることができます。逆に使うとは、以下のようなことを意味します。

普通に使うときは、分解結果と入力の誤差が0に近くなるように符号を定めつつ、を更新していきました。
- 入力が $\theta$ 、出力が $x, y$ です。
- 初期点(x, y)から $\theta$ 回転させたとき、(x, y)がどこにあるかを計算してことになります。
逆に使うときは、が0に近くなるように符号を定めつつ、を更新していきます。副産物として、分解結果が得られますが、これが求めるものです。
- 入力が $x, y$ 、出力が $\theta$ です。
- 初期点(x, y)から(x', 0)に近づけるように回転したとき、何度回転させたのかを計算していることになります。

なお、sin/cosを求めるときのCORDICを逆に使うと $\sqrt{x^2+y^2}$ を、sinh/coshを求めるときのCORDICを逆に使うと、 $\sqrt{x^2-y^2}$ を、それぞれ求めることができます。この時は、出力は $\theta$ ではなく $x$ になります（点(x, y)を通る円・双曲線のx切片を求めていることになります）。つまり、 $\theta$ は全く参照されないので、平方根を求めるだけであれば $\arctan$ や $\mathrm{arctanh}$ のテーブルは不要です。さらに、sinh/coshを求めるときのCORDICに $x = b + 1/4, y = b - 1/4$ を入力して逆に使うことで、 $\sqrt b$ を求めることもできます。