Floating-point Summation

Listing 6: Pairwise Summation

float fp_add(float * flt_arr)
{
  long i, j, limit;
  float sum[ARR_SIZE / 2 + 1];

  if (ARR_SIZE == 2)
    return flt_arr[0] + flt_arr[1];
  else if (ARR_SIZE == 1)
    return flt_arr[0];

  for (i = j = 0; i < ARR_SIZE / 2; i++, j += 2)
    sum[i] = flt_arr[j] + flt_arr[j + 1];
  if (ARR_SIZE & 1)
    sum[ARR_SIZE / 2] = flt_arr[ARR_SIZE - 1];
  limit = (ARR_SIZE + 1) / 2;

  while (limit > 2) {
    for (i = j = 0; i < limit / 2; i++, j += 2)
      sum[i] = sum[j] + sum[j + 1];

    if (limit & 1)
      sum[limit / 2] = sum[limit - 1];
    limit = (limit + 1) / 2;
  }

  return sum[0] + sum[1];
}
/* End of File */