CORDICアルゴリズムという、乗算なしにsin/cos/sinh/coshを求めるアルゴリズムがあります。 乗算器すらないような非常に小規模なプロセッサでsinやcosが必要になった場合などに有用です。
sin/cosを求める場合
戦略
加法定理を思い出します。
を因数としてくくりだすと、以下のようになります。
したがって、「と
が分かっていれば、
と
から
と
が求められる」ということになります。
が二のべき乗になっていれば、
を掛けている部分は単純なシフトで済みます。
したがって、
が二のべき乗になるような
をうまく選び、これを足したり引いたりすることで求めたい角度に近づけていくという戦略を取ります。
については、全体にかかっているので、最後にまとめて計算します。
実は、この「最後にまとめて計算」するときにかけるべき値は、求めたい角度に近づけていくときの経路によらず事前に計算した値とすることができます。
これは、
が
の符号によらず同じ値になることを利用しています。
求めたい角度に近づけていく際、「
を足す、または足さない」のように求めたい角度に近づけていくのではなく、「
を足す、または
を足す」のように求めたい角度に近づけていきます。
このようにすることで、「最後にまとめて計算」するときにかけるべき値が、どう近づけていったかによらずに決まります。
準備
CORDICアルゴリズムでは、以下の値を事前に計算して持っておきます。
計算
CORDICアルゴリズムでは、sin/cosを求めたい角度を以下のように分解します。
ここで、については、誤差
が十分小さくなるよう、うまく符号を選びます。
とはいっても、上位桁から貪欲に選んでいく、という方法が実用上十分です。
選んだ符号を
とすると、以下のように書けます。
すると、以下のように漸化式を順に計算していけば、や
が求まります。
、
コード
まずは説明のために分かりやすい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を求める場合
事前に計算する値がや
ではなく
や
になること、加法定理の符号が一部異なることを除いて、ほぼ同じ……ではありません。
はこれらを組み合わせて適当に足し引きすることで入力
を精度良く表すことができるのですが、
はそうではありません。
具体的には、例えば0.0付近の数を精度良く表すことができません。
であり、これより絶対値の小さい値を精度良く表すことができません。
この問題に対しては、の複数使用を許した分解を考えればよいです。
例えば、全部二回使って、
のようにします。
参考文献[2]によれば、重複利用するのは三回に一回程度でよいようです。
これは、
が成り立つためです。
具体的なコードは以下のようになります。
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を逆に使うと、や
を求めることができます。
逆に使うとは、以下のようなことを意味します。
- 普通に使うときは、分解結果
と入力
の誤差
が0に近くなるように符号
を定めつつ、
を更新していきました。
- 入力が
- 初期点(x, y)から
- 入力が
- 逆に使うときは、
が0に近くなるように符号
が得られますが、これが求めるものです。
- 入力が
- 初期点(x, y)から(x', 0)に近づけるように回転したとき、何度回転させたのかを計算していることになります。
- 入力が
なお、sin/cosを求めるときのCORDICを逆に使うとを、sinh/coshを求めるときのCORDICを逆に使うと、
を、それぞれ求めることができます。
この時は、出力は
ではなく
になります（点(x, y)を通る円・双曲線のx切片を求めていることになります）。
つまり、
は全く参照されないので、平方根を求めるだけであれば
や
のテーブルは不要です。
さらに、sinh/coshを求めるときのCORDICに
を入力して逆に使うことで、
を求めることもできます。
参考文献
- [1] CORDIC - You / がらくた箱
- [2] cordic methods
