This is the mail archive of the gsl-discuss@sources.redhat.com mailing list for the GSL project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

Bug fix (proposed): Incomplete gamma function, large parameters equal to each other


Hi all

  Special thanks to Gerard Jungman for his work on these functions.

  I'd like to propose a bugfix for gamma_inc_Q_asymp_unif() in
specfunc/gamma_inc.c near line 150.  Note that

  const double eps = (x-a)/a;


but soon after we set

    eta  = eps * sqrt(-2.0*ln_term.val/(eps*eps))


which of course does not work if eps is zero (or effectively zero).  It
seems to me the fix is simple, just by inserting something like:

    if (eps<GSL_SQRT_DBL_EPSILON) {
        eta = sqrt(-2.0*ln_term.val);
    }

but I must say I'm far from sure GSL_SQRT_DBL_EPSILON is the correct bound.




Brian


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]