[Bug-gsl] bug in gsl_sf_ellint_Kcomp_e due to typo

[Thies Heidecke]

i'm using version 1.11 of the library.
while testing i found that

        gsl_sf_ellint_Kcomp_e(x, GSL_PREC_DOUBLE);

gives wrong results if x approaches the singularity at 1.0 as can be  
seen with this test-code:


#include <stdio.h>
#include <math.h>
#include <gsl/gsl_sf_ellint.h>

int main (int argc, char const* argv[])
{
        int i;
        
        for (i=1;i<=16;++i) {
                double x = 1.0-pow(10.0,(double)(-i));
                double y = gsl_sf_ellint_Kcomp(sqrt(x), GSL_PREC_DOUBLE);
                printf("K(%g) = %.18e\n", 1.0-x, y);
        }
        
        return 0;
}


the value should approach infinity but stays nearly constant for  
sqrt(x) > 1 - 1e-8

looking at ellint.c the error is a simple typo in the part of the code  
which handles the computation near the singularity:

/* [Carlson, Numer. Math. 33 (1979) 1, (4.5)] */
int
gsl_sf_ellint_Kcomp_e(double k, gsl_mode_t mode, gsl_sf_result * result)
{
  if(k*k >= 1.0) {
    DOMAIN_ERROR(result);
  }
  else if(k*k >= 1.0 - GSL_SQRT_DBL_EPSILON) {
    /* [Abramowitz+Stegun, 17.3.33] */
    const double y = 1.0 - k*k;
    const double a[] = { 1.38629436112, 0.09666344259, 0.03590092383 };
    const double b[] = { 0.5, 0.12498593597, 0.06880248576 };
    const double ta = a[0] + y*(a[1] + y*a[2]);
    const double tb = -log(y) * (b[0] * y*(b[1] + y*b[2]));                     
<--------
    result->val = ta + tb;
    result->err = 2.0 * GSL_DBL_EPSILON * result->val;
    return GSL_SUCCESS;
  }
  else {
        ...
  }
}

line 489 reads:
    const double tb = -log(y) * (b[0] * y*(b[1] + y*b[2]));

but it should be
    const double tb = -log(y) * (b[0] + y*(b[1] + y*b[2]));

after recompiling the function now gives reasonable values.

Regards,
        Thies Heidecke