CORDICアルゴリズムを導き出す

CORDICアルゴリズムは組み込み系の数値計算に携わっていると随所でその名を見ます。

むろん、ワークステーション系の恵まれた環境であれば全てはライブラリ化されており、こんなものにはお目に掛かりません。手厚い支援を受けられない組み込みならではの事と言えます。

CORDICアルゴリズムは「実装が簡単でプログラムも単純である」というのがうたい文句ではあるのですが、一方で、CORDICアルゴリズムが正しく動作する理由というのは、プログラムを子細に読んでもわかりません。ある意味、FFTと同じくらい元の数学から乖離した手順を踏んでいると言ってよいでしょう。

以下では数学的背景から順を追って説明し、三角関数のCORDICアルゴリズムを導き出します。是非、このアルゴリズムがいかに巧妙であるかを味わい、楽しんでください。

複素平面と回転子

最初に複素平面の事を考えてみましょう。CORDICアルゴリズムは平面上の変形操作を利用したアルゴリズムです。ですので平面から話を始める方がすっきりするというものです。

よく使う複素平面は図1のように、横軸に実部、縦軸に虚部を配置したものです。図1はさらにこの平面の上に、点(1,0)を基点とし、(0,0)を回転中心としてコンパスを使って半径1の円弧を書き込んだものです。このとき、(0,0)から実軸に対して角θの直線を定規で引いて円弧と交わる点Pを作図します。このとき、その点Pの座標が(cos(θ), sin(θ))です。

図1

次に少し視点をひねってみましょう。この図をながめるとき、我々は上に書いたようにコンパスと定規による作図を想定しました。しかし上の図は、実軸上の点(1,0)を、回転によって点Pまで移動させる操作と考えることもできます。その場合、操作はコンパスと定規による作図ではなく、数式による計算ということになります。

さて、上の座標は複素平面ですので点(1,0)は、複素数で1+0iと表すことができます。そしてこれをPまで移動させるには、次のような演算を行います。

(1+0i)*(cosθ+i*sinθ)=cosθ+i*sinθ  ……… （式1）

ここで、cosθ+i*sinθは、複素平面上の点を原点を中心に角θだけ反時計回りに回転させる複素回転子です。

このように、複素数一つを掛けることで、複素平面上の任意の点を好きな角度だけ原点を中心に回転させてその座標を計算することができます。

回転子群からsin(x)とcos(x)を求める

与えられたxに対して、sin(x)とcos(x)を求めるにはどうすればいいでしょうか。

一つには級数展開したsin(x)、cos(x)を真面目に計算していく方法が考えられます。もう一つ考えられるのは既知の複素回転子を使用してxに等しい角の回転を合成する方法です。例えば次のような表を考えてください。

表1

表1のように、角が45度、22.5度、11.25度…と次々に1/2になっていく回転子の表をあらかじめ持っているとしましょう（θと回転子だけでいい）このような場合、表中の回転子を組み合わせて適用することで、与えられたxに最も近いθを作り出しすことができます。

例えば、xが66.5度だとしましょう。この場合、45+22.5=66.5ですので、(1,0)を回転させたときのPの座標は

(1+0i)*(0.707+0.707i)*(0.924+0.383i) ……… （式2）

を計算すれば求めることができます。すなわちθ=65度のときのPの座標が計算で分かります。これは図2のようになります。

図2

図2で分かるように、この場合の回転子は二つの回転子の合成ですが、考え方を変えると、「まず45度動かして、つぎに22度動かしている」ともとれます。その場合、点の座標を表す複素数に対して次々に回転子である複素数を掛けていくことになります。そして、式にすると一回の演算は次のようになります。

(p+qi)*(r+si) = (pr – qs) + (ps + qr)i ……… （式3）

このように、一回の複素乗算あたり4回の実数乗算が必要になります。これを、θの誤差が十分小さくなるまで繰り返していけば、与えられたxに対する十分精度の高いsin(x)、cos(x)を求めることができます。

スケーリングして乗算を減らす

