[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
## [Bug-gsl] Rounding issues in histogram tools / gsl-histogram with intege

**From**: |
Christian Leitold |

**Subject**: |
[Bug-gsl] Rounding issues in histogram tools / gsl-histogram with integer numbers |

**Date**: |
Wed, 20 Feb 2013 16:38:24 +0100 |

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

[Prev in Thread] |
**Current Thread** |
[Next in Thread] |

**[Bug-gsl] Rounding issues in histogram tools / gsl-histogram with integer numbers**,
*Christian Leitold* **<=**