Double 的 R 精度 - 为什么代码返回负值，为什么预期正结果？

问：

我正在测试计算 Prod（b-a） 的 2 种方法，其中 a 和 b 是长度为 n 的向量。 Prod（b-a）=（b1-a1）（b2-a2）（b3-a3）*...（bn-an），其中 b_i>a_i>0 表示所有 i=1,2,3，n。对于某些特殊情况，另一种计算方法（方法2）此prod（b-a）更有效。它使用以下公式，即展开项并求和：

我的问题是：当a_i非常接近 b_i 时，真正的结果可能是非常非常接近 0，比如 10^（-16）。方法 1（减法和乘法）始终返回正输出。使用公式的方法 2 有时会返回负输出（我的实验大约有 7~8% 的时间返回负值）。从数学上讲，这两个方法应该返回完全相同的输出。但在计算机语言中，它显然会产生不同的输出。

以下是我运行测试的代码。当我运行测试代码 10000 次时，大约 7~8% 的方法 2 运行返回负输出。根据官方文档，R double 的精度为 “2.225074e-308”，如 R 参数所示：“.机器$double.xmin”。为什么当差异在 10^（-16） ~ 10^（-18） 之间时它会变成负值？任何能说明这一点的帮助都将被逮捕。我还希望得到一些关于如何实际将精度提高到更高级别的建议，如R文档所示。

##########  Testing code 1.  
ftest1case<-function(a,b)  {
  n<-length(a)
  if (length(b)!=n) stop("---------  length a and b are not right.")
  if ( any(b<a) ) stop("----------  b has to be greater than a all the time.")
  out1<-prod(b-a)
  
  out2<-0
  N<-2^n
  for ( i in 1:N ) {
    tidx<-rev(as.integer(intToBits(x=i-1))[1:n])
    tsign<-ifelse( (sum(tidx)%%2)==0,1.0,-1.0)
    out2<-out2+tsign*prod(b[tidx==0])*prod(a[tidx==1])
  }
  c(out1,out2)
}


##########  Testing code 2.  
ftestManyCases<-function(N,printFreq=1000,smallNum=10^(-20))
{
  tt<-matrix(0,nrow=N,ncol=2)
  n<-12
  for ( i in 1:N) {
    a<-runif(n,0,1)
    b<-a+runif(n,0,1)*0.1
    tt[i,]<-ftest1case(a=a,b=b)
  
    if ( (i%%printFreq)==0 ) cat("----- i = ",i,"\n")
    if ( tt[i,2]< smallNum ) cat("------ i = ",i, " ---- Negative summation found.\n")
  }
  
  tout<-apply(tt,2,FUN=function(x)  { round(sum(x<smallNum)/N,6) } )
  names(tout)<-c("PerLess0_Method1","PerLee0_Method2")
  list(summary=tout, data=tt)  
}  
 

########  Step 1.  Test for 1 case.
n<-12
a<-runif(n,0,1)
b<-a+runif(n,0,1)*0.1
ftest1case(a=a,b=b)   

########  Step 2   Test Code 2 for multiple cases.
N<-300
tt<-ftestManyCases(N=N,printFreq = 100) 
tt[[1]]