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

[Alexey A. Illarionov]

-----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 <address@hidden>
To: Christian Leitold <address@hidden>

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-----