home / junkbox / cordic CORDICアルゴリズム

　COordinate Rotation DIgital Computer の略だそうです。 事実上加減算とシフトだけで三角関数を求められるアルゴリズムで、ずいぶん昔からあるらしいですが、最近まで知りませんでした。 三角関数以外にも双曲線関数をはじめとしていろいろなものが求められます。 双曲線関数CORDICについて述べられているサイトが少なかったので、まとめてみました（覚え書きとも言う）。

参考： cordic methods（英語）

home / junkbox / cordic 三角関数

　三角関数については他の方が詳しく説明されているので、ここでは流れだけを説明します。 代数的には加法定理の応用で、cosβをくくり出して、第2項のtanが2の負べき乗になるようにテーブルを変更し、くくり出したcosβを事前に計算できるよう、条件分岐を「足す・足さない」ではなく、「正・負」に変更したものです。 最後の条件分岐の変更は、除算で言う「復元法（加え戻し法）」と「非復元法（引き放し法）」の違いに相当します。

　x = cosθ, y = sinθを求めるとして、θを α₀ + S₀β₀ + S₁β₁ + ... のように分解します。 たとえば、 α₀ = 0 , β₀ = 51.2 , β₁ = 25.6 , β₂ = 12.8 , β₃ = 6.4 , β₄ = 3.2 , β₅ = 1.6 , β₆ = 0.8 , β₇ = 0.4 , β₈ = 0.2 , β₉ = 0.1 , として、θを10倍して2進数にし、各桁を S_k に当てはめれば、0.1度単位でθを表すことができます。 そして、 x₀ = cosα₀ , y₀ = sinα₀ として、次の漸化式を必要な精度まで計算することで、cosθとsinθが求まります。

　この漸化式で cos S_kβ_k をくくりだします。

　ここで、β_k による分解を、β_k を 2^m-k にするのではなく、tanβ_k = 2^m-k になるように分解すれば、つまり、β_k = tan^-1 2^m-k としてしまえば、tanβ_k に関する乗算がなくなり、単なるシフトになります。

ただし、こうするとθの2進分解ができなくなるので、漸化式を計算する都度α_k +β_k をθと比べて S_k を決定しなければなりません。 このとき、もし、S_k を「足すか足さないか」、つまり、1 か 0 ではなく、「足すか引くか」、つまり、+1 か -1 にすれば、くくり出した cos S_kβ_k は cos (±β_k ) = cosβ_k になり、S_k が影響しなくなります。 したがって、M = cosβ₀・cosβ₁・cosβ₂・・・cosβ_n をあらかじめ計算しておいて最後に乗じれば正しい答が得られます。

これがCORDICの基本的な形で、最後に M を乗じる以外、乗算がなくなります。 x と y の初期値が定数の場合、定数にあらかじめ M を乗じておけば最後の乗算もなくなります。 β_k = tan^-1 2^m-k をラジアンで求めておけば、全体の計算はラジアンで行われます。 度で求めておけば度で行われます。 2π/ 360 をかけて変換を行う必要はありません。

　初期値として α₀ = 45度 とすることが多いようですが、α₀ = 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がひとつずつ求まりますから、面積・速度とも圧倒的に有利です。

●tanθ

　CORDICでは cosθと sinθが同時に求まるので、これを使って tanθ= sinθ/ cosθで求めます。 除算が必要になってしまいますが、最後に比を求めることになるので、Mを乗じる必要はありません。

●tan^-1 t

　CORDICを逆に使います。 今まではθで S_k を決めていましたが、スタート地点を x₀ = 1 , y₀ = t として、y_k+1 が 0 に近づくように S_k を決めます。 つまり y_k > 0 なら s_k = -1 、 そうでなければ s_k = +1 です。 このとき、α₀ = 0 として積算していくと、-α_n が tan^-1 t になっています。 必要なのは角度だけなので、M を乗じる必要はありません。

