Added script front-end for primer-design code

[htsworkflow.git] / htswanalysis / MACS / lib / gsl / gsl-1.11 / BUGS

------------------------------------------------------------------------
BUG-ID:   21819   
STATUS:   Open/Confirmed  
CATEGORY: Runtime error
SUMMARY:  gsl_sf_bessel_J_CF1() crash for large arguments

Reply-To: Koichi Takahashi <ktakahashi@molsci.org>
From: Koichi Takahashi <ktakahashi.molsci@gmail.com>
To: bug-gsl@gnu.org
CC: Brian Gough <bjg@gnu.org>, Jonathan Taylor <j.m.taylor@durham.ac.uk>
Subject: gsl_sf_bessel__J_CF1 bug again
Date: Sun, 02 Dec 2007 00:11:02 -0800

Hi,

Unfortunately, I found some additional sets of numbers that still
crash gsl_sf_bessel_J_CF1() even with the cvs version.   Symptom
is exactly the same as what I reported before.  I tested on
x86_64.

main()
{
//    this used to crash, but now fixed in the current cvs.
//    double a = gsl_sf_bessel_jl( 30,  3875.6138424501978 );

//    at least the following three sets of values still crashes.
     double a = gsl_sf_bessel_jl( 49, 9912.6308956132361 );
//    double a = gsl_sf_bessel_jl( 49, 9950.3478935215589 );
//    double a = gsl_sf_bessel_jl( 52, 9930.5181281798232 );

     printf("%g\n",a);
}

Let me know if there is anything I could do to help you fixing
this issue.

thanks,
Koichi

From: Jonny Taylor <j.m.taylor@dur.ac.uk>
To: Koichi Takahashi <ktakahashi@molsci.org>
Cc: bug-gsl@gnu.org
Subject: [Bug-gsl] Re: gsl_sf_bessel__J_CF1 bug again
Date: Sun, 2 Dec 2007 10:49:33 +0000

While the symptom is the same, the cause is different. For those  
numbers it seems that 10,000 iterations is simply not enough.  
Interestingly, In fact, in all three cases it requires less than 50  
additional iterations to converge!?

The naive fix is simply to increase the maximum permitted number of  
iterations, but a more sophisticated fix would probably need to  
justify a maximum number of iterations, and propose an alternative  
method of generating the result in the cases where the number of  
iterations is prohibitive...

Jonny
From: Koichi Takahashi <ktakahashi@molsci.org>
To: Jonny Taylor <j.m.taylor@dur.ac.uk>
Cc: bug-gsl@gnu.org
Subject: [Bug-gsl] Re: gsl_sf_bessel__J_CF1 bug again
Date: Sun, 02 Dec 2007 06:06:11 -0800

At this range, x is still not large enough to use the asymptotic form?
In gsl_sf_bessel_jl_e,

   else if(x > 1000.0 && x > 100.0*l*l)
   {
     //asymptotic
   }
   else
   {
     //CF1
   }

so, for example, 100 * 50 * 50 = 250,000.
For l=50, the iteration starts to exceed 10,000 around x=9900.
If we want to stick to 10,000 max iteration, we have to switch to
the asymptotic version with something like x > 3*l*l.  Maybe this is
too small?   If so I'd consider increasing the max iteration.
Maybe we should do both.  Or there can be a better method.

koichi

------------------------------------------------------------------------
BUG-ID:   21826   
STATUS:   Open            
CATEGORY: Build
SUMMARY:  libtool problem on hp-ux

From: "Benoit, Gerard" <GERARD.BENOIT@astrium.eads.net>
To: bug-gsl@gnu.org
Subject: [Bug-gsl] gsl-1.9 : error in compilation scripts
Date: Fri, 13 Jul 2007 12:28:35 +0200

Hello,

I wish to report a bug in the compilation scripts:
GSL version 1.9  downloaded from ftp.gnu.org
computer HP 9000/785 under HP-UX  B.11.00 A 
config.guess -> hppa2.0w-hp-hpux11.00
compiler  : cc -> /opt/ansic/bin/cc   HP C  HP92453-01 A.11.01.20
compilation options : -Ae
linker  : ld -> /bin/ld

