COordinate Rotation DIgital Computer の略だそうです。 事実上加減算とシフトだけで三角関数を求められるアルゴリズムで、ずいぶん昔からあるらしいですが、最近まで知りませんでした。 三角関数以外にも双曲線関数をはじめとしていろいろなものが求められます。 双曲線関数CORDICについて述べられているサイトが少なかったので、まとめてみました(覚え書きとも言う)。
参考: cordic methods(英語、リンク切れ)
三角関数については他の方が詳しく説明されているので、ここでは流れだけを説明します。 代数的には加法定理の応用で、 \(\cos\beta\) をくくり出して、第2項の \(\tan\) が2の負べき乗になるようにテーブルを変更し、くくり出した \(\cos\beta\) を事前に計算できるよう、条件分岐を「足す・足さない」ではなく、「正・負」に変更したものです。 最後の条件分岐の変更は、除算で言う「復元法(加え戻し法)」と「非復元法(引き放し法)」の違いに相当します。
\(x = \cos\theta,\ y = \sin\theta\) を求めるとして、 \(\theta\) を \(\alpha_0+S_0\beta_0+S_1\beta_1+\dots\) のように分解します。 たとえば、 \(\alpha_0=0,\) \(\beta_0=51.2,\) \(\beta_1=25.6,\) \(\beta_2=12.8,\) \(\beta_3=6.4,\) \(\beta_4=3.2,\) \(\beta_5=1.6,\) \(\beta_6=0.8,\) \(\beta_7=0.4,\) \(\beta_8=0.2,\) \(\beta_9=0.1\) として、 \(\theta\) を10倍して2進数にし、各桁を \(S_k\) に当てはめれば、0.1度単位で \(\theta\) を表すことができます。 そして、 \(x_0=\cos\alpha_0,\) \(y_0=\sin\alpha_0\) として、次の漸化式を必要な精度まで計算することで、 \(\cos\theta\) と \(\sin\theta\) が求まります。
この漸化式で \(\cos S_k\beta_k\) をくくりだします。
ここで、 \(\beta_k\) による分解を、 \(\beta_k\) を \(2^{m-k}\) にするのではなく、 \(\tan\beta_k=2^{m-k}\) になるように分解すれば、つまり、 \(\beta_k=\tan^{-1} 2^{m-k}\) としてしまえば、 \(\tan\beta_k\) に関する乗算がなくなり、単なるシフトになります。
ただし、こうすると \(\theta\) の2進分解ができなくなるので、漸化式を計算する都度 \(\alpha_k+\beta_k\) を \(\theta\) と比べて \(S_k\) を決定しなければなりません。 このとき、もし、 \(S_k\) を「足すか足さないか」、つまり、1 か 0 ではなく、「足すか引くか」、つまり、+1 か -1 にすれば、くくり出した \(\cos S_k\beta_k\) は \(\cos(\pm\beta_k)=\) \(\cos\beta_k\) になり、 \(S_k\) が影響しなくなります。 したがって、 \(M=\cos\beta_0\cdot\:\)\(\cos\beta_1\cdot\:\)\(\cos\beta_2\dots\:\)\(\cos\beta_n\) をあらかじめ計算しておいて最後に乗じれば正しい答が得られます。
これがCORDICの基本的な形で、最後に \(M\) を乗じる以外、乗算がなくなります。 \(x\) と \(y\) の初期値が定数の場合、定数にあらかじめ \(M\) を乗じておけば最後の乗算もなくなります。 \(\beta_k=\tan^{-1}2^{m-k}\) をラジアンで求めておけば、全体の計算はラジアンで行われます。 度で求めておけば度で行われます。 \(2\pi/360\) をかけて変換を行う必要はありません。
初期値として \(\alpha_0=45\makeunit{\deg}\) 度とすることが多いようですが、 \(\alpha_0=0\) としても動作します。 前者の場合は \((\sqrt 2/2,\ \sqrt 2/2)\) が初期値ですが、 \((1,\ 1)\) を初期値として計算しても構いません。 この場合、最後に乗じる定数が変わります。 逆に、 \((M\sqrt 2/2,\ M\sqrt 2/2)\) を出発点にすれば最後の乗算を省けます。 後者の場合は \((1,\ 0)\) が出発点ですが、これも \((M,\ 0)\) を出発点にすれば最後の乗算を省けます。
級数展開と比べると乗算の回数は圧倒的に減りますが、必要な精度を得るための演算回数はそんなに減りません。 したがって、乗算器がないシステムでは圧倒的に有利ですが、乗算も加算も1クロック、といったシステムではかえって遅くなることもあります。 条件分岐も必要なので、条件分岐のコストが高いシステムも不利です。 また、実行時に計算の打ち切り位置を変えられる級数展開と違って、分解能はテーブルによって決まってしまいますから、0近くの有効桁が減ってしまいます。 このため、浮動小数には向きません。 乗算器のない8ビットCPUで、固定小数で値を求めればいい場合などに向きます。
ハードウェアで実現する場合は加算器をパイプライン接続すれば、乗算器ひとつと同程度のハードウェア規模で1クロックにつき \(\sin\) ・ \(\cos\) がひとつずつ求まりますから、面積・速度とも圧倒的に有利です。
CORDICでは \(\cos\theta\) と \(\sin\theta\) が同時に求まるので、これを使って \(\tan\theta=\sin\theta/\cos\theta\) で求めます。 除算が必要になってしまいますが、最後に比を求めることになるので、Mを乗じる必要はありません。
CORDICを逆に使います。 今までは \(\theta\) で \(S_k\) を決めていましたが、スタート地点を \(x_0=1,\ \) \(y_0=t\) として、 \(y_{k+1}\) が \(0\) に近づくように \(S_k\) を決めます。 つまり \(y_k>0\) なら \(s_k=-1\) 、 そうでなければ \(s_k=+1\) です。 このとき、 \(\alpha_0=0\) として積算していくと、 \(-\alpha_n\) が \(\tan^{-1}t\) になっています。 必要なのは角度だけなので、 \(M\) を乗じる必要はありません。
\(\tan^{-1}t\) と同じですが、初期値を \(x_0=X,\) \(y_0=Y\) とします。 \(\tan^{-1}\) だけが欲しい場合は角度だけあればいいので \(M\) の出番はありません。 こちらは \(X=0\) でも \(\tan^{-1}\) が求まります。 CORDICは幾何学的にはベクトルを \(\pm\beta_k\) 回転させて、大きさを \(1/\cos\beta_k\) 倍していることになるので、最終的に求まったベクトルの大きさに \(M\) をかけると \(\sqrt{X^2+Y^2}\) が求まります。 \(y_n\) が \(0\) になるように積算したのですから、最終的なベクトルの大きさは \(x_n\) になるはずです。 つまり、 \(M\cdot x_n\) が元のベクトル \((X,\ Y)\) の大きさ(絶対値)になります。
たとえば \(\sin^{-1}s\) ならば初期値を \(x_0=M,\ \) \(y_0=0\) として、 \(\theta\) の代わりに \(y\) と \(s\) を比較すれば計算できそうな気がしますが、これだとふたつの理由できちんと動かないと思います。
ひとつは計算しているうちに象限をまたいでしまうことがある点。 これは \(y\) の大小比較時に \(x\) もきちんと見て、隣の象限に飛び出してしまっているようなら引き戻す形にすれば対策できます。
もうひとつの理由は、ベクトルの大きさの初期値は \(M<1\) で、漸化式を計算していくうちにベクトルが大きくなっていく点で、 \(y\) が \(s\) より小さいのに、動径の角度は目的の角度よりも大きいことがあります。 この状態で分解を進めてしまうとベクトルが間違った方向に動いてしまい、それを引き戻す余裕はないため、間違った結果になります。 \(\sin\theta=0.65\) 付近を計算させてみると分かります。 この時にベクトルがどう動くかを示したのが次の図です。
円に乗っている点は最終目標地点です。 つまり、この点のY座標が0.65になっています。 求めたいのは \(\sin^{-1}0.65\) ですから、X軸とこの動径がなす角が求める結果になります。 一方、ベクトル \((x,\ y)\) の方は最初はX軸上にあります。 \(x\) の初期値は \(M<1\) なので、円の内側にあることに注意してください。
最初の分解では目標地点よりY座標が小さいので、ベクトルは角度が増える方に動きます。 動かした状態が円まで達していない黒線になります。 次の分解では角度としては目標地点を行き過ぎているので角度が減る方向に動かさなければなりませんが、Y座標は目標地点より小さい(目標地点より下にある)ので、角度を増やす方向に動いてしまいます。 この結果、間違った方向へ動いてしまい(図のY軸に近い灰色の線)、これ以後の分解でリカバリーすることができず、間違った結果が表示されます。
今までうまくいっていたのは、ベクトルの比率が大事で大きさは関係してこない「角度」や、拡大縮小しても値が変わらない「ゼロ」との比較だったからです。その他の場所での比較はうまくいかない、ということです。
\(\sin^{-1}s\) の場合、双曲線CORDICならば \(c=\sqrt{1-s^2}\) を求められますから、ここから \(\tan^{-1}(s/c)\) を求めるのがいい気がします。
双曲線関数でも加法定理が成り立つので、同じようにCORDICを使うことができます。
表面上、三角関数と違うのは x の項が減算ではなく加算になることくらいですが、問題がふたつあります。 ひとつ目の問題は分解に関する問題で、三角関数のときに \(\beta_k=\tan^{-1}2^{m-k}\) と置きなおしましたが、これが許されるのは、
の関係があるからです。 元々は \(\beta_{k+1}=\beta_k/2\) になっていたのを変えたわけですが、 \(k\) が1増えたときに、 \(\beta\) が元の半分以上ないと分解結果に隙間が開いてしまいます。 ところが、双曲線関数では
と、関係が逆になっています。 これは \(\tan\) と \(\tanh\) の曲線の曲がり方が逆だからです。 \(\tanh t\)は、原点と曲線上の点 \((\cosh t,\) \(\sinh t)\) を結ぶ直線の傾きです。 双曲線関数が作る双曲線の漸近線が \(y=\pm x\) ですから、 \(\tanh t\) は漸近線の傾き、つまり \(\pm 1\) の範囲内にあります。 したがって、 \(y=\tanh x\) の漸近線は \(y=\pm 1\) です。 これはX軸に平行な直線です。 一方で三角関数 \(y=\tan x\) の漸近線は原点付近では \(x=\pm\pi/2\) です。 こちらはY軸に平行な直線です。 つまり、 \(\tan\) の傾きはだんだん急になるのに対し、 \(\tanh\) の傾きはだんだん緩やかになります。
三角関数の場合、たとえば \(\beta_0=\tan^{-1}1,\) \(\beta_1=\tan^{-1}2^{-1},\) \(\beta_2=\tan^{-1}2^{-2},\) \(\dots\) として \(S_k=\pm 1\) のすべての組み合わせについて \(\sum S_k\beta_k\) がどの位置にくるかを示したのが下の図です。 横軸の単位は度です。
スキマなく並び、0付近では一部が重複していることが分かります。
双曲線関数で \(\beta_1=\tanh^{-1}2^{-1},\) \(\beta_2=\tanh^{-1}2^{-2},\) \(\beta_2=\tanh^{-1}2^{-3},\) \(\dots\) として、同じことをすると、以下のようになります。 横軸は双曲角。
線は128本あるはずですが、一部にスキマができ、中央のスキマが一番広くなることが分かります。 この問題を解決する必要があります。
もし、ある分解の
\(\beta_k\)
より下位の双曲角を全て足したときに
\(\beta_k\)
よりも小さいならば、
\(\beta_k\)
よりも少しだけ小さい値は分解できないでしょう。
その場合、下位の分解のどこかで2回(もしかすると3回)分解する必要がありますが、なるべくギリギリ条件を満たせる位置で分解した方が有利そうです。
以下にコードの例を示します。
実際にCORDICで双曲線関数を計算する部分(cordic関数)まで入っています。
shrは右シフトで、
\(2^i\)
で割っていますが、固定小数点で演算するなら文字通り右シフトで置き換えられます。
なお、JavaScriptの右シフトは32ビット符号つき整数で行われますので注意が必要です。
テーブルは小数点以下5桁の精度が取れるようにループしています。
このとき、.c(カウンタ)を用意して1をセットしておきます。
テーブルを求めたあと、全ての分解について、より下位の分解の和を求めます。
「より下位の分解」はArray.で簡単に取ることができます。
実際には、通常の2進分解を考えると、16の桁より下は8+4+2+1で合計は15になりますが、これで分解は足りています。
つまり、最下位の1をもう一度足して、現在の分解の値以上ならば大丈夫、ということです。
和はArray.で求めていますが、reduceはありがたいことに初期値を指定できるので、Array.で最後の要素を持ってきて、その.xを初期値にしています。
最後の分解についてはsliceの結果が空になるので、reduceは初期値を返します。
つまり、最後の要素の.xを返してくるので、その後のif文はいい具合に飛ばしてくれます。
分解が足りなくなると分かった場合は、findLastで分解を満たせる一番小さい値を見つけ、その.cをインクリメントしておきます。
分解の系列が変わるとMが変わりますので、最後にMを求め直しておきます。
v.cが1より大きい場合は1回で分解が足りると分かっても、必ずその回数だけ分解を行います。
そうしないとMがずれてしまいます。
Mは必ずv.c回乗じるのでcoshのv.c乗をMに乗じていくことになります。
このコードで
\(0.5\times 10^{-5}\)
の精度で分解できることを確認してあります。
ふたつ目の問題は定義域に関する問題で、 \(S_k\beta_k\) による分解には上限があり、これは双曲線CORDICの入力の上限になります。 上の図では \(\pm 1\) ぐらいから上は表現できないことが分かります。 これは \(\tanh^{-1}\) の定義域が1未満なためで、シフト演算を有効に活用しようとすると \(\tanh^{-1}0.5\) よりも大きな \(\beta_k\) を用意できないことになります。 三角関数でも上限はあるのですが、周期性があるのであまり困らないのです。 双曲線関数の場合は(実数については)周期がないので、XY平面上で漸近線に近い値を扱う際には注意が必要です。 \(\beta_0\) を何回も足し算するという手もないことはないですが・・・。
\(t\) が0.76を超えると答が1を超えるようになるので工夫が必要です。 \(\log\) と \(\tanh^{-1}\) は表裏一体ですが、 \(\log\) は正しく計算できる範囲を超えた場合に対数法則で細工ができるので、 一度 \(\log\) に直して、対数法則を適用すると、
となりますが、双曲線関数のページに書いてあるように、 \(\log\) は
右辺の \(\tanh^{-1}\) の分数は、 \(t\) が正で定義域の中に入っていれば、言い換えると \(0<t<1\) ならば、分母は常に正ですが、分子は \(n\) を大きくしていくとどこかで負になります。 その直前を狙って \(n\) を決めて計算するといいのかな。 この分子分母は加減算とシフトだけで計算できます。 右辺の \(\tanh^{-1}\) は次の \(\tanh^{-1}Y/X\) で計算するので、除算はいりません。 \((n\log_e 2)/2\) は \(\log_e 2\) が定数ですから、単なる乗算です。
\(t\) が-0.76以下になったときは各自の課題とします。
\(X^2-Y^2\) を \(A^2\) と置くと、 \(M\) 倍の差を考えなければ、ベクトルは双曲線 \(x^2-y^2=A^2\) に沿ってX軸上まで、つまり \(y=0\) まで動くことになります。 したがって、最終的には \(Mx=A\) になり、 \(A^2=X^2-Y^2\) より \(A=\sqrt{X^2-Y^2}\) なので、 \(Mx=A=\;\)\(\sqrt{X^2-Y^2}\) が求まったことになります。
\(Y/X\) が1に近いときは、
から、
なので、先ほどと同様、 \(x+y-\;\)\((x-y)2^n\) が負になる直前を狙って、
という変数変換をして計算するといいのかな。
双曲線CORDIC自体は \(x\) と \(y\) の符号は問わないのですが、上の式は \(x-y\) が0に近いことを前提としています。 どちらかの符号が逆になると \(x+y\) が0に近くなるので、 \(2^n\) を付ける位置が変わってきます。 \(\tanh^{-1}\) は奇関数と分かっていますので、 符号を取っ払って計算して、後で符号を調整した方が楽かもしれません。
変数変換したときに \(A^2\) がどうなるか計算すると、
と、 \(4\cdot 2^n\;\)\(=\;\)\(2^{n+2}\) 倍になることが分かります。 実際に欲しいのはこの平方根ですから、 \(n\) を偶数にして \(2m=n\) と置けば、 \(A'^2\;\)\(=\;\)\(2^{n+2}A^2\;\)\(=\;\)\(2^{2m+2}A^2\;\)\(=\;\)\(2^{2(m+1)}A^2\) より、 \(A>0,\;\)\(A'>0\) ならば \(A'=\;\)\(2^{m+1}A\) の関係があり、求まるのは \(A'\) でここから \(A\) を求めたいので、 \(A=\;\)\(2^{-(m+1)}\;\)\(A'\) です。 これは右シフトだけで計算できます。
ここから先ほど説明した \(\tanh^{-1}s/c\) で \(t=\;\) \(\sinh^{-1}s\;\) \(=\cosh^{-1}c\) が求まります、たぶん。
で求めればいいのかな。 つまり、 \(x\) が正になるまで \(\log_e2\) を足して、足した回数だけ結果を右シフトすればOK。
符号を逆さにして結果の逆数を取ることもできますが、それだと割り算が必要になってしまいます。 あ、ということはもしかして逆数もCORDICで計算できるのか。
この場合、範囲外は \(\log_e(x\cdot 2^n)\;\)\(=\;\)\(\log_e x+\;\)\(n\log_e 2\) を使って求めることになるのでしょう。 特に、 \(x\) が1より小さくなると \(\tanh^{-1}\) の入力・出力とも負になることに注意する必要があるでしょう。
\(b\) と定数の加減算だけで求めることを考えると、 \(x=b+m,\;\)\(y=b+n\) と置いて、
となるので、右辺と左辺を見比べて、
です。 これを解くと \(m=1/4,\;\)\(n=-1/4\) になります。 したがって、 \(x=b+1/4,\;\)\(y=b-1/4\) からCORDICを逆に使うと、 \(mx\) が \(\sqrt b\) になります。
\(b\) が2を超えると \(\tanh^{-1}\;\)\((b-1/4)/\;\)\((b+1/4)\) が1を超えるようになるため、 \(4^n\) で割るなどして対応する必要があります。 0に近いほうは0.03くらいまで。
05 Mar 2026: MathJax化、内容も若干見直し、特に \(\sqrt{x^2-y^2}\) は \(\tanh^{-1}(Y/X)\) と共通の項目にする
07 Oct 2012: \(\tanh^{-1}t\) で \(t\) が1に近い場合について追記
06 Oct 2012: 新規作成
ご意見・ご要望の送り先は あかもず仮店舗 の末尾をご覧ください。
Copyright (C) 2012-2026 akamoz.jp