# GNU Scientific Library - Bugs: bug #60741, Inaccurate Results for Lambert W...

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

# bug #60741: Inaccurate Results for Lambert W function

 Submitted by: Patrick Alken Submitted on: Mon 07 Jun 2021 03:49:16 AM UTC Category: Runtime error Severity: 3 - Normal Operating System: Status: None Assigned to: None Open/Closed: Open Release:
* Mandatory Fields

Comment Type & Canned Response:

## Discussion

Mon 07 Jun 2021 03:49:16 AM UTC, original submission:

from daming.zou =at= inf =dot= ethz =dot= ch

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

Patrick Alken <psa>

## Attached Files

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

Attach Files:

Comment:

No files currently attached

## Dependencies

Depends on the following items: None found

Items that depend on this one: None found

Carbon-Copy List
• -email is unavailable- added by psa (Submitted the item)
•

## History

No changes have been made to this item