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
----