fatal error during the build of the shared libraries( first for cblas/libgslcblas.sl )
command executed : $archive_cmds 
gsl-1.9/libtool ( line 214 ) : archive_cmds="\$CC -b \${wl}+h \${wl}\$soname \${wl}+b \${wl}\$install_libdir -o \$lib \$libobjs \$deplibs \$compiler_flags"
coming from gsl-1.9/aclocal.m4 ( line 6541 ) : 
   _LT_AC_TAGVAR(archive_cmds, $1)='$CC -b ${wl}+h ${wl}$soname ${wl}+b ${wl}$install_libdir -o $lib $libobjs $deplibs $compiler_flags'

the compiler cc don't recognize the "-b" option ( it is a linker option ) and throw it : the linker looks for a main module and don't find it
with the correct form ${wl}-b to pass the option to the linker, cc add a start-up module ( like /usr/ccs/lib/crt0.o or /opt/langtools/lib/crt0.o ) which is useless for a shared library and has not be compiled with a PIC option ( +z or +Z ) and the linker stops.

In the preceding version 1.6 ( I have never downloaded version 1.7 and 1.8 ) the command was :
gsl-1.6/libtool ( line 187 ) : archive_cmds="\$LD -b +h \$soname +b \$install_libdir -o \$lib \$libobjs \$deplibs \$linker_flags"
from
gsl-1.6/aclocal.m4 ( line 2441 ) :    *) archive_cmds='$LD -b +h $soname +b $install_libdir -o $lib $libobjs $deplibs $linker_flags' ;;

This command works fine

so, look for a way to return to the LD command in gsl-1.9/aclocal.m4 ( line 6541 ) :

Best regards

Gérard BENOIT

Ce courriel (incluant ses éventuelles pièces jointes) peut contenir des informations confidentielles et/ou protégées ou dont la diffusion est restreinte. Si vous avez reçu ce courriel par erreur, vous ne devez ni le copier, ni l'utiliser, ni en divulguer le contenu à quiconque. Merci d’en avertir immédiatement l’expéditeur et d’effacer ce courriel de votre système. Astrium décline toute responsabilité en cas de corruption par virus, d’altération ou de falsification de ce courriel lors de sa transmission par voie électronique.

This email (including any attachments) may contain confidential and/or privileged information or information otherwise protected from disclosure. If you are not the intended recipient, please notify the sender immediately, do not copy this message or any attachments, do not use it for any purpose or disclose its content to any person, but delete this message and any attachments from your system. Astrium disclaims any and all liability if this email transmission was virus corrupted, altered or falsified.
----------------------------------------------------------------
Astrium SAS (393 341 516 RCS Paris) - Siège social: 6 rue Laurent Pichat, 75016 Paris, France
_______________________________________________
Bug-gsl mailing list
Bug-gsl@gnu.org
http://lists.gnu.org/mailman/listinfo/bug-gsl

------------------------------------------------------------------------
BUG-ID:   21828   
STATUS:   Open/Confirmed  
CATEGORY: Performance
SUMMARY:  suboptimal performance of gsl_fdfsolver_lmsder

From: "Alexander Usov" <a.s.usov@gmail.com>
To: help-gsl@gnu.org
Subject: [Help-gsl] Strange performance of gsl_fdfsolver_lmsder
Date: Wed, 24 Oct 2007 20:45:01 +0200

Hi all,

I am currently working on the problem involving source extraction from
astronomical images, which essentially boils down to fitting a number of
2d gaussians to the image.

One of the traditionally used fitters in this field is a Levenberg-Marquardt,
which gsl_fdfsolver_lmsder is and implementation of.

At some moment I have notices that for the bigger images (about 550
pixels, 20-30 parameters) gsl's lmsder algorithm spends a large fraction
of the run-time (about 50%) doing household transform.

