From: Thies Heidecke <theidecke@astrophysik.unikiel.de>
To: email is unavailable
Subject: [Buggsl] bug in gsl_sf_ellint_Kcomp_e due to typo
Date: Mon, 25 Aug 2008 17:02:49 +0200
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 testcode:
#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.0pow(10.0,(double)(i));
double y = gsl_sf_ellint_Kcomp(sqrt(x), GSL_PREC_DOUBLE);
printf("K(%g) = %.18e\n", 1.0x, y);
}
return 0;
}
the value should approach infinity but stays nearly constant for
sqrt(x) > 1  1e8
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
