FPGAやマイコンで「角度 $\theta$ から $\sin\theta$ と $\cos\theta$ を計算したい」という場面はとても多くあります。たとえばソフトウェア無線（SDR）でIQ信号の周波数変換を行うとき、毎サンプルごとに $\cos$ と $\sin$ の値が必要になります。ところが、小さなFPGAや低消費電力のマイコンには「乗算器」が少ししか載っていなかったり、まったく無かったりします。テイラー展開で三角関数を計算しようとすると、何度も掛け算が必要になり、限られたハードウェア資源では実装が苦しくなります。
ここで登場するのが CORDIC（COordinate Rotation DIgital Computer） です。CORDICの驚くべき点は、掛け算をいっさい使わず、ビットシフトと加減算だけで $\sin$、$\cos$、$\arctan$、ベクトルの大きさ $\sqrt{x^2+y^2}$ などを計算できることです。ビットシフトはハードウェアでは配線を1本ずらすだけで実現でき、加算器も非常に小さい回路で済みます。そのため、乗算器を持たない安価なFPGAやASICでも高精度な三角関数計算が可能になります。
この技術は実際に、初期の関数電卓（HP-35が有名です）、デジタルシンセサイザのNCO（数値制御発振器）、レーダーの座標変換、ロボットの逆運動学、そして高速な複素数演算など、驚くほど広い分野で使われています。本記事では、CORDICが「なぜシフトと加算だけで動くのか」をベクトルの逐次回転という視点から一行ずつ導出し、回転モードとベクトルモードの2つの動作を完全に理解した上で、Pythonで実装して収束の速さと精度を実際に確かめます。
本記事の内容
- CORDICの直感 — 大きな回転を「決まった角度の足し引き」で近似する
- 回転行列を擬似回転の積に分解する導出（省略なし）
- スケール係数 $K = \prod \cos(\arctan 2^{-i})$ の導出
- 回転モード（角度入力 → $\sin/\cos$）の更新式
- ベクトルモード（座標入力 → 振幅・位相）の更新式
- Pythonでの固定反復回数CORDIC実装と収束・精度評価
- 乗算器レスでの三角関数計算という利点
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
回転行列、複素数の極形式、ビットシフトと2のべき乗の関係を知っていると、各ステップの意味がスムーズに頭に入ります。とはいえ、これらの知識も本文中で必要に応じて補いますので、はじめての方でも読み進められます。
CORDICとは — 回転を「決まった角度の積み木」で作る
まず、CORDICが何をやっているのかを直感的にイメージしてみましょう。あなたが「ちょうど30度だけ向きを変えたい」とします。しかし手元には「±45度」「±26.57度」「±14.04度」「±7.13度」……といった、あらかじめ決められた角度の回転だけが許されているとします。しかもそれぞれの角度は1回ずつしか使えません。このとき、これらの決まった角度を上手に足したり引いたりして、目標の30度にできるだけ近づけることを考えます。
たとえば「+45度して、次に−14.04度して、さらに−0.9度して……」というように、許された角度を順番に足し引きしていくと、合計が少しずつ30度に近づいていきます。CORDICはまさにこれをやっています。回転の角度を、あらかじめ決められた特別な角度の列で「積み木のように」組み立てていくのです。
ではなぜ、わざわざこんな回りくどいことをするのでしょうか。鍵は「どんな特別な角度を選ぶか」にあります。CORDICでは、回転角を
$$ \alpha_i = \arctan(2^{-i}) \quad (i = 0, 1, 2, \dots) $$
という角度に制限します。すなわち $\tan\alpha_i = 2^{-i}$ となる角度です。$i=0$ なら $\tan\alpha_0 = 1$ で $\alpha_0 = 45^\circ$、$i=1$ なら $\tan\alpha_1 = 0.5$ で $\alpha_1 \approx 26.57^\circ$、$i=2$ なら $\tan\alpha_2 = 0.25$ で $\alpha_2 \approx 14.04^\circ$、というように、だんだん小さくなる角度の列になります。
なぜ $\tan\alpha_i = 2^{-i}$ という一見奇妙な選び方をするのか。それは、後で見るように回転の計算式の中に $\tan\alpha_i$ という掛け算が現れるからです。$\tan\alpha_i = 2^{-i}$ にしておけば、その掛け算が「2の $-i$ 乗を掛ける」、つまりビットを $i$ 桁右にシフトするだけになります。掛け算がシフトに化けるこの一点が、CORDICがハードウェアで愛される最大の理由です。
ここで自然な疑問が生まれます。決められた角度を足し引きするだけで、本当に任意の角度に到達できるのでしょうか。そして、回転を擬似的に行うとベクトルの長さが変わってしまうのではないでしょうか。次のセクションから、回転行列の式を出発点にして、これらの疑問に数式で答えていきます。
回転行列から擬似回転への分解
出発点 — 普通の回転行列
平面上の点 $(x, y)$ を原点まわりに角度 $\theta$ だけ回転させると、新しい座標 $(x’, y’)$ は回転行列を使って
$$ \begin{pmatrix} x’ \\ y’ \end{pmatrix} = \begin{pmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix} $$
と書けます。成分で書くと
$$ \begin{align} x’ &= x\cos\theta – y\sin\theta \\ y’ &= x\sin\theta + y\cos\theta \end{align} $$
です。この式には $\sin\theta$ と $\cos\theta$ が含まれているので、そのままでは「三角関数を計算するために三角関数が必要」という堂々めぐりになってしまいます。CORDICの目標は、この式をうまく変形して、三角関数の値を知らなくても回転を実行できるようにすることです。
$\cos\theta$ をくくり出す
ここで一つ目の工夫をします。$\cos\theta$ を両式から強引にくくり出してみましょう。$\sin\theta = \cos\theta \cdot \tan\theta$ という関係を使うと、
$$ \begin{align} x’ &= \cos\theta\,(x – y\tan\theta) \\ y’ &= \cos\theta\,(y + x\tan\theta) \end{align} $$
と書き直せます。1行目は $x\cos\theta – y\sin\theta = x\cos\theta – y\cos\theta\tan\theta = \cos\theta(x – y\tan\theta)$ という変形、2行目も同様に $\sin\theta = \cos\theta\tan\theta$ を代入して $\cos\theta$ をくくり出したものです。
この形を見ると、回転は「$\tan\theta$ を使った座標の足し引き」と「$\cos\theta$ を掛けるスケーリング」の2段階に分かれていることがわかります。$\cos\theta$ のスケーリングをいったん脇に置いて、$\tan\theta$ による足し引きの部分だけに注目します。これを 擬似回転（pseudo-rotation） と呼びます。
$$ \begin{align} \tilde{x}’ &= x – y\tan\theta \\ \tilde{y}’ &= y + x\tan\theta \end{align} $$
擬似回転は、本来の回転からスケーリング $\cos\theta$ を取り除いたものです。後でまとめて補正するので、いまはこの形で進めます。
大きな角度を小さな角度の和に分解する
次に、回転したい角度 $\theta$ を、前のセクションで導入した特別な角度 $\alpha_i = \arctan(2^{-i})$ の符号付きの和で表します。
$$ \theta = \sum_{i=0}^{N-1} \sigma_i \, \alpha_i, \quad \sigma_i \in \{+1, -1\} $$
ここで $\sigma_i$ は各ステップで「足すか引くか」を決める符号です。各ステップで一つの $\alpha_i$ ずつ回転させていけば、$N$ 回の擬似回転の積み重ねで、全体として $\theta$ だけ回転したことになります。1ステップ分の擬似回転は、$\tan(\sigma_i\alpha_i) = \sigma_i\tan\alpha_i = \sigma_i 2^{-i}$ なので
$$ \begin{align} x_{i+1} &= x_i – y_i \cdot \sigma_i 2^{-i} \\ y_{i+1} &= y_i + x_i \cdot \sigma_i 2^{-i} \end{align} $$
となります。$\tan(\sigma_i\alpha_i) = \sigma_i 2^{-i}$ を擬似回転の式に代入しただけです。ここで注目してほしいのは、$2^{-i}$ を掛ける操作は $i$ ビットの右シフト にほかならないということです。$y_i \cdot 2^{-i}$ は $y_i$ を $i$ 桁右にシフトした値、$x_i \cdot 2^{-i}$ も同様です。つまり1ステップの更新は、シフトと加減算だけで完結します。掛け算は一つも要りません。
これで「掛け算をシフトに置き換える」というCORDICの心臓部が見えました。しかし、まだ擬似回転のせいでベクトルの長さが伸びてしまう問題が残っています。次のセクションで、この伸びをスケール係数 $K$ として正確に計算します。
スケール係数 $K$ の導出
擬似回転がベクトルを何倍に伸ばすか
本来の回転は長さを変えませんが、私たちが行っている擬似回転は $\cos\theta$ を掛け忘れている分だけベクトルを伸ばしてしまいます。1ステップ分の擬似回転で、ベクトルの長さがどれだけ変わるかを計算してみましょう。擬似回転後の長さの2乗は
$$ \tilde{x}’^2 + \tilde{y}’^2 = (x – y\tan\alpha_i)^2 + (y + x\tan\alpha_i)^2 $$
です。右辺を展開します。
$$ \begin{align} &= x^2 – 2xy\tan\alpha_i + y^2\tan^2\alpha_i + y^2 + 2xy\tan\alpha_i + x^2\tan^2\alpha_i \\ &= x^2 + y^2 + (x^2 + y^2)\tan^2\alpha_i \\ &= (x^2 + y^2)(1 + \tan^2\alpha_i) \end{align} $$
1行目から2行目では、$-2xy\tan\alpha_i$ と $+2xy\tan\alpha_i$ が打ち消し合い、残った項を $(x^2+y^2)$ でくくっています。2行目から3行目では $(x^2+y^2)$ を再びくくり出しました。三角関数の恒等式 $1 + \tan^2\alpha = 1/\cos^2\alpha = \sec^2\alpha$ を使うと
$$ \tilde{x}’^2 + \tilde{y}’^2 = \frac{x^2 + y^2}{\cos^2\alpha_i} $$
となります。両辺の平方根をとると、擬似回転後の長さは元の長さの $1/\cos\alpha_i$ 倍、すなわち擬似回転1回ごとにベクトルは $1/\cos\alpha_i$ 倍に伸びることがわかりました。
全ステップの伸びを掛け合わせる
擬似回転を $N$ 回繰り返すと、伸びは毎回掛け算されていくので、全体の伸び率は
$$ \prod_{i=0}^{N-1} \frac{1}{\cos\alpha_i} $$
になります。この伸びを打ち消すために、最後にまとめて掛けるべき補正係数を スケール係数（gain） $K$ と呼びます。$K$ は伸び率の逆数なので
$$ K = \prod_{i=0}^{N-1} \cos\alpha_i = \prod_{i=0}^{N-1} \cos\left(\arctan 2^{-i}\right) $$
です。ここで重要なのは、$\sigma_i$（足すか引くか）が $\cos$ の中に入っても $\cos(-\alpha_i) = \cos(\alpha_i)$ なので、スケール係数 $K$ は回転の向きによらず一定 だということです。つまり、どんな角度を回転させようとも、ステップ数 $N$ さえ決めれば $K$ はあらかじめ計算できる定数になります。これがCORDICがうまく動く理由のもう一つの柱です。
$K$ の具体的な値
$\cos(\arctan 2^{-i})$ をもう少し扱いやすい形にしましょう。$\tan\alpha_i = 2^{-i}$ のとき、直角三角形を考えると、対辺が $2^{-i}$、隣辺が $1$ なので斜辺は $\sqrt{1 + 2^{-2i}} = \sqrt{1 + 4^{-i}}$ です。したがって
$$ \cos\alpha_i = \frac{1}{\sqrt{1 + 2^{-2i}}} $$
となり、スケール係数は
$$ K = \prod_{i=0}^{N-1} \frac{1}{\sqrt{1 + 2^{-2i}}} $$
と書けます。$N$ を十分大きくすると、この無限積は
$$ K_\infty = \prod_{i=0}^{\infty} \frac{1}{\sqrt{1 + 2^{-2i}}} \approx 0.6072529350\dots $$
という定数に収束します。逆数 $1/K_\infty \approx 1.6467602581\dots$ はCORDICのゲインと呼ばれ、文献によく登場する有名な値です。実装では $K$ を定数として持っておき、最後に1回だけ掛けるか、あるいは初期値に $K$ を仕込んでおくことで、計算の途中に掛け算を一切入れずに済ませます。
これでCORDICの数学的な骨組みがすべて揃いました。あとは「各ステップで $\sigma_i$ をどう決めるか」によって、CORDICは2つの異なる使い方に分かれます。次のセクションでその2つのモードを具体的に見ていきましょう。
回転モード — 角度から $\sin$ と $\cos$ を得る
角度を消し込む戦略
CORDICには大きく2つの動作モードがあります。一つ目が 回転モード（rotation mode） です。これは「回転させたい角度 $\theta$ が与えられたとき、ベクトルを実際に $\theta$ だけ回す」モードで、$\sin\theta$ と $\cos\theta$ を計算する目的で使われます。
回転モードの戦略はこうです。「あといくつ回せばよいか」を覚えておく角度レジスタ $z_i$ を用意し、初期値を $z_0 = \theta$ とします。各ステップで、$z_i$ がまだ正なら「もっと回す必要がある」ので $+\alpha_i$ 回転（$\sigma_i = +1$）し、$z_i$ が負なら「回しすぎた」ので $-\alpha_i$ 回転（$\sigma_i = -1$）します。つまり
$$ \sigma_i = \operatorname{sign}(z_i) = \begin{cases} +1 & z_i \geq 0 \\ -1 & z_i < 0 \end{cases} $$
と決めます。回転するたびに「使った角度」を $z_i$ から差し引いていくので、$z_i$ はだんだん $0$ に近づきます。最終的に $z_N \approx 0$ になれば、ちょうど $\theta$ だけ回転し終えたことになります。
回転モードの更新式
回転モードの3本の更新式を一行ずつ書き下します。
$$ \begin{align} x_{i+1} &= x_i – \sigma_i \, 2^{-i} \, y_i \\ y_{i+1} &= y_i + \sigma_i \, 2^{-i} \, x_i \\ z_{i+1} &= z_i – \sigma_i \, \alpha_i \end{align} $$
1本目と2本目は前に導いた擬似回転そのものです。3本目が新しく加わった角度レジスタの更新で、「使った角度 $\sigma_i\alpha_i$ を残り角度 $z_i$ から引く」という意味です。ここで $\alpha_i = \arctan(2^{-i})$ の値はあらかじめテーブル（ROM）に格納しておきます。微小な角度の列を $N$ 個記憶しておくだけなので、メモリもごくわずかで済みます。
$\sin$ と $\cos$ の取り出し方
回転モードで $\sin\theta$ と $\cos\theta$ を得るには、初期ベクトルを工夫します。初期値を
$$ x_0 = K, \quad y_0 = 0, \quad z_0 = \theta $$
とおきます。$N$ ステップ後、ベクトル $(x_0, y_0) = (K, 0)$ は角度 $\theta$ だけ回転し、さらに擬似回転による伸び $1/K$ が掛かります。元の長さが $K$、伸びが $1/K$ なので、最終的なベクトルの長さはちょうど $1$ になります。回転後の長さ $1$ のベクトルが角度 $\theta$ を向くので
$$ x_N \approx \cos\theta, \quad y_N \approx \sin\theta $$
が同時に得られます。一度の計算で $\sin$ と $\cos$ の両方が手に入るのは、IQ変調のように両方を同時に必要とする用途では大きな利点です。初期値に $K$ を仕込んでおくことで、最後にスケール補正の掛け算をする必要すらありません。
ここで一点、注意があります。CORDICがそのまま収束するのは $\theta$ が約 $\pm 99.7^\circ$（すべての $\alpha_i$ の総和）の範囲に限られます。これより外側の角度は、あらかじめ象限を判定して $\pm 90^\circ$ 回転を済ませてから（これは座標の入れ替えと符号反転だけでできます）、残りをCORDICで処理します。
これで「角度から三角関数を得る」回転モードが完成しました。では逆に、座標 $(x, y)$ が与えられたときに、その大きさと角度を求めたい場合はどうすればよいでしょうか。次のセクションで、符号の決め方を変えるだけで実現できるベクトルモードを見ていきます。
ベクトルモード — 座標から振幅と位相を得る
$y$ を消し込む戦略
二つ目の動作が ベクトルモード（vectoring mode） です。これは回転モードと逆の問題を解きます。すなわち、座標 $(x, y)$ が与えられたとき、その大きさ（振幅）$\sqrt{x^2 + y^2}$ と角度（位相）$\arctan(y/x)$ を求めます。複素数 $x + jy$ を極形式に変換する操作だと思えばわかりやすいでしょう。
ベクトルモードの戦略は、「ベクトルを回転させて $x$ 軸の上にぴったり乗せる」ことです。ベクトルが $x$ 軸上に来れば $y$ 成分は $0$ になり、そのとき $x$ 成分がベクトルの長さ（に伸びを掛けたもの）、回転に使った角度の合計が元のベクトルの偏角になります。そこで、$y_i$ を $0$ に追い込むように符号を決めます。$y_i$ が正なら時計回り（$\sigma_i = -1$）に回して $y$ を減らし、$y_i$ が負なら反時計回り（$\sigma_i = +1$）に回して $y$ を増やします。式で書くと
$$ \sigma_i = -\operatorname{sign}(y_i) = \begin{cases} -1 & y_i \geq 0 \\ +1 & y_i < 0 \end{cases} $$
です。回転モードでは $z_i$ の符号を見ていましたが、ベクトルモードでは $y_i$ の符号を見る、という違いがすべてです。
ベクトルモードの更新式
更新式そのものは回転モードと同じ形ですが、$\sigma_i$ の決め方が変わります。
$$ \begin{align} x_{i+1} &= x_i – \sigma_i \, 2^{-i} \, y_i \\ y_{i+1} &= y_i + \sigma_i \, 2^{-i} \, x_i \\ z_{i+1} &= z_i – \sigma_i \, \alpha_i, \quad \sigma_i = -\operatorname{sign}(y_i) \end{align} $$
ここで $z_i$ は「これまでに回した角度の合計」を記録するアキュムレータとして働きます。初期値 $z_0 = 0$ から始めて、毎回 $\sigma_i\alpha_i$ を引いていくので、最終的に $z_N$ には「ベクトルを $x$ 軸に乗せるのに必要だった回転角の総和の符号反転」、すなわち元のベクトルの偏角が貯まります。
振幅と位相の取り出し方
初期値を
$$ x_0 = x, \quad y_0 = y, \quad z_0 = 0 $$
とすると、$N$ ステップ後には $y_N \approx 0$ となり、
$$ x_N \approx \frac{1}{K}\sqrt{x^2 + y^2}, \quad z_N \approx \arctan\frac{y}{x} $$
が得られます。$x_N$ には擬似回転による伸び $1/K$ が掛かっているので、真の振幅を得るには $K$ を掛けます（$\sqrt{x^2+y^2} = K \cdot x_N$）。$z_N$ はそのまま位相（偏角）です。
このベクトルモードは、レーダーやSDRで複素信号 $I + jQ$ の振幅と位相を求めるのに直接使えます。FFTの出力ビンごとの振幅・位相を求める処理にも応用できます。$\sqrt{\cdot}$ も $\arctan$ も、専用の関数を呼ばずにシフトと加算だけで同時に得られるのが嬉しい点です。
ここまでで、回転モードとベクトルモードの両方の更新式と、それぞれの初期値・出力の取り出し方が出揃いました。理論が正しく動くこと、そして何回反復すればどれくらいの精度が出るのかを、次のセクションでPythonを使って実際に確かめましょう。
Pythonでの実装と検証
角度テーブルとスケール係数の準備
まず、CORDICの土台となる角度テーブル $\alpha_i = \arctan(2^{-i})$ とスケール係数 $K$ を計算します。これらは反復回数 $N$ を決めれば一度だけ計算しておける定数です。
import numpy as np
def cordic_tables(n_iter):
"""CORDICの角度テーブルとスケール係数Kを生成する"""
# alpha_i = arctan(2^-i) を i=0..n_iter-1 まで
i = np.arange(n_iter)
angles = np.arctan(2.0 ** (-i.astype(float)))
# スケール係数 K = prod cos(alpha_i) = prod 1/sqrt(1+2^-2i)
K = np.prod(1.0 / np.sqrt(1.0 + 2.0 ** (-2.0 * i)))
return angles, K
for n in [4, 8, 16, 24]:
angles, K = cordic_tables(n)
total = np.degrees(np.sum(angles)) # 収束する角度範囲の上限
print(f"N={n:2d}: K={K:.10f}, 1/K={1/K:.10f}, 収束範囲=±{total:.2f}deg")
このコードを実行すると、$N$ を増やすにつれてスケール係数 $K$ が約 $0.6072529350$ に収束していく様子と、CORDICが正しく収束できる角度範囲が約 $\pm 99.88^\circ$（$N$ が大きいとき、全 $\alpha_i$ の総和）になることが確認できます。$1/K$ が約 $1.6468$ という有名な値に近づいていくのも読み取れます。$N=4$ のように反復が少ないと角度テーブルが短く、$K$ もまだ収束しきっていないことがわかります。
回転モードの実装
次に、角度 $\theta$ を入力して $\cos\theta$ と $\sin\theta$ を返す回転モードを実装します。掛け算を使わないというCORDICの精神を尊重し、$2^{-i}$ の乗算をビットシフトに相当する操作（np.ldexp で $2^{-i}$ 倍）で表現します。
import numpy as np
def cordic_rotation(theta, n_iter=16):
"""回転モード: 角度thetaからcos, sinを計算（シフトと加減算のみ）"""
angles, K = cordic_tables(n_iter)
# 初期値: x=K を仕込むことで最後のスケール補正を不要にする
x, y, z = K, 0.0, theta
for i in range(n_iter):
# z>=0 なら +回転、z<0 なら -回転
sigma = 1.0 if z >= 0 else -1.0
# 2^-i 倍はビットシフト相当（ldexp(v, -i) = v * 2^-i）
x_shift = np.ldexp(x, -i)
y_shift = np.ldexp(y, -i)
x_new = x - sigma * y_shift
y_new = y + sigma * x_shift
z = z - sigma * angles[i]
x, y = x_new, y_new
return x, y # x≈cos(theta), y≈sin(theta)
# 検証: いくつかの角度でnp.cos/np.sinと比較
for deg in [0, 30, 45, 60, 90, -45]:
th = np.radians(deg)
c, s = cordic_rotation(th, n_iter=16)
print(f"{deg:4d}deg: cos CORDIC={c:.6f} (true {np.cos(th):.6f}), "
f"sin CORDIC={s:.6f} (true {np.sin(th):.6f})")
出力を見ると、CORDICで計算した $\cos$ と $\sin$ が、np.cos / np.sin の真値と小数点以下6桁程度まで一致していることがわかります。掛け算を一切使わず、ビットシフトと加減算だけでここまでの精度が出ているのは印象的です。$\theta=45^\circ$ がぴったり合いやすいのは、$\alpha_0 = 45^\circ$ という最初のステップだけで到達できる特別な角度だからです。
反復回数と精度の関係（収束曲線）
CORDICは反復1回につき、おおよそ1ビット（約 $\log_{10}2 \approx 0.3$ 桁）ずつ精度が上がる「線形収束」のアルゴリズムです。これを実際にグラフで確かめます。
import numpy as np
import matplotlib.pyplot as plt
# 多数の角度に対する平均誤差を反復回数ごとに測る
test_deg = np.linspace(-89, 89, 179)
test_rad = np.radians(test_deg)
iters = np.arange(1, 25)
mean_err = []
for n in iters:
errs = []
for th in test_rad:
c, s = cordic_rotation(th, n_iter=n)
# cos, sin 同時の誤差をユークリッド距離で評価
err = np.hypot(c - np.cos(th), s - np.sin(th))
errs.append(err)
mean_err.append(np.mean(errs))
plt.figure(figsize=(9, 5))
plt.semilogy(iters, mean_err, 'o-', color='tab:blue')
# 1反復1ビットの理論線（参考）: 誤差 ~ 2^-n
plt.semilogy(iters, 2.0 ** (-iters.astype(float)), 'r--', label='$2^{-N}$ (1 bit/iter)')
plt.xlabel('Number of iterations N')
plt.ylabel('Mean error (log scale)')
plt.title('CORDIC Rotation Mode: Convergence')
plt.grid(True, which='both', alpha=0.3)
plt.legend()
plt.tight_layout()
plt.savefig('cordic_convergence.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフから、誤差が反復回数に対して 直線的に（対数軸で）減少 していることが読み取れます。これはまさに「1反復ごとに約1ビット精度が上がる」線形収束の特徴で、赤い破線の $2^{-N}$ とほぼ平行に走っています。$N=16$ で誤差が $10^{-5}$ 程度、$N=24$ で $10^{-7}$ 程度まで下がり、単精度浮動小数点では $N=24$ 前後で頭打ちになる（丸め誤差が支配的になる）こともわかります。必要な精度に応じて反復回数を選べばよい、というハードウェア設計上の見通しが立ちます。
ベクトルモードの実装と検証
続いて、座標 $(x, y)$ から振幅 $\sqrt{x^2+y^2}$ と位相 $\arctan(y/x)$ を求めるベクトルモードを実装します。
import numpy as np
def cordic_vectoring(x_in, y_in, n_iter=16):
"""ベクトルモード: 座標(x,y)から振幅と位相を計算"""
angles, K = cordic_tables(n_iter)
x, y, z = float(x_in), float(y_in), 0.0
for i in range(n_iter):
# y を 0 に追い込む向きに回す
sigma = -1.0 if y >= 0 else 1.0
x_shift = np.ldexp(x, -i)
y_shift = np.ldexp(y, -i)
x_new = x - sigma * y_shift
y_new = y + sigma * x_shift
z = z - sigma * angles[i]
x, y = x_new, y_new
magnitude = K * x # 伸び1/Kを補正
phase = z # 偏角
return magnitude, phase
# 検証: いくつかの複素数で確認
for (xi, yi) in [(1.0, 0.0), (1.0, 1.0), (3.0, 4.0), (-2.0, 1.5)]:
mag, ph = cordic_vectoring(xi, yi, n_iter=20)
true_mag = np.hypot(xi, yi)
true_ph = np.arctan2(yi, xi)
# 注: 素のCORDICは第1・第4象限のみ対応（x>0想定）
print(f"(x,y)=({xi},{yi}): mag={mag:.6f} (true {true_mag:.6f}), "
f"phase={np.degrees(ph):.4f}deg (true {np.degrees(true_ph):.4f}deg)")
出力から、$x>0$ の場合（第1・第4象限）に振幅と位相が真値とよく一致していることが確認できます。たとえば $(3, 4)$ では振幅が $5.000000$、位相が $53.13^\circ$ と、ピタゴラスの三角形そのものの値が得られます。$(-2, 1.5)$ のように $x<0$ の象限では素のCORDICはそのままでは対応できないため、実装では事前に $180^\circ$ 回転（符号反転）して第1・第4象限に持ち込み、最後に位相へ $\pm 180^\circ$ を足し戻す前処理が必要になります。
np.sin / np.cos との誤差の可視化
最後に、回転モードのCORDIC出力を全角度域でライブラリ関数と比較し、誤差の分布を可視化します。
import numpy as np
import matplotlib.pyplot as plt
n_iter = 16
deg = np.linspace(-90, 90, 361)
rad = np.radians(deg)
cos_cordic, sin_cordic = [], []
for th in rad:
c, s = cordic_rotation(th, n_iter=n_iter)
cos_cordic.append(c)
sin_cordic.append(s)
cos_cordic = np.array(cos_cordic)
sin_cordic = np.array(sin_cordic)
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# (1) 波形の重ね描き
ax = axes[0]
ax.plot(deg, np.cos(rad), 'k-', linewidth=2, label='cos (numpy)')
ax.plot(deg, cos_cordic, 'r--', linewidth=1.2, label='cos (CORDIC)')
ax.plot(deg, np.sin(rad), 'b-', linewidth=2, label='sin (numpy)')
ax.plot(deg, sin_cordic, 'g--', linewidth=1.2, label='sin (CORDIC)')
ax.set_xlabel('Angle [deg]')
ax.set_ylabel('Value')
ax.set_title(f'CORDIC vs numpy (N={n_iter})')
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
# (2) 誤差
ax = axes[1]
ax.plot(deg, np.abs(cos_cordic - np.cos(rad)), 'r-', label='|cos error|')
ax.plot(deg, np.abs(sin_cordic - np.sin(rad)), 'b-', label='|sin error|')
ax.set_xlabel('Angle [deg]')
ax.set_ylabel('Absolute error')
ax.set_title(f'Error vs numpy (N={n_iter})')
ax.set_yscale('log')
ax.legend(fontsize=9)
ax.grid(True, which='both', alpha=0.3)
plt.tight_layout()
plt.savefig('cordic_vs_numpy.png', dpi=150, bbox_inches='tight')
plt.show()
左のグラフでは、CORDICの破線がライブラリ関数の実線にぴったり重なり、視覚的には区別がつきません。右の誤差グラフから、誤差が全角度域でおおむね $10^{-5}$ 前後に収まっており、角度によって誤差が少し波打つことがわかります。これは、各角度を表現するのに使える $\alpha_i$ の組み合わせが角度ごとに微妙に異なり、$N$ 回で打ち切ったときの「残り角度」$z_N$ の大きさが角度依存になるためです。それでも最大誤差が $N=16$ で $10^{-5}$ オーダーに抑えられており、16ビット程度の固定小数点DSPには十分な精度です。
まとめ
本記事では、乗算器を使わずシフトと加減算だけで三角関数や座標変換を計算するCORDICアルゴリズムを、ベクトルの逐次回転という視点から導出し、Pythonで検証しました。
- 基本アイデア: 大きな回転を $\alpha_i = \arctan(2^{-i})$ という決まった角度の符号付きの和に分解する。$\tan\alpha_i = 2^{-i}$ なので、回転の式に現れる掛け算がすべて ビットシフト に化ける
- 擬似回転とスケール係数: $\cos\theta$ をくくり出した擬似回転はベクトルを $1/\cos\alpha_i$ 倍に伸ばす。全ステップの伸びを補正する係数が $K = \prod \cos(\arctan 2^{-i}) \approx 0.6073$ で、回転の向きによらず一定なので定数として扱える
- 回転モード: $z_i$ の符号で $\sigma_i$ を決め、$z_i$ を $0$ に追い込む。初期値 $(K, 0, \theta)$ から $(\cos\theta, \sin\theta)$ を同時に得る
- ベクトルモード: $y_i$ の符号で $\sigma_i$ を決め、$y_i$ を $0$ に追い込む。初期値 $(x, y, 0)$ から振幅 $K x_N = \sqrt{x^2+y^2}$ と位相 $z_N = \arctan(y/x)$ を得る
- 収束と精度: 1反復につき約1ビットの線形収束。$N=16$ で $10^{-5}$、$N=24$ で $10^{-7}$ 程度の精度が出る
- 利点: 乗算器レスのFPGA/DSPで $\sin$、$\cos$、$\arctan$、$\sqrt{x^2+y^2}$ を小さな回路で高精度に計算できる。NCO、レーダーの座標変換、ロボットの逆運動学などで広く使われる
CORDICは「掛け算をシフトに置き換える」という一つのアイデアから、これだけ多彩な関数を統一的に計算できる美しいアルゴリズムです。同じ更新式の符号判定を変えるだけで回転モードとベクトルモードが切り替わる構造は、ハードウェアでの再利用性という観点でも理想的です。
次のステップとして、以下の記事も参考にしてください。