如何精确取大型浮点数组的平均值？答案

【问题标题】：How do I take the average of a large floating point array precisely?如何精确取大型浮点数组的平均值？
【发布时间】：2020-11-03 07:34:43
【问题描述】：

如何精确取大型浮点数组（100.000+ 个值）的平均值？理想情况下使用 SIMD/AVX 指令。 rdi 中的数组指针； rsi 中数组的大小。

【问题讨论】：

现有代码..？什么是基本的和+除“不精确”？这些值是否会丢失（可能）相关信息？
目前我正在尝试找出更精确的方法 - 迭代加法或 Kahan 总和。两者似乎都是可行的方法。预计我使用的值不会处于任何极端，但我对这些数字没有“感觉”。将采用的算法将估计一个用于计算的“幻数”。
请注意，使用多个 SIMD 向量累加器展开的普通求和是朝着成对求和方向迈出的一步。我认为即使以速度为代价，您也希望最大限度地提高准确性，但也许其他未来的读者会考虑这一点（因为它也是对数组求和的最快方法，并且通常比天真的标量更准确）。 Simd matmul program gives different numerical results。 Efficient stable sum of a sorted array in AVX2 也有一些指向 SIMD Kahan 和其他求和算法的链接。
@EricPostpischil 我没有写它是“几乎不可能”，而是我在阅读后的第一刻认为题。 “几乎不可能”意味着一开始我怀疑是否存在在所有情况下都能提供正确结果的算法：在示例中，当我们按原始顺序汇总元素时，我们有 100% 的错误。示例{ +1E+100, 16, -1E+100, -12} 甚至会导致 300% 的错误。而且 300% 的错误绝对不是“精确的”。
为了澄清我所说的精确的意思：如果可以以无限精度添加所有浮点数，然后除以然后再次将结果减少为 32 位浮点数，我们将得到我定义为“精确" 所有浮点数的平均值。

标签： assembly floating-point precision simd avx

【解决方案1】：

准确的

如果精度比速度更重要：

使用浮点运算可能总是会损失精度。

但是，如果您使用定点算法，您可以计算出准确的值：

所有浮点值都可以表示为某个常数（这是所使用的数据类型的典型值）和一个大的有符号整数值的乘积。

在double 的情况下，每个值都可以表示为double 数据类型的典型常数和2102 位有符号整数的乘积。

如果您的数组有 1000 万个元素，则所有元素的总和可以表示为该常数乘以 2126 位有符号整数的乘积。（因为 1000 万适合 24 位和 2102 + 24 = 2026。）

您可以使用在 8 位 CPU 上执行 32 位整数运算的相同方法在 64 位 CPU 上执行 2126 位整数运算。

您无需将所有浮点值本身相加，而是将表示每个浮点值的 2102 位整数相加（这里 lsint 是一种可以处理 2126 位整数的有符号数据类型）：

void addNumber(lsint * sum, double d)
{
    uint64   di = *(uint64 *)&d;
    lsint    tmp;
    int      ex = (di>>52)&0x7FF;
    if(ex == 0x7FF)
    {
        /* Error: NaN or Inf found! */
    }
    else if(ex == 0)
    {
        /* Denormalized */
        tmp = di & 0xFFFFFFFFFFFFF;
    }
    else
    {
        /* Non-Denormalized */
        tmp = di & 0xFFFFFFFFFFFFF;
        tmp |= 0x10000000000000;
        tmp <<= ex-1;
    }
    if(di & 0x8000000000000000) (*sum) -= tmp;
    else (*sum) += tmp;
}

如果总和为负，则取反（计算平均值的绝对值）；在这种情况下，您必须稍后否定结果（平均值）。

对总和进行整数除法（除以元素的数量）。

现在计算得到的大整数值的（绝对值）平均值：

double lsintToDouble(lsint sum)
{
    int    ex;
    double result;
    if(sum < 0x10000000000000)
    {
        *(uint64 *)&result = (uint64)sum;
    }
    else
    {
        ex = 1;
        while(sum >= 0x20000000000000)
        {
            sum >>= 1;
            ex++;
        }
        *(uint64 *)&result = (uint64)sum & 0xFFFFFFFFFFFFF;
        *(uint64 *)&result |= ex<<52;
    }
     return result;
}

