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

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

取绝对值a和b,并在必要时进行交换,以便拥有a <= b.然后:

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

Run Code Online (Sandbox Code Playgroud)

为了直观地看到这是如何工作的,考虑到浅的斜线是一个像素的显示屏上绘制(例如,使用布氏算法)的方式.它看起来像这样:

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

Run Code Online (Sandbox Code Playgroud)

对于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%.这不是很好,但如果它对你的应用程序来说足够好,这种方法具有简单快速的优点.(乘以常数可以写成一系列的移位和加法.)

Edg*_*net 9

为了记录,这里有一些更近似的,按复杂性和准确性的大致递增顺序列出.所有这些假设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 ² + a ²)
=b√(1 +(a/b)²)
= bf(a/b)其中f(x)=√(1 + x ²)

添加约束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)²)= b + 0.428 a ²/b

Cli*_*ord 7

考虑使用CORDIC方法.Dobb博士在这里有一篇文章和相关的图书馆资料.平方根,乘法和除法在本文末尾处理.

注意:我使用过这个库,在log()函数中发现了一个错误.通过在`log_two_power_n_reversed []`数组初始化器的末尾添加一个"0x0LL"来纠正这个问题.我已经与作者确认了这一更正. (4认同)

Jer*_*fin 6

一种可能性如下:

#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;
}

Run Code Online (Sandbox Code Playgroud)

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

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