bugGNU Scientific Library - Bugs: bug #24146, bug in gsl_sf_ellint_Kcomp_e for k...

 
 

You are not allowed to post comments on this tracker with your current authentication level.

bug #24146: bug in gsl_sf_ellint_Kcomp_e for k near 1

Submitter:  -Deleted Account- <bjg>
Submitted:  Tue 26 Aug 2008 04:36:55 PM UTC
   
 
Category:  Accuracy problem Severity:  3 - Normal
Operating System:  Status:  Fixed
Assigned to:  bjg Open/Closed:  Closed
Release:  1.11

Tue 26 Aug 2008 04:39:58 PM UTC, comment #1: 

fixed by 163bd7269d3c74af6b0a81ac82d5396a2afd3b3e

-Deleted Account- <bjg>
Tue 26 Aug 2008 04:36:55 PM UTC, original submission:  

From: Thies Heidecke <theidecke@astrophysik.uni-kiel.de>
To: -email is unavailable-
Subject: [Bug-gsl] 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 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


-Deleted Account- <bjg>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

CC list is empty

 

Follow 2 latest changes.

Date Changed by Updated Field Previous Value => Replaced by
2008-08-27 bjg Open/ClosedOpen Closed
2008-08-26 bjg StatusConfirmed Fixed

Back to the top

Powered by Savane 3.13-d3ae.
Corresponding source code