●tan^-1 ( Y / X ) と ベクトルの大きさ

　tan^-1 t と同じですが、初期値を x₀ = X , y₀ = Y とします。 tan^-1 だけが欲しい場合は角度だけあればいいので M の出番はありません。 こちらは X = 0 でも tan^-1 が求まります。 CORDICは幾何学的にはベクトルを ±β_k 回転させて、大きさを 1 / cosβ_k 倍していることになるので、最終的に求まったベクトルの大きさに M をかけると sqrt ( X² + Y² ) が求まります。 y_n が 0 になるように積算したのですから、最終的なベクトルの大きさは x_n になるはずです。 つまり、M x_n が元のベクトル ( X , Y ) の大きさ（絶対値）になります。

●sin^-1 s ・ cos^-1 c

　冒頭の参考サイトでは x₀ = M 、y₀ = 0 として、θの変わりに y と t を比較することで計算しています。 これだとなんとなく動くような動かないような気がするので、c = sqrt ( 1 - s² ) を求めて、tan^-1 s/c を求めた方がいい気がします。 c は双曲線CORDICで求まります。 c から s を求めれば cos^-1 c が求められます。

home / junkbox / cordic 双曲線関数

　双曲線関数でも加法定理が成り立つので、同じようにCORDICを使うことができます。

　表面上、三角関数と違うのは x の項が減算ではなく加算になることくらいですが、問題がふたつあります。 一つ目は分解に関する問題で、三角関数のときに β_k = tan^-1 2^m-k と置きなおしましたが、これが許されるのは、

の関係があるからです。 元々は β_k+1 = β_k / 2 になっていたのを変えたわけですが、k が1増えたときに、βが元の半分以上ないと分解結果に隙間が開いてしまいます。 ところが、双曲線関数では

と、関係が逆になっています。 これは tan と tanh の曲線の曲がり方が逆だからです。 tanh t は、原点と曲線上の点 ( cosh t , sinh t ) を結ぶ直線の傾きです。 双曲線関数が作る双曲線の漸近線が y = ±x ですから、tanh t は漸近線の傾き、つまり ±1 の範囲内にあります。 したがって、y = tanh x の漸近線は y = ±1 です。 これはX軸に平行な直線です。 一方で三角関数 y = tan x の漸近線は原点付近では x =±π/ 2 です。 こちらはY軸に平行な直線です。 つまり、tan の傾きはだんだん急になるのに対し、tanh の傾きはだんだん緩やかになります。

　三角関数の場合、たとえば β₀ = tan^-1 1 , β₁ = tan^-1 2^-1 , β₂ = tan^-1 2^-2 , ・・・ として S_k = ±1 のすべての組み合わせについて ΣS_kβ_k がどの位置にくるかを示したのが下の図です。 横軸の単位は度です。

スキマなく並び、0付近では一部が重複していることが分かります。

　双曲線関数で β₀ = tanh^-1 2^-1 , β₁ = tanh^-1 2^-2 , β₂ = tanh^-1 2^-3 , ・・・ として、同じことをすると、以下のようになります。 横軸は双曲角。

線は128本あるはずですが、一部にスキマができ、中央のスキマが一番広くなることが分かります。 この問題を解決する必要があるのですが、冒頭のサイトでは3回に1回は分解を2回行うことで解決しています。

　もうひとつは定義域に関する問題で、S_kβ_k による分解には上限があり、これは双曲線CORDICの入力の上限になります。 上の図では±1ぐらいから上は表現できないことが分かります。 三角関数でも上限はあるのですが、周期性があるのであまり困らないのです。 双曲線関数の場合は（実数については）周期がないので、XY平面上で漸近線に近い値を扱う際には注意が必要です。

●cosh t ・ sinh t

　( M , 0 ) を起点とした双曲線CORDICで求まりますが、t が1以上になる場合は倍角の公式 cosh ( 2 t ) = 2 cosh² t - 1 、 sinh ( 2 t ) = 2 sint h cost h などを使って求めることになるのかな。