複素平面上で既知の角度の回転を繰り返していくことで、与えられた角xに対して十分精度の高いsin(x)とsin(y)を求めることができます。しかしながら、一回の繰り返しあたり4回の実数乗算では演算コストが大きすぎます。そこで、これを削ってみましょう。

削るための方法はスケーリングです。先に挙げたθと回転子の表を元に、回転子にスケーリング･ファクタScを掛けることにします。ここでScは1/cosθです。そうすると、以下のように修正回転子の表をえることができます。

表2

修正回転子は、スケーリングによって実部が1に固定されています。その結果、上の式の回転は修正回転子によって次のように変形されます。

(p+qi)*(1+si) = ( ( p – qs ) + ( ps + q )i ) / Sc ……… （式4）

Scによる除算が増えています！しかし、ここでは一旦そこから目をそらしましょう。そうすると、乗算の数は4から2に減っていることがわかります。Scによる除算については後でもう一度考えます。

この変形を座標上で考えてみましょう。式4からあえてScによる補正を取り除くと、座標上での点の動きは図3のようになります。これまで綺麗な円弧を描いていた回転子は、スケーリングによって三角形に置き換わっています。修正回転子は実部が常に1であるため、点(1,0)と回転した点と原点で三角形を作ると、必ず直角三角形になります。

図3

そして、図3を見れば一目瞭然ですが、原点からPへのベクトルの長さは修正回転子を適用する度に長くなっていきます。それがScの副作用です。これがScによる補正が必要な理由です。

三角形の高さを1/(2^n)にして乗算を減らす

スケーリング･ファクタの導入で乗算を減らすことができました。これをさらに減らしてみましょう。

表2で、修正回転子の虚部は1,0.414,0.198となっています。これを思い切って、1, 0.5, 0.25 … 1/2^nに変更するとどうなるでしょうか。

表3

こうすると、上の回転の式4は次のようになります。

(p+qi)*(1+(1/2^n)i) = ( ( p – q(1/2^n) ) + ( p(1/2^n) + q )i ) / Sc ……… （式5）

なにやら複雑になったように思えますが、ここで思い直してみましょう。1/2^nは右シフト演算です。

つまり、ここまでの変形によって、繰り返し一回あたり乗算4つの演算だったのが、加減算とシフト演算（およびスケーリング処理）にまで演算負荷を下げることができました。

折りたたんでスケーリングを解決する

「与えられた角xと同じ角度を、あらかじめ角が分かっている回転子の組み合わせで合成する」という考え方で、複素座標上での点を動かすことをこれまで考えてきました。

ここまでの手順の変形により、点Pを位置(1,0)から始めて既知の回転子を適用していくことにより、与えられた角xに近づけていくことができます。回転子の適用一回当たり乗算が4回必要ですが、スケーリング・ファクタの導入と、虚部の変形によって乗算4回は加減算とシフト演算だけになりました。ただし、スケーリング･ファクタを最後に取り除く必要があります。

スケーリング･ファクタの取り除きに必要な演算はどのくらいになるでしょうか。これを見積もることは簡単です。角xを作り出す（あるいは最も近い角を作り出す）ためにk個の回転子が必要な場合、必要な補正演算は次のようになります。

• k個の回転子それぞれのスケーリング・ファクタを取り出して、それら全ての積を求める
• 求めた積の逆数を求める

つまり、乗算k回と除算1回になりました。元々は乗算4k回でしたからこれでも随分軽くなりましたが、ここでもう一声がんばってみましょう。

上の方法では、必要な回転子だけを選び出して与えられた角xを合成していました。これをやめて表にある回転子全部を使用することにします。例えば、表3のインデックスが0から14までの計15個であれば、15個の回転子全部を使うことにします。一見無駄に見えるこの方法のメリットはScによる補正値の計算を事前に済ませる事ができる点です。なにしろ回転子は全部使います。ですから、コーディング時に全てのScの積とその逆数を計算しておけばいいのです。これで、演算中のScの補正は乗算一回ですむようになりました。そして、もう一度頭をひねることで、そもそもScの補正を実行中に行う必要は無いことに気がつきます。回転を始める点の初期値を(1,0)ではなく、(1/S, 0)とすればいいのです。ここでSは全てのScの積です。

