如何解释反三角(和sqrt)函数的浮点算法中的舍入误差( C)？

Question

我有一个相当复杂的函数，它有几个双值，表示形式(大小、纬度、经度)的三个空间中的两个向量，其中纬度和经度是以弧度为单位的，还有一个角度。该函数的目的是将第一个矢量按指定的角度围绕第二个矢量旋转，并返回所得到的向量。我已经验证了代码在逻辑上是正确的，并且工作正常。

函数的预期用途是用于图形，因此不需要双重精度；但是，在目标平台上，接受浮点数的trig (和sqrt)函数(特别是atan2f、asinf、acosf和sqrtf )在双倍上工作的速度比浮点数快(可能是因为计算这些值的指令实际上可能需要双倍；如果传递浮点数，则必须将值转换为双值，这需要将其复制到内存更多的区域--即开销)。因此，函数所涉及的所有变量都是双精度的。

问题是:我正在尝试优化我的函数，以便它可以被每秒调用更多次。因此，我已经将对sin、cos、sqrt等的调用替换为对这些函数的浮点版本的调用，因为它们的总体速度增长了3-4倍。这对几乎所有的输入都有效；但是，如果输入向量接近于标准单元向量(i、j或k)，那么各种函数的舍入错误就会累积起来，从而导致以后对sqrtf或逆trig函数(asinf、acosf、atan2f)的调用传递那些刚好在这些函数域之外的参数。

因此，我陷入了这样的困境:要么我只能调用双精度函数并避免这个问题(并以每秒大约1300,000次向量运算结束)，要么我可以尝试想出其他的方法。最后，我希望有一种方法来清理反三角函数的输入，以处理边缘情况(对于sqrt来说，这样做很简单:只需使用abs)。分支不是一个选项，因为即使是一个条件语句也会增加大量开销，导致任何性能增益都会丢失。

有什么想法吗？

编辑:有人对我的双打和浮点操作表示混淆。如果实际将所有值存储在双大小容器(即双类型变量)中，则函数的速度要比将其存储在浮点数容器中快得多。然而，由于明显的原因，浮点精度trig操作比双精度trig操作更快。

Answer 1

基本上，您需要找到解决问题的数值稳定算法。这类事情没有通用的解决方案，需要使用概念(如条件数 )(如果是单独的步骤)对您的具体情况进行处理。事实上，如果根本问题本身是病态的，这可能是不可能的。

Answer 2

单精度浮点本身就会带来误差。所以，你需要建立你的数学，这样所有的比较都有一定程度的“斜率”，使用epsilon因子，你需要净化输入到有限域的函数。

前者在分枝时很容易

bool IsAlmostEqual( float a, float b ) { return fabs(a-b) < 0.001f; } // or
bool IsAlmostEqual( float a, float b ) { return fabs(a-b) < (a * 0.0001f); } // for relative error

但那太乱了。夹紧域输入是有点棘手，但更好。关键是使用条件移动算子，它通常会执行以下操作

float ExampleOfConditionalMoveIntrinsic( float comparand, float a, float b ) 
{ return comparand >= 0.0f ? a : b ; }

在一次行动中，没有一根树枝。

这取决于不同的架构。在x87浮点单元上，您可以使用FCMOV条件-移动op来完成它，但是这很笨拙，因为这取决于以前设置的条件标志，所以速度很慢。另外，cmov没有一致的编译器。这就是为什么我们在可能的情况下避免使用x87浮点来支持SSE2标量数学的原因之一。

条件移动在SSE中得到更好的支持，方法是将比较算子与位和。即使对于标量数学，这也是更好的选择：

// assuming you've already used _mm_load_ss to load your floats onto registers 
__m128 fsel( __m128 comparand, __m128 a, __m128 b ) 
{
    __m128 zero = {0,0,0,0};
    // set low word of mask to all 1s if comparand > 0
    __m128 mask = _mm_cmpgt_ss( comparand, zero );  
    a = _mm_and_ss( a, mask );    // a = a & mask 
    b = _mm_andnot_ss( mask, b ); // b = ~mask & b
    return _mm_or_ss( a, b );     // return a | b
    }
}

当启用SSE2标量数学时，编译器可以更好地(但不是很好)为三元组发出这种模式。您可以使用MSVC上的编译器标志/arch:sse2或GCC上的-mfpmath=sse来实现这一点。

在PowerPC和许多其他RISC体系结构上，fsel()是一种硬件操作码，因此通常也是编译器固有的。

Answer 3

你看过图形编程黑书或者把计算交给你的GPU了吗？