如果总和是负数并且你计算了绝对值，不要忘记否定结果。

【讨论】：

Re“2102位有符号整数”：最高位为2^1024，最低位为2^(−1022-52) = 2^−1074，所以需要2099位值位，符号位多为 2100。尽管仍然必须考虑无穷大和 NaN。
@EricPostpischil 感谢您的澄清；但是，因为 2102>2100，它至少可以使用 2102 位...
1) 让我们称之为好并使用int4096_t。 ? 2) 进行精确加法的大量工作，但从 lsint 到 double 的转换会截断而不是丢失 1/2 ULP - 假设是 LAAETTR。

【解决方案2】：

为了最大程度地减少精度损失，我使用了一个由 2048 个双精度数组成的数组，由指数索引，这意味着代码是特定于实现的，并且期望双精度数是 IEEE 格式的双精度数。将数字添加到数组中，仅添加具有相同指数的数字。为了得到实际的总和，然后将数组从小到大相加。

/* clear array */
void clearsum(double asum[2048])
{
size_t i;
    for(i = 0; i < 2048; i++)
        asum[i] = 0.;
}

/* add a number into array */
void addtosum(double d, double asum[2048])
{
size_t i;
    while(1){
        /* i = exponent of d */
        i = ((size_t)((*(unsigned long long *)&d)>>52))&0x7ff;
        if(i == 0x7fe){         /* max exponent, could be overflow */
            asum[i] += d;
            return;
        }
        if(asum[i] == 0.){      /* if empty slot store d */
            asum[i] = d;
            return;
        }
        d += asum[i];           /* else add slot to d, clear slot */
        asum[i] = 0.;           /* and continue until empty slot */
    }
}

/* return sum from array */
double returnsum(double asum[2048])
{
double sum = 0.;
size_t i;
    for(i = 0; i < 2048; i++)
        sum += asum[i];
    return sum;
}

【讨论】：

好主意，虽然远非完美。例如，添加(1+eps) + (1+2eps) 将产生舍入误差（如果您稍后减去相似的值，可能会变得很重要）。一种更“直接”的方法是使用(1024)+(1077) 位定点数进行累加（可能在顶部添加一些额外位以避免溢出）。
@chtz - 我不明白你所说的 (1+eps) + (1 + 2eps) 是什么意思。
eps指的是机器精度（最小的数字，如@987654327@），1+eps和1+2*eps应该是单独的输入数字，它们具有相同的指数，并且它们的sum 将添加到2+3*eps，它被四舍五入为2+4*eps（四舍五入）。如果你稍后还添加-2.0，你的结果是4*eps，而它应该是3*eps。
如果我的观点不明确，这里有一个简单的失败示例（naive_sum 恰好可以工作，因为我以某种方式人为地对值进行排序）：godbolt.org/z/rx79x7
@chtz：“最小的数，例如1+eps > 1”是“机器精度”的错误定义，尽管它被广泛使用。对于 IEEE-754，正确的机器精度是 2^-52，但该定义给出了 2^-53+2^-105，因为 1+(0x1p-53+0x1p-105) 由于四舍五入将产生 1。正确的定义是 1 与可表示的下一个最大数之间的差。

【解决方案3】：

给定 OP：

预计我使用的值不会处于任何极端，但我对数字没有“感觉”

当值具有相同符号且彼此相差几个数量级时，采用中间方法来提高精度：

2遍，求粗略平均值，然后求平均值与平均值的偏差。

double average(size_t rsi, const double *rdi) {
   double sum = 0.0;
   for (size_t i=0; i<rsi; i++) {
     sum += rdi[i];
   }
   double course_average = sum/rsi;

   sum = 0.0;
   for (size_t i=0; i<rsi; i++) {
     sum += rdi[i] - course_average;
   }
   double differnce_average = sum/rsi;

   return course_average + differnce_average;
}

【讨论】：