我*不*想要函数 exp 的正确舍入
I do *not* want correct rounding for function exp
GCC C 数学库在 Debian 系统上的实现显然有一个 (符合 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）。特别是，当 x 的数量级为 2-44 时，pow (1. + x, 0.5) 涉及的计算倾向于使用这样的参数调用 exp .
由于在某些情况下正确舍入的实现可能非常耗时，我猜开发人员也会设计出一种方法来获得稍微不太精确的结果（比如，最多只能达到 0.6 ULP 或类似的东西this) 在给定范围内参数的 每个 值（大致）有界的时间内…… (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.: 我希望我的问题不是重复的：我彻底搜索了一个类似的问题但没有成功，但也许我确实注意到使用了相关的关键字...... :-/
回答有关为什么需要库函数才能提供正确舍入结果的一般问题：
浮点数很难，而且常常违反直觉。不是每个程序员都读过 what they should have。当库过去允许一些稍微不准确的舍入时，人们抱怨库函数的精度，因为他们不准确的计算不可避免地出错并产生废话。作为回应，图书馆的作者把他们的图书馆做得很圆，所以现在人们不能把责任推给他们。
在许多情况下，有关浮点算法的特定知识可以显着提高准确性 and/or 性能，如测试用例：
在浮点数中取非常接近 0 的数字的 exp() 是有问题的，因为结果是接近 1 的数字，而所有精度都在与一的差异，因此丢失了最重要的数字。通过 C 数学库函数 expm1(x) 计算 exp(x) - 1 更精确（并且在此测试用例中明显更快）。如果 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 等库的不精确舍入的专业函数，放宽编译器的浮点标准合规性，将精度降低到 ieee754 binary32 (float)，或这些的组合。
请注意，对问题的更好描述会带来更好的答案。
关于您对@EOF 的回答的评论，@NominalAnimal 的"write your own" 评论在这里似乎很简单，甚至微不足道，如下所示。
你上面的原始代码似乎有一个 exp() 的最大可能参数 a=(1+2*0x400)/0x2000...=4.55e-13 （这实际上应该是 2*0x3FF，我在你的 0x2000... 之后数了 13 个零，这使得它 2x16^13)。所以 4.55e-13 最大参数非常非常小。
然后平凡的泰勒展开式是 exp(a)=1+a+(a^2)/2+(a^3)/6+... 其中已经为您提供了如此小的参数的所有双精度。现在，您必须丢弃 1 部分，如上所述，然后它会减少到 expm1(a)=a*(1.+a* (1.+a/3.)/2.) 这应该会很快！只要确保 a 保持较小即可。如果它变得更大一点，只需添加下一项，a^4/24（你知道怎么做了吗？）。
>>编辑<<
我修改了 OP 的测试程序如下，以测试更多的东西（讨论遵循代码）
/*
i-do-not-want-correct-rounding-for-function-exp/44397261 */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#define BASE 16 /*denominator will be (multiplier)xBASE^EXPON*/
#define EXPON 13
#define taylorm1(a) (a*(1.+a*(1.+a/3.)/2.)) /*expm1() approx for small args*/
int main (int argc, char *argv[]) {
int N = (argc>1?atoi(argv[1]):1e6),
multiplier = (argc>2?atoi(argv[2]):2),
isexp = (argc>3?atoi(argv[3]):1); /* flags to turn on/off exp() */
int isexpm1 = 1; /* and expm1() for timing tests*/
int i, n=0;
double denom = ((double)multiplier)*pow((double)BASE,(double)EXPON);
double a, c=0.0, cm1=0.0, tm1=0.0;
clock_t start = clock();
n=0; c=cm1=tm1=0.0;
/* --- to smooth random fluctuations, do the same type of computation
a large number of (N) times with different values --- */
for (i=0; i<N; i++) {
n++;
a = (double)(1 + 2*(rand()%0x400)) / denom; /* "a" has only a few
significant digits, and its last non-zero
digit is at (fixed-point) position 53. */
if ( isexp ) c += exp(a); /* turn this off to time expm1() alone */
if ( isexpm1 ) { /* you can turn this off to time exp() alone, */
cm1 += expm1(a); /* but difference is negligible */
tm1 += taylorm1(a); }
} /* --- end-of-for(i) --- */
int nticks = (int)(clock()-start);
printf ("N=%d, denom=%dx%d^%d, Clock time: %d (%.2f secs)\n",
n, multiplier,BASE,EXPON,
nticks, ((double)nticks)/((double)CLOCKS_PER_SEC));
printf ("\t c=%.20e,\n\t c-n=%e, cm1=%e, tm1=%e\n",
c,c-(double)n,cm1,tm1);
return 0;
} /* --- end-of-function main() --- */
编译并运行它作为测试重现OP的0x2000...场景，或运行它带有（最多三个）可选参数 test #trials multiplier timeexp 其中 #trials 默认为 OP 的 1000000，并且 multipler 默认为 2 用于 OP 的 2x16^13（将其更改为 4，等等，用于她的其他测试）。对于最后一个参数，timeexp，输入 0 只执行 expm1()（和我的不必要的泰勒式）计算。这样做的目的是表明 OP 显示的错误时间案例随着 expm1() 消失，无论 multiplier.
所以默认运行s，test和test 1000000 4，产生（好吧，我调用程序四舍五入)...
bash-4.3$ ./rounding
N=1000000, denom=2x16^13, Clock time: 11155070 (11.16 secs)
c=1.00000000000000023283e+06,
c-n=2.328306e-10, cm1=1.136017e-07, tm1=1.136017e-07
bash-4.3$ ./rounding 1000000 4
N=1000000, denom=4x16^13, Clock time: 200211 (0.20 secs)
c=1.00000000000000011642e+06,
c-n=1.164153e-10, cm1=5.680083e-08, tm1=5.680083e-08
所以你首先要注意的是，OP 的 c-n 使用 exp() 与 cm1==tm1 使用 expm1() 和我的泰勒大约。如果你减少N他们达成一致，如下...
N=10, denom=2x16^13, Clock time: 941 (0.00 secs)
c=1.00000000000007140954e+01,
c-n=7.140954e-13, cm1=7.127632e-13, tm1=7.127632e-13
bash-4.3$ ./rounding 100
N=100, denom=2x16^13, Clock time: 5506 (0.01 secs)
c=1.00000000000010103918e+02,
c-n=1.010392e-11, cm1=1.008393e-11, tm1=1.008393e-11
bash-4.3$ ./rounding 1000
N=1000, denom=2x16^13, Clock time: 44196 (0.04 secs)
c=1.00000000000011345946e+03,
c-n=1.134595e-10, cm1=1.140730e-10, tm1=1.140730e-10
bash-4.3$ ./rounding 10000
N=10000, denom=2x16^13, Clock time: 227215 (0.23 secs)
c=1.00000000000002328306e+04,
c-n=2.328306e-10, cm1=1.131288e-09, tm1=1.131288e-09
bash-4.3$ ./rounding 100000
N=100000, denom=2x16^13, Clock time: 1206348 (1.21 secs)
c=1.00000000000000232831e+05,
c-n=2.328306e-10, cm1=1.133611e-08, tm1=1.133611e-08
关于 exp() 与 expm1() 的时间关系，请自行查看...
bash-4.3$ ./rounding 1000000 2
N=1000000, denom=2x16^13, Clock time: 11168388 (11.17 secs)
c=1.00000000000000023283e+06,
c-n=2.328306e-10, cm1=1.136017e-07, tm1=1.136017e-07
bash-4.3$ ./rounding 1000000 2 0
N=1000000, denom=2x16^13, Clock time: 24064 (0.02 secs)
c=0.00000000000000000000e+00,
c-n=-1.000000e+06, cm1=1.136017e-07, tm1=1.136017e-07
问题：你会注意到一旦 exp() 计算达到 N=10000 次试验，无论较大 N。不确定为什么会这样。
>>__SECOND EDIT__<<
好的，@EOF，"you made me look" 接受你的 "heirarchical accumulation" 评论。这确实可以使 exp() 总和更接近（更接近）（可能是正确的）expm1() 总和。修改后的代码紧随其后进行讨论。但是这里有一个讨论记录：回忆上面的 multiplier。那不见了，在同一个地方是 expon 所以分母现在是 2^expon ，默认是 53，匹配 OP 的默认值（我相信更好地匹配她的想法）。好的，这是代码...
/*
i-do-not-want-correct-rounding-for-function-exp/44397261 */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#define BASE 2 /*denominator=2^EXPON, 2^53=2x16^13 default */
#define EXPON 53
#define taylorm1(a) (a*(1.+a*(1.+a/3.)/2.)) /*expm1() approx for small args*/
int main (int argc, char *argv[]) {
int N = (argc>1?atoi(argv[1]):1e6),
expon = (argc>2?atoi(argv[2]):EXPON),
isexp = (argc>3?atoi(argv[3]):1), /* flags to turn on/off exp() */
ncparts = (argc>4?atoi(argv[4]):1), /* #partial sums for c */
binsize = (argc>5?atoi(argv[5]):10);/* #doubles to sum in each bin */
int isexpm1 = 1; /* and expm1() for timing tests*/
int i, n=0;
double denom = pow((double)BASE,(double)expon);
double a, c=0.0, cm1=0.0, tm1=0.0;
double csums[10], cbins[10][65537]; /* c partial sums and heirarchy */
int nbins[10], ibin=0; /* start at lowest level */
clock_t start = clock();
n=0; c=cm1=tm1=0.0;
if ( ncparts > 65536 ) ncparts=65536; /* array size check */
if ( ncparts > 1 ) for(i=0;i<ncparts;i++) cbins[0][i]=0.0; /*init bin#0*/
/* --- to smooth random fluctuations, do the same type of computation
a large number of (N) times with different values --- */
for (i=0; i<N; i++) {
n++;
a = (double)(1 + 2*(rand()%0x400)) / denom; /* "a" has only a few
significant digits, and its last non-zero
digit is at (fixed-point) position 53. */
if ( isexp ) { /* turn this off to time expm1() alone */
double expa = exp(a); /* exp(a) */
c += expa; /* just accumulate in a single "bin" */
if ( ncparts > 1 ) cbins[0][n%ncparts] += expa; } /* accum in ncparts */
if ( isexpm1 ) { /* you can turn this off to time exp() alone, */
cm1 += expm1(a); /* but difference is negligible */
tm1 += taylorm1(a); }
} /* --- end-of-for(i) --- */
int nticks = (int)(clock()-start);
if ( ncparts > 1 ) { /* need to sum the partial-sum bins */
nbins[ibin=0] = ncparts; /* lowest-level has everything */
while ( nbins[ibin] > binsize ) { /* need another heirarchy level */
if ( ibin >= 9 ) break; /* no more bins */
ibin++; /* next available heirarchy bin level */
nbins[ibin] = (nbins[ibin-1]+(binsize-1))/binsize; /*#bins this level*/
for(i=0;i<nbins[ibin];i++) cbins[ibin][i]=0.0; /* init bins */
for(i=0;i<nbins[ibin-1];i++) {
cbins[ibin][(i+1)%nbins[ibin]] += cbins[ibin-1][i]; /*accum in nbins*/
csums[ibin-1] += cbins[ibin-1][i]; } /* accumulate in "one bin" */
} /* --- end-of-while(nprevbins>binsize) --- */
for(i=0;i<nbins[ibin];i++) csums[ibin] += cbins[ibin][i]; /*highest level*/
} /* --- end-of-if(ncparts>1) --- */
printf ("N=%d, denom=%d^%d, Clock time: %d (%.2f secs)\n", n, BASE,expon,
nticks, ((double)nticks)/((double)CLOCKS_PER_SEC));
printf ("\t c=%.20e,\n\t c-n=%e, cm1=%e, tm1=%e\n",
c,c-(double)n,cm1,tm1);
if ( ncparts > 1 ) { printf("\t binsize=%d...\n",binsize);
for (i=0;i<=ibin;i++) /* display heirarchy */
printf("\t level#%d: #bins=%5d, c-n=%e\n",
i,nbins[i],csums[i]-(double)n); }
return 0;
} /* --- end-of-function main() --- */
好的，现在您可以注意到旧的 timeexp 后面有两个额外的命令行参数。它们是 ncparts，表示整个 #trials 将分配到的初始垃圾箱数量。因此，在层次结构的最低级别，每个 bin 应该（模数 bugs:) 的总和为 #trials/ncparts 双倍。之后的参数是 binsize，这将是每个连续级别的每个 bin 中总和的双倍数，直到最后一级具有更少（或等于）#bins 为 binsize。所以这是一个将 1000000 次试验分成 50000 个箱子的示例，这意味着最低级别为 20doubles/bin，此后为 5doubles/bin...
bash-4.3$ ./rounding 1000000 53 1 50000 5
N=1000000, denom=2^53, Clock time: 11129803 (11.13 secs)
c=1.00000000000000465661e+06,
c-n=4.656613e-09, cm1=1.136017e-07, tm1=1.136017e-07
binsize=5...
level#0: #bins=50000, c-n=4.656613e-09
level#1: #bins=10002, c-n=1.734588e-08
level#2: #bins= 2002, c-n=7.974450e-08
level#3: #bins= 402, c-n=1.059379e-07
level#4: #bins= 82, c-n=1.133885e-07
level#5: #bins= 18, c-n=1.136214e-07
level#6: #bins= 5, c-n=1.138542e-07
请注意 c-n 对于 exp() 如何很好地收敛到 expm1() 价值。但请注意它在第 5 级时是如何最好的，并且根本没有统一收敛。请注意，如果您将 #trials 分成 5000 个初始分箱，您会得到同样好的结果，
bash-4.3$ ./rounding 1000000 53 1 5000 5
N=1000000, denom=2^53, Clock time: 11165924 (11.17 secs)
c=1.00000000000003527384e+06,
c-n=3.527384e-08, cm1=1.136017e-07, tm1=1.136017e-07
binsize=5...
level#0: #bins= 5000, c-n=3.527384e-08
level#1: #bins= 1002, c-n=1.164153e-07
level#2: #bins= 202, c-n=1.158332e-07
level#3: #bins= 42, c-n=1.136214e-07
level#4: #bins= 10, c-n=1.137378e-07
level#5: #bins= 4, c-n=1.136214e-07
其实用ncparts和binsize玩玩好像灵敏度不高，也不一定总是"more is better"（即 binsize 的更少）。所以我不确定到底发生了什么。可能是一个（或两个）错误，或者可能是@EOF 的另一个问题......???
>>编辑——显示对添加的示例 "binary tree" 层次结构<<
根据@EOF 的评论添加了以下示例
（注意：重新复制前面的代码。我必须将每个下一级的 nbins[ibin] 计算编辑为 nbins[ibin]=(nbins[ibin-1]+(binsize-1))/binsize； 来自 nbins[ibin]=(nbins[ibin-1]+2*binsize)/binsize; 这是 "too conservative" 创建 ...16,8,4,2 序列)
bash-4.3$ ./rounding 1024 53 1 512 2
N=1024, denom=2^53, Clock time: 36750 (0.04 secs)
c=1.02400000000011573320e+03,
c-n=1.157332e-10, cm1=1.164226e-10, tm1=1.164226e-10
binsize=2...
level#0: #bins= 512, c-n=1.159606e-10
level#1: #bins= 256, c-n=1.166427e-10
level#2: #bins= 128, c-n=1.166427e-10
level#3: #bins= 64, c-n=1.161879e-10
level#4: #bins= 32, c-n=1.166427e-10
level#5: #bins= 16, c-n=1.166427e-10
level#6: #bins= 8, c-n=1.166427e-10
level#7: #bins= 4, c-n=1.166427e-10
level#8: #bins= 2, c-n=1.164153e-10
>>编辑——在下面的评论中展示@EOF 的优雅解决方案<<
"Pair addition" 可以优雅地递归完成，根据下面@EOF 的评论，我在这里复制。 （注意递归结束时的情况 0/1 以处理 n even/odd。）
/* Quoting from EOF's comment...
What I (EOF) proposed is effectively a binary tree of additions:
a+b+c+d+e+f+g+h as ((a+b)+(c+d))+((e+f)+(g+h)).
Like this: Add adjacent pairs of elements, this produces
a new sequence of n/2 elements.
Recurse until only one element is left.
(Note that this will require n/2 elements of storage,
rather than a fixed number of bins like your implementation) */
double trecu(double *vals, double sum, int n) {
int midn = n/2;
switch (n) {
case 0: break;
case 1: sum += *vals; break;
default: sum = trecu(vals+midn, trecu(vals,sum,midn), n-midn); break; }
return(sum);
}
这是 EOF 之前关于他的 trecu() 算法和代码的评论的 "answer"/跟进 "binary tree summation" 建议。 "Prerequisites" 在阅读本文之前正在阅读该讨论。把所有这些都收集在一个有组织的地方会很好，但我还没有这样做...
...我所做的是将 EOF 的 trecu() 构建到我通过修改 OP 的原始测试程序编写的前面答案中的测试程序中。但后来我发现 trecu() 生成的答案与 "plain sum" c 使用 exp()，而不是使用 expm1() 的总和 cm1，我们希望从更准确的二进制文件中得到树求和。
但是那个测试程序有点（可能是两位:）"convoluted"（或者，如 EOF 所说，"unreadable"），所以我写了一个单独的较小的测试程序，如下所示（带有示例运行s 和下面的讨论），单独 test/exercise trecu()。此外，我还将函数 bintreesum() 写入下面的代码中，其中 abstracts/encapsulates 我在前面的测试程序中嵌入的二叉树求和的迭代代码。在前面的例子中，我的迭代代码确实接近 cm1 答案，这就是为什么我期望 EOF 的递归 trecu() 做同样的事情。总而言之，下面发生了同样的事情——bintreesum() 仍然接近正确答案，而 trecu() 离正确答案更远，准确地再现了 "plain sum".
我们下面求和的就是sum(i),i=1...n，也就是众所周知的n(n+1)/2。但这不太正确——为了重现 OP 的问题，被加数不是单独的 sum(i)，而是 sum(1+i*10^(-e))，其中 e 可以在命令行上给出。所以，比如说，n=5，你得到的不是 15，而是 5.000...00015，或者对于 n=6，你得到 6.000...00021，等等。为了避免很长很长的格式，我 printf( ) sum-n 删除整数部分。好的？？？所以这是代码...
/* Quoting from EOF's comment...
What I (EOF) proposed is effectively a binary tree of additions:
a+b+c+d+e+f+g+h as ((a+b)+(c+d))+((e+f)+(g+h)).
Like this: Add adjacent pairs of elements, this produces
a new sequence of n/2 elements.
Recurse until only one element is left. */
#include <stdio.h>
#include <stdlib.h>
double trecu(double *vals, double sum, int n) {
int midn = n/2;
switch (n) {
case 0: break;
case 1: sum += *vals; break;
default: sum = trecu(vals+midn, trecu(vals,sum,midn), n-midn); break; }
return(sum);
} /* --- end-of-function trecu() --- */
double bintreesum(double *vals, int n, int binsize) {
double binsum = 0.0;
int nbin0 = (n+(binsize-1))/binsize,
nbin1 = (nbin0+(binsize-1))/binsize,
nbins[2] = { nbin0, nbin1 };
double *vbins[2] = {
(double *)malloc(nbin0*sizeof(double)),
(double *)malloc(nbin1*sizeof(double)) },
*vbin0=vbins[0], *vbin1=vbins[1];
int ibin=0, i;
for ( i=0; i<nbin0; i++ ) vbin0[i] = 0.0;
for ( i=0; i<n; i++ ) vbin0[i%nbin0] += vals[i];
while ( nbins[ibin] > 1 ) {
int jbin = 1-ibin; /* other bin, 0<-->1 */
nbins[jbin] = (nbins[ibin]+(binsize-1))/binsize;
for ( i=0; i<nbins[jbin]; i++ ) vbins[jbin][i] = 0.0;
for ( i=0; i<nbins[ibin]; i++ )
vbins[jbin][i%nbins[jbin]] += vbins[ibin][i];
ibin = jbin; /* swap bins for next pass */
} /* --- end-of-while(nbins[ibin]>0) --- */
binsum = vbins[ibin][0];
free((void *)vbins[0]); free((void *)vbins[1]);
return ( binsum );
} /* --- end-of-function bintreesum() --- */
#if defined(TESTTRECU)
#include <math.h>
#define MAXN (2000000)
int main(int argc, char *argv[]) {
int N = (argc>1? atoi(argv[1]) : 1000000 ),
e = (argc>2? atoi(argv[2]) : -10 ),
binsize = (argc>3? atoi(argv[3]) : 2 );
double tens = pow(10.0,(double)e);
double *vals = (double *)malloc(sizeof(double)*MAXN),
sum = 0.0;
double trecu(), bintreesum();
int i;
if ( N > MAXN ) N=MAXN;
for ( i=0; i<N; i++ ) vals[i] = 1.0 + tens*(double)(i+1);
for ( i=0; i<N; i++ ) sum += vals[i];
printf(" N=%d, Sum_i=1^N {1.0 + i*%.1e} - N = %.8e,\n"
"\t plain_sum-N = %.8e,\n"
"\t trecu-N = %.8e,\n"
"\t bintreesum-N = %.8e \n",
N, tens, tens*((double)N)*((double)(N+1))/2.0,
sum-(double)N,
trecu(vals,0.0,N)-(double)N,
bintreesum(vals,N,binsize)-(double)N );
} /* --- end-of-function main() --- */
#endif
因此，如果您将其保存为 trecu.c，然后将其编译为 cc DTESTTRECU trecu.c lm o trecu 然后 运行 具有零到三个可选的命令行参数，如 trecu #trials e binsize 默认值是 #trials=1000000（如 OP 的程序）、e=10 和 binsize=2（对于我的bintreesum() 函数执行二叉树求和而不是更大尺寸的 bins）。
下面是一些说明上述问题的测试结果，
bash-4.3$ ./trecu
N=1000000, Sum_i=1^N {1.0 + i*1.0e-10} - N = 5.00000500e+01,
plain_sum-N = 5.00000500e+01,
trecu-N = 5.00000500e+01,
bintreesum-N = 5.00000500e+01
bash-4.3$ ./trecu 1000000 -15
N=1000000, Sum_i=1^N {1.0 + i*1.0e-15} - N = 5.00000500e-04,
plain_sum-N = 5.01087168e-04,
trecu-N = 5.01087168e-04,
bintreesum-N = 5.00000548e-04
bash-4.3$
bash-4.3$ ./trecu 1000000 -16
N=1000000, Sum_i=1^N {1.0 + i*1.0e-16} - N = 5.00000500e-05,
plain_sum-N = 6.67552231e-05,
trecu-N = 6.67552231e-05,
bintreesum-N = 5.00001479e-05
bash-4.3$
bash-4.3$ ./trecu 1000000 -17
N=1000000, Sum_i=1^N {1.0 + i*1.0e-17} - N = 5.00000500e-06,
plain_sum-N = 0.00000000e+00,
trecu-N = 0.00000000e+00,
bintreesum-N = 4.99992166e-06
因此您可以看到，对于默认值 运行，e=10，每个人都做得很好。也就是说，上面写着 "Sum" 的行只是做了 n(n+1)/2 的事情，所以大概会显示正确的答案。下面的每个人都同意默认的 e=10 测试用例。但是对于 e=15 和 e=16 下面的情况，trecu() 与 plain_sum 完全一致，而 bintreesum 非常接近正确答案。最后，对于 e=17，plain_sum 和 trecu() 有 "disappeared"，而 bintreesum() 仍然很好地挂在那里。
所以 trecu() 正确地进行了求和，但是它的递归显然没有做 "binary tree" 我更直接的迭代 bintreesum() 显然做正确的事情。这确实表明，对于这些 1+epsilon 类型的情况，EOF 对 "binary tree summation" 的建议比 plain_sum 实现了相当大的改进。所以我们真的很想看看他的 trecu() 递归工作！！！当我最初看它时，我认为它确实有效。但是在他的 default: 案例中，那个双重递归（它有一个特殊的名称吗？）显然比我想象的更令人困惑（至少对我而言：）。就像我说的，它 是 求和，而不是 "binary tree"。
好的，那么谁愿意接受挑战并解释 trecu() 递归中发生了什么？而且，也许更重要的是，对其进行修复，使其按预期运行。谢谢
GCC C 数学库在 Debian 系统上的实现显然有一个 (符合 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）。特别是，当 x 的数量级为 2-44 时，pow (1. + x, 0.5) 涉及的计算倾向于使用这样的参数调用 exp .
由于在某些情况下正确舍入的实现可能非常耗时，我猜开发人员也会设计出一种方法来获得稍微不太精确的结果（比如，最多只能达到 0.6 ULP 或类似的东西this) 在给定范围内参数的 每个 值（大致）有界的时间内…… (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.: 我希望我的问题不是重复的：我彻底搜索了一个类似的问题但没有成功，但也许我确实注意到使用了相关的关键字...... :-/
回答有关为什么需要库函数才能提供正确舍入结果的一般问题：
浮点数很难，而且常常违反直觉。不是每个程序员都读过 what they should have。当库过去允许一些稍微不准确的舍入时，人们抱怨库函数的精度，因为他们不准确的计算不可避免地出错并产生废话。作为回应，图书馆的作者把他们的图书馆做得很圆，所以现在人们不能把责任推给他们。
在许多情况下，有关浮点算法的特定知识可以显着提高准确性 and/or 性能，如测试用例：
在浮点数中取非常接近 0 的数字的 exp() 是有问题的，因为结果是接近 1 的数字，而所有精度都在与一的差异，因此丢失了最重要的数字。通过 C 数学库函数 expm1(x) 计算 exp(x) - 1 更精确（并且在此测试用例中明显更快）。如果 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 等库的不精确舍入的专业函数，放宽编译器的浮点标准合规性，将精度降低到 ieee754 binary32 (float)，或这些的组合。
请注意，对问题的更好描述会带来更好的答案。
关于您对@EOF 的回答的评论，@NominalAnimal 的"write your own" 评论在这里似乎很简单，甚至微不足道，如下所示。
你上面的原始代码似乎有一个 exp() 的最大可能参数 a=(1+2*0x400)/0x2000...=4.55e-13 （这实际上应该是 2*0x3FF，我在你的 0x2000... 之后数了 13 个零，这使得它 2x16^13)。所以 4.55e-13 最大参数非常非常小。
然后平凡的泰勒展开式是 exp(a)=1+a+(a^2)/2+(a^3)/6+... 其中已经为您提供了如此小的参数的所有双精度。现在，您必须丢弃 1 部分，如上所述，然后它会减少到 expm1(a)=a*(1.+a* (1.+a/3.)/2.) 这应该会很快！只要确保 a 保持较小即可。如果它变得更大一点，只需添加下一项，a^4/24（你知道怎么做了吗？）。
>>编辑<<
我修改了 OP 的测试程序如下，以测试更多的东西（讨论遵循代码）
/*
i-do-not-want-correct-rounding-for-function-exp/44397261 */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#define BASE 16 /*denominator will be (multiplier)xBASE^EXPON*/
#define EXPON 13
#define taylorm1(a) (a*(1.+a*(1.+a/3.)/2.)) /*expm1() approx for small args*/
int main (int argc, char *argv[]) {
int N = (argc>1?atoi(argv[1]):1e6),
multiplier = (argc>2?atoi(argv[2]):2),
isexp = (argc>3?atoi(argv[3]):1); /* flags to turn on/off exp() */
int isexpm1 = 1; /* and expm1() for timing tests*/
int i, n=0;
double denom = ((double)multiplier)*pow((double)BASE,(double)EXPON);
double a, c=0.0, cm1=0.0, tm1=0.0;
clock_t start = clock();
n=0; c=cm1=tm1=0.0;
/* --- to smooth random fluctuations, do the same type of computation
a large number of (N) times with different values --- */
for (i=0; i<N; i++) {
n++;
a = (double)(1 + 2*(rand()%0x400)) / denom; /* "a" has only a few
significant digits, and its last non-zero
digit is at (fixed-point) position 53. */
if ( isexp ) c += exp(a); /* turn this off to time expm1() alone */
if ( isexpm1 ) { /* you can turn this off to time exp() alone, */
cm1 += expm1(a); /* but difference is negligible */
tm1 += taylorm1(a); }
} /* --- end-of-for(i) --- */
int nticks = (int)(clock()-start);
printf ("N=%d, denom=%dx%d^%d, Clock time: %d (%.2f secs)\n",
n, multiplier,BASE,EXPON,
nticks, ((double)nticks)/((double)CLOCKS_PER_SEC));
printf ("\t c=%.20e,\n\t c-n=%e, cm1=%e, tm1=%e\n",
c,c-(double)n,cm1,tm1);
return 0;
} /* --- end-of-function main() --- */
编译并运行它作为测试重现OP的0x2000...场景，或运行它带有（最多三个）可选参数 test #trials multiplier timeexp 其中 #trials 默认为 OP 的 1000000，并且 multipler 默认为 2 用于 OP 的 2x16^13（将其更改为 4，等等，用于她的其他测试）。对于最后一个参数，timeexp，输入 0 只执行 expm1()（和我的不必要的泰勒式）计算。这样做的目的是表明 OP 显示的错误时间案例随着 expm1() 消失，无论 multiplier.
所以默认运行s，test和test 1000000 4，产生（好吧，我调用程序四舍五入)...
bash-4.3$ ./rounding
N=1000000, denom=2x16^13, Clock time: 11155070 (11.16 secs)
c=1.00000000000000023283e+06,
c-n=2.328306e-10, cm1=1.136017e-07, tm1=1.136017e-07
bash-4.3$ ./rounding 1000000 4
N=1000000, denom=4x16^13, Clock time: 200211 (0.20 secs)
c=1.00000000000000011642e+06,
c-n=1.164153e-10, cm1=5.680083e-08, tm1=5.680083e-08
所以你首先要注意的是，OP 的 c-n 使用 exp() 与 cm1==tm1 使用 expm1() 和我的泰勒大约。如果你减少N他们达成一致，如下...
N=10, denom=2x16^13, Clock time: 941 (0.00 secs)
c=1.00000000000007140954e+01,
c-n=7.140954e-13, cm1=7.127632e-13, tm1=7.127632e-13
bash-4.3$ ./rounding 100
N=100, denom=2x16^13, Clock time: 5506 (0.01 secs)
c=1.00000000000010103918e+02,
c-n=1.010392e-11, cm1=1.008393e-11, tm1=1.008393e-11
bash-4.3$ ./rounding 1000
N=1000, denom=2x16^13, Clock time: 44196 (0.04 secs)
c=1.00000000000011345946e+03,
c-n=1.134595e-10, cm1=1.140730e-10, tm1=1.140730e-10
bash-4.3$ ./rounding 10000
N=10000, denom=2x16^13, Clock time: 227215 (0.23 secs)
c=1.00000000000002328306e+04,
c-n=2.328306e-10, cm1=1.131288e-09, tm1=1.131288e-09
bash-4.3$ ./rounding 100000
N=100000, denom=2x16^13, Clock time: 1206348 (1.21 secs)
c=1.00000000000000232831e+05,
c-n=2.328306e-10, cm1=1.133611e-08, tm1=1.133611e-08
关于 exp() 与 expm1() 的时间关系，请自行查看...
bash-4.3$ ./rounding 1000000 2
N=1000000, denom=2x16^13, Clock time: 11168388 (11.17 secs)
c=1.00000000000000023283e+06,
c-n=2.328306e-10, cm1=1.136017e-07, tm1=1.136017e-07
bash-4.3$ ./rounding 1000000 2 0
N=1000000, denom=2x16^13, Clock time: 24064 (0.02 secs)
c=0.00000000000000000000e+00,
c-n=-1.000000e+06, cm1=1.136017e-07, tm1=1.136017e-07
问题：你会注意到一旦 exp() 计算达到 N=10000 次试验，无论较大 N。不确定为什么会这样。
>>__SECOND EDIT__<<
好的，@EOF，"you made me look" 接受你的 "heirarchical accumulation" 评论。这确实可以使 exp() 总和更接近（更接近）（可能是正确的）expm1() 总和。修改后的代码紧随其后进行讨论。但是这里有一个讨论记录：回忆上面的 multiplier。那不见了，在同一个地方是 expon 所以分母现在是 2^expon ，默认是 53，匹配 OP 的默认值（我相信更好地匹配她的想法）。好的，这是代码...
/*
i-do-not-want-correct-rounding-for-function-exp/44397261 */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#define BASE 2 /*denominator=2^EXPON, 2^53=2x16^13 default */
#define EXPON 53
#define taylorm1(a) (a*(1.+a*(1.+a/3.)/2.)) /*expm1() approx for small args*/
int main (int argc, char *argv[]) {
int N = (argc>1?atoi(argv[1]):1e6),
expon = (argc>2?atoi(argv[2]):EXPON),
isexp = (argc>3?atoi(argv[3]):1), /* flags to turn on/off exp() */
ncparts = (argc>4?atoi(argv[4]):1), /* #partial sums for c */
binsize = (argc>5?atoi(argv[5]):10);/* #doubles to sum in each bin */
int isexpm1 = 1; /* and expm1() for timing tests*/
int i, n=0;
double denom = pow((double)BASE,(double)expon);
double a, c=0.0, cm1=0.0, tm1=0.0;
double csums[10], cbins[10][65537]; /* c partial sums and heirarchy */
int nbins[10], ibin=0; /* start at lowest level */
clock_t start = clock();
n=0; c=cm1=tm1=0.0;
if ( ncparts > 65536 ) ncparts=65536; /* array size check */
if ( ncparts > 1 ) for(i=0;i<ncparts;i++) cbins[0][i]=0.0; /*init bin#0*/
/* --- to smooth random fluctuations, do the same type of computation
a large number of (N) times with different values --- */
for (i=0; i<N; i++) {
n++;
a = (double)(1 + 2*(rand()%0x400)) / denom; /* "a" has only a few
significant digits, and its last non-zero
digit is at (fixed-point) position 53. */
if ( isexp ) { /* turn this off to time expm1() alone */
double expa = exp(a); /* exp(a) */
c += expa; /* just accumulate in a single "bin" */
if ( ncparts > 1 ) cbins[0][n%ncparts] += expa; } /* accum in ncparts */
if ( isexpm1 ) { /* you can turn this off to time exp() alone, */
cm1 += expm1(a); /* but difference is negligible */
tm1 += taylorm1(a); }
} /* --- end-of-for(i) --- */
int nticks = (int)(clock()-start);
if ( ncparts > 1 ) { /* need to sum the partial-sum bins */
nbins[ibin=0] = ncparts; /* lowest-level has everything */
while ( nbins[ibin] > binsize ) { /* need another heirarchy level */
if ( ibin >= 9 ) break; /* no more bins */
ibin++; /* next available heirarchy bin level */
nbins[ibin] = (nbins[ibin-1]+(binsize-1))/binsize; /*#bins this level*/
for(i=0;i<nbins[ibin];i++) cbins[ibin][i]=0.0; /* init bins */
for(i=0;i<nbins[ibin-1];i++) {
cbins[ibin][(i+1)%nbins[ibin]] += cbins[ibin-1][i]; /*accum in nbins*/
csums[ibin-1] += cbins[ibin-1][i]; } /* accumulate in "one bin" */
} /* --- end-of-while(nprevbins>binsize) --- */
for(i=0;i<nbins[ibin];i++) csums[ibin] += cbins[ibin][i]; /*highest level*/
} /* --- end-of-if(ncparts>1) --- */
printf ("N=%d, denom=%d^%d, Clock time: %d (%.2f secs)\n", n, BASE,expon,
nticks, ((double)nticks)/((double)CLOCKS_PER_SEC));
printf ("\t c=%.20e,\n\t c-n=%e, cm1=%e, tm1=%e\n",
c,c-(double)n,cm1,tm1);
if ( ncparts > 1 ) { printf("\t binsize=%d...\n",binsize);
for (i=0;i<=ibin;i++) /* display heirarchy */
printf("\t level#%d: #bins=%5d, c-n=%e\n",
i,nbins[i],csums[i]-(double)n); }
return 0;
} /* --- end-of-function main() --- */
好的，现在您可以注意到旧的 timeexp 后面有两个额外的命令行参数。它们是 ncparts，表示整个 #trials 将分配到的初始垃圾箱数量。因此，在层次结构的最低级别，每个 bin 应该（模数 bugs:) 的总和为 #trials/ncparts 双倍。之后的参数是 binsize，这将是每个连续级别的每个 bin 中总和的双倍数，直到最后一级具有更少（或等于）#bins 为 binsize。所以这是一个将 1000000 次试验分成 50000 个箱子的示例，这意味着最低级别为 20doubles/bin，此后为 5doubles/bin...
bash-4.3$ ./rounding 1000000 53 1 50000 5
N=1000000, denom=2^53, Clock time: 11129803 (11.13 secs)
c=1.00000000000000465661e+06,
c-n=4.656613e-09, cm1=1.136017e-07, tm1=1.136017e-07
binsize=5...
level#0: #bins=50000, c-n=4.656613e-09
level#1: #bins=10002, c-n=1.734588e-08
level#2: #bins= 2002, c-n=7.974450e-08
level#3: #bins= 402, c-n=1.059379e-07
level#4: #bins= 82, c-n=1.133885e-07
level#5: #bins= 18, c-n=1.136214e-07
level#6: #bins= 5, c-n=1.138542e-07
请注意 c-n 对于 exp() 如何很好地收敛到 expm1() 价值。但请注意它在第 5 级时是如何最好的，并且根本没有统一收敛。请注意，如果您将 #trials 分成 5000 个初始分箱，您会得到同样好的结果，
bash-4.3$ ./rounding 1000000 53 1 5000 5
N=1000000, denom=2^53, Clock time: 11165924 (11.17 secs)
c=1.00000000000003527384e+06,
c-n=3.527384e-08, cm1=1.136017e-07, tm1=1.136017e-07
binsize=5...
level#0: #bins= 5000, c-n=3.527384e-08
level#1: #bins= 1002, c-n=1.164153e-07
level#2: #bins= 202, c-n=1.158332e-07
level#3: #bins= 42, c-n=1.136214e-07
level#4: #bins= 10, c-n=1.137378e-07
level#5: #bins= 4, c-n=1.136214e-07
其实用ncparts和binsize玩玩好像灵敏度不高，也不一定总是"more is better"（即 binsize 的更少）。所以我不确定到底发生了什么。可能是一个（或两个）错误，或者可能是@EOF 的另一个问题......???
>>编辑——显示对添加的示例 "binary tree" 层次结构<<
根据@EOF 的评论添加了以下示例 （注意：重新复制前面的代码。我必须将每个下一级的 nbins[ibin] 计算编辑为 nbins[ibin]=(nbins[ibin-1]+(binsize-1))/binsize； 来自 nbins[ibin]=(nbins[ibin-1]+2*binsize)/binsize; 这是 "too conservative" 创建 ...16,8,4,2 序列)
bash-4.3$ ./rounding 1024 53 1 512 2
N=1024, denom=2^53, Clock time: 36750 (0.04 secs)
c=1.02400000000011573320e+03,
c-n=1.157332e-10, cm1=1.164226e-10, tm1=1.164226e-10
binsize=2...
level#0: #bins= 512, c-n=1.159606e-10
level#1: #bins= 256, c-n=1.166427e-10
level#2: #bins= 128, c-n=1.166427e-10
level#3: #bins= 64, c-n=1.161879e-10
level#4: #bins= 32, c-n=1.166427e-10
level#5: #bins= 16, c-n=1.166427e-10
level#6: #bins= 8, c-n=1.166427e-10
level#7: #bins= 4, c-n=1.166427e-10
level#8: #bins= 2, c-n=1.164153e-10
>>编辑——在下面的评论中展示@EOF 的优雅解决方案<<
"Pair addition" 可以优雅地递归完成，根据下面@EOF 的评论，我在这里复制。 （注意递归结束时的情况 0/1 以处理 n even/odd。）
/* Quoting from EOF's comment...
What I (EOF) proposed is effectively a binary tree of additions:
a+b+c+d+e+f+g+h as ((a+b)+(c+d))+((e+f)+(g+h)).
Like this: Add adjacent pairs of elements, this produces
a new sequence of n/2 elements.
Recurse until only one element is left.
(Note that this will require n/2 elements of storage,
rather than a fixed number of bins like your implementation) */
double trecu(double *vals, double sum, int n) {
int midn = n/2;
switch (n) {
case 0: break;
case 1: sum += *vals; break;
default: sum = trecu(vals+midn, trecu(vals,sum,midn), n-midn); break; }
return(sum);
}
这是 EOF 之前关于他的 trecu() 算法和代码的评论的 "answer"/跟进 "binary tree summation" 建议。 "Prerequisites" 在阅读本文之前正在阅读该讨论。把所有这些都收集在一个有组织的地方会很好，但我还没有这样做...
...我所做的是将 EOF 的 trecu() 构建到我通过修改 OP 的原始测试程序编写的前面答案中的测试程序中。但后来我发现 trecu() 生成的答案与 "plain sum" c 使用 exp()，而不是使用 expm1() 的总和 cm1，我们希望从更准确的二进制文件中得到树求和。
但是那个测试程序有点（可能是两位:）"convoluted"（或者，如 EOF 所说，"unreadable"），所以我写了一个单独的较小的测试程序，如下所示（带有示例运行s 和下面的讨论），单独 test/exercise trecu()。此外，我还将函数 bintreesum() 写入下面的代码中，其中 abstracts/encapsulates 我在前面的测试程序中嵌入的二叉树求和的迭代代码。在前面的例子中，我的迭代代码确实接近 cm1 答案，这就是为什么我期望 EOF 的递归 trecu() 做同样的事情。总而言之，下面发生了同样的事情——bintreesum() 仍然接近正确答案，而 trecu() 离正确答案更远，准确地再现了 "plain sum".
我们下面求和的就是sum(i),i=1...n，也就是众所周知的n(n+1)/2。但这不太正确——为了重现 OP 的问题，被加数不是单独的 sum(i)，而是 sum(1+i*10^(-e))，其中 e 可以在命令行上给出。所以，比如说，n=5，你得到的不是 15，而是 5.000...00015，或者对于 n=6，你得到 6.000...00021，等等。为了避免很长很长的格式，我 printf( ) sum-n 删除整数部分。好的？？？所以这是代码...
/* Quoting from EOF's comment...
What I (EOF) proposed is effectively a binary tree of additions:
a+b+c+d+e+f+g+h as ((a+b)+(c+d))+((e+f)+(g+h)).
Like this: Add adjacent pairs of elements, this produces
a new sequence of n/2 elements.
Recurse until only one element is left. */
#include <stdio.h>
#include <stdlib.h>
double trecu(double *vals, double sum, int n) {
int midn = n/2;
switch (n) {
case 0: break;
case 1: sum += *vals; break;
default: sum = trecu(vals+midn, trecu(vals,sum,midn), n-midn); break; }
return(sum);
} /* --- end-of-function trecu() --- */
double bintreesum(double *vals, int n, int binsize) {
double binsum = 0.0;
int nbin0 = (n+(binsize-1))/binsize,
nbin1 = (nbin0+(binsize-1))/binsize,
nbins[2] = { nbin0, nbin1 };
double *vbins[2] = {
(double *)malloc(nbin0*sizeof(double)),
(double *)malloc(nbin1*sizeof(double)) },
*vbin0=vbins[0], *vbin1=vbins[1];
int ibin=0, i;
for ( i=0; i<nbin0; i++ ) vbin0[i] = 0.0;
for ( i=0; i<n; i++ ) vbin0[i%nbin0] += vals[i];
while ( nbins[ibin] > 1 ) {
int jbin = 1-ibin; /* other bin, 0<-->1 */
nbins[jbin] = (nbins[ibin]+(binsize-1))/binsize;
for ( i=0; i<nbins[jbin]; i++ ) vbins[jbin][i] = 0.0;
for ( i=0; i<nbins[ibin]; i++ )
vbins[jbin][i%nbins[jbin]] += vbins[ibin][i];
ibin = jbin; /* swap bins for next pass */
} /* --- end-of-while(nbins[ibin]>0) --- */
binsum = vbins[ibin][0];
free((void *)vbins[0]); free((void *)vbins[1]);
return ( binsum );
} /* --- end-of-function bintreesum() --- */
#if defined(TESTTRECU)
#include <math.h>
#define MAXN (2000000)
int main(int argc, char *argv[]) {
int N = (argc>1? atoi(argv[1]) : 1000000 ),
e = (argc>2? atoi(argv[2]) : -10 ),
binsize = (argc>3? atoi(argv[3]) : 2 );
double tens = pow(10.0,(double)e);
double *vals = (double *)malloc(sizeof(double)*MAXN),
sum = 0.0;
double trecu(), bintreesum();
int i;
if ( N > MAXN ) N=MAXN;
for ( i=0; i<N; i++ ) vals[i] = 1.0 + tens*(double)(i+1);
for ( i=0; i<N; i++ ) sum += vals[i];
printf(" N=%d, Sum_i=1^N {1.0 + i*%.1e} - N = %.8e,\n"
"\t plain_sum-N = %.8e,\n"
"\t trecu-N = %.8e,\n"
"\t bintreesum-N = %.8e \n",
N, tens, tens*((double)N)*((double)(N+1))/2.0,
sum-(double)N,
trecu(vals,0.0,N)-(double)N,
bintreesum(vals,N,binsize)-(double)N );
} /* --- end-of-function main() --- */
#endif
因此，如果您将其保存为 trecu.c，然后将其编译为 cc DTESTTRECU trecu.c lm o trecu 然后 运行 具有零到三个可选的命令行参数，如 trecu #trials e binsize 默认值是 #trials=1000000（如 OP 的程序）、e=10 和 binsize=2（对于我的bintreesum() 函数执行二叉树求和而不是更大尺寸的 bins）。
下面是一些说明上述问题的测试结果，
bash-4.3$ ./trecu
N=1000000, Sum_i=1^N {1.0 + i*1.0e-10} - N = 5.00000500e+01,
plain_sum-N = 5.00000500e+01,
trecu-N = 5.00000500e+01,
bintreesum-N = 5.00000500e+01
bash-4.3$ ./trecu 1000000 -15
N=1000000, Sum_i=1^N {1.0 + i*1.0e-15} - N = 5.00000500e-04,
plain_sum-N = 5.01087168e-04,
trecu-N = 5.01087168e-04,
bintreesum-N = 5.00000548e-04
bash-4.3$
bash-4.3$ ./trecu 1000000 -16
N=1000000, Sum_i=1^N {1.0 + i*1.0e-16} - N = 5.00000500e-05,
plain_sum-N = 6.67552231e-05,
trecu-N = 6.67552231e-05,
bintreesum-N = 5.00001479e-05
bash-4.3$
bash-4.3$ ./trecu 1000000 -17
N=1000000, Sum_i=1^N {1.0 + i*1.0e-17} - N = 5.00000500e-06,
plain_sum-N = 0.00000000e+00,
trecu-N = 0.00000000e+00,
bintreesum-N = 4.99992166e-06
因此您可以看到，对于默认值 运行，e=10，每个人都做得很好。也就是说，上面写着 "Sum" 的行只是做了 n(n+1)/2 的事情，所以大概会显示正确的答案。下面的每个人都同意默认的 e=10 测试用例。但是对于 e=15 和 e=16 下面的情况，trecu() 与 plain_sum 完全一致，而 bintreesum 非常接近正确答案。最后，对于 e=17，plain_sum 和 trecu() 有 "disappeared"，而 bintreesum() 仍然很好地挂在那里。
所以 trecu() 正确地进行了求和，但是它的递归显然没有做 "binary tree" 我更直接的迭代 bintreesum() 显然做正确的事情。这确实表明，对于这些 1+epsilon 类型的情况，EOF 对 "binary tree summation" 的建议比 plain_sum 实现了相当大的改进。所以我们真的很想看看他的 trecu() 递归工作！！！当我最初看它时，我认为它确实有效。但是在他的 default: 案例中，那个双重递归（它有一个特殊的名称吗？）显然比我想象的更令人困惑（至少对我而言：）。就像我说的，它 是 求和，而不是 "binary tree"。
好的，那么谁愿意接受挑战并解释 trecu() 递归中发生了什么？而且，也许更重要的是，对其进行修复，使其按预期运行。谢谢