c - 我*不*想要函数 exp 的正确舍入

Question

在 Debian 系统上 C 数学库的 GCC 实现显然具有符合 (IEEE 754-2008) 的函数实现exp ，暗示四舍五入应始终正确:

(from Wikipedia) The IEEE floating point standard guarantees that add, subtract, multiply, divide, fused multiply–add, square root, and floating point remainder will give the correctly rounded result of the infinite precision operation. No such guarantee was given in the 1985 standard for more complex functions and they are typically only accurate to within the last bit at best. However, the 2008 standard guarantees that conforming implementations will give correctly rounded results which respect the active rounding mode; implementation of the functions, however, is optional.

事实证明，我遇到了这个功能实际上有阻碍的情况，因为 exp 的确切结果函数通常几乎恰好位于两个连续 double 之间的中间值 (1)，然后程序进行大量的进一步计算，速度损失高达 400 (!) 倍:这实际上是对我的解释(病态的:-S)Question #43530011 .

(1) 更准确地说，这发生在 exp 的参数时结果是 (2 k + 1) × 2^-53 的形式，其中 k 是一个相当小的整数(例如 242) .特别是，pow (1. + x, 0.5) 涉及的计算倾向于打电话 exp当 x 时有这样的论点数量级为2^-44。

由于在某些情况下正确舍入的实现可能非常耗时，我猜开发人员也会设计出一种方法来获得稍微不太精确的结果(比如，最多只能达到 0.6 ULP 或类似的东西)在给定范围内参数的每个值(大致)有界的时间内......(2)

……但是怎么做呢？？

(2) 我的意思是我只是不希望参数的某些异常值像 (2 k + 1) × 2^-53 会是比大多数相同数量级的值更耗时；但我当然不介意参数的某些特殊值是否变得更快，或者大参数(绝对值)是否需要更长的计算时间。

这是一个显示该现象的最小程序:

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>

int main (void)
 {
  int i;
  double a, c;
  c = 0;
  clock_t start = clock ();
  for (i = 0; i < 1e6; ++i) // Doing a large number of times the same type of computation with different values, to smoothen random fluctuations.
   {
    a = (double) (1 + 2 * (rand () % 0x400)) / 0x20000000000000; // "a" has only a few significant digits, and its last non-zero digit is at (fixed-point) position 53.
    c += exp (a); // Just to be sure that the compiler will actually perform the computation of exp (a).
   }
  clock_t stop = clock ();
  printf ("%e\n", c); // Just to be sure that the compiler will actually perform the computation.
  printf ("Clock time spent: %d\n", stop - start);
  return 0;
 }

现在在 gcc -std=c99 program53.c -lm -o program53 之后:

$ ./program53
1.000000e+06
Clock time spent: 13470008
$ ./program53 
1.000000e+06
Clock time spent: 13292721
$ ./program53 
1.000000e+06
Clock time spent: 13201616

另一方面，program52和 program54 (通过将 0x20000000000000 分别替换为 0x10000000000000 和 0x40000000000000 得到):

$ ./program52
1.000000e+06
Clock time spent: 83594
$ ./program52
1.000000e+06
Clock time spent: 69095
$ ./program52
1.000000e+06
Clock time spent: 54694
$ ./program54
1.000000e+06
Clock time spent: 86151
$ ./program54
1.000000e+06
Clock time spent: 74209
$ ./program54
1.000000e+06
Clock time spent: 78612

注意，这种现象是依赖于实现的!显然，在常见的实现中，只有 Debian 系统(包括 Ubuntu)显示这种现象。

P.-S.:我希望我的问题不是重复的:我彻底搜索了一个类似的问题但没有成功，但也许我确实注意到使用了相关的关键字...... :-/

Answer 1

回答关于为什么需要库函数才能给出正确舍入结果的一般问题:

float 很难，而且常常违反直觉。不是每个程序员都读过 what they should have .当库过去允许一些稍微不准确的舍入时，人们提示库函数的精度，因为他们不准确的计算不可避免地出错并产生废话。作为回应，图书馆的作者们把他们的图书馆做得恰到好处，所以现在人们不能把责任推给他们。

在许多情况下，有关浮点算法的特定知识可以显着提高准确性和/或性能，例如在测试用例中:

以exp()非常接近 0 的数字在 float 中是有问题的，因为结果是接近 1 的数字虽然所有精度都在与一的差中，因此大多数有效数字都丢失了。计算 exp(x) - 1 更精确(并且在此测试用例中明显更快)通过 C 数学库函数 expm1(x) .如果exp()本身是真的需要的，它仍然要快得多 expm1(x) + 1 .

计算 log(1 + x) 也存在类似的问题, 其中有函数 log1p(x) .

加速提供的测试用例的快速修复:

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>

int main (void)
{
  int i;
  double a, c;
  c = 0;
  clock_t start = clock ();
  for (i = 0; i < 1e6; ++i) // Doing a large number of times the same type of computation with different values, to smoothen random fluctuations.
    {
      a = (double) (1 + 2 * (rand () % 0x400)) / 0x20000000000000; // "a" has only a few significant digits, and its last non-zero digit is at (fixed-point) position 53.
      c += expm1 (a) + 1; // replace exp() with expm1() + 1
    }
  clock_t stop = clock ();
  printf ("%e\n", c); // Just to be sure that the compiler will actually perform the computation.
  printf ("Clock time spent: %d\n", stop - start);
  return 0;
}

对于这个案例，我机器上的时间是这样的:

原始代码

1.000000e+06

Clock time spent: 21543338

修改后的代码

1.000000e+06

Clock time spent: 55076

对伴随权衡具有高级知识的程序员有时可能会考虑在精度不重要的情况下使用近似结果

对于有经验的程序员，可以使用 Newton-Raphson、Taylor 或 Maclaurin 多项式等方法编写慢函数的近似实现，特别是来自 Intel 的 MKL、AMD 的 AMCL 等库的不精确舍入的特殊函数，放宽 float -点编译器的标准合规性，将精度降低到 ieee754 binary32 ( float )，或这些的组合。

请注意，更好地描述问题会带来更好的答案。