●tanh t

　双曲線CORDICでは cosh t と sinh t が同時に求まるので、tanh t = sinh t / cosh t で求めます。 M を乗じる必要はありません。

●tanh^-1 t

　x = 1 , y = t として双曲線CORDICを逆に使います。 必要なのは双曲角なので、Mは不要です。 t が 0.76 を超えると答が1を超えるようになるので、工夫が必要です。 log と tanh^-1 は表裏一体ですが、 log は正しく計算できる範囲を超えた場合に対数法則で細工ができるので、 一度 log に直して、対数法則を適用してから log を tan^-1 を使って計算しなおすと、

tanh^-1 t = 0.5 log { ( 1 + t ) / ( 1 - t ) }

= 0.5 { log ( 1 + t ) - log ( 1 - t ) }

= 0.5 [ log ( 1 + t ) - log { ( 1 - t ) 2ⁿ / 2ⁿ } ]

= 0.5 [ log ( 1 + t ) - log { ( 1 - t ) 2ⁿ } + log 2ⁿ ]

= 0.5 [ log ( 1 + t ) - log { ( 1 - t ) 2ⁿ } + n log 2 ]

= 0.5 log [ ( 1 + t ) / { ( 1 - t ) 2ⁿ } ] + 0.5 n log 2

= tanh-1	( 1 + t ) / { ( 1 - t ) 2n } - 1	+ 0.5 n log 2
	( 1 + t ) / { ( 1 - t ) 2n } + 1

= tanh-1	{ ( 1 + t ) - ( 1 - t ) 2n } / { ( 1 - t ) 2n }	+ 0.5 n log 2
	{ ( 1 + t ) + ( 1 - t ) 2n } / { ( 1 - t ) 2n }

= tanh-1	( 1 + t ) - ( 1 - t ) 2n	+ 0.5 n log 2
	( 1 + t ) + ( 1 - t ) 2n

という関係があるので、( 1 - t ) 2ⁿ　が1〜4くらいになるように n を決めて計算するといいのかな。 tanh^-1の分子分母は加減算とシフトだけで計算できます。 tanh^-1自体は次の tanh^-1 Y / X で計算するので、除算はいりません。 0.5 n log 2 は log 2 が定数ですから、単なる乗算です。

　t が-0.76以下になったときは各自の課題とします。

●tanh^-1 Y / X

　x = X , y = Y として双曲線CORDICを逆に使います。 Y / X が 1 に近いときは、

tanh-1 t = tanh-1	( 1 + t ) - ( 1 - t ) 2n	+ 0.5 n log 2
	( 1 + t ) + ( 1 - t ) 2n

から、

tanh-1 y/x = tanh-1	( 1 + y/x ) - ( 1 - y/x ) 2n	+ 0.5 n log 2
	( 1 + y/x ) + ( 1 - y/x ) 2n

= tanh-1	x { ( 1 + y/x ) - ( 1 - y/x ) 2n }	+ 0.5 n log 2
	x { ( 1 + y/x ) + ( 1 - y/x ) 2n }

= tanh-1	x ( 1 + y/x ) - x ( 1 - y/x ) 2n	+ 0.5 n log 2
	x ( 1 + y/x ) + x ( 1 - y/x ) 2n

= tanh-1	( x + y ) - ( x - y ) 2n	+ 0.5 n log 2
	( x + y ) + ( x - y ) 2n

なので、( x - y ) 2ⁿ が x + y の 0.5〜2倍程度になるように n を決めればいいのかな。

●sinh^-1 s , cosh^-1 c

　cosh² c - sinh² s = 1 より cosh c = sqrt ( 1 + sinh² s ) 、これは三角関数CORDICで求まります。 sinh s = sqrt ( cosh² c - 1 ) 、これは双曲線CORDICで求まります。 ここから tanh^-1 s / c で sinh^-1 s , cosh^-1 c が求まります、たぶん。

