如何计算反向触发（和sqrt）函数（C语言）的浮点运算中的舍入错误？

Question

我有一个相当复杂的函数，它需要几个double值，它们表示以纬度和经度为弧度的形式（大小，纬度，经度）的3维空间中的两个向量以及一个角度。函数的作用是将第一个向量绕第二个向量旋转指定的角度并返回合成向量。我已经验证该代码在逻辑上是正确的并且可行。

函数的预期用途是用于图形，所以不需要双精度;然而，在目标平台上，使用浮点数（sinf，cosf，atan2f，asinf，acosf和sqrtf）的trig（和sqrt）函数在double上的运行速度要快于浮点数（可能是因为计算这些值的指令实际上可能需要double;如果传递一个float值，则该值必须转换为double值，这要求将其复制到具有更多内存的区域 - 即开销）。因此，函数中涉及的所有变量都是双精度的。

这是问题：我试图优化我的功能，以便它可以被称为每秒更多次。因此，我将这些调用替换为sin，cos，sqrt等调用这些函数的浮点版本，因为它们导致总体速度提高3-4倍。这适用于几乎所有的输入;然而，如果输入矢量接近于与标准单位矢量（i，j或k）平行，则各种函数的舍入误差会足以导致稍后调用sqrtf或反函数trig函数（asinf，acosf， atan2f）在这些函数的域之外传递参数仅仅是。因此，我留下了这个困境：要么我只能调用双精度函数，并避免这个问题（并且每秒最多可以有大约1,300,000个向量操作），或者我可以试着想出一些东西其他。最终，我想要一种方法来消除对反向trig函数的输入以处理边缘情况（对于sqrt，这样做很简单：只需使用abs）。分支不是一种选择，因为即使是单个条件语句也会增加很多开销，以至于性能增益都会丢失。

那么，有什么想法？

编辑：有人对我使用双打和浮点操作表示困惑。如果我实际上将所有值存储在双倍大小的容器中（I.E. double-type变量），则将函数存储在浮动大小容器中的速度要快得多。但是，由于显而易见的原因，浮点精度触发操作比双精度触发操作要快。

Answer 1

基本上，你需要找到一个解决你的问题的numerically stable算法。这种事情没有通用的解决方案，需要根据您的具体情况完成，如果单个步骤使用的概念如condition number。如果潜在的问题本身是病态的，事实上可能是不可能的。

Answer 2

你看过Graphics Programming Black Book还是将计算交给GPU？

Answer 3

单精度浮点固有地引入错误。所以，你需要通过使用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 

是很容易的，但是这是混乱的。限制域输入有点棘手，但更好。关键是要使用conditional move operators，这在一般做这样的事情

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

在单个运算

，而不会产生一个分支。

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

条件此举是更好的配对与位与一个comparison operator在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标量运算发射这种用于ternaries格局。您可以使用MSVC上的编译器标记/arch:sse2或GCC上的-mfpmath=sse来完成此操作。

在PowerPC和许多其他的RISC体系结构，fsel()是硬件操作码和因此通常编译器本征为好。