[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[bug #60741] Inaccurate Results for Lambert W function
From: |
Patrick Alken |
Subject: |
[bug #60741] Inaccurate Results for Lambert W function |
Date: |
Sun, 6 Jun 2021 23:49:18 -0400 (EDT) |
User-agent: |
Mozilla/5.0 (X11; Ubuntu; Linux x86_64; rv:89.0) Gecko/20100101 Firefox/89.0 |
URL:
<https://savannah.gnu.org/bugs/?60741>
Summary: Inaccurate Results for Lambert W function
Project: GNU Scientific Library
Submitted by: psa
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:
Discussion Lock: Any
_______________________________________________________
Details:
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
----
_______________________________________________________
Reply to this item at:
<https://savannah.gnu.org/bugs/?60741>
_______________________________________________
Message sent via Savannah
https://savannah.gnu.org/
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [bug #60741] Inaccurate Results for Lambert W function,
Patrick Alken <=