[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Bug-gsl] Fwd: Re: Rounding issues in histogram tools / gsl-histogram wi

From: Alexey A. Illarionov
Subject: [Bug-gsl] Fwd: Re: Rounding issues in histogram tools / gsl-histogram with integer numbers
Date: Thu, 21 Feb 2013 15:14:49 -0500
User-agent: Mozilla/5.0 (X11; Linux x86_64; rv:18.0) Gecko/20100101 Firefox/18.0 SeaMonkey/2.15.1

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>


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

Version: GnuPG v2.0.19 (GNU/Linux)


reply via email to

[Prev in Thread] Current Thread [Next in Thread]