While looking around for are different minimization algorithms I have made
a surprising finding that original netlib/minpack/lmder is almost twice faster
that that of gsl.

Could anyone explain such a big difference in performace?

-- 
Best regards,
  Alexander.

_______________________________________________
Help-gsl mailing list
Help-gsl@gnu.org
http://lists.gnu.org/mailman/listinfo/help-gsl

Reply-To: help-gsl@gnu.org
From: Brian Gough <bjg@network-theory.co.uk>
To: "Alexander Usov" <a.s.usov@gmail.com>
Cc: help-gsl@gnu.org
Subject: Re: [Help-gsl] Strange performance of gsl_fdfsolver_lmsder
Date: Thu, 25 Oct 2007 21:57:08 +0100

At Wed, 24 Oct 2007 20:45:01 +0200,
Alexander Usov wrote:
> At some moment I have notices that for the bigger images (about 550
> pixels, 20-30 parameters) gsl's lmsder algorithm spends a large fraction
> of the run-time (about 50%) doing household transform.
> 
> While looking around for are different minimization algorithms I have made
> a surprising finding that original netlib/minpack/lmder is almost twice faster
> that that of gsl.
> 
> Could anyone explain such a big difference in performace?

I have a vague memory that there was some quantity (Jacobian?) that
MINPACK only computes fully at the end, but in GSL it is accessible to
the user at each step so I felt I had to update it on each iteration
in the absence of some alternate scheme.  Sorry this is not a great
answer but I am not able to look at it in detail now.

-- 
Brian Gough

_______________________________________________
Help-gsl mailing list
Help-gsl@gnu.org
http://lists.gnu.org/mailman/listinfo/help-gsl

------------------------------------------------------------------------
BUG-ID:   21830   
STATUS:   Open/Confirmed  
CATEGORY: None
SUMMARY:  provide function to retrieve size of chebyshev data

Reply-To: help-gsl@gnu.org
From: Brian Gough <bjg@network-theory.co.uk>
To: help-gsl@gnu.org
Subject: Re: [Help-gsl] Writing Chebychev coefficients to disk
Date: Mon, 12 Mar 2007 21:40:21 +0000

At Mon, 12 Mar 2007 10:28:10 -0600,
Patrick Alken wrote:
> 
> It seems to me the simplest way to do this is just:
> 
> fwrite(gsl_cheb_ptr, sizeof(gsl_cheb_series), 1, fp);
> fwrite(gsl_cheb_ptr->c, sizeof(double), n, fp);
> 
> and later:
> 
> fread(gsl_cheb_ptr, sizeof(gsl_cheb_series), 1, fp);
> fread(gsl_cheb_ptr->c, sizeof(double), n, fp);
> 

Yes, the library should help in getting the correct value of "n"
(which is c->order + 1 in this case) to avoid subtle off-by-one
errors.

It's certainly more useful to expose the memory than provide the read
and write functions directly though.

-- 
Brian Gough

Network Theory Ltd,
Publishing the GSL Manual - http://www.network-theory.co.uk/gsl/manual/

_______________________________________________
Help-gsl mailing list
Help-gsl@gnu.org
http://lists.gnu.org/mailman/listinfo/help-gsl

None

------------------------------------------------------------------------
BUG-ID:   21831   
STATUS:   Open            
CATEGORY: Accuracy problem
SUMMARY:  Levý random number generator for alpha < 1

From: rafael@fis.unb.br
To: bug-gsl@gnu.org
Subject: [Bug-gsl] Levý random number generator
Date: Mon, 26 Mar 2007 19:48:01 -0300

