bug in gsl_rstat_workspace

[Wilhelm Braunschober]

Yours Wilhelm Braunschober

/*!lcc
 * gcc -O3 %* -lgsl
 * a.exe
 * del a.exe
 */
/* gsl_rstat_median has 2 bugs:
 * 1. input vector size 5: result wrong unless sorted
 * 2. input vector size > 5 and even: always wrong
 * */
#include <stdio.h>
#include <gsl/gsl_rstat.h>

int main(void) {
    double d1[] { 1 };
    double d2[] { 2, 1 };
    double d3[] { 2, 3, 1 };
    double d4[] { 2, 4, 1, 3 };
    double d5[] { 3, 2, 5, 1, 4 };
    double d6[] { 3, 5, 1, 2, 4, 6 };
    double d7[] { 4, 5, 7, 1, 2, 6, 3 };
    double d8[] { 4, 6, 2, 1, 8, 3, 7, 5 };
    double* data[8] { d1, d2, d3, d4, d5, d6, d7, d8 };
    for (int i=0; i<8; ++i) {
        gsl_rstat_workspace *rstat_p = gsl_rstat_alloc();
        /* add data to rstat accumulator */
        for (int j = 0; j < i+1; ++j) {
            printf("%g ", data[i][j]);
            gsl_rstat_add(data[i][j], rstat_p);
        }
        double median = gsl_rstat_median(rstat_p);
        double should_be = (i+2.0)/2.0;
        printf(" median %g %s %g\n", median, median == should_be ? "==" : "!=", 
should_be);
    }
}