一方、全ての回転子をつかって反時計回りに回転を行えば、求める角xを追い越してしまうことは自明です。そこで、次のような方策を採ります。

• インデックス0から始める
• 合成した角がxより小さければ、回転子をつかって反時計回りに回転させる。
• 合成した角がxより大きければ、回転子をつかって時計回りに回転させる。
• これを全てのインデックスに昇順に適用する

この方法を使うと、回転による点の動きは図3のように折り返しを含むものになります。

図4

また、今回は時計回りの回転が入ってきます。この計算は(式5)を書き換えて次の式で求めることができます。

(p+qi)*(1-(1/2^n)i) = ( ( p + q(1/2^n) ) + ( – p(1/2^n) + q )i ) / Sc ……… （式6）

sin(x)、cos(x)を求めるCORDICアルゴリズム

以上の考察を元に導き出されたCORDICアルゴリズムはおおむね以下のようになります。

sincos( x )
  p = 1/S
  q = 0
  angle = 0

  for ( i = 0; i< MAXINDEX; i++ )

    if x > angle

      p_next = p – q >> i

      q_next = q + p >> i

      angle += theta[i]

    else

      p_next = p + q >> i

      q_next = q – p >> i

      angle -= theta[i]

    end if

    p = p_next

    q = q_next

  end for

  cos = p

  sin = q

よく言われているとおり、乗除算を用いずに加減算とシフト演算のみでsin(x), cos(x)を求めることができます。

収束するのか？

さて、こうしてCORDICアルゴリズムを求めることができました。求めたアルゴリズムは乗算を含まず、加減算とシフト演算だけで実装できます。また事前に必要な定数も1/Sと、theta[n]だけで、非常にコンパクトです。

しかし、そもそもこのアルゴリズムは収束するのでしょうか。

ポイントは与えられたxに対して十分近いthetaを合成できるか否かです。thetaさえ合成できれば回転アルゴリズムには問題ありません。従って、「CORDICアルゴリズムは収束するのか」という問いは、「xに十分近いthetaを常に合成できるか」という問いに変換できます。

表1に基づく回転だけを用いるのであれば、直感的にこのアルゴリズムは収束すると言えます。コンピュータ技術者であれば、2の巾乗だけを足し合わせて任意の整数を作ることができると知っています。二進数の基本です。

しかし、乗算を削るために導き出した表3は、θの数列が等比数列になっていません。従って、収束するかどうかは自明ではありません。

詳しい考察はめんどくさいのでしていませんが、直感的に「任意のtheta[k]を選んだときに、theta[k+1]からtheta[n]の総和がtheta[k]より大きければ、xに十分近いthetaを合成できる」と言えるようです。

なんにせよ、世間では広くCORDICアルゴリズムが使われています。精度に関して言えば、N回繰り返しを行うCORDICアルゴリズムは、Nが大きな領域ではNを1増やす毎にアルゴリズム固有の雑音が1/2になっていきます。

atan(x,y)

ここまで、与えられた角xからsin(x),cos(x)を求めるCORDICアルゴリズムについて解説しました。

では、与えられた座標(X,Y)から角xを求めるatan2(Y,X)アルゴリズムはどうでしょうか。

結論からいえば、sin,cosを求めるアルゴリズムをほんの少し修正すればそれで実現できます。ポイントとなるのは、ループの中のif文の条件式です。sin,cosを求めるときには「合成したthetaがxより大きいか否か」を判別しています。atan2(Y,2)アルゴリズムでは、ここを「qの値が正か負か」で判別するようにします。qを0に近づくように追い込んでいけば、ループから出たときには-atan2(Y,X)の値がthetaに格納されています。

終わりに

CORDICアルゴリズムについて解説しました。

このアルゴリズムは実装が容易で必要な資源も軽いことで有名です。しかし、それを導き出す過程は最適化という観点から見たときに非常に興味深いものです。信号処理の最適化では時たまアルゴリズムから数学まで引き返して最適化を行う必要がありますが、CORDICはその好例と言えるでしょう。