Re: [Bug-gsl] lambert W function

[Brian Gough]

At Thu, 26 Apr 2007 11:27:16 +0200,
Katrin Wolff wrote:
> An example program to reproduce the error on my installation of apple darwin 
> is:
> 
> #include <gsl/gsl_sf_lambert.h>
> #include <stdio.h>
> #include <math.h>
> #include <stdlib.h>
> 
> int main(int argc, char *argv[])
> { 
>       double x = 1.6849341956993197921832461361191235482692718505859;
>       double w = gsl_sf_lambert_W0(x);
> 
>       printf("x= %.50g\n", x);
>       printf("w= %.50g\n", w);
>    
>       return 0;
>   
> }
> 
> On Linux suse 9.3 it ran without error, however.


Thanks for the bug report.  I was able to reproduce the problem.
I've committed a fix to the development CVS repository.

-- 
Brian Gough

Network Theory Ltd,
Publishing Free Software Manuals --- http://www.network-theory.co.uk/

--- gsl/specfunc/lambert.c:1.6  Sun Jun 26 14:27:09 2005
+++ gsl/specfunc/lambert.c      Fri Apr 27 15:47:34 2007
@@ -59,7 +59,7 @@
 
     w -= t;
 
-    tol = GSL_DBL_EPSILON * GSL_MAX_DBL(fabs(w), 1.0/(fabs(p)*e));
+    tol = 10 * GSL_DBL_EPSILON * GSL_MAX_DBL(fabs(w), 1.0/(fabs(p)*e));
 
     if(fabs(t) < tol)
     {