[Bug-gsl] Bug in gsl_cdf_gamma_Q (x86-only?)

[Sean O'Rourke]

I stumbled across an apparently x86-only (or only since 1.6?) bug in
gsl_cdf_gamma_P:

#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>
#include <stdio.h>

int
main(int argc, char * argv[])
{
    double a, b, x, r;
    a = strtod(argv[1], NULL);
    b = strtod(argv[2], NULL);
    printf("a=%g,b=%g\n", a, b);
    while (scanf("%lg", &x) > 0)
        printf("%lg: P = %lg, Q = %lg\n", x, gsl_cdf_gamma_P(x, a,
b),
               gsl_cdf_gamma_Q(x, a, b));
}

######## x86 machine ###########

puzzle:~% gsl-config --version
1.6
puzzle:~% uname -m
i686
puzzle:~% ./a.out 5670 0.00240449
a=5670,b=0.00240449
10.9
10.9: P = 1.33244e-59, Q = 1
11
11: P = 3.69123e-55, Q = -3.69123e-55

######## PPC ###########

80363334:data/0626% gsl-config --version
1.5
80363334:data/0626% uname -m
Power Macintosh
80363334:data/0626% /tmp/a.out 5670 0.00240449   
a=5670,b=0.00240449
10.9
10.9: P = 1.33244e-59, Q = 1
11
11: P = 3.69123e-55, Q = 1

Let me know if you have any questions,

Sean


                
____________________________________________________
Sell on Yahoo! Auctions – no fees. Bid on great items.  
http://auctions.yahoo.com/