使用正弦函数示例可以快速准确地计算浮点数。第3部分：定点

我们继续讲课的周期（第1部分和第2部分）。在第2部分中，我们研究了libm库中的内容，在这项工作中，我们将尝试稍微更改do_sin函数以提高其准确性和速度。我将再次引用此函数do_sin）：







如上一篇文章132-145部分所示。在[0.126，0.855469]范围内对x执行。好。让我们尝试编写一个函数，该函数在给定的范围内会更准确，甚至可能更快。



我们将使用它的方式非常明显。有必要扩大计算的准确性，以包括更多的小数位。显而易见的解决方案是选择long double类型，对其进行计数，然后转换回去。就准确性而言，解决方案应该不错，但就性能而言，可能会有问题。毕竟，long double是一种非常奇特的数据类型，它在现代处理器中的支持也不是优先事项。在x86_64上，此数据类型的SSE / AVX指令不起作用。数学协处理器将被“吹走”。



那你应该选择什么呢？让我们仔细看一下参数和函数限制。



它们在1.0区域中。那些。实际上，我们不需要浮点数。让我们在计算函数时使用64位整数。这将使我们额外获得10-11位的原始精度。让我们弄清楚如何使用这些数字。用这种格式的数字表示为a / d，其中a是整数，并且d是除数，我们为所有变量选择常数，然后“存储在我们的内存中”，而不是计算机的内存中。以下是此类号码的一些操作：

$$

$$\frac{c}{d}=\frac{a}{d}\pm\frac{b}{d}=\frac{a\pm b}{d}\\ \frac{c}{d}=\frac{a}{d}\cdot\frac{b}{d}=\frac{a \cdot b}{d^2}\\ \frac{c}{d}=\frac{a}{d}\cdot x=\frac{a \cdot x}{d}$$

如您所见，没有什么复杂的。最后一个公式显示乘以任何整数。还要注意的一件很明显的事是，将两个大小为N的无符号整数变量相乘的结果通常是大小最大为2 * N（含）的数。加法可能导致最多1个额外的位溢出。



让我们尝试选择除数d。显然，在二进制世界中，最好将其选择为2的幂，以免除法，而只是移动寄存器。您应该选择二的幂？在乘法器指令中找到提示。例如，x86系统中的标准MUL指令将2个寄存器相乘并将结果也写入2个寄存器，其中寄存器之一是结果的“上部”，第二个寄存器是下部。



例如，如果我们有2个64位数字，那么结果将是将128位数字写入两个64位寄存器中。让我们将RH称为“大写”情况，并将RL称为“小写”情况¹。然后，在数学上，结果可以写成$$$R=R_H \cdot 2^{64}+R_L$... 现在我们使用上面的公式并为$$$d=2^{-64}$

$$

$$\frac{c}{d}=\frac{a}{2^{64}}\cdot\frac{b}{2^{64}}=\frac{a \cdot b}{2^{128}}=\frac{R_H \cdot 2^{64} + R_L}{2^{128}}=\frac{R_H + R_L \cdot 2^{-64}}{2^{64}}$$

结果表明，将这两个定点数相乘的结果是寄存器 $$$R=R_H$...



对于Aarch64系统，它甚至更容易。“ UMULH”指令将两个寄存器相乘，并将相乘的“上”部分写入第三个寄存器。



好吧。我们已经指定了一个定点数，但是仍然存在一个问题。负数。在泰勒级数中，展开带有可变符号。为了解决这个问题，我们根据Goner方法将用于计算多项式的公式转换为以下形式：

$$

$$\sin(x)\approx x(1-x^2(1/3!-x^2(1/5!-x^2(1/7!-x^2\cdot1/9!))))$$

检查它在数学上是否与原始公式完全相同。但是每个括号中都有许多形式$$$1/(2n + 1)! - x^2\cdot(\cdots)$总是积极的。那些。此转换允许将表达式评估为无符号整数。

constexpr mynumber toint    = {{0x00000000, 0x43F00000}};  /*  18446744073709551616 = 2^64     */
constexpr mynumber todouble = {{0x00000000, 0x3BF00000}};  /*  ~5.42101086242752217003726400434E-20 = 2^-64     */

double sin_e7(double xd) {
  uint64_t x = xd * toint.x;
  uint64_t xx = mul2(x, x);
  uint64_t res = tsx[19]; 
  for(int i = 17; i >= 3; i -= 2) {
    res = tsx[i] - mul2(res, xx);
  }
  res = mul2(res, xx);
  res = x - mul2(x, res);
  return res * todouble.x;
}

哪里 $$$tsx[i]=1/i!$以定点格式。这次，为方便起见，我将所有代码发布在fast_sine github上，摆脱了四进制，以便与clang和arm兼容。而且我改变了一些计算误差的方法。



下表两张表给出了标准正弦函数和定点函数的比较。第一张表显示了计算精度（对于x86_64和ARM完全相同）。第二张表是性能比较。







在测试过程中，使用MPFR库计算了“真”正弦值...最大ULP值被视为与“真”值的最大偏差。错误百分比-我们或libm二进制文件的正弦函数的计算值与四舍五入到正弦值不相等的情况下的数目。偏差的平均值表示计算误差的“方向”：对该值的高估或低估。从表中可以看出，我们的函数倾向​​于高估正弦值。这可以解决！谁说tsx值应该恰好等于泰勒级数的系数。一个相当明显的想法建议自己改变系数的值，以提高近似的准确性并消除误差的常数分量。正确地做出这样的变化是非常困难的。但是我们可以尝试。让我们举个例子tsx系数数组中的第4个值（tsx [3]），并将最后一个数字a更改为f。让我们重新启动程序并查看精度（sin_e7a）。看，有点，但是增加了！我们将此方法添加到我们的存钱罐。



现在，让我们看看性能。为了进行测试，我使用了当时可用的i5移动版和稍微超频的第四个树莓派（Raspberry PI 4 8GB），这两个系统都来自Ubuntu 20.04 x64发行版的GCC10。







我不假装在这些测量中更准确。根据处理器负载，可能会有百分之几十的变化。可以这样得出主要结论。在现代处理器²上，切换为整数运算不会提高性能。现代处理器中数量惊人的晶体管使快速执行复杂的计算成为可能。但是，我认为在Intel Atom等处理器以及弱控制器上，这种方法可以显着提高性能。你怎么看？



尽管此方法提高了准确性，但这种准确性的提高似乎更有趣而不是有用。在性能方面，这种方法可以在IoT中找到自己。但是对于高性能计算，它不再是主流。在当今世界，SSE / AVX / CUDA更喜欢使用并行函数计算。并在浮点运算中。没有MUL功能的并行类似物。该功能本身是对传统的致敬。



在下一章中，我将描述如何有效地使用AVX进行计算。再次，让我们进入libm代码并尝试对其进行改进。



¹在我所知的处理器中，没有这样名称的寄存器。例如，已经选择了名称。

²这里应该注意，我的ARM配备了最新版本的数学协处理器。如果处理器模拟浮点计算，则结果可能会大不相同。