Editing Lambert W function (section)

== Numerical evaluation ==
The {{mvar|W}} function may be approximated using [[Newton's method]], with successive approximations to {{math|1=''w'' = ''W''(''z'')}} (so {{math|1=''z'' = ''we<sup>w</sup>''}}) being
: <math>w_{j+1}=w_j-\frac{w_j e^{w_j}-z}{e^{w_j}+w_j e^{w_j}}.</math>

The {{mvar|W}} function may also be approximated using [[Halley's method]],
: <math>
w_{j+1}=w_j-\frac{w_j e^{w_j}-z}{e^{w_j}\left(w_j+1\right)-\dfrac{\left(w_j+2\right)\left(w_je^{w_j}-z\right)}{2w_j+2}}
</math>
given in Corless et al.<ref name="Corless" /> to compute {{mvar|W}}.

For real <math>x \ge -1/e</math>, it may be approximated by the quadratic-rate recursive formula of R. Iacono and J.P. Boyd:<ref name="doi.org"/>
: <math>w_{n+1} (x) = \frac{w_{n} (x)}{1 + w_{n} (x)} \left( 1 + \log \left(\frac{x}{w_{n} (x)} \right) \right).</math>

Lajos Lóczi proves<ref>{{Cite journal |last=Lóczi |first=Lajos |date=2022-11-15 |title=Guaranteed- and high-precision evaluation of the Lambert W function |journal=Applied Mathematics and Computation |language=en |volume=433 |pages=127406 |doi=10.1016/j.amc.2022.127406 |doi-access=free |issn=0096-3003 |hdl=10831/89771 |hdl-access=free |url=https://www.researchgate.net/publication/362219191}}</ref> that by using this iteration with an appropriate starting value <math>w_0 (x)</math>,
* For the principal branch <math>W_0:</math>
** if <math>x \in (e,\infty)</math>: <math>w_0 (x) = \log(x) - \log(\log(x)),</math>
** if <math>x \in (0, e):</math> <math>w_0 (x) = x/e,</math>
** if <math>x \in (-1/e, 0):</math> <math>w_0 (x) = \frac{ ex \log(1+\sqrt{1+ex}) }{ 1+ ex + \sqrt{1+ex} },</math> 
* For the branch <math>W_{-1}:</math>
** if <math>x \in (-1/4, 0):</math> <math>w_0 (x) = \log(-x) - \log(-\log(-x)),</math>
** if <math>x \in (-1/e, -1/4]:</math> <math>w_0 (x) = -1 - \sqrt{2}\sqrt{1+ex},</math>
one can determine the maximum number of iteration steps in advance for any precision:
* if <math>x \in (e,\infty)</math> (Theorem 2.4): <math>0 < W_0 (x) - w_n(x) < \left( \log(1+1/e) \right)^{2^n},</math>
* if <math>x \in (0, e)</math> (Theorem 2.9): <math>0 < W_0 (x) - w_n(x) < \frac{\left( 1 - 1/e \right)^{2^n-1}}{5},</math>
* if <math>x \in (-1/e, 0):</math>
** for the principal branch <math>W_0</math> (Theorem 2.17): <math>0 < w_n(x) - W_0 (x) < \left( 1/10 \right)^{2^n},</math>
** for the branch <math>W_{-1}</math>(Theorem 2.23): <math>0 < W_{-1} (x) - w_n(x) < \left( 1/2 \right)^{2^n}.</math>


Toshio Fukushima has presented a fast method for approximating the real valued parts of the principal and secondary branches of the {{mvar|W}} function without using any iteration.<ref>{{Cite web|last=Fukushima |first=Toshio |date=2020-11-25 |title=Precise and fast computation of Lambert W function by piecewise minimax rational function approximation with variable transformation |doi=10.13140/RG.2.2.30264.37128 |doi-access=free |url=https://www.researchgate.net/publication/346309410}}</ref> In this method the {{mvar|W}} function is evaluated as a conditional switch of [[rational functions]] on transformed variables:
<math display="block">W_0(z) = \begin{cases}
  X_k(x), & (z_{k-1}<=z<z_k, \quad k=1,2,\ldots,17), \\
  U_k(u), & (z_{k-1}<=z<z_k, \quad k=18,19),
\end{cases}
</math>
<math display="block">W_{-1}(z) = \begin{cases}
  Y_k(y), & (z_{k-1}<=z<z_k, \quad k=-1,-2,\ldots,-7), \\
  V_k(u), & (z_{k-1}<=z<z_k, \quad k=-8,-9,-10),
\end{cases}
</math>
where {{mvar|x}}, {{mvar|u}}, {{mvar|y}} and {{mvar|v}} are transformations of {{mvar|z}}:
: <math>x=\sqrt{z+1/e}, \quad u=\ln{z}, \quad y=-z/(x+1/\sqrt{e}), \quad v=\ln(-z)</math>.

Here <math>X_k(x)</math>, <math>U_k(u)</math>, <math>Y_k(y)</math>, and <math>V_k(v)</math> are rational functions whose coefficients for different {{mvar|k}}-values are listed in the referenced paper together with the <math>z_k</math> values that determine the subdomains. With higher degree polynomials in these rational functions the method can approximate the {{mvar|W}} function more accurately.

For example, when <math>-1/e\leq z\leq2.0082178115844727</math>, <math>W_0(z)</math> can be approximated to 24 bits of accuracy on 64-bit floating point values as <math>W_0(z)\approx X_1(x)=\frac{\sum_i^4P_ix^i}{\sum_i^3Q_ix^i}</math> where {{mvar|x}} is defined with the transformation above and the coefficients <math>P_i</math> and <math>Q_i</math> are given in the table below.

{| class="wikitable"
|+ Coefficients
|-
! <math>i</math> !! <math>P_i</math> !! <math>Q_i</math>
|-
| 0 || {{val|-0.9999999403954019}} || 1
|-
| 1 || {{val|0.0557300521617778}} || {{val|2.275906559863465}}
|-
| 2 || {{val|2.1269732491053173}} || {{val|1.367597013868904}}
|-
| 3 || {{val|0.8135112367835288}} || {{val|0.18615823452831623}}
|-
| 4 || {{val|0.01632488014607016}} || 0
|}
Fukushima also offers an approximation with 50 bits of accuracy on 64-bit floats that uses 8th- and 7th-degree polynomials.