The Levý skew random number generator (gsl_ran_levy_skew) does not  
procuce a Levý random number when beta=0 (symmetric case), and the  
gsl_ran_levy function does not work as stated in the docs. I made some  
histograms from 10^6 samples to check the accuracy of the algorithms,  
by comparison agaisnt the numerical integration of the equation of  
Levý's PDF. For the gsl_ran_levy function there is a good precison for  
alpha [1,2], for alpha (0.3,1) you must sum a series of random numbers  
to get the same precision (tipicaly 100 or more gsl_ran_levy numbers).  
For alpha<=0.3 the algorithm does not work properly, even worse, the  
error increases as you add more random numbers. This contradicts the  
manual that says "the algoritm only works for alpha (0,2]". The  
function gsl_ran_levy_skew does not produce levy random numbers when  
beta=0, instead the pdf of the random numbers is a linear (?!?!) one.

----------------------------------------------------------------
This message was sent using IMP, the Internet Messaging Program.

_______________________________________________
Bug-gsl mailing list
Bug-gsl@gnu.org
http://lists.gnu.org/mailman/listinfo/bug-gsl

From: rafael@fis.unb.br
To: Brian Gough <bjg@network-theory.co.uk>
Cc: 
Subject: Re: [Bug-gsl] Lev� random number generator
Date: Tue, 27 Mar 2007 09:35:15 -0300

Thanks for your quick answer, and sorry about my poor english, it is  
not my natural language.

The code below generates 10^6 random numbers, and makes a normalized  
histogram wich is compared to the levy pdf. To get the levy pdf, it  
numericaly integrates the characteristic function for levy process  
(the function f in the code). The n parameter just adds a series of  
levy numbers to get better precision. The code saves 2 files:  
lhist-$alpha-$n (the normalized histogram) and lpdf-$alpha (the pdf  
for the levy process). It also prints to the stdout the absolute error  
(square of the difference) between the histogram and the pdf.

The function levy skew shows problems for alpha<1.

With this code, you can also check the problems related to the  
gsl_ran_levy function (just change gsl_ran_levy_skew by gsl_ran_levy,  
cutting the last paramenter).

I am using the pre-compiled gsl that comes with debian etch (gsl  
version 1.8.2).

If you are interested, I also encoded a routine to generate levy skew  
random numbers, it is not fully tested, but it works for beta=0 and  
alpha<1 (it suffers from the same precision problem as gsl_ran_levy  
function for small alpha)

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <sys/time.h>
#include <unistd.h>
#include <string.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_statistics_double.h>
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_integration.h>

double f (double x, void *params)
{
        double alpha = *(double *) params;
        double f = exp(-pow(x,alpha))/M_PI;
        return f;
}

double *levy_pdf(double alpha,int hist_size,double a,double b)
{
        double abserr,*lpdf,dx;
        gsl_function F;
        int i;
        gsl_integration_workspace *w=gsl_integration_workspace_alloc(1000);
        gsl_integration_workspace *cw=gsl_integration_workspace_alloc(1000);
        gsl_integration_qawo_table  
*wf=gsl_integration_qawo_table_alloc(0.0,1.0,GSL_INTEG_COSINE,200);
        lpdf=(double*)calloc(hist_size,sizeof(double));
        F.function=&f;
        F.params=&alpha;
        dx=(double)(b-a)/(double)hist_size;
        for (i=0;i<hist_size;i++)
        {
                gsl_integration_qawo_table_set(wf,i*dx+a,1.0,GSL_INTEG_COSINE);
                gsl_integration_qawf (&F,0.0,1e-10,1000,w,cw,wf,&lpdf[i],&abserr);
        }
        gsl_integration_qawo_table_free(wf);
        gsl_integration_workspace_free(w);
        gsl_integration_workspace_free(cw);
        return (lpdf);
}

