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

我继续关于浮点运算的系列文章。在第一篇文章中，我给出了一个小的数学介绍，并通​​过具有不同“陷阱”的程序示例展示了最简单，最明显的正弦计算方法。今天的文章将在样式上稍有不同。这里不会进行任何练习，但是我们将更深入地研究数学并深入研究神圣的圣物-标准库代码。在第一篇文章的结尾，我还将回答这个问题。所以走吧

再次数学

显然，泰勒级数的绝对值减小得越快，获得所需精度所需的项就越少。因此，结果似乎会更加准确（下面将对此进行更详细的讨论）。为了进行比较，例如，以泰勒级数的第七次项（$$$\frac{x^7}{7!}$） 在 $$$x=1.0$ 和 $$$x=0.1$... 表达式值将是$$$1.198\cdot 10^{-4}$ 和 $$$1.198\cdot 10^{-11}$分别。有很大的不同，不是吗？因此，让我们尝试找到一种减小正弦函数计算间隔上限的方法。

围绕给定值进行级数展开

为了理解这种方法，我们需要回到研究所的第一年，并回顾泰勒级数（wiki）的定义。简而言之：知道某个点的函数及其派生类，您可以通过扩展为泰勒级数来找到该点附近函数的值。对于正弦函数，这意味着

$$

$$\sin(x_0+\Delta x) \approx \sin(x_0)+\sin'(x_0)\frac{\Delta x}{1!}+\sin''(x_0)\frac{\Delta x^2}{2!}+\sin'''(x_0)\frac{\Delta x^3}{3!}+\cdots$$

从实践的角度来看，这种方法给我们带来了什么？想象一下，我们与$$$0$ 之前 $$$\pi/2$... 让我们在此间隔上选择10个线性分布的点（选择不是最优的）：$$$x_0=0$， $$$x_1=\pi/20$， $$$x_2=2\pi/20$， $$$x_i=i\cdot \pi/20$... 对于每个点，在该点上使用正弦及其导数计算板。现在，您可以修改函数，以便在获取值时$$$x$ 该函数取最接近的值 $$$x_i$ 并围绕该点连续布置 $$$x_i$，而不是零（$$$\Delta x = x - x_i$）。

使用三角变换

如果您再回到高年级，那么您会记住一个非常重要的公式：

$$

$$\sin(x_0+\Delta x)=\sin(x_0)\cos(\Delta x)+\cos(x_0)\sin(\Delta x)$$

然后一切都与上一段相同。我们选择区间内的点，为其计算正弦和余弦，并在调用正弦函数时，寻找最接近的点，并使用上述公式，使用较小的值计算正弦$$$\Delta x$...

考虑一下这两种方法中的哪一种是更好的选择，但是现在我们将从数学转向实用计算。

浮点世界中乘法的分布特性

我不得不向互联网寻求建议，这叫什么 $$$a\cdot(b+c)=a\cdot b + a\cdot c$... 事实证明这是一种分配属性。让我们回到我在第一部分结尾处提出的问题。即，为什么数学上等价的表达式$$$a\cdot(b+c)$ 和 $$$a\cdot b + a\cdot c$可以在浮点计算中给出不同的结果？举例说明最简单的方法。让我们假设一个使用十进制格式的浮点数且精度为4位的系统。让我们假装$$$b=1.000\text{E}0$， $$$c=1.234\text{E-}2$和 $$$a=5.678\text{E-}2$... 首先，让我们使用带括号的表达式并逐步进行计算，并记住在每一步都进行四舍五入：

1）$$$b+c=1.000\text{E}0 + 1.234\text{E-}2=1.012\text{E}0$

2） $$$a\cdot(b+c)=5.678\text{E-}2\times 1.012\text{E}0=5.746\text{E-}2$

收到回应 $$$y=5.746\text{E-}2$

现在，让我们以相同的方式逐步计算第二个表达式：

$$$a\cdot b + a\cdot c=5.678\text{E-}2\times 1.000\text{E}0 + 5.678\text{E-}2\times 1.234\text{E-}2=5.678\text{E-}2 + 7.007\text{E-}4=5.748\text{E-}2$

