Inaccurate Results for Lambert W function

[Daming Zou]

Hi,

The gsl_sf_lambert_W0(x) function gives inaccurate results for inputs 
near 0. I attached a code snippet to demonstrate it at the end of this 
mail. The details are as following:

[Error Examples]
- For very small input, e.g., 1e-30, 1e-40, the returned results are 
inaccurate;
- For smaller ones, e.g., 1e-50, 1e-60, the returned results are 0.

[Accurate Result]
- The series of lambertW(x) near 0 is W(x) = x - x^2 + 3/2 x^3 - 8/3 x^4 
+ 125/24 x^5 + O(x^6), for more terms, please refer the Table 4 of this 
paper[1].
- For inputs |x|<1e-20 (more precisely, the double epsilon 
GSL_DBL_EPSILON), and for double precision, we can consider W(x) = x, 
since the -x^2 term can be neglected. Thus, this applies to 1e-30, 
1e-40, etc.

[Reason for Inaccurate Results]
- Currently, for input values larger than -1/E+1e-3, 
gsl_sf_lambert_W0_e(x) uses *halley_iteration* for solving the result. 
However, for extremely small values near 0, there is severe cancellation 
at line 61 of lambert.c : `w -= t;`, so the results are inaccurate or 
even 0.

[Possible Solution]
- Just return x for |x| < 1e-20.
- More aggressively, using polynomial to approximate the function for 
|x| < 1e-4. (A 5-terms polynomial should be sufficient for double 
precision). The interval and the number of terms could be further 
investigated.

[Reference]
[1] Fukushima, Toshio. "Precise and fast computation of Lambert 
W-functions without transcendental function evaluations." 
https://www.sciencedirect.com/science/article/pii/S0377042712005213

Not sure if this is defined as a bug in GSL, since from the aspect of 
absolute error, it may be acceptable; but from the aspect of relative 
error, it is significant. However, I believe the result W(1e-50) = 1e-50 
is better than W(1e-50) = 0 in most scenarios.

Thanks,
Daming


[Attached Code Example]

#include <gsl/gsl_sf.h>
#include <stdio.h>

const double c[5] = {
     1.0,
    -1.0,
     3.0/2.0,
    -8.0/3.0,
     125.0/24.0
};

double horner_eval(double x) {
    double y = x*(c[0] + x*(c[1] + x*(c[2] + x*(c[3] + x*c[4]))));
    return y;
}

double inputs[7] = {
    1e-5,
    1e-10,
    1e-20,
    1e-30,
    1e-40,
    1e-50,
    1e-60
};

int main() {
    double x, y, y_poly;
    for (int i = 0; i < 7; i++) {
        x = inputs[i];
        y = gsl_sf_lambert_W0(x);
        y_poly = horner_eval(x);
        printf("input:              %.20e\n", x);
        printf("gsl_sf_lambertW:    %.20e\n", y);
        printf("poly_output:        %.20e\n", y_poly);
        printf("----\n");
    }
    return 0;
}


$ ./gsl-lambert.out
input:              1.00000000000000008180e-05
gsl_sf_lambertW:    9.99990000149997397560e-06
poly_output:        9.99990000149997397560e-06
----
input:              1.00000000000000003643e-10
gsl_sf_lambertW:    9.99999999899999974982e-11
poly_output:        9.99999999899999974982e-11
----
input:              9.99999999999999945153e-21
gsl_sf_lambertW:    9.99999999999999945153e-21
poly_output:        9.99999999999999945153e-21
----
input:              1.00000000000000008334e-30
gsl_sf_lambertW:    1.00000000000275223352e-30
poly_output:        1.00000000000000008334e-30
----
input:              9.99999999999999929293e-41
gsl_sf_lambertW:    1.03314933177740113005e-40
poly_output:        9.99999999999999929293e-41
----
input:              1.00000000000000000762e-50
gsl_sf_lambertW:    0.00000000000000000000e+00
poly_output:        1.00000000000000000762e-50
----
input:              9.99999999999999970433e-61
gsl_sf_lambertW:    0.00000000000000000000e+00
poly_output:        9.99999999999999970433e-61
----