Source-linked AI summary
Lambert W Function for Applications in Physics
Darko Veberic
TL;DR
The paper addresses the need for fast and accurate numerical evaluation of Lambert W across its real branches. It develops a C++ implementation combining Halley’s and Fritsch’s iterations with range-specific initial approximations. The implementation achieves machine-size accuracy and substantial speedups over GSL in the reported tests.
Problem
The paper addresses the need for a fast and accurate numerical implementation of the Lambert W function.
Method
The implementation combines Halley’s and Fritsch’s iterations with branch-point expansions, asymptotic series, rational fits, and continued-logarithm recursion for initial approximations.
Results
The implementation is at least 2× faster than GSL broadly and up to 5× faster where rational fits are used, while the reported branch approximations reach machine-size accuracy.
Takeaways & Limitations
The open-source C++ implementation is suitable for inclusion in analysis packages used across physics.
Takeaways & Limitations
The recursion-based form described is restricted to x > 0 and W0(x) > 1, corresponding to x > e.
Abstract
from arXiv · showhide
The Lambert W(x) function and its possible applications in physics are presented. The actual numerical implementation in C++ consists of Halley's and Fritsch's iterations with initial approximations based on branch-point expansion, asymptotic series, rational fits, and continued-logarithm recursion.
Program summary
The program provides a fast, accurate C++ implementation of the Lambert W function using iterative solvers supplied with range-specific initial approximations.
- The C++ program can be used on Unix and Windows systems with a C++ compiler and has a memory footprint below 1 MB.
- The implementation targets fast and accurate numerical evaluation of the Lambert W function.
- Halley’s and Fritsch’s iterations are combined with branch-point expansions, asymptotic series, rational fits, and continued-logarithm recursion.
- The distribution includes the lambert-w command-line utility, source documentation, and a Makefile.
1. Introduction
The introduction defines Lambert W as the inverse of a simple exponential mapping, describes its two real branches and broad applications, and motivates numerical implementation through astrophysical examples.
- Lambert W is the inverse function associated with the exponential mapping W(x)e^W(x)=x.
- The two real branches meet at (−e−1, −1), with W−1 defined on [−e−1, 0] and W0 on [−e−1, ∞].
- Lambert W has applications across mathematics, numerics, computer science, engineering, and physics, including quantum mechanics, statistics, relativity, fluid dynamics, and optics.
- The paper highlights cosmic-ray physics as a motivation because Lambert W has already been used there and new applications continue to appear.
- The inverse Moyal function and inverse Gaisser-Hillas function can be expressed using Lambert W branches, with the latter’s branch choice selecting the side of its maximum.
2. Numerical methods
The paper reviews recursive logarithmic approaches and develops Halley’s and Fritsch’s iterations for evaluating both real Lambert W branches. Fritsch’s fourth-order scheme can reach machine-size precision in one iteration from sufficiently accurate initial approximations.
- Recursion: Continued-logarithm recursions provide successive approximations for W0(x) and W−1(x), with the principal-branch form applicable under specified positive-domain conditions.For W0(x), the recursive form is suitable when x > e; analogous recursion is derived for W−1(x).
- Halley’s iteration: Halley’s iteration is third order, reducing an O(ε) error to O(ε^3) per step.An initial approximation of order O(10^-4) is expected to reach O(10^-16) machine precision in at least two iterations.
- Fritsch’s iteration: Fritsch’s iteration is fourth order, reducing an O(ε_n) error to O(ε_n^4).With an initial guess accurate to O(10^-4), it delivers O(10^-16) precision in one iteration and O(10^-64) in two.
- Initial approximations: The higher-order iteration requires reliable first approximations across the distinct behavioral ranges of both Lambert W branches.The following section therefore constructs initial approximations over the complete definition ranges.
3. Initial approximations
Initial values are constructed across the two real branches using branch-point expansions, asymptotic series, rational fits, and continued-logarithm recursion. These approximations are selected and combined over intervals to support accurate numerical iteration.
- Branch-point expansion: Near x = −e^-1, branch-point expansion captures the square-root behavior of both real Lambert W branches.The expansion is developed around the minimum and plotted at successive orders for W−1(x) and W0(x).
- Asymptotic series: Asymptotic expansions use logarithmic variables defined differently for W0 and W−1.For W0, a = ln x and b = ln ln x; for W−1, a = ln(−x) and b = ln(−ln(−x)).
- Rational fits: Rational approximations are fitted from sampled pairs {w_i e^w_i, w_i} and selected by minimizing maximal absolute residual over an interval.Separate fits are constructed for W0 ranges and for the W−1 branch, with polynomial orders varied during selection.
4. Implementation
The implementation combines branch-specific approximations into piecewise functions, then refines them with Halley’s or Fritsch’s iteration. For W0(x), one Fritsch step is selected for the C++ implementation because it provides machine-size accuracy across the definition interval without Halley’s extra step in a restricted range.
- W0(x) construction: For W0(x), branch-point, rational-fit, and asymptotic approximations are assigned optimal intervals and combined into a piecewise approximation over the full domain.The intervals are evaluated on linear [−e−1, 0.3] and logarithmic [0.3, 105] ranges.
- W0(x) construction: The W0(x) piecewise approximation is accurate to at least 5 decimal places on [−e−1, 7] and at least 3 decimal places across the whole definition range.The construction uses transition parameters a = −0.323581′, b = 0.145469′, and c = 8.706658′.
- Iteration choice: Both iterations produce machine-size accurate W0(x) floating-point results across the definition interval, except Halley’s method needs another step between 6.5 and 190.The implementation therefore uses one step of Fritsch’s iteration for the Lambert W function.
- W−1(x) construction: For W−1(x), Fritsch’s iteration yields machine-size accuracy across the definition range, while Halley’s iteration reaches at least 13 decimal places after one step.The combined W−1(x) approximation is refined after selecting branch-specific approximations.
- Component approximations: The W0(x) branch-point approximation B[9] is accurate to at least 5 decimal places throughout [−e−1, 0].Its coefficients are a = −0.302985′ and b = −0.051012′.
5. Source availability, installation and usage
The open-source C++ implementation provides Lambert W functions, initial approximations, command-line access, and testing support across common systems and scripting environments.
- Availability: The implementation is freely available in C++ and is described as the only known implementation using Fritsch’s iteration.It is released under dual GPL/BSD licensing.
- API: The supplied code implements approximation and evaluation functions for compile-time branches 0 and −1, plus a runtime branch-selection interface.The functions are provided in LambertW.h and LambertW.cc.
- Installation: The source files can be compiled directly into a project with a suitable C++ compiler, without special installation procedures.The stated target environments include Unix flavors and Windows, with a memory footprint below 1 MB.
- Usage: The command-line utility accepts an optional branch number and x, returning the corresponding Lambert W value for use in shell scripts and other languages.Branch 0 is the default.
- Testing: The distribution includes a test suite that reports discrepancies larger than approximately 10−14 when comparing obtained and expected return values.Tests are invoked with make tests.
6. Timing
The authors benchmark their implementation against GSL using overhead-corrected timing and repeated function calls. It is generally faster, with the largest gains occurring where rational fits are used.
- Method: Timing compares pure numerical execution by subtracting surrounding driver and function-call overhead from measured runtimes.The corrected ratio is t′gsl/t′.
- Method: The benchmark uses 3 000 000 function calls and modifies x between calls to prevent compiler optimization from removing the computation.An identity-function control estimates timing overhead.
- Interpretation: GSL performs better only in small regions where branch-point and asymptotic expansions are used without subsequent Halley iteration refinements.The authors attribute the broader speed advantage to accurate initial approximations and one-step Fritsch iteration.
7. Conclusions
The paper concludes that combining accurate initial approximations with Fritsch’s iteration improves the efficiency of evaluating the real Lambert W branches. The open-source C++ implementation is intended for physics analysis packages.
- Conclusion: Accurate initial approximations coupled with Fritsch’s iteration deliver significant efficiency gains for numerical evaluation of the real Lambert W branches.The conclusion presents this combination as the central computational result.
- Conclusion: The released C++ implementation is suitable for inclusion in analysis packages used across physics fields.The stated scope is the open-source implementation and its application in physics software.