bugGNU Scientific Library - Bugs: bug #22293, integer overflow in...

 
 

You are not allowed to post comments on this tracker with your current authentication level.

bug #22293: integer overflow in gsl_cdf_hypergeometric_P/gsl_cdf_hypergeometric_Q

Submitted by:  -Deleted Account- <bjg>
Submitted on:  Tue 12 Feb 2008 01:35:31 PM UTC  
 
Category: Runtime errorSeverity: 4 - Important
Operating System: Status: Fixed
Assigned to: -Deleted Account- <bjg>Open/Closed: Closed
Release: 1.10

Tue 12 Feb 2008 01:35:31 PM UTC, original submission:

Reply-To: Stijn van Dongen <svd@sanger.ac.uk>
From: Stijn van Dongen <svd@sanger.ac.uk>
To: -unavailable-
Subject: [Bug-gsl] bugs in gsl_cdf_hypergeometric_P/gsl_cdf_hypergeometric_Q
Date: Sat, 9 Feb 2008 13:57:55 +0000

[1 <text/plain; us-ascii (7bit)>]
Hi,

I believe there is a bug in the functions gsl_cdf_hypergeometric_P and
gsl_cdf_hypergeometric_Q in gsl-1.10 when computing the midpoint. This can
result in severely wrong probabilities (as a result of computing the
non-optimal tail), including negative probabilities. A demonstration program
is included below.

Explanation: The relevant statement is

double midpoint = (int) (t * n1 / (n1 + n2));

The integer multiplication t * n1 may overflow. With certain
input parameters it does; this causes the wrong tail to be computed,
leading to accumulation of rounding errors, and erroneous results
including negative probabilities. Running the example program yields

midpoint should be: 3899.52, but it is: 9
midpoint should be: 3899.52, but it is: 9
hyper [ > 3511 1.000000006 ] [ <= 3511 -6.46390208e-09 ]

The midpoint is then compared with the variable k which has value 3511,
leading to a choice of the wrong tail to compute.
(the first two lines are from a patched version of cdf/hypergeometric.c that
prints out the midpoint). I believe that the correct computation of midpoint
should be:

double midpoint = (1.0 * t * n1) / (n1 + n2);

A (trivial) diff to this effect is attached.

On a side note, is there any knowledge/documentation available on the range
of parameters for which these functions are supposed to give good/reasonable
precision?

regards,
Stijn

Example program (run on a system where integers are 4 bytes):

----------------------------------------
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <gsl/gsl_cdf.h>

int main
( int argc
, char* argv[]
)
{ int k = 3511
; int n1 = 76200
; int bg = 13250474
; int t = 678090
; double f = gsl_cdf_hypergeometric_Q(k, n1, bg-n1, t)
; double g = gsl_cdf_hypergeometric_P(k, n1, bg-n1, t)
; fprintf(stdout, "hyper [ > %d %.10g ] [ <= %d %.10g ]\n", k, f, k, g)
; return 0
; }
----------------------------------------

--
Stijn van Dongen >8< -o) O< forename pronunciation: [Stan]
Wellcome Trust Sanger Institute, /\\ Tel: +44-(0)1223-494785
Hinxton, Cambridge, CB10 1SA, UK \/ http://www.sanger.ac.uk/Users/svd/

--
The Wellcome Trust Sanger Institute is operated by Genome Research
Limited, a charity registered in England with number 1021457 and a
company registered in England with number 2742969, whose registered
office is 215 Euston Road, London, NW1 2BE.
[2 hypergeometric.c.diff <text/plain; us-ascii (7bit)>]
131c131
< double midpoint = (int) (t * n1 / (n1 + n2));
---

> double midpoint = (1.0 * t * n1) / (n1 + n2);

148a149

>

170c171
< double midpoint = (int) (t * n1 / (n1 + n2));
---

> double midpoint = (1.0 * t * n1) / (n1 + n2);

[3 <text/plain; us-ascii (7bit)>]
_____________________________________________
Bug-gsl mailing list
-unavailable-
http://lists.gnu.org/mailman/listinfo/bug-gsl

-Deleted Account- <bjg>In charge of this item.

 

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

CC list is empty

 

Please enter the title of George Orwell's famous dystopian book (it's a date):

 

 

Follow 3 latest changes.

Date Changed By Updated Field Previous Value => Replaced By
Tue 12 Feb 2008 01:41:09 PM UTCbjgStatusConfirmed=>Fixed
  Open/ClosedOpen=>Closed
  Summarybugs in gsl_cdf_hypergeometric_P/gsl_cdf_hypergeometric_Q=>integer overflow in gsl_cdf_hypergeometric_P/gsl_cdf_hypergeometric_Q

Back to the top


Powered by Savane 3.1-cleanup1