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 authentification level.

bug #24146: bug in gsl_sf_ellint_Kcomp_e for k near 1

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

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

fixed by 163bd7269d3c74af6b0a81ac82d5396a2afd3b3e

-Deleted Account- <bjg>In charge of this item.
Tue 26 Aug 2008 04:36:55 PM UTC, original submission:

From: Thies Heidecke <theidecke@astrophysik.uni-kiel.de>
To: -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] + ya[2]);
const double tb = -log(y) * (b[0] * y(b[1] + yb[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] + yb[2]));

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

after recompiling the function now gives reasonable values.

Regards,
Thies Heidecke

-Deleted Account- <bjg>In charge of this item.

 

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

CC list is empty

 

Please enter the title of George Orwell's famous dystopian book (it's a date):

 

 

Follow 2 latest changes.

Date Changed By Updated Field Previous Value => Replaced By
Wed 27 Aug 2008 07:39:39 AM UTCbjgOpen/ClosedOpen=>Closed
Tue 26 Aug 2008 04:39:58 PM UTCbjgStatusConfirmed=>Fixed

Back to the top


Powered by Savane 3.1-cleanup