[Bug-gsl] [bug #38548] Rounding issues in gsl-histogram with integer numbers

[Christian Leitold]

URL:
  <http://savannah.gnu.org/bugs/?38548>

                 Summary: Rounding issues in gsl-histogram with integer
numbers
                 Project: GNU Scientific Library
            Submitted by: cleitold
            Submitted on: Di 19 Mär 2013 15:53:57 GMT
                Category: Accuracy problem
                Severity: 3 - Normal
        Operating System: Debian 7.0 amd64
                  Status: None
             Assigned to: None
             Open/Closed: Open
                 Release: 1.15
         Discussion Lock: Any

    _______________________________________________________

Details:

A rounding problem occurs in gsl-histogram 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. I have attached a sample
list of integer data where the problem occurs (varlist.dat). As discussed in
the corresponding thread [1] on the mailing list, a proposed solution is

double dx = (xmax-xmin) / (double) n;
double f1 = xmin + i*dx;
double f2 = xmax - (n-i) * dx;
range[i] = (f1 + f2) * 0.5;

As I understand it, in particular after reading [2], there are certain special
cases where floating-point arithmetic is just as accurate and error-free as
integer arithmetic. This should be fulfilled here, when xmin, xmax, n are real
integer values (even though represented as doubles), and thus dx (=1), f1 and
f2 are integer values as well.

[1] http://lists.gnu.org/archive/html/bug-gsl/2013-02/msg00006.html
[2] http://en.wikipedia.org/wiki/Double-precision_floating-point_format



    _______________________________________________________

File Attachments:


-------------------------------------------------------
Date: Di 19 Mär 2013 15:53:57 GMT  Name: varlist.dat  Size: 1kB   By:
cleitold
Sample data, use e. g. xmin = 0, xmax = n = 60
<http://savannah.gnu.org/bugs/download.php?file_id=27634>

    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?38548>

_______________________________________________
  Nachricht gesendet von/durch Savannah
  http://savannah.gnu.org/