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

precisely

If precision is more important than speed:

Using floating-point arithmetic you will probably always have a loss of precision.

However, you can calculate the exact value if you use fixed-point arithmetic:

All floating-point values can be expressed as the product of some constant (which is typical for the data type used) and a large signed integer value.

In the case of double, each value can be expressed as product of a constant typical for the double data type and a 2102-bit signed integer.

If your array has 10 million elements, the sum of all elements can be expressed as product of that constant multiplied with a 2126-bit signed integer. (Because 10 million fits into 24 bits and 2102 + 24 = 2026.)

You can use the same methods that are used to do 32-bit integer arithmetic on an 8-bit CPU to perform 2126-bit integer arithmetic on a 64-bit CPU.

Instead of adding up all floating-point values itself, you add up the 2102-bit integers representing each floating point value (here lsint is a signed data type that can handle 2126-bit integers):

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

If the sum is negative, negate it (calculate the absolute value of the average); in this case, you have to negate the result (average) later.

Perform an integer division of the sum (divide it by the number of elements).

Now calculate the (absolute value of) the average from the resulting large integer value:

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

If the sum was negative and you calculate the absolute value, don't forget to negate the result.

To minimize loss in precision, I use an array of 2048 doubles, indexed by the exponent, which means the code is implementation specific and expects doubles to be IEEE formatted doubles. Numbers are added into the array, only adding numbers with identical exponents. To get the actual sum, the array is then added from smallest to largest.

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

Given OP's:

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

A middle of the road approach for increased precision when the values have the same sign and within a few magnitudes of each other:

2 passes, find coarse average and then find average deviation from the average.

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