Based on discussion on the newsgroup today about the errors in summing large arrays, I implemented the Kahan summation algorithm as a routine in a DLM. Download, unzip, and build it with the following command:

IDL> mg_make_dll, 'mg_analysis.c', /show_all_output

Don’t forget to put the directory with the DLM in your IDL_DLM_PATH.

I thought it would be interesting because of the discussion in the newsgroup and the fact that it was pretty simple to make it work for all the different numeric data types IDL offers (for an example of when it is not so simple, check out the implementation of MG_ARRAY_EQUAL in the source code).

The main algorithm is very simple, but complicated a bit by the fact that it is a C macro:

#define IDL_MG_TOTAL_TYPE(TYPE, IDL_TYPE, INIT)        \
IDL_VPTR IDL_mg_total_ ## IDL_TYPE (IDL_VPTR arg) {    \
  TYPE sum = INIT, c = INIT, y = INIT, t = INIT;       \
  int i;                                               \
                                                       \
  TYPE *arr = (TYPE *)arg->value.arr->data;            \
                                                       \
  for (i = 0; i < arg->value.arr->n_elts; i++) {       \
    y = arr[i] - c;                                    \
    t = sum + y;                                       \
    c = (t - sum) - y;                                 \
    sum = t;                                           \
  }                                                    \
  return(IDL_Gettmp ## IDL_TYPE (sum));                \
}

Ugly yes, but better than repeating the code nine times.

So this macro needs to be called for each of the supported data types to make the routines that handle each type:

IDL_MG_TOTAL_TYPE(UCHAR, Byte, 0)
IDL_MG_TOTAL_TYPE(short int, Int, 0)
IDL_MG_TOTAL_TYPE(int, Long, 0)
IDL_MG_TOTAL_TYPE(float, Float, 0.0)
IDL_MG_TOTAL_TYPE(double, Double, 0.0)
IDL_MG_TOTAL_TYPE(IDL_UINT, UInt, 0)
IDL_MG_TOTAL_TYPE(IDL_ULONG, ULong, 0)
IDL_MG_TOTAL_TYPE(IDL_LONG64, Long64, 0)
IDL_MG_TOTAL_TYPE(IDL_ULONG64, ULong64, 0)

And then the master routine needs to call the appropriate routine for the input type (it always computes a result that is the same type as the input):

IDL_VPTR IDL_mg_total(int argc, IDL_VPTR *argv, char *argk) {
  IDL_ENSURE_ARRAY(argv[0]);
  switch (argv[0]->type) {
    case IDL_TYP_BYTE:
      return(IDL_mg_total_Byte(argv[0]));
    case IDL_TYP_INT:
      return(IDL_mg_total_Int(argv[0]));
    case IDL_TYP_LONG:
      return(IDL_mg_total_Long(argv[0]));
    case IDL_TYP_FLOAT:
      return(IDL_mg_total_Float(argv[0]));
    case IDL_TYP_DOUBLE:
      return(IDL_mg_total_Double(argv[0]));
    case IDL_TYP_UINT:
      return(IDL_mg_total_UInt(argv[0]));
    case IDL_TYP_ULONG:
      return(IDL_mg_total_ULong(argv[0]));
    case IDL_TYP_LONG64:
      return(IDL_mg_total_Long64(argv[0]));
    case IDL_TYP_ULONG64:
      return(IDL_mg_total_ULong64(argv[0]));
    default:
      IDL_Message(IDL_M_NAMED_GENERIC, IDL_MSG_LONGJMP, "unknown type");    
  }
}