我惊讶地发现很难用俄语找到有关此问题的信息，好像很少有人关心现代编译器中使用的数学库有时不能给出正确的舍入结果。我担心这种情况，因为我只是在从事此类数学库的开发工作。在外国文学中，这个问题已得到很好的解决，因此我决定依靠西方资源和一些个人经验，以一种流行的科学形式用俄语来介绍这个问题。

---

朋友，为方便起见，本文还提供视频演示格式（将近34分钟），此格式更适合那些难以在脑海中构建必要的数学图像的读者，因为演示中有很多说明性材料。视频中的信息与文章内容完全相同。请在方便时采取行动。

---

我再说一遍，这不是一本科学的文章，而是一本通俗的科学文章，阅读后，您将简短地了解这一点。

• 使用浮点运算的先验基本函数（exp，sin，log，cosh等）被错误地舍入，有时会在最后一位出错。
• 错误的原因并不总是在于开发人员的懒惰或资格低下，而是在一种基本情况下，这是现代科学尚未能够克服的。
• «», - .
• , , , , exp2(x) pow(2.0, x).

要了解本文的内容，您需要熟悉IEEE-754浮点格式。例如，您至少应该了解以下内容就足够了：0x400921FB54442D18-双精度格式的数字pi（binary64或double），也就是说，您仅了解此记录的含义；我不要求能够即时进行这种转换。我会提醒您有关本文中的舍入模式的知识，这是故事的重要组成部分。还希望知道“程序员”英语，因为会有西方文学的术语和引语，但是您可以通过在线翻译获得帮助。



首先使用示例，以便您立即了解对话的主题。现在，我将使用C ++给出代码，但是如果这不是您的语言，那么我相信您仍然会很容易理解所写内容。请看下面的代码：

#include <stdio.h>
#include <cmath>

int main() {
  float x = 0.00296957581304013729095458984375f;  // ,  .
  float z;
  z = exp2f(x);  // z = 2**x  .
  printf ("%.8f\n", z);  //      8   .
  z = powf(2.0f, x);  // z = 2**x  
  printf ("%.8f\n", z);  //   .
  return 0;
}

数字x刻意写有这么多的有效数字，以便可以在浮点类型中准确表示它，也就是说，编译器无需舍入即可将其转换为二进制代码。毕竟，您很清楚某些编译器无法四舍五入而没有错误（如果您不知道，请在评论中指出，我将另写一篇带有示例的文章）。在程序的下一步中，我们需要计算2 ^x，但是我们可以通过两种方式进行计算：函数exp2f（x）和两个powf（2.0f，x）的显式求幂。结果当然会有所不同，因为我在上面说基本功能不能在所有情况下都能正常工作，因此我特意选择了一个示例来说明这一点。这是输出结果：

1.00206053
1.00206041

四个编译器为我提供了这些值：Microsoft C ++（19.00.23026），Intel C ++ 15.0，GCC（6.3.0）和Clang（3.7.0）。它们的最低有效位不同。这是这些数字的十六进制代码：

0x3F804385  // 
0x3F804384  // 

请记住这个例子，稍后我们将讨论问题的实质，但是现在，为了让您有一个更清晰的印象，请参阅具有其他一些基本功能的双精度数据类型（double，binary64）的示例。我在表中显示结果。正确答案（如果有）末尾带有*。







希望您不要对我故意进行一些几乎找不到的完全独特的测试有印象吗？如果是这样，让我们​​跪下来，对float数据类型的2 ^x函数的所有可能的小数参数进行完整的枚举。显然，我们只对0到1之间的x值感兴趣，因为其他参数将给出仅在指数字段中的值不同且不感兴趣的结果。您自己了解：

$$

$$2^x = 2^{[x]}\cdot2^{\{x\}}.$$

编写了这样的程序（隐藏的文字将在下面）后，我检查了exp2f函数以及它在从0到1的间隔x上产生了多少个错误值。







从下面的程序中可以明显看出，测试的参数x的数量为197612997。事实证明，例如，Microsoft C ++错误地计算了其中几乎百分之一的2 ^x函数。亲爱的GCC和Clang粉丝们，不要为它们高兴而已，只是在这些编译器中正确实现了此功能，而在其他编译器中却充满了错误。





