[Bug-gsl] Rounding issues in histogram tools / gsl-histogram with integer numbers

[Christian Leitold]

Hello,

I have recently discovered a bug (or at least what is a bug for me) in
gsl-histogram respectively one of the helper function it calls. It
seems to occur in all versions, in any case in 1.15. The problem
occurs when you try to make a histogram from pure integer data, where
it is appropriate to choose a bin size of 1. I have tracked this down
to the file histogram/init.c, where the array range is filled with the
bin limits. There, in make_uniform one has

double f1 = ((double) (n-i) / (double) n);
double f2 = ((double) i / (double) n);
range[i] = f1 * xmin +  f2 * xmax;

In the case of e. g. xmin = 0.0, xmax = 60.0, n = 60, this (on my
system, gcc 4.6, standard make command from the file archive) leads to

range[31] = 31.000000000000003552713679

where it should be, without rounding error,

range[31] = 31.000000000000000000000000

This then leads to the number 31 being put into bin number 30, instead
of 31. A bunch of uniformely distributed integer data would then lead
to a peak of double the average height at 30, and no entry at 31. So
my solution would be

double dx = (xmax-xmin) / (double) n;
range[i] = xmin + i*dx;

which is apart from rounding error equivalent and as I understand it
does not not involve any rounding error at all as long as xmin, xmax
are true integer values (represented as doubles) and furthermore, n =
xmax - xmin and thus dx = 1.0, which again can be represented as
double without any rounding error. Is there any (presumably numerical)
reason I did overlook why one would still prefer the original
solution?

 - Christian Leitold