当给出相同的输入时，R的sum()函数和RcppArmadillo的accu()函数的结果存在微小差异。例如，以下代码：
R：
vec <- runif(100, 0, 0.00001)
accu(vec)
sum(vec)
C ++：
// [[Rcpp::depends("RcppArmadillo")]]
// [[Rcpp::export]]
double accu(arma::vec& obj)
{
return arma::accu(obj);
}
给出结果：
0.00047941851844312633（C ++）
0.00047941851844312628（R）
据http://keisan.casio.com/calculator说，真正的答案是：
4,794,185,444,316,270,948 J-4
这些小的差异在我的算法中加起来并且显着影响它的执行方式。有没有办法在C ++中更准确地总结向量？或者至少在不必调用R代码的情况下获得与R相同的结果？
更新：根据其他人在源代码中发现的内容，我错了 - sum()没有排序。我在下面找到的一致性模式源于以下事实：排序（如下面某些情况所做）和使用扩展精度中间值（如sum()中所做的那样）可以对精度产生类似的影响......
@ user2357112评论如下：
src/main/summary.c ...不进行任何排序。 （添加到求和运算需要很多费用。）它甚至不使用成对或补偿求和;它只是天真地在LDOUBLE中添加从左到右的所有内容（长双或双，取决于HAVE_LONG_DOUBLE）。
我已经筋疲力尽在R源代码中寻找这个（没有成功 - sum很难搜索），但是
我可以通过实验证明，在执行sum()时，R将输入向量从最小到最大排序，以便最大限度地提高准确性
; sum()和Reduce()结果之间的差异是由于使用了扩展精度。我不知道accu做了什么......
set.seed(101)
vec <- runif(100, 0, 0.00001)
options(digits=20)
(s1 <- sum(vec))
## [1] 0.00052502325481269514554
使用Reduce("+",...)只是按顺序添加元素。
(s2 <- Reduce("+",sort(vec)))
## [1] 0.00052502325481269514554
(s3 <- Reduce("+",vec))
## [1] 0.00052502325481269503712
identical(s1,s2) ## TRUE
?sum()也说
在可能的情况下使用扩展精度累加器，但这取决于平台。
在RcppArmadillo中对有序向量执行此操作会得到与R中相同的答案;在原始顺序中对矢量进行处理给出了一个不同的答案（我不知道为什么;我的猜测是前面提到的扩展精度累加器，当数据未排序时会更多地影响数值结果）。
suppressMessages(require(inline))
code <- '
arma::vec ax = Rcpp::as<arma::vec>(x);
return Rcpp::wrap(arma::accu(ax));
'
## create the compiled function
armasum <- cxxfunction(signature(x="numeric"),
code,plugin="RcppArmadillo")
(s4 <- armasum(vec))
## [1] 0.00052502325481269525396
(s5 <- armasum(sort(vec)))
## [1] 0.00052502325481269514554
identical(s1,s5) ## TRUE
但正如评论中所指出的，这并不适用于所有种子：在这种情况下，Reduce()结果更接近于sum()的结果
set.seed(123)
vec2 <- runif(50000,0,0.000001)
s4 <- sum(vec2); s5 <- Reduce("+",sort(vec2))
s6 <- Reduce("+",vec2); s7 <- armasum(sort(vec2))
rbind(s4,s5,s6,s7)
## [,1]
## s4 0.024869900535651481843
## s5 0.024869900535651658785
## s6 0.024869900535651523477
## s7 0.024869900535651343065
我在这里难过。我原以为至少s6和s7是相同的......
我将指出，一般情况下，当你的算法依赖于这些微小的数字差异时，你可能会非常沮丧，因为结果可能会因许多小而且可能不合适的情况而有所不同。控制因素，如您使用的特定操作系统，编译器等。
我发现了什么：
我成功地设法编写了一个能够模仿R和函数的函数。看起来R使用更高精度的变量来存储每个加法运算的结果。
我写的：
// [[Rcpp::depends("RcppArmadillo")]]
// [[Rcpp::export]]
double accu2(arma::vec& obj)
{
long double result = 0;
for (auto iter = obj.begin(); iter != obj.end(); ++iter)
{
result += *iter;
}
return result;
}
速度如何比较：
set.seed(123)
vec <- runif(50000, 0, 0.000001)
microbenchmark(
sum(vec),
accu(vec),
accu2(vec)
)
expr min lq mean median uq max neval
sum(vec) 72.155 72.351 72.61018 72.6755 72.7485 75.068 100
accu(vec) 48.275 48.545 48.84046 48.7675 48.9975 52.128 100
accu2(vec) 69.087 69.409 70.80095 69.6275 69.8275 182.955 100
所以，我的c ++解决方案仍然比R的总和更快，但它明显慢于犰狳的精度（）
你可以使用mpfr包（多精度浮点可靠）并指定小数点
library("Rmpfr")
set.seed(1)
vec <- runif(100, 0, 0.00001)
# [1] 2.655087e-06 3.721239e-06 5.728534e-06 9.082078e-06 2.016819e-06 8.983897e-06 9.446753e-06 6.607978e-06 6.291140e-06 6.178627e-07 2.059746e-06
# [12] 1.765568e-06 6.870228e-06 3.841037e-06 7.698414e-06 4.976992e-06 7.176185e-06 9.919061e-06 3.800352e-06 7.774452e-06 9.347052e-06 2.121425e-06
# [23] 6.516738e-06 1.255551e-06 2.672207e-06 3.861141e-06 1.339033e-07 3.823880e-06 8.696908e-06 3.403490e-06 4.820801e-06 5.995658e-06 4.935413e-06
# [34] 1.862176e-06 8.273733e-06 6.684667e-06 7.942399e-06 1.079436e-06 7.237109e-06 4.112744e-06 8.209463e-06 6.470602e-06 7.829328e-06 5.530363e-06
# [45] 5.297196e-06 7.893562e-06 2.333120e-07 4.772301e-06 7.323137e-06 6.927316e-06 4.776196e-06 8.612095e-06 4.380971e-06 2.447973e-06 7.067905e-07
# [56] 9.946616e-07 3.162717e-06 5.186343e-06 6.620051e-06 4.068302e-06 9.128759e-06 2.936034e-06 4.590657e-06 3.323947e-06 6.508705e-06 2.580168e-06
# [67] 4.785452e-06 7.663107e-06 8.424691e-07 8.753213e-06 3.390729e-06 8.394404e-06 3.466835e-06 3.337749e-06 4.763512e-06 8.921983e-06 8.643395e-06
# [78] 3.899895e-06 7.773207e-06 9.606180e-06 4.346595e-06 7.125147e-06 3.999944e-06 3.253522e-06 7.570871e-06 2.026923e-06 7.111212e-06 1.216919e-06
# [89] 2.454885e-06 1.433044e-06 2.396294e-06 5.893438e-07 6.422883e-06 8.762692e-06 7.789147e-06 7.973088e-06 4.552745e-06 4.100841e-06 8.108702e-06
# [100] 6.049333e-06
sum(mpfr(vec,10))
# 1 'mpfr' number of precision 53 bits
# [1] 0.00051783234812319279