●e^x

　cosh x + sinh x = e^x から求めます。

●log x

　log x = 2 tanh^-1 { ( x - 1 ) / (x + 1 ) } から求めます。 どうしてこうなるかは 双曲線関数 に書いてあります。 tanh^-1 Y / X を使えば割り算する必要はありません。 上限は8くらい、下限は0.12くらい。 範囲外は log ( x 2ⁿ ) = log x + n log 2 を使って求めることになるのでしょう。

●sqrt b

　x² - y² = a² を満たす x , y を選んで y が 0 になるように双曲線CORDICを逆に使うと、 双曲線と座標系 に書いたように、双曲線 x² - y² = a² を y = 0 に向かって動くことになるので、最後にMを乗じると x = a , y = 0 になります。 a² を b で置き換えると、a = sqrt ( b ) ですから、b の 平方根が求まったことになるのですが、問題はこのような x , y を簡単に決められるか、です。

　bと定数の加減算だけで求めることを考えると、x = b + M , y = b + N と置いて、

( b + M )² - ( b + N )² = b

b² + M² + 2 M b - b² - N² - 2 N b = b

M² + 2 M b - N² - 2 N b = b

2 b ( M - N ) + M² - N² = b

となるので、右辺と左辺を見比べて、

2 ( M - N ) = 1

M² = N²

です。 これを解くと M = 1 / 4 , N = - 1 / 4 になります。 したがって、 x = b + 1 / 4 , y = b - 1 / 4 からCORDICを逆に使うと、M x が sqrt ( b ) になります。

bが2を超えると tanh^-1 ( b - 1/4 ) / ( b + 1/4 ) が 1 を超えるようになるため、4ⁿ で割るなどして対応する必要があります。 0に近いほうは0.03くらいまで。

●sqrt ( x² - y² )

　( x , y ) から双曲線CORDICを逆に使うと、M x が sqrt ( x² - y² ) になります。 たとえば sin²θ = 1 - cos²θ ですから、x = 1 , y = cosθ からはじめれば cosθ から sinθ を求めることができます。 同様に sinθ から cosθ を求めることもできますし、 cosh² u - sinh² u = 1 より sinh² u = cosh² u - 1 ですから、( x = cosh u , y = 1 ) から始めれば、sinh u が求まります。 いずれも結果は必ず正になりますので、符号については別に考える必要があります。

　結果が0に近くなる条件では tanh が1を超えるようになり、正しい結果が得られなくなります。 そこで、tanh^-1 のときの変形

tanh-1 y/x = tanh-1	( x + y ) - ( x - y ) 2n	+ 0.5 n log 2
	( x + y ) + ( x - y ) 2n

から、( u = ( x + y ) + ( x - y ) 2ⁿ , v = ( x + y ) - ( x - y ) 2ⁿ ) を開始点として選ぶと、求まるのは w = sqrt ( u² - v² ) ですから、

w² = u² - v²

= { ( x + y ) + ( x - y ) 2ⁿ }² - { ( x + y ) - ( x - y ) 2ⁿ }²

= ( x + y )² + 2 ( x + y ) ( x - y ) 2ⁿ + { ( x - y ) 2ⁿ }² - ( x + y )² + 2 ( x + y ) ( x - y ) 2ⁿ - { ( x - y ) 2ⁿ }²

= 4 ( x + y ) ( x - y ) 2ⁿ

= 4・2ⁿ ( x² - y² )

となって、 a² = x² - y² と比べると w² = 4・2ⁿ a² となっているので、 n を偶数 2m にしておけば w から簡単に a = w / 2^m+1 と a が求まることが分かります。

home / junkbox / cordic

Copyright (C) 2012 You SUZUKI

$Id: cordic.htm,v 1.4 2014/03/29 10:43:07 you Exp $