我不会用这些示例来让读者感到厌烦，这里的主要目的是表明先验函数的现代实现可能会错误地舍入最后一点，并且不同的编译器会在不同的地方犯错误，但是它们都无法正常工作。顺便说一句，IEEE-754标准在最后一位允许该错误（我将在下面讨论），但对我来说还是很奇怪：好的，double，这是一个大数据类型，但是浮点可以通过蛮力检查！难做吗？一点也不难，我已经展示了一个例子。



我们的枚举代码包含正确计算的“自写”功能2 ^x使用十分之一级的近似多项式，并且在几分钟之内完成编写，因为这样的多项式是自动得出的，例如在Maple计算机代数系统中。设置多项式以提供54位精度的精度就足够了（此函数为2 ^x）。为什么是54？但是您很快就会发现问题，就在我告诉您问题的本质并告诉您原理上为什么现在无法为四精度数据类型（binary128）创建快速且正确的先验函数，尽管理论上已经有尝试解决此问题的方法。

默认舍入及其问题

如果您不沉迷于数学库的开发，那么忘记根据IEEE-754标准的浮点数的默认舍入规则没有什么不妥。因此，我会提醒您。如果您对所有内容都记忆犹新，无论如何至少要看一下本节的结尾，您会感到很惊讶：我将向您展示一个很难取整数字的情况。



您可以通过名称轻松回忆起“向上舍入”（至正无穷大），“向下舍入”（至负无穷大）或“四舍五入”（如果有的话，请参阅维基百科））。程序员面临的主要困难是“四舍五入到四舍五入，但是如果从最接近的数到最后一位数字相等，则四舍五入”。是的，这就是这种舍入模式的翻译方式，西方文献简称为“舍入为偶”。



默认情况下使用此舍入模式，其工作方式如下。计算得出的结果是，如果尾数的长度大于结果数据类型可以容纳的长度，则将舍入到两个可能值中的最接近值。但是，当原始数字恰好位于两个最接近的数字之间的中间时，可能会出现这样的情况，然后选择结果，最后一位（舍入后）变成偶数，即等于零。考虑四个示例，您需要将二进制小数点后的位数四舍五入为两位：

1. 圆形1.00 1 001小数点后第三位是1，但随后还有另外一个第6位，这是1，这意味着舍必达，因为原来的数字接近1.01比1.00。
2. 1,001000. , 1,00 1,01, .
3. 1,011000. 1,01 1,10. , .
4. 1,010111. , 1,01, 1,10.

从这些示例中，似乎一切都很简单，但事实并非如此。事实是，有时我们无法确定我们是否真的处于两个值之间的中间位置。看一个例子。假设我们再次想舍入到小数点后两位：



1.00 1 0000000000000000000000000000000000001



现在对您来说很明显，舍入应该向上，即舍入到1.01。但是，您正在看的是小数点后40位的数字。如果您的算法不能提供40位精度而只能达到30位，该怎么办？然后它将给出另一个数字：



1.00 1 000000000000000000000000000



不知道在第40个位置（算法无法计算）时，您会珍惜一个，将这个数字四舍五入并得到1.00，这是错误的。您不正确地舍入了最后一位-这是我们讨论的主题。从上面的结果可以看出，为了仅使第二位正确，您必须计算最多40位的函数！哇！并且如果零的“机车”竟然还要更长？这就是我们将在下一节中讨论的内容。



顺便说一句，这是许多编译器在将浮点数的十进制表示法转换为结果二进制格式时所犯的错误。如果程序代码中的原始十进制数太靠近两个可精确表示的二进制值之间的中间值，则将无法正确舍入。但这不是本文的主题，而是另外一个故事的原因。

四舍五入的问题的实质

该问题表现出来的原因有两个。首先是故意拒绝耗时的计算，而倾向于速度。在这种情况下，只要遵守指定的精度，并且响应中将包含哪些位是第二要务。第二个原因是制表者的困境，这是我们谈话的主要主题。让我们更详细地考虑这两个原因。

第一个原因

