- CORDICは、FPUや大きなルックアップテーブルなしで
sin、cos、tan のような三角関数を計算するために、複雑な演算を加算とビットシフト中心に置き換えるアルゴリズム
- この方式は高性能システムよりも組み込み環境、特に性能の低いマイクロコントローラやFPGAで有用であり、速度だけで価値を判断するのは難しい
- 浮動小数点の代わりに固定小数点を使えば、
int32_t の上位16ビットを整数部、下位16ビットを小数部として分け、約 -32768.99997 から 32767.99997 まで表現できる
- ベクトルを目標角度に合わせて徐々に小さい角度で回転させ、
atan(2**-i) テーブル16個と x=39796 の初期値を使えば、反復ごとに乗算をビットシフトで置き換えられる
- 例の角度
0.9152 を16回反復すると、sin(0.9152) の絶対誤差は 0.00000956、cos(0.9152) の絶対誤差は 0.0000434 の水準まで小さくなる
CORDICが適している計算環境
- CORDICは、
sin、cos、tan のような三角関数を低消費電力ハードウェアで計算するためのアルゴリズム
- FPU、つまり浮動小数点演算装置がない、あるいは大きなルックアップテーブルを使いにくい環境でも動作する
- 実際の演算は単純な加算とビットシフトを中心に構成される
- ベクトル数学、三角法、収束、計算機科学的なアイデアを組み合わせて、複雑な関数を単純な演算で近似する
- 高性能ハードウェアでは、この手法が必須でない場合もある
- 主な適用対象は組み込み環境
- 特に性能の低いマイクロコントローラやFPGAに適している
- より高速なハードウェアや周辺装置があるかもしれないが、速度だけが有用性の唯一の尺度ではない
浮動小数点を避ける固定小数点表現
sin(x) のように -1.0 から 1.0 の値を返す関数も、必ずしも浮動小数点で表現する必要はない- 固定小数点は、整数型の中で小数点位置を固定して有理数を表す
- 例では
int32_t を上位16ビットの整数部と下位16ビットの小数部に分ける
- この場合の範囲はおよそ
-32768.99997 から 32767.99997 まで
- 小数点位置をどこに置くかによって、整数部の範囲と小数部の精度をトレードオフできる
- 値そのものは依然として
int32_t であり、プログラマがビット列に追加の意味を与えている
固定小数点変換と基本演算
- 小数精度が16ビットなら、
42.01 のような float 値に (1 << 16) を掛けて固定小数点値にできる
42.01 * (1 << 16) は int32_t にキャストすると 2753167 になる- 再び float に戻すには
2753167 / (1 << 16) を計算し、約 42.0099945 を得る
- 浮動小数点をまったく使わずに、
1.5 のような値を直接エンコードすることもできる
- 整数部
1 は (1 << 16) で持ち上げる
- 小数部の半分は
0x0000 と 0xffff の中間値である 0x7fff にできる
- この方式の結果は10進数で
98303 になる
- 同じスケーリング係数を使う値同士なら、加算と減算はそのまま機能する
- 乗算は2つの固定小数点値を掛けた後、結果をスケーリング係数ぶんだけ再び右シフトする
- 除算は被除数をあらかじめスケーリング係数ぶんだけ左シフトしてから除数で割ると、追加の精度を得られる
ベクトル回転で三角関数を近似する
- CORDICは “co-ordinate rotation digital computer” の略で、1950年代半ばに作られた
- 核となるアイデアは、単位円上でベクトルを徐々に小さい角度で回転させ、目標角度に到達したときにベクトル成分がサインとコサインの値になるようにすること
- この過程は二分探索に似た形で進む
- 目標角度に向かって大きな角度で移動する
- 目標を行き過ぎたかを確認する
- その後はより小さな角度で時計回りまたは反時計回りの回転を繰り返す
- 例として
sin(0.7) を求めるときは、初期ベクトル (1, 0) と目標 0.7 ラジアンから始める
- まず
0.7853 ラジアン、すなわち 45˚ を反時計回りに回転する
- 残りの目標は
0.7 - 0.7853 = -0.0853 になる
- 値が負なので、次は
0.3926 ラジアン、すなわち 22.5˚ を時計回りに回転する
- 以後は残り目標の符号に応じて
0.1963 ラジアンなど、さらに小さい角度へと方向を切り替えながら回転する
- 16回反復すると、ベクトルは元の目標角度にほぼ一致し、
y は sin(a)、x は cos(a) の近似になる
回転行列で高価な演算を減らす
- 一般的なベクトル回転は、サインとコサインを含む行列積を使う
- CORDICは三角恒等式を使って、回転行列を
tan(a) 中心の形に変換する
- 最初は
45˚、22.5˚、11.25˚ のような決まった回転角を使うため、tan(a) の値を事前計算したテーブルとして持てる
- このテーブルは16個の
uint32_t だけで済み、64バイトになる
- 比較対象として、
-1 から 1 までの4096個の値を持つ最適化されていない sin(x) テーブルは16KiB必要で、精度も低いとみなせる
- 回転ごとに先頭に付く
cos(a) 項は毎回現れるが、それらをすべて掛け合わせた値は定数に収束する
45˚、22.5˚、11.25˚ といった角度を使う場合、この積は約 0.6366 になる- この定数は、すべての反復の後に1回だけ掛ければよい
シフトと加算だけを残す角度選択
- 乗算をなくすには、
tan(a) の結果が常に2の負の冪になるように角度を選ぶ
- そのため、各反復
i=0 から 15 までについて atan(2**-i) の値を持つ16エントリのテーブルを作る
- 実際の回転角は
45˚、26.565˚、14.036˚、7.125˚ などになる
- 角度は正確に半分ずつ減っていくわけではないが、これらの角度を使っても過程は正しい結果へ収束する
tan(a) の乗算は、反復番号 i に応じたビットシフトに置き換わるcos(a) 項の積も、新しい角度選択に合わせて再計算される
- 値は約
0.60725
- 16ビット固定小数点では
39796 になる
- 最後に掛ける代わりに、初期ベクトルの
x を 1 ではなく 39796 に設定すればよい
アルゴリズム手順
- 事前計算段階では、各エントリが
atan(2**-i) のテーブルを作り、各値を固定小数点に変換する
- 変換式は
atan(2**-i) * (1 << 16)
sin や cos を求めるときは、入力角度も固定小数点に変換する
- 例の
0.9152 は 0.9152 * (1 << 16) = 59978 になる
- 初期状態は次のとおり
x = 39796
y = 0
z = 59978
z はベクトルの一部ではなく、残っている目標角度を追跡する値z の符号が回転方向を決める
z >= 0 なら反時計回りに回転し、z -= table[i] を行うz < 0 なら時計回りに回転し、z += table[i] を行う
- 各反復では
x と y に対して、加算、減算、>> i シフトだけを使う
if z >= 0:
x_next = x - (y >> i)
y_next = y + (x >> i)
z -= table[i]
else:
x_next = x + (y >> i)
y_next = y - (x >> i)
z += table[i]
x = x_next
y = y_next
収束結果の例と残るテーマ
0.9152 ラジアンの例では、最初の反復で z が正なので約 0.785 ラジアンぶん反時計回りに回転する- 2回目の反復でも
z が正なので約 0.436 ラジアン反時計回りに回転するが、目標を行き過ぎる
- 3回目の反復では
z が負になり、約 0.244 ラジアンぶん時計回りに回転する
- 4回目の反復でも
z が負なので約 0.124 ラジアン時計回りに回転する
- 角度変化が小さくなるほど、ベクトルは実際の結果の近くで前後しながら収束していく
- 16回の反復後、
y は sin(0.9152) の非常に近い近似になる
- サインの絶対誤差は
0.00000956
x のコサイン絶対誤差は 0.0000434
- まだ扱っていないテーマもある
- 関心のある角度が単位円の第1象限や第4象限の外にある場合に必要な特別処理
- CORDICの派生形で計算できる
tan、atan、asin、acos、sinh、cosh、tanh、sqrt、ln、e^x
- 対数や指数の計算向けに設計された関連アルゴリズム BKM
- 関連内容は Low Byte Productions YouTube channel でさらに詳しく扱う予定
1件のコメント
Hacker News のコメント
著者は主に FPGA のような場所に適用されると言っていたが、ゲーム開発や分散物理シミュレーションにも使える
浮動小数点計算はプラットフォーム間の決定性をそろえるのが難しく、1つの解決策は浮動小数点をそもそも避けて固定小数点の物理エンジンを実装すること
三角関数機能を実装するには CORDIC のようなものが必要になる
数年前、面白半分でこういうものを作り始めたが完成はできず、いつかまたやってみたい
https://randomascii.wordpress.com/2013/07/16/floating-point-...
要するに x87 には奇妙な点があり、丸めモードやゼロへのフラッシュといった設定を一貫して合わせる必要があり、古いプロセッサには FMA がなく、
mmsqrtpsのような近似命令には一貫した仕様がなく、コンパイラが式を再結合することもある小さなルーチンや自作ライブラリなら、苦痛ではあっても、こうした問題を避けることを保証するのは可能
IEEE-754 2008 は仕様をより明確にし、事実上 x87 の死を前提にしており、2024年には x87 を確実に避けられる
FMA も IEEE-754 2008 仕様の一部で、Intel Haswell 以降を含む現代のプロセッサには搭載されている
それでも 8-wide AVX2 と 4-wide NEON のようなアーキテクチャ差が足を引っ張ることはあるが、アセンブリやイントリンシック、あるいは Compiler Explorer や objdump で確認する C を使えば、出力を見て「これは一貫しているはずだ」と判断できる
「実際、IEEE 754 が今日のように普及した標準になる前は、固定小数点が常に使われていた。1980年から2000年ごろにかけて仕事をしていたゲーム開発者に聞けば、詳しく話してくれるだろう」
nphysics を書き直した新しいライブラリ Rapier は、代わりに IEEE-754 2008 の保証に依拠してプラットフォーム間の決定性を提供している
そのため古いプラットフォームでは動作しないが、wasm を含む現代的なプラットフォームでは決定的である
もちろん、各プラットフォームが提供する sin、cos のような超越関数ルーチンには依存できず、どこでも同じように動作するよう自前で実装する必要がある
しかし非準拠プラットフォームで実行しないなら、可能なアプローチだ
https://www.rustsim.org/blog/2020/06/01/this-month-in-rustsi...
https://rapier.rs/docs/user_guides/rust/determinism/
CORDIC はサインとコサインの計算・生成だけでなく、対数、指数、平方根、ベクトルの大きさ、極座標-直交座標変換、ベクトル回転のようなさまざまな演算にも使える
著者も結論でこうした可能性を示唆している
既存の正規直交行列の代わりにクォータニオンを使えば、CORDIC ベースの演算をより効率的に、つまり計算サイクルとメモリを少なくし、誤差も減らして実行できそうな気がする
https://core.ac.uk/works/8439118
高校のプレ微積分でテイラー級数を学び、先生が電卓の三角関数は実際にそう実装されていると言っていた
調べてみると実際には CORDIC で、TI Basic で実装して楽しく遊んだ
CORDIC ではなかったが、アルゴリズムには似た点がある
http://files.righto.com/calculator/sinclair_scientific_simul...
ハードウェア実装に関する記事:
https://arxiv.org/pdf/2211.04053
https://hal.science/hal-01327460/document
https://archive.ll.mit.edu/HPEC/agendas/proc05/Day_1/Abstrac...
時代ごとのさまざまなハードウェアで、一般的なソフトウェア／ハードウェアによる三角関数実装とどう比較されるのか見てみたい
IoTや機械間通信が成長しており、CORDICの実装と演算効率を考えると利用はおそらく大きく増えるはずなので、正しく最適化された実装のための良い参考文献が必要
例外的に、Omondi教授とDeschamps教授の本がある
https://www.worldscientific.com/worldscibooks/10.1142/p1054
http://www.arithmetic-circuits.org/guide2fpga/vhdl_codes.htm
sinとcosはベクトル回転によく使われる
この場合のCORDICのコツは、従来の sin/cos/乗算 の計算を避け、回転させるベクトルそのものをCORDICの入力として与えること
そうするとCORDICは、sin/cosを計算したり複素数乗算をしたりせずに、回転後のベクトルを直接生成する
CORDICはレイテンシがそれほど重要でないときに特に真価を発揮する
計算の各段をパイプライン化すれば大きなスループットが得られるため、無線システムのデジタルミキシングによく合う
2023年時点では、一部の現代的な MCU は安価でありながらFPUを搭載している
STM32G4が良い例で、M0 MCUのような場合と違い、固定小数点を使いたくなければ
f32を自由に使えるこうしたチップはMCUあたりおおよそ1〜2ドルで入手できる
ただしG4には、固定小数点用途向けにこのアルゴリズムを実装したハードウェアCORDICペリフェラルもある
これが主に浮動小数点の精度損失を避けるためのものなのか気になる
レジスタでプログラムするが、CPUで直接CORDICを実装するのではなく、IC内部の専用ハードウェアが処理する
最安のSTM32G4はSTM32G441KBT6で、丸めると4ドル https://www.digikey.com/en/products/detail/microchip-technol...
2ドル未満ではどこで入手できるのか気になる
Digi-KeyではNuvotonのチップが500個数量でようやく2ドル未満になる程度
高速で、64ビットの中間積を扱うため、除算と三角関数の精度はほとんどの用途に十分
必要ならソフトウェアでさらに精度を上げることもできる
CORDICを知ったのは遅かったが、それ以前は性能と決定性のために8ビット／16ビットアセンブリの世界で固定小数点を多用していた
知ってからは驚いた
高速で、有用に使うために必要な数学の知識も基本的なものだけだった
以前関わった、なかなか可愛らしいコード片を思い出す
単位円の弧が作る角の二等分線の座標を求める必要があり、両腕の
(x,y)座標はすでに分かっていた既存実装は、
(x,y)座標を極座標(r,θ)に変換し、計算されたθが正しい象限にあるか確認してから、θを半分に割り、再び(x,y)に変換するという三角関数だらけのものだった結果として三角関数と逆関数を大量に呼び出していた
Pythonでは複素数を第一級の値として扱えるので、
(x1,y1)からz1、(x2,y2)からz2という2つの複素数を定義し、積の幾何平均√(z1*z2)を取るだけで済んだ新しいコードには、明示的な三角関数も、明示的な変換と逆変換もなかった
https://fgiesen.wordpress.com/2010/10/21/finish-your-derivat...
「22.75˚だけ回転するのは、45˚回転してから-22.5˚回転するのと同じだというのはかなり明白だ」とあるが、それなら 22.5° 回転ではないのか？
記事の誤りなのか、自分が誤解しているのか気になる
Meagher のオクトリーシステムは、整数の乗算・除算を使わず、整数演算だけを用いることで有名
「ブール演算（和集合、共通集合、差集合）、幾何演算（移動、サイズ変更、回転）、N次元の干渉検出、空間内の任意地点での隠面消去を含む表示のための、効率的な線形時間アルゴリズムが開発された。このアルゴリズムは浮動小数点演算、整数乗算、整数除算を必要としない」
https://doi.org/10.1016/0146-664X(82)90104-6
そのおかげで、オクトリー表現向けの高速なカスタムVLSIグラフィックスアクセラレーションハードウェアを作りやすかった
CORDIC が、小さなテーブルを使う3次補間や他の多項式補間と比べて、どの程度の性能を出すのか気になる
リソースの限られたシンセサイザーが時々3次補間を使うと教わったが、おそらく CORDIC が比較的新しかった時期の話だと思う
ざっくり見ると、CORDIC は反復ごとに精度を1ビット得るので計算はより高コストだが、空間は多項式より少なくて済みそう
ただし空間の面では、記事で
sin(x)用に提示されていた4096項目のルックアップテーブルよりも安くできる点を強調すべき対称性のおかげで、円全体の1/4だけで足りる
バイトサイズの角度を使うと自動的に循環して便利で、2Dゲームの回転には
2^8でもかなり十分だったただし滑らかな動きを求めるなら、3Dではそれほど遠くまでは行けない