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

[Christian Leitold]

Hello,

On 21 February 2013 21:14, Alexey A. Illarionov
<address@hidden> 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 <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-----

varlist.dat
Description: Binary data