您当然可以理解，先验函数的计算是通过一些近似方法实现的，例如，通过近似多项式的方法或什至（很少）通过级数展开的方法。为了使计算尽可能快地进行，开发人员同意对数值方法进行尽可能少的迭代（或者采用最小可能次数的多项式），只要算法允许误差不超过尾数最后一位的一半即可。在文献中，这写成0.5ulp（ulp =最后一位的单位）。



例如，如果我们在区间（0.5; 1）中讨论float类型的数字x，则值ulp = 2 ^-23。在间隔（1; 2）上ulp = 2 ^-22。换句话说，如果x在间隔（0; 1）上，则2 ^x将位于间隔（1,2）上，并且为了确保精度为0.5ulp，您需要大致选择EPS = 2 ^-23（因为在普通人中，我们将常数“ epsilon”表示为“错误”或“准确性”，根据您的喜好，请不要发现错误）。



对于应用计算，这已经足够了，但是对于几乎100％的程序员来说，最后一位可能与绝对结果不一致可能并不重要，因为对于他们而言，重要的不是位将是多少，而是准确度是什么。



对于那些不了解的人，我将以十进制数字系统为例。这是两个数字：1.999999和2.0。假设第一个是程序员收到的，第二个是如果我们有无限可能的话应该发生的标准。它们之间的差异仅为百万分之一，也就是说，答案的误差为EPS = 10 ^-6。但是，此答案中没有一个正确的数字。不好吗不，从应用程序的角度来看，这是紫色的，程序员会将答案四舍五入到小数点后两位，并且会收到2.00（例如，大约是货币2.00美元），他不需要其他任何东西，但事实上他将EPS = 10 ^-6放入我的程序中，然后做得好，为中间计算的误差留出余量，并正确解决了问题。



换句话说，不要混淆：正确位数（或位数）的精度和数量是两个不同的东西。那些需要准确性的人（这几乎是程序员的100％），所讨论的问题根本与那些人无关。任何需要位序列来匹配正确舍入的引用的人都非常担心此问题，例如，基本函数库的开发人员。尽管如此，对于每个人来说，了解这一点对于总体开发还是很有用的。



让我提醒您，这是问题的第一个方向：答案的最后部分可能是错误的，因为这是有意解决的。最主要的是要保持0.5ulp（或更高）的精度。因此，仅在极快的情况下才从此条件中选择数值算法。在这种情况下，该标准允许实现基本功能，而无需正确舍入最后一位。我引用[1，第12.1节]（英语）：

1985年版本的IEEE 754浮点算法标准未指定任何与基本函数有关的内容。这是因为多年以来一直认为，至少对于某些输入参数而言，正确舍入的函数会太慢。此后情况发生了变化，该标准的2008版建议（但不要求）对某些函数进行正确取整。

以下是推荐的功能，但不需要正确舍入：

第二个原因

最后，我们讨论了话题：制表者的困境（缩写为TMD）。我不能充分地将其名称翻译成俄语，它是由威廉·卡汉（William Kahan）（IEEE-754的创始人）在文章[2]中引入的。也许，如果您阅读该文章，您将理解为什么名称就是这个名称。简而言之，难题的实质是我们需要对函数z = f（x）进行绝对精确的舍入，就好像我们拥有可以完美处理的结果z的无限位记录一样。但是每个人都清楚，我们无法获得无限的序列。那要花几位呢？上面，我展示了一个示例，当我们需要查看结果的40位以便四舍五入后至少获得2位正确的位时。 TMD问题的实质是我们事先不知道，最多可以计算多少位以计算z的值，以便在舍入后获得所需的正确位数。如果有一百或一千怎么办？我们事先不知道！



例如，正如我所说，对于函数2 ^x，对于数据类型float（尾数的小数部分只有23位），我们需要以2 ^-54的精度执行计算，以便对所有可能的x参数无条件正确地进行舍入。通过穷举搜索获得此估计并不难，但是对于大多数其他功能，尤其是对于double或long double类型（如果知道的话，请输入“ class”），此类估计是未知的。



让我们已经了解了为什么会这样。我特意在本文中使用浮点数据类型给出了第一个示例，并请您记住它，因为在这种类型中只有32位，并且在其他数据类型中情况类似，因此更容易查看。



