Hello, On 21 February 2013 21:14, Alexey A. Illarionov <[email protected]> wrote:
> 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; I have tried this out and in my testcase it performs in the case of the integer numbers just as well as my solution. Also, I have attached my list of integer used to compare the different versions. > Dealing with doubles is always tricky. The rounding error is > inevitable in floating point computations. True, and you are also right when you say that one should actually treat the integer case completely separately. However, as I understand it, in particular after reading [1], there are certain special cases where floating-point arithmetic is just as accurate and error-free as integer arithmetic. My hope is that we could stay within the limits of these special cases in our problem and thus avoid the extra integer-only routines and still achieving maximal precision. - Christian [1] http://en.wikipedia.org/wiki/Double-precision_floating-point_format > - -------- 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-----
varlist.dat
Description: Binary data