收到回应 $$$y=5.748\text{E-}2$

真正的答案是0.0574806652。



如您所见，第二种情况下获得的答案比第一种情况下更接近真实答案。如果我们用手指来解释这一点，那么想象一下，在第一种情况下，我们将数字加到1.0$$$c=1.234\text{E-}2$我们只丢掉最后两位数字。他们没有了。在第二种情况下，丢弃发生在乘法之后的最后。那些。在第二种情况下，乘法运算更为准确。



看来您可以完成此操作，但是请仔细研究第一种方法，并告诉我计算结果将是$$$b+c-b$... 而且...我们有一种舍入浮点数的方法！不要错过这个例子。给自己时间解决这个问题。我们将在本文和以下文章中进一步大量使用数字舍入。



让我们注意到该表达式的另一个功能。想象一下，变量的4位精度对我们来说还不够。该怎么办？在这里，我们已经有了答案-以表格形式表示数字$$$b+c$并将其作为两位数的总和存储在内存中。并且，因此，分别对两个术语执行运算（例如乘法）。在增加两个浮点数而又不损失精度的文章中对此技术进行了更详细的描述。



在上一篇文章中，我还写道$$$a\cdot(b+c)$有一个令人不快的功能。如下。数$$$c$ 总是在数字的最后一个有效数字处被截断 $$$b$... 这意味着无论数量$$$c$， 如果一个 $$$b+c \neq b$，那么即使对于较小的符号，也总是可能出现最后一个符号的错误 $$$c$... 下一章的方法中不允许这样做。

以GNU库为例，它如何工作

如何？您已经选择了本文开头介绍的两种方法中的哪一种来精确计算正弦？无论选择哪种方法，都是正确的。而且，它们是绝对相同的。相信我，看看吧。下面我将使用学校公式。它们更容易解释。



有了上一篇文章和本文中获得的知识，您可以轻松地了解标准库的代码。让我们打开s_sin.c文件并在其中找到__sin函数：



其代码非常简单。很容易理解，它会根据输入变量的限制调用一组不同的函数。在本文中，我们将讨论角度2 ^ -26 <| x | <0.855469的218-224代码部分。您可以看到在代码的这一部分中，调用了do_sin（x，0）函数。我们将更详细地介绍此功能：

1. , dx=0 .
2. 129-130 , abs(x)<0.126, .. x , . , , , .
3. 136-137. , . x 2 . u x. , 0.345678. u=0.34, 0.005678.
4. 140-142. ( s ) ( c ) x . , cos(x)=1-c, 1.0, (. ), .
5. 143. u. , u=0.34 34. sin(u)=sn+ssn, cos(u)=cs+ccs. sn cs — «» u, ssn ccs — .
6. 144-145. sin(u+x)=(sn+ssn)*(1-c)+(cs+ccs)*s. , , 144-145. — .

实际上，我仅描述了以这种方式计算正弦的最简单部分。剩下的数学很多。例如，如何计算表及其中元素的大小？魔法数字0.126和0.855469从何而来？什么时候用泰勒数来计算呢？对泰勒级数的系数进行校正以完善结果。



当然，所有这些都是有趣的，但是，客观地讲，所提出的方法有许多缺点：必须同时计算正弦和余弦（c），这需要泰勒级数^1的两倍计算。如我们所见，通过表格值进行乘法也不是免费的。同样，将3520字节的表存储在RAM中当然不是问题，但是访问它（即使在高速缓存中）也可能很昂贵。



因此，在下一部分中，我们将尝试摆脱极板并直接在区间[0.126，0.855469]中计算正弦值，但比第一章更准确。



在结束之前-快速机智的问题。在此示例中，big的数字是52776558133248 = 3 * 2 ⁴⁴。这个数字是从哪里来的，例如不是2 ⁴⁵？我将更精确地提出这个问题。为什么在对数字进行四舍五入时，数字3 * 2 ^N是最优的，而不是例如2 ^{N +1}？另一个问题，您应该选择哪个N将数字四舍五入为整数？



¹值得注意的是，当从相同角度同时计算正弦和余弦时，这种方法的显着优势将显现出来。可以免费计算第二个函数。