我们从数字x = 0.00296957581304013729095458984375开始，这是浮点数据类型中可精确表示的数字，也就是说，它被编写为可以直接转换为二进制浮点系统。我们计算2 ^x，如果我们有一个无限精度的计算器，那么我们应该得到（因此您可以检查一下，该计算是在在线系统WolframAlpha中完成的）：



1.0020604729652405753669743044108123031635398201893943954577320057 ...



让我们将此数字转换为二进制，例如64位就足够了：



1.00000000100001110000100 1 000000000000000000000000000000000001101111101



下划线的是舍入位（小数点后第24位）。问题：在哪里取整？上或下？显然，您知道这一点是因为您看到足够多的位并且可以做出决定。但是仔细看...



在舍入位之后，我们有29个零。这意味着我们非常接近两个最接近的数字之间的中间，并且只要向下舍入一点就足够了，因为四舍五入的方向会发生变化。但是问题是：这种转变在哪里？数值算法可以逐步，逐步地从各个侧面逼近精确值，直到在这“机车”中通过所有这29个零并且达到超过最后一个零的值的精度之前，我们不会知道舍入方向...如果实际上正确的答案应该是：



1.00000000100001110000100 0 11111111111111111111111111111，该怎么办？



然后舍入将向下。



直到精度达到小数点后第54位时，我们才知道这一点。当确切知道第54位时，我们将确切知道实际上更接近两个最接近的数字中的哪个。这样的数字称为最难舍入点[1，第12.3节]（舍入的临界点），而数字54称为硬度舍入，并在引用的书中用字母m表示。



舍入（m）的复杂度是位数，它是确保对于某个函数f（x）的所有自变量和预选范围，将函数f（x）正确舍入到最后一位所必需的最少位数（对于不同的舍入模式，可能会有所不同值m）。换句话说，对于数据类型float和参数x，取整模式的范围（0; 1）到“最接近的偶数”取整时间m = 54。这意味着对于间隔（0; 1）中的所有x，我们都可以将相同的精度ESP = 2 ^-54放入算法中，并且所有结果将正确舍入为二进制小数点后的23位。



实际上，某些算法能够提供准确的结果，并且基于53位甚至52位，蛮力表明了这一点，但是从理论上讲，您恰好需要54位。如果不是为了消除蛮力，我们就无法“欺骗”就像我在上面的蛮力程序中所做的那样，并节省了一些时间。我选了一个比应有的次数低的多项式，但是它仍然有效，只是因为我很幸运。



因此，无论舍入模式如何，我们都有两种可能的情况：在舍入区域中出现零的“蒸汽机车”，或在其附近出现“蒸汽机车”。用于计算先验函数f（x）的正确算法的任务是细化该函数的值，直到精度超过该“机车”的最后一位的值，并且直到由于计算f的数值算法的后续波动而变得明显为止为止（x）零不会变成一，反之亦然。一旦一切稳定下来，并且算法达到了超出“蒸汽机车”极限的精度，那么我们就可以舍入为无数位。此舍入将以正确的最后一位进行。但是，如何实现呢？

“拐杖”

如前所述，主要问题是获得一种算法来克服零位或舍入位后立即出现的机车。当机车克服并且我们将其视为一个整体时，这等效于以下事实：这些零或一已经被精确地计算出，并且我们已经确切地知道了现在朝哪个方向进行舍入。但是，如果我们不知道机车的长度，那么如何设计算法？

第一个“拐杖”

在读者看来，答案很明显：无限精确地进行算术，并故意放入过多的位数，如果不够用，则放入另一个并重新计算。一般来说，这是正确的。当计算机的速度和资源没有特别作用时，可以执行此操作。这种方法有一个名称：Ziv的多级策略[1，第12.3节]。其本质非常简单。该算法应支持多个级别的计算：快速的初步计算（在大多数情况下它是最终的），较慢但更准确的计算（在大多数紧急情况下为节省），甚至更慢但更准确的计算（当绝对“不好”时） “必须）等等。