int main (int argc,char *argv[])
{
        double *l,*lpdf,a=-20,b=20,alpha,dx,n,errabs=0.0;
        unsigned long int rnd_seed;
        int i,j,rand_numbers=1e6,hist_size=400;
        gsl_histogram *h;
        gsl_rng *r;
        struct timeval *tv;
        struct timezone *tz;
        char filename[50];
        FILE *f1,*f2;
        if(argc!=3)
        {
                printf("\nThe program must be called with 2 parameters: alpha and n\n \n");
                exit(1);
        }
        dx=(double)(b-a)/(double)hist_size;
        h=gsl_histogram_alloc(hist_size);
        alpha=atof(argv[1]);
        n=atof(argv[2]);
        strcpy(filename,"lhist-");
        strcat(filename,argv[1]);
        strcat(filename,"-");
        strcat(filename,argv[2]);
        f1=fopen(filename,"w+");
        strcpy(filename,"lpdf-");
        strcat(filename,argv[1]);
        f2=fopen(filename,"w+");
        l=(double*)calloc(rand_numbers,sizeof(double));
        lpdf=(double*)calloc(hist_size,sizeof(double));
        r=gsl_rng_alloc (gsl_rng_mt19937);
        gettimeofday(tv,tz);
        rnd_seed=(unsigned long int)tv->tv_usec;
        gsl_rng_set(r,rnd_seed);
        i=0;
        do
        {
                l[i]=0.0;
                for (j=1;j<n;j++) l[i]+=gsl_ran_levy_skew(r,1.0,alpha,0.0);
                l[i]=l[i]/pow(n,1.0/alpha);
                if (abs(l[i])<=20) i++;//picks only random numbers in the interval  
[a,b] to get good precision in the histogram
        }while(i<rand_numbers);
        gsl_histogram_set_ranges_uniform(h,a,b);
        for(i=0;i<rand_numbers;i++) gsl_histogram_increment(h,l[i]);
        gsl_histogram_scale(h,(double)hist_size/((b-a)*gsl_histogram_sum(h)));
        gsl_histogram_fprintf(f1,h,"%g","%g");
        lpdf=levy_pdf(alpha,hist_size,a,b);
        for (i=0;i<hist_size;i++) fprintf(f2,"%e\t%e\n",i*dx+a,lpdf[i]);
        for (i=0;i<hist_size;i++) errabs+=pow((gsl_histogram_get(h,i)-lpdf[i]),2.0);
        printf("%e\n",errabs/(double)hist_size);
        gsl_histogram_free(h);
        gsl_rng_free(r);
        fclose(f1);
        exit (0);
}

> At Mon, 26 Mar 2007 19:48:01 -0300,
> rafael@fis.unb.br wrote:
>>
>> The Lev� skew random number generator (gsl_ran_levy_skew) does not
>> procuce a Lev� random number when beta=0 (symmetric case), and the
>> gsl_ran_levy function does not work as stated in the docs. I made some
>> histograms from 10^6 samples to check the accuracy of the algorithms,
>> by comparison agaisnt the numerical integration of the equation of
>> Lev�'s PDF. For the gsl_ran_levy function there is a good precison for
>> alpha [1,2], for alpha (0.3,1) you must sum a series of random numbers
>> to get the same precision (tipicaly 100 or more gsl_ran_levy numbers).
>> For alpha<=0.3 the algorithm does not work properly, even worse, the
>> error increases as you add more random numbers. This contradicts the
>> manual that says "the algoritm only works for alpha (0,2]". The
>> function gsl_ran_levy_skew does not produce levy random numbers when
>> beta=0, instead the pdf of the random numbers is a linear (?!?!) one.
>
> Thanks for your email.  Please can you send a small example program
> which demonstrates the problem.
>
> Note that the Levy skew generator is tested in the GSL test suite for
> several cases where beta=0 -- if you have not done so, can you run
> "make check" and confirm that it works for these cases.
>
> --
> Brian Gough
> (GSL Maintainer)
>
> Network Theory Ltd,
> Publishing the GSL Manual - http://www.network-theory.co.uk/gsl/manual/
>

----------------------------------------------------------------
This message was sent using IMP, the Internet Messaging Program.

_______________________________________________
Bug-gsl mailing list
Bug-gsl@gnu.org
http://lists.gnu.org/mailman/listinfo/bug-gsl

-

------------------------------------------------------------------------
BUG-ID:   21832   
STATUS:   Open/Fixed      
CATEGORY: None
SUMMARY:  gsl infnan.c configure/build bug on solaris

