1

我正在处理一系列 Y 值相同但 X 值不同的点。我通过将 X 加一来遍历这些点。例如，我可能有 Y = 50，X 是从 -30 到 30 的整数。我的算法的一部分涉及找到从每个点到原点的距离，然后进行进一步处理。

分析后，我发现距离计算中的 sqrt 调用占用了我大量的时间。有没有迭代的方法来计算距离？

换句话说：

我想有效地计算： r[n] = sqrt(x[n]*x[n] + y*y)). 我可以保存上一次迭代的信息。每次迭代都通过增加 x 来改变，所以x[n] = x[n-1] + 1. 我不能使用 sqrt 或 trig 函数，因为它们太慢了，除了在每个扫描线的开头。

我可以使用近似值，只要它们足够好（小于 0.l% 的误差）并且引入的误差是平滑的（我不能将它们放入预先计算的近似值表中）。

附加信息：x 和 y 始终是介于 -150 和 150 之间的整数

明天我将尝试几个想法，并根据哪个最快来标记最佳答案。

结果

我做了一些计时

  • 距离公式：16 ms/迭代
  • Pete 的 interperlating 解决方案：8 ms / 迭代
  • wrang-wrang 预计算解：8ms/迭代

我希望测试能在两者之间做出决定，因为我喜欢这两个答案。我将选择Pete's，因为它使用的内存更少。

4

6 回答 6

4

只是为了感受一下，对于您的范围 y = 50，x = 0 给出 r = 50 和 y = 50，x = +/- 30 给出 r ~= 58.3。您需要 +/- 0.1% 或 +/- 0.05 绝对值的近似值。这比大多数图书馆 sqrts 的准确性要低得多。

两种近似方法 - 您可以根据前一个值的插值计算 r，或者使用合适系列的一些项。

从前一个 r 插值

r = ( x 2 + y 2 ) 1/2

博士/dx = 1/2。2倍。( x 2 + y 2 ) -1/2 = x/r

    double r = 50;
    
    for ( int x = 0; x <= 30; ++x ) {
        
        double r_true = Math.sqrt ( 50*50 + x*x );
        
        System.out.printf ( "x: %d r_true: %f r_approx: %f error: %f%%\n", x, r, r_true, 100 * Math.abs ( r_true - r ) / r );
        
        r = r + ( x + 0.5 ) / r; 
    }

给出：

x: 0 r_true: 50.000000 r_approx: 50.000000 error: 0.000000%
x: 1 r_true: 50.010000 r_approx: 50.009999 error: 0.000002%
....
x: 29 r_true: 57.825065 r_approx: 57.801384 error: 0.040953%
x: 30 r_true: 58.335225 r_approx: 58.309519 error: 0.044065%

这似乎满足 0.1% 误差的要求，所以我没有费心编写下一个，因为它需要更多的计算步骤。

截断系列

对于 x 接近零的 sqrt ( 1 + x ) 的泰勒级数是

sqrt ( 1 + x ) = 1 + 1/2 x - 1/8 x 2 ... + ( - 1 / 2 ) n+1 x n

使用 r = y sqrt ( 1 + (x/y) 2 ) 然后你正在寻找一个术语 t = ( - 1 / 2 ) n+1 0.36 n，其幅度小于 0.001, log ( 0.002 ) > n log ( 0.18 ) 或 n > 3.6，所以对 x^4 取项应该没问题。

于 2009-09-14T21:41:54.823 回答
2 
Y=10000
Y2=Y*Y
for x=0..Y2 do
  D[x]=sqrt(Y2+x*x)

norm(x,y)=
  if (y==0) x
  else if (x>y) norm(y,x) 
  else {
     s=Y/y
     D[round(x*s)]/s
  }

如果你的坐标是平滑的，那么这个想法可以通过线性插值来扩展。要获得更高的精度，请增加 Y。

这个想法是 s*(x,y) 位于 y=Y 线上，您已经为其预先计算了距离。得到距离，然后除以 s。

我假设您确实需要距离而不是平方。

您也许还可以找到一个通用的 sqrt 实现，它会牺牲一些准确性来提高速度，但我很难想象它会击败 FPU 能做什么。

通过线性插值，我的意思是更改D[round(x)]为：

f=floor(x)
a=x-f
D[f]*(1-a)+D[f+1]*a
于 2009-09-14T21:21:43.623 回答
1

这并不能真正回答您的问题，但可能会有所帮助...

我要问的第一个问题是：

  • “我需要 sqrt 吗？”。
  • “如果没有，我怎样才能减少sqrts的数量？”
  • 然后你的：“我可以用聪明的计算替换剩余的 sqrts 吗？”

所以我会开始：

  • 您是否需要确切的半径，或者半径平方是否可以接受？sqrt 有快速的近似值，但对于您的规范可能不够准确。
  • 您可以使用镜像象限或八分之一处理图像吗？通过批量处理相同半径值的所有像素，可以将计算次数减少 8 倍。
  • 你能预先计算半径值吗？您只需要一个大小为您正在处理的图像大小的四分之一（或八分之一）的表格，并且该表格只需要预先计算一次，然后重新用于算法的多次运行。

所以聪明的数学可能不是最快的解决方案。

于 2009-09-14T21:24:52.050 回答
1

好吧，总是在尝试优化您的 sqrt，我见过的最快的是旧的 carmack quake 3 sqrt：

http://betterexplained.com/articles/understanding-quakes-fast-inverse-square-root/

也就是说，由于 sqrt 是非线性的，因此您将无法沿直线进行简单的线性插值以获得结果。最好的办法是使用表查找，因为它可以让您快速访问数据。而且，由于您似乎是按整数进行迭代，因此表查找应该非常准确。

于 2009-09-14T21:35:18.457 回答
0

好吧，您可以从 x=0 开始镜像（您只需要计算 n>=0，并将这些结果复制到相应的 n<0）。之后，我会看看使用 sqrt(a^2+b^2) （或相应的 sin）上的导数来利用常数 dx。

如果这还不够准确，我可以指出这对于 SIMD 来说是一项非常好的工作，它将为您提供 SSE 和 VMX（以及着色器模型 2）上的倒数平方根运算。

于 2009-09-14T21:34:41.513 回答
0

这有点与HAKMEM 项目有关：

第 149 项（明斯基）：圆算法 这是在点绘图显示器上绘制几乎圆的一种优雅方法：

NEW X = OLD X - epsilon * OLD Y
NEW Y = OLD Y + epsilon * NEW(!) X

这会生成一个以原点为中心的非常圆的椭圆，其大小由初始点决定。epsilon 决定了循环点的角速度，对偏心率有轻微的影响。如果 epsilon 是 2 的幂，那么我们甚至不需要乘法，更不用说平方根、正弦和余弦了！“圆”将是完全稳定的，因为这些点很快就会变成周期性的。

当我试图在显示黑客中保存一个寄存器时，我错误地发明了圆形算法！Ben Gurley 只用了大约六到七条指令就完成了一个惊人的显示技巧，这真是一个奇迹。但它基本上是面向线的。我突然想到有曲线会很令人兴奋，我试图用最少的指令来获得曲线显示技巧。

于 2009-09-15T13:40:44.850 回答