在绝大多数情况下，足以使精度略高于0.5ulp，但是如果出现“机车”，则将其提高。只要保留“蒸汽机车”，我们就可以提高精度，直到严格地确定数值方法的进一步波动不会影响该“蒸汽机车”为止。因此，例如，在我们的情况下，如果达到ESP = 2 ^-54，则在第54位会出现一个单位，该单位“保护”机车不受零影响，并确保不会发生大于或等于2 ^-53的减法。零不会变成1，而是将舍入位拖到零。



这是一个流行的科学演示，与Ziv的四舍五入测试相同，可以在[1，第12章]或[3，第3节]中阅读，它一步一步显示了检查我们是否达到所需精度的速度。 10.5]。



这种方法的问题很明显。有必要设计一种算法来计算每个先验函数f（x），以便在计算的过程中可以提高计算的准确性。对于软件实现而言，这还不是那么令人恐惧，例如，牛顿的方法大致允许每次迭代将小数点后的精确位数加倍。您可以加倍，直到变得“足够”为止，尽管这是一个相当耗时的过程，但我必须承认牛顿法并不总是合理的，因为它需要计算反函数f ^-1（x），在某些情况下可能并不比计算f（x）本身简单。对于硬件实施，“ Ziva策略”是完全不合适的。硬编码到处理器中的算法必须使用已经预先设置的位数执行一系列操作，如果我们事先不知道该位数，则实现起来会很麻烦。拿股票？多少？



解决问题的概率方法[1，第12.6节]使我们能够估计m的值（请记住，这是位数，足够用于正确舍入）。事实证明，在概率意义上的“机车”的长度略大于该数字的尾数的长度。因此，在大多数情况下，取m等于尾数的两倍就足够了，只有在极少数情况下才需要取更多。我引用这项工作的作者的话：“我们推论，实际上，m必须略大于2p”（它们具有p-尾数的长度以及整数部分，即对于浮点数，p = 24）。在本文的进一步内容中，他们表明采用这种策略的错误概率接近于零，但仍为正值，这已通过实验得到证实。



然而，在某些情况下，必须取更大的m值，而最坏的情况事先也不知道。存在对最坏情况的理论估计[1，第12.7.2节]，但它们会产生不可思议的数百万个比特，这是不好的。这是引用的工作的表格（用于从-ln（2）到ln（2）的时间间隔上的函数exp（x））：

第二“拐杖”

实际上，m不会那么大。为了确定最坏的情况，应用了第二个“拐杖”，这被称为“详尽的预先计算”。对于数据类型float（32位），如果函数f具有一个参数（x），那么我们可以轻松地``运行''x的所有可能值。只有具有多个参数（在pow（x，y）中）的函数才会出现问题，对此我们无法想到。在检查了x的所有可能值之后，我们为每个函数f（x）和每个舍入模式计算常数m。然后，设计需要在硬件中实现的计算算法，以提供2 ^-m的精度。然后，在所有情况下都保证舍入f（x）是正确的。



对于双精度类型（64位），几乎不可能进行简单枚举。但是，他们正在整理！但是如何？答案在[4]中给出。我会简短地告诉你。



函数f（x）的域分为非常小的段，因此可以在每个段内用形式为b-ax的线性函数替换f（x）（当然，对于不同的段，系数a和b是不同的）。这些段的大小是通过分析计算得出的，因此，这样的线性函数实际上与每个段中的原始线性函数几乎没有区别。



然后，在执行了一些缩放和平移操作后，我们遇到以下问题：直线b轴能否“足够接近”到整数点？



事实证明，给出是或否的答案相对容易。也就是说，如果可能的危险点接近直线，则为“是”；如果原则上没有此类危险点，则为“否”。该方法的优点在于，在绝大多数情况下，实际上都会获得答案“否”，而很少获得答案“是”，这会迫使您仔细搜索该段以确定哪些特定点是关键的。



但是，对f（x）的参数进行迭代可以减少很多次，并且可以检测出double（binary64）和long double（80 bits！）之类的数字的断点。这是在超级计算机上完成的，当然还有视频卡，这些都是在采矿的空闲时间内完成的。但是，还没有人知道如何处理binary128数据类型。让我提醒您，此类数字的尾数的小数部分为112位。因此，在有关该主题的外国文献中，您仍然只能找到半哲学推理，从“我们希望...”（“我们希望...”）开始。