From: Richard Smith <Richard.Smith@Sun.COM>
To: bug-gsl@gnu.org
Subject: [Bug-gsl] gsl infnan.c configure/build bug
Date: Tue, 24 Jul 2007 10:24:22 +1000

infnan.c currently fails to compile on Solaris with Sun Studio 12 compilers
when using default compiler options as generated by running configure with
no options. Here is an extract of the output when running make:

/bin/bash ../libtool --tag=CC --mode=compile cc -DHAVE_CONFIG_H  -I. 
-I../../sys
 -I.. -I..     -g -c -o infnan.lo ../../sys/infnan.c
 cc -DHAVE_CONFIG_H -I. -I../../sys -I.. -I.. -g -c ../../sys/infnan.c  
-KPIC -D
PIC -o .libs/infnan.o
"/usr/include/ieeefp.h", line 74: syntax error before or at: 
__builtin_isfinite
cc: acomp failed for ../../sys/infnan.c

The problem seems to occur as a result of the following combination of 
things:
1. "finite" function is used in various source files
2. configure.ac checks for <ieeefp.h> and finds it
3. On Solaris, "finite" is defined in <ieeefp.h>
4. configure.ac only checks <math.h> for "finite", and doesn't find it.
5. Since its a c99 environment, "isfinite" is found.
6. config.h ends up containing amongst other things the following lines:
#define HAVE_DECL_FINITE 0
#define HAVE_DECL_ISFINITE 1
#if !HAVE_DECL_FINITE
#if HAVE_DECL_ISFINITE
#define finite isfinite
#else
#define finite gsl_finite
#endif
#endif
   This means that "finite" macro has the value "isfinite"
7. "isfinite" is a macro, ultimately defined in <iso/math_c99.h>:
#define        isfinite(x)     __builtin_isfinite(x)
8. <ieeefp.h> declares "finite" function:
extern int      finite(double);
9. After macro expansion the compiler sees
extern int      __builtin_isfinite ( double ) ;
and complains with an error.

