-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1 Sorry, I forgot to cc to mailing list.
Also I think you could try to avoid linear rounding error with something like this. double dx = (xmax-xmin) / (double) n; double f1 = xmin + i*dx; double f2 = xmax - (n-i) * dx; range[i] = (f1 + f2) * 0.5; But, again there are no good solutions, especially if we take into account the compiler optimizers. - -------- Original Message -------- Subject: Re: [Bug-gsl] Rounding issues in histogram tools / gsl-histogram with integer numbers Date: Thu, 21 Feb 2013 14:58:22 -0500 From: Alexey A. Illarionov <[email protected]> To: Christian Leitold <[email protected]> Hi, Dealing with doubles is always tricky. The rounding error is inevitable in floating point computations. Your code is also non-universal. I believe that the original code is written to have homogeneous rounding error from both ends. As a downside this error is a little big bigger (31.000000000000003552713679 - 31.)/ 31. = 2.220446049250313e-16 which is almost DBL_EPSILON, The only 'correct' solution would be to write histogram routines for integers separately. P.S.: Could you write the failing test case? Christian Leitold wrote: > 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 > -----BEGIN PGP SIGNATURE----- Version: GnuPG v2.0.19 (GNU/Linux) iQEcBAEBAgAGBQJRJoA5AAoJEEBWYSFzoNKerDAIAInuzs4r1r/uXWAHdnsac28G 4MEGaeSikl0/PWZKQLsFpwtYPYHf3XV7P17j+Uzli7hcSuO2T/6ukcUM5dmPGo2/ UrZiYbBmRDkPCcBzywFd8lzH8ezQvx1iq8T3TxDVpLFnLoWrjqauptWb5V45XWZj EERkRTnBbITMN8MhXec5na3ATKM+AqJO0nZW6VyuBX2nIDtlCBKw0T50oaAxbV6Y YWiOzlfT0uFLj4rKXX3Y3LaUf2/iGqPA/E9TEH8UsEaLr6XxmEXLBTB0BWsImFyA QN0te45HJ7OfUlNGYZUM30ngY4TPFwaHrNReTshfaq69fNBP3oBG6UryiXgfD4o= =w2QG -----END PGP SIGNATURE-----
