Как получить точное среднее значение для большого массива с плавающей запятой?

Question

Как получить точное среднее значение для большого массива с плавающей запятой (более 100 000 значений)? В идеале с использованием инструкций SIMD / AVX. Указатель на массив в rdi; размер массива в RSI.

Answer 1

точно

Если точность важнее скорости:

Используя арифметику с плавающей запятой c, вы, вероятно, всегда будете иметь потерю точности.

Однако вы можете рассчитать точное значение, если используете арифметические значения с фиксированной запятой c:

Все значения с плавающей запятой могут быть выражены как произведение некоторой константы ( что типично для используемого типа данных) и большое целое число со знаком.

В случае double каждое значение может быть выражено как произведение константы, типичной для типа данных double, и 2102-битное целое число со знаком.

Если в вашем массиве 10 миллионов элементов, сумма всех элементов может быть выражена как произведение этой константы на 2126-битное целое число со знаком. (Поскольку 10 миллионов умещаются в 24 бита и 2102 + 24 = 2026.)

Вы можете использовать те же методы, которые используются для 32-битной целочисленной арифметики c на 8-битном процессоре для выполнения 2126 -разрядная целочисленная арифметика c на 64-битном ЦП.

Вместо того, чтобы складывать все значения с плавающей запятой, вы складываете 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;
}

Если сумма отрицательна, отмените ее (вычислите абсолютное значение среднего); в этом случае вы должны позже отрицать результат (среднее значение).

Выполните целочисленное деление суммы (разделите ее на количество элементов).

Теперь вычислите (абсолютное значение of) среднее из полученного большого целочисленного значения:

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 2

Чтобы минимизировать потерю точности, я использую массив из 2048 чисел двойной точности, индексированных показателем степени, что означает, что код является спецификацией реализации c и ожидает, что двойники будут двойниками в формате 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 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;
}