There's probably multiple ways of overcoming the problem. A workaround is
to force the compiler not to use a C99 environment e.g. -xc99=%none. However
since it would be desirable to have configure work well by default, maybe
the test in configure.ac for "finite" should be something like:
AC_CHECK_DECLS(finite,,,[#include <math.h>
#if HAVE_IEEEFP_H
# include <ieeefp.h>
#endif])

Alternatively all uses of "finite" could be changed to use "isfinite", and
the configuration/build be based around assuming a C99 environment, with
substitute functions and macros being provided where the necessary C99
facilities are missing. On Linux, the BSD floating point classification
functions are documented as being obsolete, so their use should probably
be avoided.

-- 
============================================================================
   ,-_|\   Richard Smith - Technical Specialist
  /     \  Sun Microsystems Australia         Phone : +61 3 9869 6200
richard.smith@Sun.COM                        Direct : +61 3 9869 6224
  \_,-._/  476 St Kilda Road                    Fax : +61 3 9869 6290
       v   Melbourne Vic 3004 Australia
=========================================================================== 

_______________________________________________
Bug-gsl mailing list
Bug-gsl@gnu.org
http://lists.gnu.org/mailman/listinfo/bug-gsl

------------------------------------------------------------------------
BUG-ID:   21833   
STATUS:   Open            
CATEGORY: Performance
SUMMARY:  suboptimal performance of gsl permutation?

From: "djamel anonymous" <djam8193ah@hotmail.com>
To: bjg@network-theory.co.uk
Subject: gsl permutation
Date: Wed, 24 Jan 2007 07:42:06 +0000

Hi.
i am sending you this email about a possible issue in glibc permutation 
algorithm.i think it has qudratic worst case running time.for example if we 
have the permutation
1,2,3,4,,,,,n,0
we will permute all elements in first iteration then for each iteration we 
will traverse on average half the cycle (which is of size n) before we know 
that the elements of cycle have already been permuted.there is two possible 
solutions to make the algorithm linear:
1-at each step we traverse the full cycle and permute all elements in the 
cycle.for each permuted element i  we assign p[i]=i.the disavantage of this 
is that it will destroy the original permutation.
2-use a bit array of size n bits.each time we permute an element we set the 
relevant bit.if at iteration i we find that bit i is already set we skip to 
next iteration.the disavantage of this is that it equires additional storage 
allocation.
best regards.

_________________________________________________________________
MSN Hotmail sur i-mode\99 : envoyez et recevez des e-mails depuis votre 
téléphone portable ! http://www.msn.fr/hotmailimode/

3 - Normal

------------------------------------------------------------------------
BUG-ID:   21835   
STATUS:   Open            
CATEGORY: Accuracy problem
SUMMARY:  gsl_sf_hyperg_2F1 problematic arguments

BUG#1 -- gsl_sf_hyperg_2F1_e fails for some arguments 

From: keith.briggs@bt.com
Subject: gsl_sf_hyperg_2F1 bug report
Date: Thu, 31 Jan 2002 12:30:04 -0000

gsl_sf_hyperg_2F1_e fails with arguments (1,13,14,0.999227196008978,&r).
It should return 53.4645... .

#include <gsl/gsl_sf.h>
#include <stdio.h>

int main (void)
{
  gsl_sf_result r;
  gsl_sf_hyperg_2F1_e (1,13,14,0.999227196008978,&r);
  printf("r = %g %g\n", r.val, r.err);
}

NOTES: The program overflows the maximum number of iterations in
gsl_sf_hyperg_2F1, due to the presence of a nearby singularity at
(c=a+b,x=1) so the sum is slowly convergent.

The exact result is 53.46451441879150950530608621 as calculated by
gp-pari using sumpos(k=0,gamma(a+k)*gamma(b+k)*gamma(c)*gamma(1)/
(gamma(c+k)*gamma(1+k)*gamma(a)*gamma(b))*x^k)

The code needs to be extended to handle the case c=a+b. This is the
main problem. The case c=a+b is special and needs to be computed
differently.  There is a special formula given for it in Abramowitz &
Stegun 15.3.10

As reported by Lee Warren <warren@atom.chem.utk.edu> another set of
arguments which fail are:

#include <gsl/gsl_sf.h>
#include <stdio.h>

int main (void)
{
  gsl_sf_result r;
  gsl_sf_hyperg_2F1_e (-1, -1, -0.5, 1.5, &r);
  printf("r = %g %g\n", r.val, r.err);
}

The correct value is -2.

See also, 

From: Olaf Wucknitz <wucknitz@jive.nl>
To: bug-gsl@gnu.org
Subject: [Bug-gsl] gsl_sf_hyperg_2F1

Hi,

I am having a problem with gsl_sf_hyperg_2F1.
With the parameters (-0.5, 0.5, 1, x) the returned values show a jump at 
x=0.5. For x<0.5 the results seem to be correct, while for x>0.5 they 
aren't.
The function gsl_sf_hyperg_2F1_e calls hyperg_2F1_series for x<0.5, but
hyperg_2F1_reflect for x>0.5. The latter function checks for c-a-b being 
an integer (which it is in my case). If I change one of the parameters 
a,b,c by a small amount, the other branch of the function is taken and the 
results are correct again.
Unfortunately I know too little about the numerics of 2F1 to suggest a 
patch at the moment.

Regards,
Olaf Wucknitz
-- 
Joint Institute for VLBI in Europe                wucknitz@jive.nl

------------------------------------------------------------------------
BUG-ID:   21836   
STATUS:   Open/Confirmed  
CATEGORY: Accuracy problem
SUMMARY:  gamma_inc_P and gamma_inc_Q only satisfy P+Q=1 within errors

BUG#44 -- gamma_inc_P and gamma_inc_Q only satisfy P+Q=1 within errors

The sum of gamma_inc_P and gamma_inc_Q doesn't always satisfy the
identity P+Q=1 exactly (although it is correct within errors), due the
slightly different branch conditions for the series and continued
fraction expansions.  These could be made identical so that P+Q=1 exactly.

#include <stdio.h>
#include <gsl/gsl_sf_gamma.h>

int
main (void)
{
  gsl_sf_result r1, r2;
  double a = 0.3, x = 1.0;
  gsl_sf_gamma_inc_P_e (a, x, &r1);
  gsl_sf_gamma_inc_Q_e (a, x, &r2);
  printf("%.18e\n", r1.val);
  printf("%.18e\n", r2.val);
  printf("%.18e\n", r1.val + r2.val);
}

$ ./a.out
9.156741562411074842e-01
8.432584375889111417e-02
9.999999999999985567e-01

3 - NormalNone

------------------------------------------------------------------------
BUG-ID:   21837   
STATUS:   Open/Confirmed  
CATEGORY: Runtime error
SUMMARY:  gsl_linalg_solve_symm_tridiag requires positive definite matrix

A zero on the diagonal will cause NaNs even though a reasonable
solution could be computed in principle. 

#include <gsl/gsl_linalg.h>

int main (void)
{
  double d[] = { 0.00, 1.21, 0.80, 1.55, 0.76 } ;
  double e[] = { 0.82, 0.39, 0.09, 0.68 } ;
  double b[] = { 0.07, 0.62, 0.81, 0.11, 0.65} ;
  double x[] = { 0.00, 0.00, 0.00, 0.00, 0.00} ;

  gsl_vector_view dv = gsl_vector_view_array(d, 5);
  gsl_vector_view ev = gsl_vector_view_array(e, 4);
  gsl_vector_view bv = gsl_vector_view_array(b, 5);
  gsl_vector_view xv = gsl_vector_view_array(x, 5);

  gsl_linalg_solve_symm_tridiag(&dv.vector, &ev.vector, &bv.vector, &xv.vector);
  gsl_vector_fprintf(stdout, &xv.vector, "% .5f");

  d[0] += 1e-5;
  gsl_linalg_solve_symm_tridiag(&dv.vector, &ev.vector, &bv.vector, &xv.vector);
  gsl_vector_fprintf(stdout, &xv.vector, "% .5f");
}

$ ./a.out
 nan
 nan
 nan
 nan
 nan
 0.13626
 0.08536
 1.03840
-0.60009
 1.39219

AUG 2007: We now return an error code for this case.  To return a solution
we would need to do a permutation, see slatec/dgtsl.f

3 - NormalNone

------------------------------------------------------------------------
BUG-ID:   22478   
STATUS:   Open/Confirmed  
CATEGORY: None
SUMMARY:  Missing functions for complex vectors

From: Federico Zenith <federico.zenith@member.fsf.org>
To: help-gsl@gnu.org
Subject: [Help-gsl] Missing functions for complex vectors
Date: Mon, 03 Mar 2008 13:50:43 +0100

Hi,
I am working on a C++ GSL wrapper (yes, there are not enough of them
already around :-), and I have been working on vectors last week; while
writing the wrapper, I noticed a couple of oddities I'd like to point
out before reporting them as a bug (I am using the GSL version 1.10):

1) complex vectors have defined all properties mentioned on this page:
http://www.gnu.org/software/gsl/manual/html_node/Vector-properties.html
except isnonneg(). I am not sure what's the meaning of ispos() or
isneg() for a vector of complex numbers, I suppose it means "both real
AND imaginary part are positive", but this should be documented;
furthermore, I do not understand why, if ispos() and isneg() could be
defined, why could not isnonneg() be as well.

2) Vector operations, defined on this page:
http://www.gnu.org/software/gsl/manual/html_node/Vector-operations.html
are not defined for complex types. As far as I can guess, all these
operations would still make sense for complex numbers. What is really
odd is, the same operations are defined for complex matrices.

Did I miss something?

Cheers,
-Federico

------------------------------------------------------------------------