C 中三角函数的单精度参数约简

Question

我已经在 C 中实现了用单精度（32 位浮点）计算的三角函数（sin、cos、arctan）的一些近似值。它们精确到大约 +/- 2 ulp。

我的目标设备不支持任何<cmath>方法<math.h>。它不提供FMA，而是提供MAC ALU。ALU 和 LU 以 32 位格式进行计算。

我的反正切近似实际上是N.juffa 近似的修改版本，它在整个范围内近似反正切。正弦和余弦函数在 [-pi,pi] 范围内精确度高达 2 ulp。

我现在的目标是为正弦和余弦提供更大的输入范围（尽可能大，最好是 [FLT_MIN,FLT_MAX]），这使我能够减少参数。

我目前正在阅读不同的论文，例如KCNg 的 A RGUMENT REDUCTION FOR HUGE ARGUMENTS: Good to the Last Bit或有关此新参数缩减算法的论文，但我无法从中导出实现。

另外，我想提一下两个涉及相关问题的 stackoverflow 问题：有一种使用 matlab 和 c++ 的方法，它基于我链接的第一篇论文。它实际上使用 matlab、cmath 方法，并将输入限制为 [0,20.000]。另一种已经在评论中提到了。这是一种在 C 中实现 sin 和 cos 的方法，使用了我无法使用的各种 c 库。由于这两篇文章已经有几年的历史了，可能会有一些新的发现。

看起来这种情况下最常用的算法是将 2/pi 的数量精确存储到所需的位数，以便能够准确地进行模计算，同时避免取消。我的设备不提供大型 DMEM，这意味着无法查找具有数百位的大型查找表。该过程实际上在本参考文献的第 70 页上进行了描述，顺便说一句，它提供了许多有关浮点数学的有用信息。

所以我的问题是：是否有另一种有效的方法来减少正弦和余弦的参数以获得单精度，避免使用大的 LUT？上面提到的论文实际上专注于双精度并使用最多 1000 位数字，这不适合我的用例。

实际上我还没有找到任何 C 语言的实现，也没有找到针对单精度计算的实现，我将不胜感激任何类型的提示/链接/示例...

Answer 1

以下代码基于之前的答案，其中我演示了如何通过对数值较小的参数使用分割常量的 Cody-Waite 方法，以及对数值较大的参数使用 Payne-Hanek 方法，对三角函数执行相当准确的参数约简。震级。有关 Payne-Hanek 算法的详细信息，请参阅此处，有关 Cody-Waite 算法的详细信息，请参阅我之前的回答。

\n

在这里，我进行了必要的调整以适应询问者平台的限制，因为不支持 64 位类型，不支持融合乘加，并且辅助函数不可math.h用。我假设float映射到 IEEE-754binary32格式，并且有一种方法可以将这样的 32 位浮点数重新解释为 32 位无符号整数，反之亦然。我已经通过标准可移植习惯用法（即使用 ）实现了这种重新解释，memcpy()但可以选择适合未指定目标平台的其他方法，例如内联汇编、特定于机器的内在函数或易失性联合。

\n

由于此代码基本上是我以前的代码到更严格的环境的端口，因此它可能缺乏专门针对该环境的从头设计的优雅性。我基本上用一些位调整替换了frexp()辅助函数math.h，用 32 位整数对模拟 64 位整数计算，用 32 位定点计算替换双精度计算（效果比我预期的要好得多） ），并将所有 FMA 替换为未融合的等效项。

\n

重新处理论证还原的科迪-韦特部分需要做相当多的工作。显然，在没有 FMA 的情况下，我们需要确保常量 \xcf\x80/2 的组成部分中有足够数量的尾随零位（最低有效位除外），以确保乘积准确。我花了几个小时实验性地找出一个特定的分割，该分割可以提供准确的结果，但也将切换点尽可能推高到 Payne-Hanek 方法。

\n

当USE_FMA = 1指定时，测试应用程序的输出在使用高质量数学库编译时应类似于以下内容：

\n

Testing sinf ...  PASSED. max ulp err = 1.493253  diffsum = 337633490\nTesting cosf ...  PASSED. max ulp err = 1.495098  diffsum = 342020968\n

\n

随着USE_FMA = 0准确度稍微变差：

\n

Testing sinf ...  PASSED. max ulp err = 1.498012  diffsum = 359702532\nTesting cosf ...  PASSED. max ulp err = 1.504061  diffsum = 364682650\n

\n

输出diffsum是总体精度的粗略指标，此处显示所有输入中约 90% 会产生正确舍入的单精度响应。

\n

请注意，使用最严格的浮点设置和编译器提供的最高程度的 IEEE-754 遵守来编译代码非常重要。对于我用来开发和测试此代码的英特尔编译器，可以通过使用/fp:strict. 此外，用于参考的数学库的质量对于准确评估该单精度代码的 ulp 误差至关重要。英特尔编译器附带一个数学库，提供双精度基本数学函数，HA（高精度）变体的误差略高于 0.5 ulp。使用多精度参考库可能更好，但会减慢我的速度。

\n

Testing sinf ...  PASSED. max ulp err = 1.493253  diffsum = 337633490\nTesting cosf ...  PASSED. max ulp err = 1.495098  diffsum = 342020968\n

\n