是否有一种聪明/有效的算法来确定角度的斜边(即sqrt(a² + b²)
),在没有硬件乘法的嵌入式处理器上使用定点数学?
7 回答
如果结果不必特别准确,您可以很简单地得到一个粗略的近似值:
取 和 的绝对值,a
并b
在必要时交换,以便获得a <= b
. 然后:
h = ((sqrt(2) - 1) * a) + b
要直观地了解其工作原理,请考虑在像素显示器上绘制浅角线的方式(例如,使用 Bresenham 算法)。它看起来像这样:
+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
| | | | | | | | | | | | | | | | |*|*|*| ^
+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ |
| | | | | | | | | | | | |*|*|*|*| | | | |
+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ |
| | | | | | | | |*|*|*|*| | | | | | | | a pixels
+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ |
| | | | |*|*|*|*| | | | | | | | | | | | |
+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ |
|*|*|*|*| | | | | | | | | | | | | | | | v
+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
<-------------- b pixels ----------->
对于该b
方向上的每一步,要绘制的下一个像素要么紧邻右侧,要么向右上方一个像素。
从一端到另一端的理想线可以通过将每个像素的中心连接到相邻像素的中心的路径来近似。这是一系列长度为 的a
片段sqrt(2)
,以及b-a
长度为 1 的片段(以像素为测量单位)。于是有了上面的公式。
这清楚地给出了a == 0
和的准确答案a == b
;但高估了两者之间的值。
误差取决于比率b/a
;最大误差发生在b = (1 + sqrt(2)) * a
并且结果为2/sqrt(2+sqrt(2))
或超过真实值的 8.24% 时。这不是很好,但如果它对您的应用程序来说足够好,那么这种方法的优点是简单快速。(乘以一个常数可以写成一系列的移位和加法。)
作为记录,这里还有一些近似值,按照复杂性和准确性的大致递增顺序列出。所有这些都假设 0 ≤ a ≤ b。
h = b + 0.337 * a // max error ≈ 5.5 %
h = max(b, 0.918 * (b + (a>>1))) // max error ≈ 2.6 %
h = b + 0.428 * a * a / b // max error ≈ 1.04 %
编辑:回答 Ecir Hana 的问题,这是我得出这些近似值的方法。
第一步。逼近两个变量的函数可能是一个复杂的问题。因此,我首先将其转化为逼近一个变量的函数的问题。这可以通过选择最长边作为“比例”因子来完成,如下所示:
h = √(b 2 + a 2 )
= b √(1 + (a/b) 2 )
= bf(a/b) 其中 f(x) = √(1+x 2 )
添加约束 0 ≤ a ≤ b 意味着我们只关心在区间 [0, 1] 中逼近 f(x)。
下面是相关区间中 f(x) 的图,以及 Matthew Slattery 给出的近似值(即 (√2−1)x + 1)。
第二步。下一步是盯着这个图,同时问自己一个问题“我怎样才能便宜地近似这个函数?”。由于曲线看起来大致为抛物线,我的第一个想法是使用二次函数(第三近似)。但由于这仍然相对昂贵,我还研究了线性和分段线性近似。这是我的三个解决方案:
数值常数(0.337、0.918 和 0.428)最初是自由参数。选择特定值是为了使近似值的最大绝对误差最小化。最小化当然可以通过某种算法完成,但我只是“手动”完成,绘制绝对误差并调整常数直到最小化。在实践中,这工作得非常快。编写代码来自动化这将花费更长的时间。
第三步是回到近似两个变量的函数的初始问题:
- h ≈ b (1 + 0.337 (a/b)) = b + 0.337 a
- h ≈ b max(1, 0.918 (1 + (a/b)/2)) = max(b, 0.918 (b + a/2))
- h ≈ b (1 + 0.428 (a/b) 2 ) = b + 0.428 a 2 /b
考虑使用 CORDIC 方法。Dobb 博士在这里有一篇文章和相关的图书馆资源。平方根,乘法和除法在文章的最后处理。
一种可能性如下所示:
#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;
}
这仍然会在每次迭代中进行几次除法和 4 次乘法运算,但每次输入很少需要超过 3 次迭代(通常两次就足够了)。至少对于我见过的大多数处理器来说,这通常会比单独的处理器要快sqrt
。
目前它是为double
s 编写的,但假设您已经实现了基本操作,将其转换为与定点一起工作应该不会很困难。
关于“相当稳健”的评论引发了一些质疑。至少正如最初所写的那样,这基本上是一种相当反手的说法,“它可能并不完美,但至少比直接实现勾股定理要好很多。”
特别是,当对每个输入求平方时,表示平方结果所需的位数大约是表示输入值所需位数的两倍。添加后(只需要一个额外的位),您取平方根,这使您回到需要与输入大致相同的位数。除非您的类型比输入的精度高得多,否则很容易产生非常糟糕的结果。
该算法不直接平方任一输入。中间结果仍然可能下溢,但它的设计是这样的,当它这样做时,结果仍然会出现以及使用的格式支持。基本上,发生这种情况的情况是您有一个非常锐角的三角形(例如,90 度、0.000001 度和 89.99999 度)。如果它足够接近 90、0、90,我们可能无法表示两个长边之间的差异,因此它将计算斜边与另一个长边的长度相同。
相比之下,当勾股定理失败时,结果通常是 NaN(即,什么都不告诉我们),或者根据使用的浮点格式,很可能看起来是一个合理的答案,但实际上是非常不正确的。
如果你需要的话,你可以从重新评估开始sqrt
。很多时候,您计算斜边只是为了将其与另一个值进行比较 - 如果您将要比较的值平方,则可以完全消除平方根。
除非您以 >1kHz 的频率执行此操作,否则即使在没有硬件的 MCU 上进行乘法也不会MUL
很糟糕。更糟糕的是sqrt
. 我会尝试修改我的应用程序,因此它根本不需要计算它。
如果您确实需要标准库,它可能是最好的,但您可以考虑使用牛顿方法作为可能的替代方案。然而,它需要几个乘法/除法周期才能执行。
AVR 资源
- Atmel 应用笔记 AVR200:乘法和除法例程 (pdf)
sqrt
AVR Freaks 论坛上的此功能- 另一个AVR Freaks 帖子
也许您可以使用一些 Elm Chans汇编程序库并将 ihypot 函数调整为您的 ATtiny。您需要更换 MUL 并且可能(我还没有检查过)一些其他说明。