sqrt(a² + b²)ハードウェア乗算を行わずに組み込みプロセッサで固定小数点演算を使用して、角度の斜辺を決定するための巧妙で効率的なアルゴリズムはありますか?
7 に答える
結果が特に正確である必要がない場合は、非常に簡単に大まかな近似を得ることができます。
との絶対値を取りa、b必要に応じて交換して、 を取得しますa <= b。それで:
h = ((sqrt(2) - 1) * a) + b
これがどのように機能するかを直感的に理解するために、浅い角度の線がピクセル ディスプレイにプロットされる方法を考えてみましょう (例: Bresenham のアルゴリズムを使用)。次のようになります。
+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
| | | | | | | | | | | | | | | | |*|*|*| ^
+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ |
| | | | | | | | | | | | |*|*|*|*| | | | |
+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ |
| | | | | | | | |*|*|*|*| | | | | | | | a pixels
+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ |
| | | | |*|*|*|*| | | | | | | | | | | | |
+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ |
|*|*|*|*| | | | | | | | | | | | | | | | v
+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
<-------------- b pixels ----------->
方向の各ステップで、bプロットされる次のピクセルは、すぐ右または 1 ピクセル上と右のいずれかになります。
一方の端から他方の端までの理想的な線は、各ピクセルの中心を隣接するピクセルの中心に接続するパスによって近似できます。これは、長さ の一連のaセグメントとsqrt(2)、b-a長さ 1 のセグメントです (ピクセルを測定単位とします)。したがって、上記の式。
a == 0これにより、との正確な答えが明確に得られa == bます。しかし、その間の値を過大評価します。
エラーは比率に依存しb/aます。最大誤差は、 が である場合に発生b = (1 + sqrt(2)) * aし2/sqrt(2+sqrt(2))、真の値を約 8.24% 上回っています。それは素晴らしいことではありませんが、アプリケーションにとって十分である場合、この方法には単純で高速であるという利点があります。(定数による乗算は、一連のシフトと加算として記述できます。)
CORDIC 法の使用を検討してください。Dr. Dobb's には記事と関連するライブラリ ソースがあります。平方根、乗算、除算については、記事の最後で説明します。
1 つの可能性は次のようになります。
#include <math.h>
/* Iterations Accuracy
* 2 6.5 digits
* 3 20 digits
* 4 62 digits
* assuming a numeric type able to maintain that degree of accuracy in
* the individual operations.
*/
#define ITER 3
double dist(double P, double Q) {
/* A reasonably robust method of calculating `sqrt(P*P + Q*Q)'
*
* Transliterated from _More Programming Pearls, Confessions of a Coder_
* by Jon Bentley, pg. 156.
*/
double R;
int i;
P = fabs(P);
Q = fabs(Q);
if (P<Q) {
R = P;
P = Q;
Q = R;
}
/* The book has this as:
* if P = 0.0 return Q; # in AWK
* However, this makes no sense to me - we've just insured that P>=Q, so
* P==0 only if Q==0; OTOH, if Q==0, then distance == P...
*/
if ( Q == 0.0 )
return P;
for (i=0;i<ITER;i++) {
R = Q / P;
R = R * R;
R = R / (4.0 + R);
P = P + 2.0 * R * P;
Q = Q * R;
}
return P;
}
これでも反復ごとに 2 回の除算と 4 回の倍数が行われますが、入力ごとに 3 回以上の反復が必要になることはほとんどありません (通常は 2 回で十分です)。少なくとも私が見たほとんどのプロセッサでは、それは通常、sqrtそれ自体よりも高速です。
今のところdoubles 用に書かれていますが、基本的な操作を実装していると仮定すると、固定小数点で動作するように変換することはそれほど難しくありません。
「適度に堅牢」についてのコメントによって、いくつかの疑問が提起されました。少なくとも最初に書かれたように、これは基本的に、「完璧ではないかもしれませんが、ピタゴラスの定理を直接実装するよりも少なくともかなり優れている」というかなり裏返しの言い方でした。
特に、各入力を 2 乗する場合、2 乗した結果を表すには、入力値を表すのに必要なビット数の約 2 倍が必要です。追加した後 (余分なビットが 1 つだけ必要です)、平方根を取得します。これにより、入力とほぼ同じ数のビットが必要になります。入力よりも大幅に精度の高い型を使用しない限り、非常に悪い結果が簡単に得られます。
このアルゴリズムは、どちらの入力も直接二乗しません。中間結果がアンダーフローする可能性は依然としてありますが、アンダーフローした場合でも、使用中のフォーマットがサポートするのと同じように結果が出力されるように設計されています。基本的に、それが発生する状況は、非常に鋭い三角形 (たとえば、90 度、0.000001 度、および 89.99999 度など) がある場合です。90、0、90 に十分近い場合、2 つの長い辺の差を表すことができない可能性があるため、斜辺はもう一方の長辺と同じ長さとして計算されます。
対照的に、ピタゴラスの定理が失敗した場合、結果は多くの場合 NaN (つまり、何も教えてくれません) になるか、使用する浮動小数点形式によっては、妥当な答えのように見えても、実際には非常に間違っている可能性が非常に高くなります。
必要な場合は、再評価することから始めることができますsqrt。別の値と比較するためだけに斜辺を計算することがよくあります。比較する値を 2 乗すると、平方根を完全に削除できます。
>1kHz でこれを行っていない限り、ハードウェアのMULない MCU でも乗算はひどいものではありません。さらに悪いのは、sqrt. まったく計算する必要がないように、アプリケーションを変更しようとします。
実際に必要な場合は、おそらく標準ライブラリが最適ですが、可能な代替手段としてニュートンの方法を使用することを検討できます。ただし、実行するには数回の乗算/除算サイクルが必要です。
AVR リソース
- Atmel アプリケーション ノート AVR200: 乗算および除算ルーチン (pdf)
sqrtAVR Freaks フォーラムのこの機能- 別のAVR フリークスの投稿