この記事は Gemini に記事として作らせ、わずかに修正したものです。(初出:2017)
はじめに
組み込み開発で角度計算が必要になった際、math.h の atan() が重すぎたり、浮動小数点演算が使えない環境(テーブル参照だとメモリを食う)で困ったことはありませんか?
一般に近似といえばマクローリン展開が有名ですが、$arctan$ の級数展開は収束が絶望的に遅く、実用的ではありません。
本記事では、私が過去にエクセルのゴールシークを駆使して導き出した、計算負荷が極めて低く、かつ誤差0.03度以内を叩き出す魔法の近似式を紹介します。
導き出された近似式
$$atan(x) \approx \frac{69x^2 + 5x + 286}{3x^2 + 5} \cdot x \quad [deg]$$
※対象範囲:$0 \le x \le 1$ ($0^\circ \le \theta \le 45^\circ$)
この式のすごいところ
- 整数係数のみ: $69, 5, 286, 3, 5$ というシンプルな整数のみで構成。
- 端点の一致: $x=0$ で $0^\circ$、$x=1$ で正確に $45^\circ$ を返す設計。
- 高精度: 全域での最大誤差が約 0.022度。
- 爆速: $x^2$ を共通化すれば、乗算3回・除算1回・加算3回で完了。
開発の背景(マクローリン展開との比較)
当初、標準的な級数展開を試みましたが、納得のいく精度を得るには項数を増やす必要があり、計算負荷が増すばかりでした。
そこで、有理関数形式の近似式をベースに、エクセルのゴールシークを用いて係数を微調整しました。最初は1次の項がない形からスタートしましたが、試行錯誤の末に $5x$ という項を加えたところ、誤差が劇的に改善。この「黄金の係数」に辿り着きました。
精度検証
実際に真値と比較した結果がこちらです。
| $x$ | 近似値 [$deg$] | 真値 [$deg$] | 誤差 [$deg$] |
|---|---|---|---|
| 0.0 | 0.0000 | 0.0000 | +0.0000 |
| 0.1 | 5.7095 | 5.7106 | -0.0011 |
| 0.2 | 11.3188 | 11.3099 | +0.0089 |
| 0.3 | 16.7197 | 16.6992 | +0.0205 |
| 0.4 | 21.8277 | 21.8014 | +0.0263 |
| 0.5 | 26.5870 | 26.5651 | +0.0219 |
| 0.6 | 30.9711 | 30.9638 | +0.0073 |
| 0.7 | 34.9794 | 34.9920 | -0.0126 |
| 0.8 | 38.6312 | 38.6598 | -0.0286 |
| 0.9 | 41.9584 | 41.9872 | -0.0288 |
| 1.0 | 45.0000 | 45.0000 | +0.0000 |
実用上の懸念点と検証
実務(特に制御系)に投入する際に気になるポイントを検証しました。
- ゼロ除算のリスク: 分母 $3x^2 + 5$ の最小値は $5$ です。入力に関わらず 0 になることはありません。
- 単調増加性の保証: $0 \le x \le 1$ の全域で微係数は正。角度が逆流する心配はありません。
- 不連続性(微係数): $x=1$ 付近で $atan2$ 等へ切り替える際、値は連続(45.0)していますが傾きには微小な差があります。角速度を算出するような超高精度制御では切り替え時のスパイクに留意してください。
組み込み用実装例(C言語)
float fast_atan_deg(float x) {
// |x| <= 1.0 の範囲を想定
float x2 = x * x;
return ((69.0f * x2 + 5.0f * x + 286.0f) / (3.0f * x2 + 5.0f)) * x;
}
組み込み用実装例(C言語 / atan2形式)
全周($-180^\circ$ ~ $180^\circ$)対応版の実装例です。
float fast_atan2_deg(float y, float x) {
if (x == 0.0f && y == 0.0f) return 0.0f;
float abs_y = (y < 0) ? -y : y;
float abs_x = (x < 0) ? -x : x;
float angle;
if (abs_x >= abs_y) {
// 0 ~ 45度の範囲
float z = abs_y / abs_x;
float z2 = z * z;
angle = ((69.0f * z2 + 5.0f * z + 286.0f) / (3.0f * z2 + 5.0f)) * z;
} else {
// 45 ~ 90度の範囲 ( atan(z) = 90 - atan(1/z) )
float z = abs_x / abs_y;
float z2 = z * z;
angle = 90.0f - ((69.0f * z2 + 5.0f * z + 286.0f) / (3.0f * z2 + 5.0f)) * z;
}
// 象限に応じた符号調整
if (x < 0) angle = 180.0f - angle;
if (y < 0) angle = (y < 0) ? -angle : angle; // yが負なら角度も負に
return angle;
}
おわりに
この式は、理論的な美しさよりも「現場での速度と精度」を追求した結果生まれたものです。浮動小数点の計算コストを削りたい、あるいは固定小数点演算への移植を考えている方の助けになれば幸いです。