嵌入式处理器的快速斜边算法?

是否有一个聪明/有效的算法来确定角度的斜边(即sqrt(a² + b²) ),在没有硬件乘法的嵌入式处理器上使用定点数学运算?

除非你在大于1kHz的频率下进行此操作,否则即使在没有硬件的情况下,也可以在没有硬件的情况下倍增。 更糟糕的是sqrt 。 我会尝试修改我的应用程序,因此根本不需要计算它。

如果你确实需要它,标准库可能是最好的,但你可以看看使用牛顿方法作为一种可能的选择。 然而,它需要几个乘法/除法循环才能执行。

AVR资源

  • Atmel应用笔记AVR200:乘法和除法例程(pdf)
  • 这个sqrt函数在AVR Freaks论坛上
  • 另一个 AVR Freakspost

如果结果不一定非常准确,您可以非常简单地得到粗略的近似值:

ab绝对值,并在必要时进行交换,以便a <= b 。 然后:

 h = ((sqrt(2) - 1) * a) + b 

为了直观地了解其工作原理,请考虑在像素显示器上绘制浅角度线的方式(例如,使用Bresenham算法)。 它看起来像这样:

 +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ | | | | | | | | | | | | | | | | |*|*|*| ^ +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ | | | | | | | | | | | | | |*|*|*|*| | | | | +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ | | | | | | | | | |*|*|*|*| | | | | | | | a pixels +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ | | | | | |*|*|*|*| | | | | | | | | | | | | +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ | |*|*|*|*| | | | | | | | | | | | | | | | v +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ <-------------- b pixels -----------> 

对于b方向上的每个步骤,要绘制的下一个像素要么直接向右,要么向上和向右一个像素。

从一端到另一端的理想线可以通过将每个像素的中心连接到相邻像素的中心的路径来近似。 这是一系列长度为sqrt(2)段,以及长度为1的ba段(以像素为单位)。 因此上面的公式。

这清楚地给出a == 0a == 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≈bmax(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  /* 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, so * P==0 only if Q==0; OTOH, if Q==0, then distance == P... */ if ( Q == 0.0 ) return P; for (i=0;i 

这仍然会在每次迭代时进行几次除法和四次乘法运算,但每次输入很少需要三次以上的迭代(两次通常就足够了)。 至少对于我见过的大多数处理器来说,这通常会比sqrt本身更快。

目前它是为double s编写的,但假设你已经实现了基本操作,将其转换为使用固定点工作应该不是非常困难。

如果您需要sqrt ,可以从重新评估开始。 很多时候你只是计算斜边以将它与另一个值进行比较 – 如果你将你所比较的值平方,你可以完全消除平方根。

也许您可以使用一些Elm Chans 汇编程序库并将ihypotfunction调整到您的ATtiny。 您需要更换MUL,并且可能(我还没有检查过)其他一些说明。