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");
}
}
```

July 9th, 2011 at 4:29 am

Mike,

I think you’ve linked to the mgunit-src… file instead of the DLM

Mort

July 9th, 2011 at 7:40 am

Sorry, the link should be fixed now!

October 25th, 2011 at 2:20 pm

Kahan summation is for floating point. No point having that for integer types.

October 25th, 2011 at 3:09 pm

Once you have the mechanism for float and double, it’s easy to get the others. Implementing all the types let’s you just use it without having to worry what type of data you have.

I should implement complex and double complex, though.