允许您快速确定直线在整数点附近通过的方法的详细信息在此处不合适。对于那些希望学习该过程的人，我建议更仔细地寻找寻找直线与Z ²之间距离的问题，例如，在文章[5]中。它描述了一种改进的算法，该算法在构建过程中类似于著名的Euclid算法，用于找到最大公约数。我将从[4]和[5]中得到相同的图片，它显示了问题的进一步转变：







有大量表格，其中包含每个先验函数在不同时间间隔舍入的最坏情况。它们可以在[1第12.8.4节]和[3，第10.5.3.2节]中找到，也可以在单独的文章中找到，例如，在[6]中。



我将通过从此类表中随机抽取一些示例。我强调这些并不是所有x的最坏情况，而只是在很小的间隔内，如果您有兴趣，请参阅源。







如何阅读表格？值x以十六进制浮点双符号指定。首先，正如预期的那样，存在一个前导，然后是尾数的小数部分的52位和字母P。此字母的意思是“乘以2乘以幂”，再加上一个度数。例如，P-23表示指定的尾数需要乘以2 ^-23。



进一步，假设函数f（x）的计算精度是无限的，并且前53位从中被切除（没有舍入！）。 f（x）列中指示的是这53个位（其中之一直到逗号）。后续位在最后一栏中指示。最后一列中的位序列的“度”符号表示位重复的次数，例如10 ⁵³1011表示首先到达等于1的位，然后是53个零然后是1011。然后是省略号，这意味着我们通常不需要所有其余的位。



此外，这是技术问题-我们知道单独采取功能的每个间隔的最坏情况，我们可以为该间隔选择一个近似值，以便它以其准确性覆盖最坏情况。仅用了几年的超级计算机计算，就可以创建基本功能的快速而准确的硬件实现。问题很小：至少要教编译器开发人员使用这些表。

为什么需要这个？

好问题！毕竟，我已经在上面反复说过，几乎100％的程序员不需要了解精确到正确舍入的最后一位的基本函数（通常他们甚至不需要一半的位），为什么科学家要驱动超级计算机并编译表来解决“无用”的问题？



首先，挑战是根本。并非为了精确舍入而进行精确舍入是很有趣的，但是从原理上理解如何解决这个有趣的问题，计算数学的解决方案将向我们揭示什么秘密？这些机密如何用于其他任务？基础科学-就像这样，您可以在数十年内完成某种“废话”，然后一百年后，由于这种“废话”，在其他领域发生了科学突破。



第二，代码可移植性问题。如果函数能够以所需的任何方式处理结果的最后几位，则意味着在不同的平台和不同的编译器上，可以获得略有不同的结果（即使在指定的错误范围内）。在某些情况下，这并不重要，但在某些情况下则可能很重要，尤其是当程序的错误出现在一个平台上，而又由于结果的不同而没有出现在另一个平台上时。但是，为什么我要向您描述与不同程序行为相关的众所周知的头痛呢？没有我，你知道这一切。拥有一个在所有平台上都可以完全相同运行的数学系统，这将是很棒的，无论它如何编译。那就是你需要正确做的 四舍五入最后一位。

来源清单

[1] Jean-Michel Muller，“基本函数：算法和实现”，2016



[2] William Kahan，“对数太聪明了”，2004



[3] Jean-Michel Muller，“浮点算术手册” ，2018



[4] VincentLefèvre，让-米歇尔·穆勒（Jean-Michel Muller），“迈向正确的舍入超越”，IEEE计算机上的交易，第一卷。47号 1998年11月11日。1235-1243



[5] VincentLefèvre。“关于线段和Z ²之间距离的新结果”。应用于精确舍入。第17届IEEE计算机算法研讨会-Arith'17，2005年6月，

美国马萨诸塞州科德角。pp.68-75



[6] VincentLefèvre，Jean-Michel Muller，“以双精度对基本函数进行正确舍入的最坏情况”，Rapport de recherche（信息自然语言和自然语言）n4044-2000年11月- 19页。