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

Question

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

Answer 1

为了尽量减少精度损失，我使用了一个由 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;
}

Answer 2

precisely

如果精度比速度更重要：

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

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

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

在 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;
}

如果总和为负，而您计算的是绝对值，请不要忘记对结果取反。

Answer 3

给定 OP：

The values I work with are not expected to be on any extreme side, but I do not have a "feel" for the numbers

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

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;
}

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

How do I take the average of a large floating point array precisely?

floating-point

precision

assembly

simd

avx