1<html> 2<head> 3<meta http-equiv="Content-Type" content="text/html; charset=UTF-8"> 4<title>Lambert W function</title> 5<link rel="stylesheet" href="../math.css" type="text/css"> 6<meta name="generator" content="DocBook XSL Stylesheets V1.79.1"> 7<link rel="home" href="../index.html" title="Math Toolkit 2.12.0"> 8<link rel="up" href="../special.html" title="Chapter 8. Special Functions"> 9<link rel="prev" href="jacobi/jacobi_sn.html" title="Jacobi Elliptic Function sn"> 10<link rel="next" href="zetas.html" title="Zeta Functions"> 11</head> 12<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF"> 13<table cellpadding="2" width="100%"><tr> 14<td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../boost.png"></td> 15<td align="center"><a href="../../../../../index.html">Home</a></td> 16<td align="center"><a href="../../../../../libs/libraries.htm">Libraries</a></td> 17<td align="center"><a href="http://www.boost.org/users/people.html">People</a></td> 18<td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td> 19<td align="center"><a href="../../../../../more/index.htm">More</a></td> 20</tr></table> 21<hr> 22<div class="spirit-nav"> 23<a accesskey="p" href="jacobi/jacobi_sn.html"><img src="../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../special.html"><img src="../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../index.html"><img src="../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="zetas.html"><img src="../../../../../doc/src/images/next.png" alt="Next"></a> 24</div> 25<div class="section"> 26<div class="titlepage"><div><div><h2 class="title" style="clear: both"> 27<a name="math_toolkit.lambert_w"></a><a class="link" href="lambert_w.html" title="Lambert W function">Lambert <span class="emphasis"><em>W</em></span> 28 function</a> 29</h2></div></div></div> 30<h5> 31<a name="math_toolkit.lambert_w.h0"></a> 32 <span class="phrase"><a name="math_toolkit.lambert_w.synopsis"></a></span><a class="link" href="lambert_w.html#math_toolkit.lambert_w.synopsis">Synopsis</a> 33 </h5> 34<pre class="programlisting"><span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">boost</span><span class="special">/</span><span class="identifier">math</span><span class="special">/</span><span class="identifier">special_functions</span><span class="special">/</span><span class="identifier">lambert_w</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span> 35</pre> 36<pre class="programlisting"><span class="keyword">namespace</span> <span class="identifier">boost</span> <span class="special">{</span> <span class="keyword">namespace</span> <span class="identifier">math</span> <span class="special">{</span> 37 38 <span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">></span> 39 <a class="link" href="result_type.html" title="Calculation of the Type of the Result"><span class="emphasis"><em>calculated-result-type</em></span></a> <span class="identifier">lambert_w0</span><span class="special">(</span><span class="identifier">T</span> <span class="identifier">z</span><span class="special">);</span> <span class="comment">// W0 branch, default policy.</span> 40 <span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">></span> 41 <a class="link" href="result_type.html" title="Calculation of the Type of the Result"><span class="emphasis"><em>calculated-result-type</em></span></a> <span class="identifier">lambert_wm1</span><span class="special">(</span><span class="identifier">T</span> <span class="identifier">z</span><span class="special">);</span> <span class="comment">// W-1 branch, default policy.</span> 42 <span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">></span> 43 <a class="link" href="result_type.html" title="Calculation of the Type of the Result"><span class="emphasis"><em>calculated-result-type</em></span></a> <span class="identifier">lambert_w0_prime</span><span class="special">(</span><span class="identifier">T</span> <span class="identifier">z</span><span class="special">);</span> <span class="comment">// W0 branch 1st derivative.</span> 44 <span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">></span> 45 <a class="link" href="result_type.html" title="Calculation of the Type of the Result"><span class="emphasis"><em>calculated-result-type</em></span></a> <span class="identifier">lambert_wm1_prime</span><span class="special">(</span><span class="identifier">T</span> <span class="identifier">z</span><span class="special">);</span> <span class="comment">// W-1 branch 1st derivative.</span> 46 47 <span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">,</span> <span class="keyword">class</span> <a class="link" href="../policy.html" title="Chapter 21. Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">></span> 48 <a class="link" href="result_type.html" title="Calculation of the Type of the Result"><span class="emphasis"><em>calculated-result-type</em></span></a> <span class="identifier">lambert_w0</span><span class="special">(</span><span class="identifier">T</span> <span class="identifier">z</span><span class="special">,</span> <span class="keyword">const</span> <a class="link" href="../policy.html" title="Chapter 21. Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">&);</span> <span class="comment">// W0 with policy.</span> 49 <span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">,</span> <span class="keyword">class</span> <a class="link" href="../policy.html" title="Chapter 21. Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">></span> 50 <a class="link" href="result_type.html" title="Calculation of the Type of the Result"><span class="emphasis"><em>calculated-result-type</em></span></a> <span class="identifier">lambert_wm1</span><span class="special">(</span><span class="identifier">T</span> <span class="identifier">z</span><span class="special">,</span> <span class="keyword">const</span> <a class="link" href="../policy.html" title="Chapter 21. Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">&);</span> <span class="comment">// W-1 with policy.</span> 51 <span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">,</span> <span class="keyword">class</span> <a class="link" href="../policy.html" title="Chapter 21. Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">></span> 52 <a class="link" href="result_type.html" title="Calculation of the Type of the Result"><span class="emphasis"><em>calculated-result-type</em></span></a> <span class="identifier">lambert_w0_prime</span><span class="special">(</span><span class="identifier">T</span> <span class="identifier">z</span><span class="special">,</span> <span class="keyword">const</span> <a class="link" href="../policy.html" title="Chapter 21. Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">&);</span> <span class="comment">// W0 derivative with policy.</span> 53 <span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">,</span> <span class="keyword">class</span> <a class="link" href="../policy.html" title="Chapter 21. Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">></span> 54 <a class="link" href="result_type.html" title="Calculation of the Type of the Result"><span class="emphasis"><em>calculated-result-type</em></span></a> <span class="identifier">lambert_wm1_prime</span><span class="special">(</span><span class="identifier">T</span> <span class="identifier">z</span><span class="special">,</span> <span class="keyword">const</span> <a class="link" href="../policy.html" title="Chapter 21. Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">&);</span> <span class="comment">// W-1 derivative with policy.</span> 55 56 <span class="special">}</span> <span class="comment">// namespace boost</span> 57 <span class="special">}</span> <span class="comment">// namespace math</span> 58</pre> 59<h5> 60<a name="math_toolkit.lambert_w.h1"></a> 61 <span class="phrase"><a name="math_toolkit.lambert_w.description"></a></span><a class="link" href="lambert_w.html#math_toolkit.lambert_w.description">Description</a> 62 </h5> 63<p> 64 The <a href="http://en.wikipedia.org/wiki/Lambert_W_function" target="_top">Lambert W 65 function</a> is the solution of the equation <span class="emphasis"><em>W</em></span>(<span class="emphasis"><em>z</em></span>)<span class="emphasis"><em>e</em></span><sup><span class="emphasis"><em>W</em></span>(<span class="emphasis"><em>z</em></span>)</sup> = 66 <span class="emphasis"><em>z</em></span>. It is also called the Omega function, the inverse of 67 <span class="emphasis"><em>f</em></span>(<span class="emphasis"><em>W</em></span>) = <span class="emphasis"><em>We</em></span><sup><span class="emphasis"><em>W</em></span></sup>. 68 </p> 69<p> 70 On the interval [0, ∞), there is just one real solution. On the interval (-<span class="emphasis"><em>e</em></span><sup>-1</sup>, 71 0), there are two real solutions, generating two branches which we will denote 72 by <span class="emphasis"><em>W</em></span><sub>0</sub> and <span class="emphasis"><em>W</em></span><sub>-1</sub>. In Boost.Math, we call 73 these principal branches <code class="computeroutput"><span class="identifier">lambert_w0</span></code> 74 and <code class="computeroutput"><span class="identifier">lambert_wm1</span></code>; their derivatives 75 are labelled <code class="computeroutput"><span class="identifier">lambert_w0_prime</span></code> 76 and <code class="computeroutput"><span class="identifier">lambert_wm1_prime</span></code>. 77 </p> 78<div class="blockquote"><blockquote class="blockquote"><p> 79 <span class="inlinemediaobject"><img src="../../graphs/lambert_w_graph.svg" align="middle"></span> 80 81 </p></blockquote></div> 82<div class="blockquote"><blockquote class="blockquote"><p> 83 <span class="inlinemediaobject"><img src="../../graphs/lambert_w_graph_big_w.svg" align="middle"></span> 84 85 </p></blockquote></div> 86<div class="blockquote"><blockquote class="blockquote"><p> 87 <span class="inlinemediaobject"><img src="../../graphs/lambert_w0_prime_graph.svg" align="middle"></span> 88 89 </p></blockquote></div> 90<div class="blockquote"><blockquote class="blockquote"><p> 91 <span class="inlinemediaobject"><img src="../../graphs/lambert_wm1_prime_graph.svg" align="middle"></span> 92 93 </p></blockquote></div> 94<p> 95 There is a singularity where the branches meet at <span class="emphasis"><em>e</em></span><sup>-1</sup> ≅ <code class="literal">-0.367879</code>. 96 Approaching this point, the condition number of function evaluation tends to 97 infinity, and the only method of recovering high accuracy is use of higher 98 precision. 99 </p> 100<p> 101 This implementation computes the two real branches <span class="emphasis"><em>W</em></span><sub>0</sub> and 102 <span class="emphasis"><em>W</em></span><sub>-1</sub> 103with the functions <code class="computeroutput"><span class="identifier">lambert_w0</span></code> 104 and <code class="computeroutput"><span class="identifier">lambert_wm1</span></code>, and their 105 derivatives, <code class="computeroutput"><span class="identifier">lambert_w0_prime</span></code> 106 and <code class="computeroutput"><span class="identifier">lambert_wm1_prime</span></code>. Complex 107 arguments are not supported. 108 </p> 109<p> 110 The final <a class="link" href="../policy.html" title="Chapter 21. Policies: Controlling Precision, Error Handling etc">Policy</a> argument is optional and can 111 be used to control how the function deals with errors. Refer to <a class="link" href="../policy.html" title="Chapter 21. Policies: Controlling Precision, Error Handling etc">Policies</a> 112 for more details and see examples below. 113 </p> 114<h6> 115<a name="math_toolkit.lambert_w.h2"></a> 116 <span class="phrase"><a name="math_toolkit.lambert_w.applications"></a></span><a class="link" href="lambert_w.html#math_toolkit.lambert_w.applications">Applications 117 of the Lambert <span class="emphasis"><em>W</em></span> function</a> 118 </h6> 119<p> 120 The Lambert <span class="emphasis"><em>W</em></span> function has a myriad of applications. 121 <a href="http://www.apmaths.uwo.ca/~djeffrey/Offprints/W-adv-cm.pdf" target="_top">Corless 122 et al.</a> provide a summary of applications, from the mathematical, like 123 iterated exponentiation and asymptotic roots of trinomials, to the real-world, 124 such as the range of a jet plane, enzyme kinetics, water movement in soil, 125 epidemics, and diode current (an example replicated <a href="../../../example/lambert_w_diode.cpp" target="_top">here</a>). 126 Since the publication of their landmark paper, there have been many more applications, 127 and also many new implementations of the function, upon which this implementation 128 builds. 129 </p> 130<h5> 131<a name="math_toolkit.lambert_w.h3"></a> 132 <span class="phrase"><a name="math_toolkit.lambert_w.examples"></a></span><a class="link" href="lambert_w.html#math_toolkit.lambert_w.examples">Examples</a> 133 </h5> 134<p> 135 The most basic usage of the Lambert-<span class="emphasis"><em>W</em></span> function is demonstrated 136 below: 137 </p> 138<pre class="programlisting"><span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">boost</span><span class="special">/</span><span class="identifier">math</span><span class="special">/</span><span class="identifier">special_functions</span><span class="special">/</span><span class="identifier">lambert_w</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span> <span class="comment">// For lambert_w function.</span> 139 140<span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">lambert_w0</span><span class="special">;</span> 141<span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">lambert_wm1</span><span class="special">;</span> 142</pre> 143<pre class="programlisting"><span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span><span class="special">.</span><span class="identifier">precision</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="keyword">double</span><span class="special">>::</span><span class="identifier">max_digits10</span><span class="special">);</span> 144<span class="comment">// Show all potentially significant decimal digits,</span> 145<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">showpoint</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 146<span class="comment">// and show significant trailing zeros too.</span> 147 148<span class="keyword">double</span> <span class="identifier">z</span> <span class="special">=</span> <span class="number">10.</span><span class="special">;</span> 149<span class="keyword">double</span> <span class="identifier">r</span> <span class="special">=</span> <span class="identifier">lambert_w0</span><span class="special">(</span><span class="identifier">z</span><span class="special">);</span> <span class="comment">// Default policy for double.</span> 150<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"lambert_w0(z) = "</span> <span class="special"><<</span> <span class="identifier">r</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 151<span class="comment">// lambert_w0(z) = 1.7455280027406994</span> 152</pre> 153<p> 154 Other floating-point types can be used too, here <code class="computeroutput"><span class="keyword">float</span></code>, 155 including user-defined types like <a href="../../../../../libs/multiprecision/doc/html/index.html" target="_top">Boost.Multiprecision</a>. 156 It is convenient to use a function like <code class="computeroutput"><span class="identifier">show_value</span></code> 157 to display all (and only) potentially significant decimal digits, including 158 any significant trailing zeros, (<code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="identifier">T</span><span class="special">>::</span><span class="identifier">max_digits10</span></code>) for the type <code class="computeroutput"><span class="identifier">T</span></code>. 159 </p> 160<pre class="programlisting"><span class="keyword">float</span> <span class="identifier">z</span> <span class="special">=</span> <span class="number">10.F</span><span class="special">;</span> 161<span class="keyword">float</span> <span class="identifier">r</span><span class="special">;</span> 162<span class="identifier">r</span> <span class="special">=</span> <span class="identifier">lambert_w0</span><span class="special">(</span><span class="identifier">z</span><span class="special">);</span> <span class="comment">// Default policy digits10 = 7, digits2 = 24</span> 163<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"lambert_w0("</span><span class="special">;</span> 164<span class="identifier">show_value</span><span class="special">(</span><span class="identifier">z</span><span class="special">);</span> 165<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">") = "</span><span class="special">;</span> 166<span class="identifier">show_value</span><span class="special">(</span><span class="identifier">r</span><span class="special">);</span> 167<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> <span class="comment">// lambert_w0(10.0000000) = 1.74552798</span> 168</pre> 169<p> 170 Example of an integer argument to <code class="computeroutput"><span class="identifier">lambert_w0</span></code>, 171 showing that an <code class="computeroutput"><span class="keyword">int</span></code> literal is 172 correctly promoted to a <code class="computeroutput"><span class="keyword">double</span></code>. 173 </p> 174<pre class="programlisting"><span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span><span class="special">.</span><span class="identifier">precision</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="keyword">double</span><span class="special">>::</span><span class="identifier">max_digits10</span><span class="special">);</span> 175<span class="keyword">double</span> <span class="identifier">r</span> <span class="special">=</span> <span class="identifier">lambert_w0</span><span class="special">(</span><span class="number">10</span><span class="special">);</span> <span class="comment">// Pass an int argument "10" that should be promoted to double argument.</span> 176<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"lambert_w0(10) = "</span> <span class="special"><<</span> <span class="identifier">r</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> <span class="comment">// lambert_w0(10) = 1.7455280027406994</span> 177<span class="keyword">double</span> <span class="identifier">rp</span> <span class="special">=</span> <span class="identifier">lambert_w0</span><span class="special">(</span><span class="number">10</span><span class="special">);</span> 178<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"lambert_w0(10) = "</span> <span class="special"><<</span> <span class="identifier">rp</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 179<span class="comment">// lambert_w0(10) = 1.7455280027406994</span> 180<span class="keyword">auto</span> <span class="identifier">rr</span> <span class="special">=</span> <span class="identifier">lambert_w0</span><span class="special">(</span><span class="number">10</span><span class="special">);</span> <span class="comment">// C++11 needed.</span> 181<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"lambert_w0(10) = "</span> <span class="special"><<</span> <span class="identifier">rr</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 182<span class="comment">// lambert_w0(10) = 1.7455280027406994 too, showing that rr has been promoted to double.</span> 183</pre> 184<p> 185 Using <a href="../../../../../libs/multiprecision/doc/html/index.html" target="_top">Boost.Multiprecision</a> 186 types to get much higher precision is painless. 187 </p> 188<pre class="programlisting"><span class="identifier">cpp_dec_float_50</span> <span class="identifier">z</span><span class="special">(</span><span class="string">"10"</span><span class="special">);</span> 189<span class="comment">// Note construction using a decimal digit string "10",</span> 190<span class="comment">// NOT a floating-point double literal 10.</span> 191<span class="identifier">cpp_dec_float_50</span> <span class="identifier">r</span><span class="special">;</span> 192<span class="identifier">r</span> <span class="special">=</span> <span class="identifier">lambert_w0</span><span class="special">(</span><span class="identifier">z</span><span class="special">);</span> 193<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"lambert_w0("</span><span class="special">;</span> <span class="identifier">show_value</span><span class="special">(</span><span class="identifier">z</span><span class="special">);</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">") = "</span><span class="special">;</span> 194<span class="identifier">show_value</span><span class="special">(</span><span class="identifier">r</span><span class="special">);</span> 195<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 196<span class="comment">// lambert_w0(10.000000000000000000000000000000000000000000000000000000000000000000000000000000) =</span> 197<span class="comment">// 1.7455280027406993830743012648753899115352881290809413313533156980404446940000000</span> 198</pre> 199<div class="warning"><table border="0" summary="Warning"> 200<tr> 201<td rowspan="2" align="center" valign="top" width="25"><img alt="[Warning]" src="../../../../../doc/src/images/warning.png"></td> 202<th align="left">Warning</th> 203</tr> 204<tr><td align="left" valign="top"><p> 205 When using multiprecision, take very great care not to construct or assign 206 non-integers from <code class="computeroutput"><span class="keyword">double</span></code>, <code class="computeroutput"><span class="keyword">float</span></code> ... silently losing precision. Use 207 <code class="computeroutput"><span class="string">"1.2345678901234567890123456789"</span></code> 208 rather than <code class="computeroutput"><span class="number">1.2345678901234567890123456789</span></code>. 209 </p></td></tr> 210</table></div> 211<p> 212 Using multiprecision types, it is all too easy to get multiprecision precision 213 wrong! 214 </p> 215<pre class="programlisting"><span class="identifier">cpp_dec_float_50</span> <span class="identifier">z</span><span class="special">(</span><span class="number">0.7777777777777777777777777777777777777777777777777777777777777777777777777</span><span class="special">);</span> 216<span class="comment">// Compiler evaluates the nearest double-precision binary representation,</span> 217<span class="comment">// from the max_digits10 of the floating_point literal double 0.7777777777777777777777777777...,</span> 218<span class="comment">// so any extra digits in the multiprecision type</span> 219<span class="comment">// beyond max_digits10 (usually 17) are random and meaningless.</span> 220<span class="identifier">cpp_dec_float_50</span> <span class="identifier">r</span><span class="special">;</span> 221<span class="identifier">r</span> <span class="special">=</span> <span class="identifier">lambert_w0</span><span class="special">(</span><span class="identifier">z</span><span class="special">);</span> 222<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"lambert_w0("</span><span class="special">;</span> 223<span class="identifier">show_value</span><span class="special">(</span><span class="identifier">z</span><span class="special">);</span> 224<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">") = "</span><span class="special">;</span> <span class="identifier">show_value</span><span class="special">(</span><span class="identifier">r</span><span class="special">);</span> 225<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 226<span class="comment">// lambert_w0(0.77777777777777779011358916250173933804035186767578125000000000000000000000000000)</span> 227<span class="comment">// = 0.48086152073210493501934682309060873341910109230469724725005039758139532631901386</span> 228</pre> 229<div class="note"><table border="0" summary="Note"> 230<tr> 231<td rowspan="2" align="center" valign="top" width="25"><img alt="[Note]" src="../../../../../doc/src/images/note.png"></td> 232<th align="left">Note</th> 233</tr> 234<tr><td align="left" valign="top"><p> 235 See spurious non-seven decimal digits appearing after digit #17 in the argument 236 0.7777777777777777...! 237 </p></td></tr> 238</table></div> 239<p> 240 And similarly constructing from a literal <code class="computeroutput"><span class="keyword">double</span> 241 <span class="number">0.9</span></code>, with more random digits after digit 242 number 17. 243 </p> 244<pre class="programlisting"><span class="identifier">cpp_dec_float_50</span> <span class="identifier">z</span><span class="special">(</span><span class="number">0.9</span><span class="special">);</span> <span class="comment">// Construct from floating_point literal double 0.9.</span> 245<span class="identifier">cpp_dec_float_50</span> <span class="identifier">r</span><span class="special">;</span> 246<span class="identifier">r</span> <span class="special">=</span> <span class="identifier">lambert_w0</span><span class="special">(</span><span class="number">0.9</span><span class="special">);</span> 247<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"lambert_w0("</span><span class="special">;</span> 248<span class="identifier">show_value</span><span class="special">(</span><span class="identifier">z</span><span class="special">);</span> 249<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">") = "</span><span class="special">;</span> <span class="identifier">show_value</span><span class="special">(</span><span class="identifier">r</span><span class="special">);</span> 250<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 251<span class="comment">// lambert_w0(0.90000000000000002220446049250313080847263336181640625000000000000000000000000000)</span> 252<span class="comment">// = 0.52983296563343440510607251781038939952850341796875000000000000000000000000000000</span> 253<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"lambert_w0(0.9) = "</span> <span class="special"><<</span> <span class="identifier">lambert_w0</span><span class="special">(</span><span class="keyword">static_cast</span><span class="special"><</span><span class="keyword">double</span><span class="special">>(</span><span class="number">0.9</span><span class="special">))</span> 254<span class="comment">// lambert_w0(0.9)</span> 255<span class="comment">// = 0.52983296563343441</span> 256 <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 257</pre> 258<p> 259 Note how the <code class="computeroutput"><span class="identifier">cpp_float_dec_50</span></code> 260 result is only as correct as from a <code class="computeroutput"><span class="keyword">double</span> 261 <span class="special">=</span> <span class="number">0.9</span></code>. 262 </p> 263<p> 264 Now see the correct result for all 50 decimal digits constructing from a decimal 265 digit string "0.9": 266 </p> 267<pre class="programlisting"><span class="identifier">cpp_dec_float_50</span> <span class="identifier">z</span><span class="special">(</span><span class="string">"0.9"</span><span class="special">);</span> <span class="comment">// Construct from decimal digit string.</span> 268<span class="identifier">cpp_dec_float_50</span> <span class="identifier">r</span><span class="special">;</span> 269<span class="identifier">r</span> <span class="special">=</span> <span class="identifier">lambert_w0</span><span class="special">(</span><span class="identifier">z</span><span class="special">);</span> 270<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"lambert_w0("</span><span class="special">;</span> 271<span class="identifier">show_value</span><span class="special">(</span><span class="identifier">z</span><span class="special">);</span> 272<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">") = "</span><span class="special">;</span> <span class="identifier">show_value</span><span class="special">(</span><span class="identifier">r</span><span class="special">);</span> 273<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 274<span class="comment">// 0.90000000000000000000000000000000000000000000000000000000000000000000000000000000)</span> 275<span class="comment">// = 0.52983296563343441213336643954546304857788132269804249284012528304239956413801252</span> 276</pre> 277<p> 278 Note the expected zeros for all places up to 50 - and the correct Lambert 279 <span class="emphasis"><em>W</em></span> result! 280 </p> 281<p> 282 (It is just as easy to compute even higher precisions, at least to thousands 283 of decimal digits, but not shown here for brevity. See <a href="../../../example/lambert_w_simple_examples.cpp" target="_top">lambert_w_simple_examples.cpp</a> 284 for comparison of an evaluation at 1000 decimal digit precision with <a href="http://www.wolframalpha.com/" target="_top">Wolfram Alpha</a>). 285 </p> 286<p> 287 Policies can be used to control what action to take on errors: 288 </p> 289<pre class="programlisting"><span class="comment">// Define an error handling policy:</span> 290<span class="keyword">typedef</span> <span class="identifier">policy</span><span class="special"><</span> 291 <span class="identifier">domain_error</span><span class="special"><</span><span class="identifier">throw_on_error</span><span class="special">>,</span> 292 <span class="identifier">overflow_error</span><span class="special"><</span><span class="identifier">ignore_error</span><span class="special">></span> <span class="comment">// possibly unwise?</span> 293<span class="special">></span> <span class="identifier">my_throw_policy</span><span class="special">;</span> 294 295<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span><span class="special">.</span><span class="identifier">precision</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="keyword">double</span><span class="special">>::</span><span class="identifier">max_digits10</span><span class="special">);</span> 296<span class="comment">// Show all potentially significant decimal digits,</span> 297<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">showpoint</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 298<span class="comment">// and show significant trailing zeros too.</span> 299<span class="keyword">double</span> <span class="identifier">z</span> <span class="special">=</span> <span class="special">+</span><span class="number">1</span><span class="special">;</span> 300<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Lambert W ("</span> <span class="special"><<</span> <span class="identifier">z</span> <span class="special"><<</span> <span class="string">") = "</span> <span class="special"><<</span> <span class="identifier">lambert_w0</span><span class="special">(</span><span class="identifier">z</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 301<span class="comment">// Lambert W (1.0000000000000000) = 0.56714329040978384</span> 302<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"\nLambert W ("</span> <span class="special"><<</span> <span class="identifier">z</span> <span class="special"><<</span> <span class="string">", my_throw_policy()) = "</span> 303 <span class="special"><<</span> <span class="identifier">lambert_w0</span><span class="special">(</span><span class="identifier">z</span><span class="special">,</span> <span class="identifier">my_throw_policy</span><span class="special">())</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 304<span class="comment">// Lambert W (1.0000000000000000, my_throw_policy()) = 0.56714329040978384</span> 305</pre> 306<p> 307 An example error message: 308 </p> 309<pre class="programlisting"><span class="identifier">Error</span> <span class="identifier">in</span> <span class="identifier">function</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">lambert_wm1</span><span class="special"><</span><span class="identifier">RealType</span><span class="special">>(<</span><span class="identifier">RealType</span><span class="special">>):</span> 310<span class="identifier">Argument</span> <span class="identifier">z</span> <span class="special">=</span> <span class="number">1</span> <span class="identifier">is</span> <span class="identifier">out</span> <span class="identifier">of</span> <span class="identifier">range</span> <span class="special">(</span><span class="identifier">z</span> <span class="special"><=</span> <span class="number">0</span><span class="special">)</span> <span class="keyword">for</span> <span class="identifier">Lambert</span> <span class="identifier">W</span><span class="special">-</span><span class="number">1</span> <span class="identifier">branch</span><span class="special">!</span> <span class="special">(</span><span class="identifier">Try</span> <span class="identifier">Lambert</span> <span class="identifier">W0</span> <span class="identifier">branch</span><span class="special">?)</span> 311</pre> 312<p> 313 Showing an error reported if a value is passed to <code class="computeroutput"><span class="identifier">lambert_w0</span></code> 314 that is out of range, (and was probably meant to be passed to <code class="computeroutput"><span class="identifier">lambert_wm1</span></code> instead). 315 </p> 316<pre class="programlisting"><span class="keyword">double</span> <span class="identifier">z</span> <span class="special">=</span> <span class="special">+</span><span class="number">1.</span><span class="special">;</span> 317<span class="keyword">double</span> <span class="identifier">r</span> <span class="special">=</span> <span class="identifier">lambert_wm1</span><span class="special">(</span><span class="identifier">z</span><span class="special">);</span> 318<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"lambert_wm1(+1.) = "</span> <span class="special"><<</span> <span class="identifier">r</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 319</pre> 320<p> 321 The full source of these examples is at <a href="../../../example/lambert_w_simple_examples.cpp" target="_top">lambert_w_simple_examples.cpp</a> 322 </p> 323<h6> 324<a name="math_toolkit.lambert_w.h4"></a> 325 <span class="phrase"><a name="math_toolkit.lambert_w.diode_resistance"></a></span><a class="link" href="lambert_w.html#math_toolkit.lambert_w.diode_resistance">Diode 326 Resistance Example</a> 327 </h6> 328<p> 329 A typical example of a practical application is estimating the current flow 330 through a diode with series resistance from a paper by Banwell and Jayakumar. 331 </p> 332<p> 333 Having the Lambert <span class="emphasis"><em>W</em></span> function available makes it simple 334 to reproduce the plot in their paper (Fig 2) comparing estimates using with 335 Lambert <span class="emphasis"><em>W</em></span> function and some actual measurements. The colored 336 curves show the effect of various series resistance on the current compared 337 to an extrapolated line in grey with no internal (or external) resistance. 338 </p> 339<p> 340 Two formulae relating the diode current and effect of series resistance can 341 be combined, but yield an otherwise intractable equation relating the current 342 versus voltage with a varying series resistance. This was reformulated as a 343 generalized equation in terms of the Lambert W function: 344 </p> 345<p> 346 Banwell and Jakaumar equation 5 347 </p> 348<div class="blockquote"><blockquote class="blockquote"><p> 349 <span class="serif_italic">I(V) = μ V<sub>T</sub>/ R <sub>S</sub> ․ W<sub>0</sub>(I<sub>0</sub> R<sub>S</sub> / (μ V<sub>T</sub>))</span> 350 </p></blockquote></div> 351<p> 352 Using these variables 353 </p> 354<pre class="programlisting"><span class="keyword">double</span> <span class="identifier">nu</span> <span class="special">=</span> <span class="number">1.0</span><span class="special">;</span> <span class="comment">// Assumed ideal.</span> 355<span class="keyword">double</span> <span class="identifier">vt</span> <span class="special">=</span> <span class="identifier">v_thermal</span><span class="special">(</span><span class="number">25</span><span class="special">);</span> <span class="comment">// v thermal, Shockley equation, expect about 25 mV at room temperature.</span> 356<span class="keyword">double</span> <span class="identifier">boltzmann_k</span> <span class="special">=</span> <span class="number">1.38e-23</span><span class="special">;</span> <span class="comment">// joules/kelvin</span> 357<span class="keyword">double</span> <span class="identifier">temp</span> <span class="special">=</span> <span class="number">273</span> <span class="special">+</span> <span class="number">25</span><span class="special">;</span> 358<span class="keyword">double</span> <span class="identifier">charge_q</span> <span class="special">=</span> <span class="number">1.6e-19</span><span class="special">;</span> <span class="comment">// column</span> 359<span class="identifier">vt</span> <span class="special">=</span> <span class="identifier">boltzmann_k</span> <span class="special">*</span> <span class="identifier">temp</span> <span class="special">/</span> <span class="identifier">charge_q</span><span class="special">;</span> 360<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"V thermal "</span> <span class="special"><<</span> <span class="identifier">vt</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> <span class="comment">// V thermal 0.0257025 = 25 mV</span> 361<span class="keyword">double</span> <span class="identifier">rsat</span> <span class="special">=</span> <span class="number">0.</span><span class="special">;</span> 362<span class="keyword">double</span> <span class="identifier">isat</span> <span class="special">=</span> <span class="number">25.e-15</span><span class="special">;</span> <span class="comment">// 25 fA;</span> 363<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Isat = "</span> <span class="special"><<</span> <span class="identifier">isat</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 364<span class="keyword">double</span> <span class="identifier">re</span> <span class="special">=</span> <span class="number">0.3</span><span class="special">;</span> <span class="comment">// Estimated from slope of straight section of graph (equation 6).</span> 365<span class="keyword">double</span> <span class="identifier">v</span> <span class="special">=</span> <span class="number">0.9</span><span class="special">;</span> 366<span class="keyword">double</span> <span class="identifier">icalc</span> <span class="special">=</span> <span class="identifier">iv</span><span class="special">(</span><span class="identifier">v</span><span class="special">,</span> <span class="identifier">vt</span><span class="special">,</span> <span class="number">249.</span><span class="special">,</span> <span class="identifier">re</span><span class="special">,</span> <span class="identifier">isat</span><span class="special">);</span> 367<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"voltage = "</span> <span class="special"><<</span> <span class="identifier">v</span> <span class="special"><<</span> <span class="string">", current = "</span> <span class="special"><<</span> <span class="identifier">icalc</span> <span class="special"><<</span> <span class="string">", "</span> <span class="special"><<</span> <span class="identifier">log</span><span class="special">(</span><span class="identifier">icalc</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> <span class="comment">// voltage = 0.9, current = 0.00108485, -6.82631</span> 368</pre> 369<p> 370 the formulas can be rendered in C++ 371 </p> 372<pre class="programlisting"><span class="keyword">double</span> <span class="identifier">iv</span><span class="special">(</span><span class="keyword">double</span> <span class="identifier">v</span><span class="special">,</span> <span class="keyword">double</span> <span class="identifier">vt</span><span class="special">,</span> <span class="keyword">double</span> <span class="identifier">rsat</span><span class="special">,</span> <span class="keyword">double</span> <span class="identifier">re</span><span class="special">,</span> <span class="keyword">double</span> <span class="identifier">isat</span><span class="special">,</span> <span class="keyword">double</span> <span class="identifier">nu</span> <span class="special">=</span> <span class="number">1.</span><span class="special">)</span> 373<span class="special">{</span> 374 <span class="comment">// V thermal 0.0257025 = 25 mV</span> 375 <span class="comment">// was double i = (nu * vt/r) * lambert_w((i0 * r) / (nu * vt)); equ 5.</span> 376 377 <span class="identifier">rsat</span> <span class="special">=</span> <span class="identifier">rsat</span> <span class="special">+</span> <span class="identifier">re</span><span class="special">;</span> 378 <span class="keyword">double</span> <span class="identifier">i</span> <span class="special">=</span> <span class="identifier">nu</span> <span class="special">*</span> <span class="identifier">vt</span> <span class="special">/</span> <span class="identifier">rsat</span><span class="special">;</span> 379 <span class="comment">// std::cout << "nu * vt / rsat = " << i << std::endl; // 0.000103223</span> 380 381 <span class="keyword">double</span> <span class="identifier">x</span> <span class="special">=</span> <span class="identifier">isat</span> <span class="special">*</span> <span class="identifier">rsat</span> <span class="special">/</span> <span class="special">(</span><span class="identifier">nu</span> <span class="special">*</span> <span class="identifier">vt</span><span class="special">);</span> 382<span class="comment">// std::cout << "isat * rsat / (nu * vt) = " << x << std::endl;</span> 383 384 <span class="keyword">double</span> <span class="identifier">eterm</span> <span class="special">=</span> <span class="special">(</span><span class="identifier">v</span> <span class="special">+</span> <span class="identifier">isat</span> <span class="special">*</span> <span class="identifier">rsat</span><span class="special">)</span> <span class="special">/</span> <span class="special">(</span><span class="identifier">nu</span> <span class="special">*</span> <span class="identifier">vt</span><span class="special">);</span> 385 <span class="comment">// std::cout << "(v + isat * rsat) / (nu * vt) = " << eterm << std::endl;</span> 386 387 <span class="keyword">double</span> <span class="identifier">e</span> <span class="special">=</span> <span class="identifier">exp</span><span class="special">(</span><span class="identifier">eterm</span><span class="special">);</span> 388<span class="comment">// std::cout << "exp(eterm) = " << e << std::endl;</span> 389 390 <span class="keyword">double</span> <span class="identifier">w0</span> <span class="special">=</span> <span class="identifier">lambert_w0</span><span class="special">(</span><span class="identifier">x</span> <span class="special">*</span> <span class="identifier">e</span><span class="special">);</span> 391<span class="comment">// std::cout << "w0 = " << w0 << std::endl;</span> 392 <span class="keyword">return</span> <span class="identifier">i</span> <span class="special">*</span> <span class="identifier">w0</span> <span class="special">-</span> <span class="identifier">isat</span><span class="special">;</span> 393<span class="special">}</span> <span class="comment">// double iv</span> 394</pre> 395<p> 396 to reproduce their Fig 2: 397 </p> 398<div class="blockquote"><blockquote class="blockquote"><p> 399 <span class="inlinemediaobject"><img src="../../graphs/diode_iv_plot.svg" align="middle"></span> 400 401 </p></blockquote></div> 402<p> 403 The plotted points for no external series resistance (derived from their published 404 plot as the raw data are not publicly available) are used to extrapolate back 405 to estimate the intrinsic emitter resistance as 0.3 ohm. The effect of external 406 series resistance is visible when the colored lines start to curve away from 407 the straight line as voltage increases. 408 </p> 409<p> 410 See <a href="../../../example/lambert_w_diode.cpp" target="_top">lambert_w_diode.cpp</a> 411 and <a href="../../../example/lambert_w_diode_graph.cpp" target="_top">lambert_w_diode_graph.cpp</a> 412 for details of the calculation. 413 </p> 414<h6> 415<a name="math_toolkit.lambert_w.h5"></a> 416 <span class="phrase"><a name="math_toolkit.lambert_w.implementations"></a></span><a class="link" href="lambert_w.html#math_toolkit.lambert_w.implementations">Existing 417 implementations</a> 418 </h6> 419<p> 420 The principal value of the Lambert <span class="emphasis"><em>W</em></span> function is implemented 421 in the <a href="http://mathworld.wolfram.com/LambertW-Function.html" target="_top">Wolfram 422 Language</a> as <code class="computeroutput"><span class="identifier">ProductLog</span><span class="special">[</span><span class="identifier">k</span><span class="special">,</span> 423 <span class="identifier">z</span><span class="special">]</span></code>, 424 where <code class="computeroutput"><span class="identifier">k</span></code> is the branch. 425 </p> 426<p> 427 The symbolic algebra program <a href="https://www.maplesoft.com" target="_top">Maple</a> 428 also computes Lambert <span class="emphasis"><em>W</em></span> to an arbitrary precision. 429 </p> 430<h5> 431<a name="math_toolkit.lambert_w.h6"></a> 432 <span class="phrase"><a name="math_toolkit.lambert_w.precision"></a></span><a class="link" href="lambert_w.html#math_toolkit.lambert_w.precision">Controlling 433 the compromise between Precision and Speed</a> 434 </h5> 435<h6> 436<a name="math_toolkit.lambert_w.h7"></a> 437 <span class="phrase"><a name="math_toolkit.lambert_w.small_floats"></a></span><a class="link" href="lambert_w.html#math_toolkit.lambert_w.small_floats">Floating-point 438 types <code class="computeroutput"><span class="keyword">double</span></code> and <code class="computeroutput"><span class="keyword">float</span></code></a> 439 </h6> 440<p> 441 This implementation provides good precision and excellent speed for __fundamental 442 <code class="computeroutput"><span class="keyword">float</span></code> and <code class="computeroutput"><span class="keyword">double</span></code>. 443 </p> 444<p> 445 All the functions usually return values within a few <a href="http://en.wikipedia.org/wiki/Unit_in_the_last_place" target="_top">Unit 446 in the last place (ULP)</a> for the floating-point type, except for very 447 small arguments very near zero, and for arguments very close to the singularity 448 at the branch point. 449 </p> 450<p> 451 By default, this implementation provides the best possible speed. Very slightly 452 average higher precision and less bias might be obtained by adding a <a class="link" href="roots_deriv.html#math_toolkit.roots_deriv.halley">Halley</a> step refinement, but 453 at the cost of more than doubling the runtime. 454 </p> 455<h6> 456<a name="math_toolkit.lambert_w.h8"></a> 457 <span class="phrase"><a name="math_toolkit.lambert_w.big_floats"></a></span><a class="link" href="lambert_w.html#math_toolkit.lambert_w.big_floats">Floating-point 458 types larger than double</a> 459 </h6> 460<p> 461 For floating-point types with precision greater than <code class="computeroutput"><span class="keyword">double</span></code> 462 and <code class="computeroutput"><span class="keyword">float</span></code> <a href="http://en.cppreference.com/w/cpp/language/types" target="_top">fundamental 463 (built-in) types</a>, a <code class="computeroutput"><span class="keyword">double</span></code> 464 evaluation is used as a first approximation followed by Halley refinement, 465 using a single step where it can be predicted that this will be sufficient, 466 and only using <a class="link" href="roots_deriv.html#math_toolkit.roots_deriv.halley">Halley</a> 467 iteration when necessary. Higher precision types are always going to be <span class="bold"><strong>very, very much slower</strong></span>. 468 </p> 469<p> 470 The 'best' evaluation (the nearest <a href="http://en.wikipedia.org/wiki/Floating_point#Representable_numbers.2C_conversion_and_rounding" target="_top">representable</a>) 471 can be achieved by <code class="computeroutput"><span class="keyword">static_cast</span></code>ing 472 from a higher precision type, typically a <a href="../../../../../libs/multiprecision/doc/html/index.html" target="_top">Boost.Multiprecision</a> 473 type like <code class="computeroutput"><span class="identifier">cpp_bin_float_50</span></code>, 474 but at the cost of increasing run-time 100-fold; this has been used here to 475 provide some of our reference values for testing. 476 </p> 477<p> 478 For example, we get a reference value using a high precision type, for example; 479 </p> 480<pre class="programlisting"><span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">multiprecision</span><span class="special">::</span><span class="identifier">cpp_bin_float_50</span><span class="special">;</span> 481</pre> 482<p> 483 that uses Halley iteration to refine until it is as precise as possible for 484 this <code class="computeroutput"><span class="identifier">cpp_bin_float_50</span></code> type. 485 </p> 486<p> 487 As a further check we can compare this with a <a href="http://www.wolframalpha.com/" target="_top">Wolfram 488 Alpha</a> computation using command <code class="literal">N[ProductLog[10.], 50]</code> 489 to get 50 decimal digits and similarly <code class="literal">N[ProductLog[10.], 17]</code> 490 to get the nearest representable for 64-bit <code class="computeroutput"><span class="keyword">double</span></code> 491 precision. 492 </p> 493<pre class="programlisting"> <span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">multiprecision</span><span class="special">::</span><span class="identifier">cpp_bin_float_50</span><span class="special">;</span> 494 <span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">float_distance</span><span class="special">;</span> 495 496 <span class="identifier">cpp_bin_float_50</span> <span class="identifier">z</span><span class="special">(</span><span class="string">"10."</span><span class="special">);</span> <span class="comment">// Note use a decimal digit string, not a double 10.</span> 497 <span class="identifier">cpp_bin_float_50</span> <span class="identifier">r</span><span class="special">;</span> 498 <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span><span class="special">.</span><span class="identifier">precision</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="identifier">cpp_bin_float_50</span><span class="special">>::</span><span class="identifier">digits10</span><span class="special">);</span> 499 500 <span class="identifier">r</span> <span class="special">=</span> <span class="identifier">lambert_w0</span><span class="special">(</span><span class="identifier">z</span><span class="special">);</span> <span class="comment">// Default policy.</span> 501 <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"lambert_w0(z) cpp_bin_float_50 = "</span> <span class="special"><<</span> <span class="identifier">r</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 502 <span class="comment">//lambert_w0(z) cpp_bin_float_50 = 1.7455280027406993830743012648753899115352881290809</span> 503 <span class="comment">// [N[productlog[10], 50]] == 1.7455280027406993830743012648753899115352881290809</span> 504 <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span><span class="special">.</span><span class="identifier">precision</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="keyword">double</span><span class="special">>::</span><span class="identifier">max_digits10</span><span class="special">);</span> 505 <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"lambert_w0(z) static_cast from cpp_bin_float_50 = "</span> 506 <span class="special"><<</span> <span class="keyword">static_cast</span><span class="special"><</span><span class="keyword">double</span><span class="special">>(</span><span class="identifier">r</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 507 <span class="comment">// double lambert_w0(z) static_cast from cpp_bin_float_50 = 1.7455280027406994</span> 508 <span class="comment">// [N[productlog[10], 17]] == 1.7455280027406994</span> 509<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"bits different from Wolfram = "</span> 510 <span class="special"><<</span> <span class="keyword">static_cast</span><span class="special"><</span><span class="keyword">int</span><span class="special">>(</span><span class="identifier">float_distance</span><span class="special">(</span><span class="keyword">static_cast</span><span class="special"><</span><span class="keyword">double</span><span class="special">>(</span><span class="identifier">r</span><span class="special">),</span> <span class="number">1.7455280027406994</span><span class="special">))</span> 511 <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> <span class="comment">// 0</span> 512</pre> 513<p> 514 giving us the same nearest representable using 64-bit <code class="computeroutput"><span class="keyword">double</span></code> 515 as <code class="computeroutput"><span class="number">1.7455280027406994</span></code>. 516 </p> 517<p> 518 However, the rational polynomial and Fukushima Schroder approximations are 519 so good for type <code class="computeroutput"><span class="keyword">float</span></code> and <code class="computeroutput"><span class="keyword">double</span></code> that negligible improvement is gained 520 from a <code class="computeroutput"><span class="keyword">double</span></code> Halley step. 521 </p> 522<p> 523 This is shown with <a href="../../../example/lambert_w_precision_example.cpp" target="_top">lambert_w_precision_example.cpp</a> 524 for Lambert <span class="emphasis"><em>W</em></span><sub>0</sub>: 525 </p> 526<pre class="programlisting"><span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">lambert_w_detail</span><span class="special">::</span><span class="identifier">lambert_w_halley_step</span><span class="special">;</span> 527<span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">epsilon_difference</span><span class="special">;</span> 528<span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">relative_difference</span><span class="special">;</span> 529 530<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">showpoint</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> <span class="comment">// and show any significant trailing zeros too.</span> 531<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span><span class="special">.</span><span class="identifier">precision</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="keyword">double</span><span class="special">>::</span><span class="identifier">max_digits10</span><span class="special">);</span> <span class="comment">// 17 decimal digits for double.</span> 532 533<span class="identifier">cpp_bin_float_50</span> <span class="identifier">z50</span><span class="special">(</span><span class="string">"1.23"</span><span class="special">);</span> <span class="comment">// Note: use a decimal digit string, not a double 1.23!</span> 534<span class="keyword">double</span> <span class="identifier">z</span> <span class="special">=</span> <span class="keyword">static_cast</span><span class="special"><</span><span class="keyword">double</span><span class="special">>(</span><span class="identifier">z50</span><span class="special">);</span> 535<span class="identifier">cpp_bin_float_50</span> <span class="identifier">w50</span><span class="special">;</span> 536<span class="identifier">w50</span> <span class="special">=</span> <span class="identifier">lambert_w0</span><span class="special">(</span><span class="identifier">z50</span><span class="special">);</span> 537<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span><span class="special">.</span><span class="identifier">precision</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="identifier">cpp_bin_float_50</span><span class="special">>::</span><span class="identifier">max_digits10</span><span class="special">);</span> <span class="comment">// 50 decimal digits.</span> 538<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Reference Lambert W ("</span> <span class="special"><<</span> <span class="identifier">z</span> <span class="special"><<</span> <span class="string">") =\n "</span> 539 <span class="special"><<</span> <span class="identifier">w50</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 540<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span><span class="special">.</span><span class="identifier">precision</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="keyword">double</span><span class="special">>::</span><span class="identifier">max_digits10</span><span class="special">);</span> <span class="comment">// 17 decimal digits for double.</span> 541<span class="keyword">double</span> <span class="identifier">wr</span> <span class="special">=</span> <span class="keyword">static_cast</span><span class="special"><</span><span class="keyword">double</span><span class="special">>(</span><span class="identifier">w50</span><span class="special">);</span> 542<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Reference Lambert W ("</span> <span class="special"><<</span> <span class="identifier">z</span> <span class="special"><<</span> <span class="string">") = "</span> <span class="special"><<</span> <span class="identifier">wr</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 543 544<span class="keyword">double</span> <span class="identifier">w</span> <span class="special">=</span> <span class="identifier">lambert_w0</span><span class="special">(</span><span class="identifier">z</span><span class="special">);</span> 545<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Rat/poly Lambert W ("</span> <span class="special"><<</span> <span class="identifier">z</span> <span class="special"><<</span> <span class="string">") = "</span> <span class="special"><<</span> <span class="identifier">lambert_w0</span><span class="special">(</span><span class="identifier">z</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 546<span class="comment">// Add a Halley step to the value obtained from rational polynomial approximation.</span> 547<span class="keyword">double</span> <span class="identifier">ww</span> <span class="special">=</span> <span class="identifier">lambert_w_halley_step</span><span class="special">(</span><span class="identifier">lambert_w0</span><span class="special">(</span><span class="identifier">z</span><span class="special">),</span> <span class="identifier">z</span><span class="special">);</span> 548<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Halley Step Lambert W ("</span> <span class="special"><<</span> <span class="identifier">z</span> <span class="special"><<</span> <span class="string">") = "</span> <span class="special"><<</span> <span class="identifier">lambert_w_halley_step</span><span class="special">(</span><span class="identifier">lambert_w0</span><span class="special">(</span><span class="identifier">z</span><span class="special">),</span> <span class="identifier">z</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 549 550<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"absolute difference from Halley step = "</span> <span class="special"><<</span> <span class="identifier">w</span> <span class="special">-</span> <span class="identifier">ww</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 551<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"relative difference from Halley step = "</span> <span class="special"><<</span> <span class="identifier">relative_difference</span><span class="special">(</span><span class="identifier">w</span><span class="special">,</span> <span class="identifier">ww</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 552<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"epsilon difference from Halley step = "</span> <span class="special"><<</span> <span class="identifier">epsilon_difference</span><span class="special">(</span><span class="identifier">w</span><span class="special">,</span> <span class="identifier">ww</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 553<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"epsilon for float = "</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="keyword">double</span><span class="special">>::</span><span class="identifier">epsilon</span><span class="special">()</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 554<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"bits different from Halley step = "</span> <span class="special"><<</span> <span class="keyword">static_cast</span><span class="special"><</span><span class="keyword">int</span><span class="special">>(</span><span class="identifier">float_distance</span><span class="special">(</span><span class="identifier">w</span><span class="special">,</span> <span class="identifier">ww</span><span class="special">))</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 555</pre> 556<p> 557 with this output: 558 </p> 559<pre class="programlisting"><span class="identifier">Reference</span> <span class="identifier">Lambert</span> <span class="identifier">W</span> <span class="special">(</span><span class="number">1.2299999999999999822364316059974953532218933105468750</span><span class="special">)</span> <span class="special">=</span> 560<span class="number">0.64520356959320237759035605255334853830173300262666480</span> 561<span class="identifier">Reference</span> <span class="identifier">Lambert</span> <span class="identifier">W</span> <span class="special">(</span><span class="number">1.2300000000000000</span><span class="special">)</span> <span class="special">=</span> <span class="number">0.64520356959320235</span> 562<span class="identifier">Rat</span><span class="special">/</span><span class="identifier">poly</span> <span class="identifier">Lambert</span> <span class="identifier">W</span> <span class="special">(</span><span class="number">1.2300000000000000</span><span class="special">)</span> <span class="special">=</span> <span class="number">0.64520356959320224</span> 563<span class="identifier">Halley</span> <span class="identifier">Step</span> <span class="identifier">Lambert</span> <span class="identifier">W</span> <span class="special">(</span><span class="number">1.2300000000000000</span><span class="special">)</span> <span class="special">=</span> <span class="number">0.64520356959320235</span> 564<span class="identifier">absolute</span> <span class="identifier">difference</span> <span class="identifier">from</span> <span class="identifier">Halley</span> <span class="identifier">step</span> <span class="special">=</span> <span class="special">-</span><span class="number">1.1102230246251565e-16</span> 565<span class="identifier">relative</span> <span class="identifier">difference</span> <span class="identifier">from</span> <span class="identifier">Halley</span> <span class="identifier">step</span> <span class="special">=</span> <span class="number">1.7207329236029286e-16</span> 566<span class="identifier">epsilon</span> <span class="identifier">difference</span> <span class="identifier">from</span> <span class="identifier">Halley</span> <span class="identifier">step</span> <span class="special">=</span> <span class="number">0.77494921535422934</span> 567<span class="identifier">epsilon</span> <span class="keyword">for</span> <span class="keyword">float</span> <span class="special">=</span> <span class="number">2.2204460492503131e-16</span> 568<span class="identifier">bits</span> <span class="identifier">different</span> <span class="identifier">from</span> <span class="identifier">Halley</span> <span class="identifier">step</span> <span class="special">=</span> <span class="number">1</span> 569</pre> 570<p> 571 and then for <span class="emphasis"><em>W</em></span><sub>-1</sub>: 572 </p> 573<pre class="programlisting"><span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">lambert_w_detail</span><span class="special">::</span><span class="identifier">lambert_w_halley_step</span><span class="special">;</span> 574<span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">epsilon_difference</span><span class="special">;</span> 575<span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">relative_difference</span><span class="special">;</span> 576 577<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">showpoint</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> <span class="comment">// and show any significant trailing zeros too.</span> 578<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span><span class="special">.</span><span class="identifier">precision</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="keyword">double</span><span class="special">>::</span><span class="identifier">max_digits10</span><span class="special">);</span> <span class="comment">// 17 decimal digits for double.</span> 579 580<span class="identifier">cpp_bin_float_50</span> <span class="identifier">z50</span><span class="special">(</span><span class="string">"-0.123"</span><span class="special">);</span> <span class="comment">// Note: use a decimal digit string, not a double -1.234!</span> 581<span class="keyword">double</span> <span class="identifier">z</span> <span class="special">=</span> <span class="keyword">static_cast</span><span class="special"><</span><span class="keyword">double</span><span class="special">>(</span><span class="identifier">z50</span><span class="special">);</span> 582<span class="identifier">cpp_bin_float_50</span> <span class="identifier">wm1_50</span><span class="special">;</span> 583<span class="identifier">wm1_50</span> <span class="special">=</span> <span class="identifier">lambert_wm1</span><span class="special">(</span><span class="identifier">z50</span><span class="special">);</span> 584<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span><span class="special">.</span><span class="identifier">precision</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="identifier">cpp_bin_float_50</span><span class="special">>::</span><span class="identifier">max_digits10</span><span class="special">);</span> <span class="comment">// 50 decimal digits.</span> 585<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Reference Lambert W-1 ("</span> <span class="special"><<</span> <span class="identifier">z</span> <span class="special"><<</span> <span class="string">") =\n "</span> 586 <span class="special"><<</span> <span class="identifier">wm1_50</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 587<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span><span class="special">.</span><span class="identifier">precision</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="keyword">double</span><span class="special">>::</span><span class="identifier">max_digits10</span><span class="special">);</span> <span class="comment">// 17 decimal digits for double.</span> 588<span class="keyword">double</span> <span class="identifier">wr</span> <span class="special">=</span> <span class="keyword">static_cast</span><span class="special"><</span><span class="keyword">double</span><span class="special">>(</span><span class="identifier">wm1_50</span><span class="special">);</span> 589<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Reference Lambert W-1 ("</span> <span class="special"><<</span> <span class="identifier">z</span> <span class="special"><<</span> <span class="string">") = "</span> <span class="special"><<</span> <span class="identifier">wr</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 590 591<span class="keyword">double</span> <span class="identifier">w</span> <span class="special">=</span> <span class="identifier">lambert_wm1</span><span class="special">(</span><span class="identifier">z</span><span class="special">);</span> 592<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Rat/poly Lambert W-1 ("</span> <span class="special"><<</span> <span class="identifier">z</span> <span class="special"><<</span> <span class="string">") = "</span> <span class="special"><<</span> <span class="identifier">lambert_wm1</span><span class="special">(</span><span class="identifier">z</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 593<span class="comment">// Add a Halley step to the value obtained from rational polynomial approximation.</span> 594<span class="keyword">double</span> <span class="identifier">ww</span> <span class="special">=</span> <span class="identifier">lambert_w_halley_step</span><span class="special">(</span><span class="identifier">lambert_wm1</span><span class="special">(</span><span class="identifier">z</span><span class="special">),</span> <span class="identifier">z</span><span class="special">);</span> 595<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Halley Step Lambert W ("</span> <span class="special"><<</span> <span class="identifier">z</span> <span class="special"><<</span> <span class="string">") = "</span> <span class="special"><<</span> <span class="identifier">lambert_w_halley_step</span><span class="special">(</span><span class="identifier">lambert_wm1</span><span class="special">(</span><span class="identifier">z</span><span class="special">),</span> <span class="identifier">z</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 596 597<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"absolute difference from Halley step = "</span> <span class="special"><<</span> <span class="identifier">w</span> <span class="special">-</span> <span class="identifier">ww</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 598<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"relative difference from Halley step = "</span> <span class="special"><<</span> <span class="identifier">relative_difference</span><span class="special">(</span><span class="identifier">w</span><span class="special">,</span> <span class="identifier">ww</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 599<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"epsilon difference from Halley step = "</span> <span class="special"><<</span> <span class="identifier">epsilon_difference</span><span class="special">(</span><span class="identifier">w</span><span class="special">,</span> <span class="identifier">ww</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 600<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"epsilon for float = "</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="keyword">double</span><span class="special">>::</span><span class="identifier">epsilon</span><span class="special">()</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 601<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"bits different from Halley step = "</span> <span class="special"><<</span> <span class="keyword">static_cast</span><span class="special"><</span><span class="keyword">int</span><span class="special">>(</span><span class="identifier">float_distance</span><span class="special">(</span><span class="identifier">w</span><span class="special">,</span> <span class="identifier">ww</span><span class="special">))</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 602</pre> 603<p> 604 with this output: 605 </p> 606<pre class="programlisting"><span class="identifier">Reference</span> <span class="identifier">Lambert</span> <span class="identifier">W</span><span class="special">-</span><span class="number">1</span> <span class="special">(-</span><span class="number">0.12299999999999999822364316059974953532218933105468750</span><span class="special">)</span> <span class="special">=</span> 607<span class="special">-</span><span class="number">3.2849102557740360179084675531714935199110302996513384</span> 608<span class="identifier">Reference</span> <span class="identifier">Lambert</span> <span class="identifier">W</span><span class="special">-</span><span class="number">1</span> <span class="special">(-</span><span class="number">0.12300000000000000</span><span class="special">)</span> <span class="special">=</span> <span class="special">-</span><span class="number">3.2849102557740362</span> 609<span class="identifier">Rat</span><span class="special">/</span><span class="identifier">poly</span> <span class="identifier">Lambert</span> <span class="identifier">W</span><span class="special">-</span><span class="number">1</span> <span class="special">(-</span><span class="number">0.12300000000000000</span><span class="special">)</span> <span class="special">=</span> <span class="special">-</span><span class="number">3.2849102557740357</span> 610<span class="identifier">Halley</span> <span class="identifier">Step</span> <span class="identifier">Lambert</span> <span class="identifier">W</span> <span class="special">(-</span><span class="number">0.12300000000000000</span><span class="special">)</span> <span class="special">=</span> <span class="special">-</span><span class="number">3.2849102557740362</span> 611<span class="identifier">absolute</span> <span class="identifier">difference</span> <span class="identifier">from</span> <span class="identifier">Halley</span> <span class="identifier">step</span> <span class="special">=</span> <span class="number">4.4408920985006262e-16</span> 612<span class="identifier">relative</span> <span class="identifier">difference</span> <span class="identifier">from</span> <span class="identifier">Halley</span> <span class="identifier">step</span> <span class="special">=</span> <span class="number">1.3519066740696092e-16</span> 613<span class="identifier">epsilon</span> <span class="identifier">difference</span> <span class="identifier">from</span> <span class="identifier">Halley</span> <span class="identifier">step</span> <span class="special">=</span> <span class="number">0.60884463935795785</span> 614<span class="identifier">epsilon</span> <span class="keyword">for</span> <span class="keyword">float</span> <span class="special">=</span> <span class="number">2.2204460492503131e-16</span> 615<span class="identifier">bits</span> <span class="identifier">different</span> <span class="identifier">from</span> <span class="identifier">Halley</span> <span class="identifier">step</span> <span class="special">=</span> <span class="special">-</span><span class="number">1</span> 616</pre> 617<h6> 618<a name="math_toolkit.lambert_w.h9"></a> 619 <span class="phrase"><a name="math_toolkit.lambert_w.differences_distribution"></a></span><a class="link" href="lambert_w.html#math_toolkit.lambert_w.differences_distribution">Distribution 620 of differences from 'best' <code class="computeroutput"><span class="keyword">double</span></code> 621 evaluations</a> 622 </h6> 623<p> 624 The distribution of differences from 'best' are shown in these graphs comparing 625 <code class="computeroutput"><span class="keyword">double</span></code> precision evaluations with 626 reference 'best' z50 evaluations using <code class="computeroutput"><span class="identifier">cpp_bin_float_50</span></code> 627 type reduced to <code class="computeroutput"><span class="keyword">double</span></code> with <code class="computeroutput"><span class="keyword">static_cast</span><span class="special"><</span><span class="keyword">double</span><span class="special">(</span><span class="identifier">z50</span><span class="special">)</span></code> : 628 </p> 629<div class="blockquote"><blockquote class="blockquote"><p> 630 <span class="inlinemediaobject"><img src="../../graphs/lambert_w0_errors_graph.svg" align="middle"></span> 631 632 </p></blockquote></div> 633<div class="blockquote"><blockquote class="blockquote"><p> 634 <span class="inlinemediaobject"><img src="../../graphs/lambert_wm1_errors_graph.svg" align="middle"></span> 635 636 </p></blockquote></div> 637<p> 638 As noted in the implementation section, the distribution of these differences 639 is somewhat biased for Lambert <span class="emphasis"><em>W</em></span><sub>-1</sub> and this might be reduced 640 using a <code class="computeroutput"><span class="keyword">double</span></code> Halley step at 641 small runtime cost. But if you are seriously concerned to get really precise 642 computations, the only way is using a higher precision type and then reduce 643 to the desired type. Fortunately, <a href="../../../../../libs/multiprecision/doc/html/index.html" target="_top">Boost.Multiprecision</a> 644 makes this very easy to program, if much slower. 645 </p> 646<h5> 647<a name="math_toolkit.lambert_w.h10"></a> 648 <span class="phrase"><a name="math_toolkit.lambert_w.edge_cases"></a></span><a class="link" href="lambert_w.html#math_toolkit.lambert_w.edge_cases">Edge 649 and Corner cases</a> 650 </h5> 651<h6> 652<a name="math_toolkit.lambert_w.h11"></a> 653 <span class="phrase"><a name="math_toolkit.lambert_w.w0_edges"></a></span><a class="link" href="lambert_w.html#math_toolkit.lambert_w.w0_edges">The 654 <span class="emphasis"><em>W</em></span><sub>0</sub> Branch</a> 655 </h6> 656<p> 657 The domain of <span class="emphasis"><em>W</em></span><sub>0</sub> is [-<span class="emphasis"><em>e</em></span><sup>-1</sup>, ∞). Numerically, 658 </p> 659<div class="itemizedlist"><ul class="itemizedlist" style="list-style-type: disc; "> 660<li class="listitem"> 661 <code class="computeroutput"><span class="identifier">lambert_w0</span><span class="special">(-</span><span class="number">1</span><span class="special">/</span><span class="identifier">e</span><span class="special">)</span></code> is exactly -1. 662 </li> 663<li class="listitem"> 664 <code class="computeroutput"><span class="identifier">lambert_w0</span><span class="special">(</span><span class="identifier">z</span><span class="special">)</span></code> for 665 <code class="computeroutput"><span class="identifier">z</span> <span class="special"><</span> 666 <span class="special">-</span><span class="number">1</span><span class="special">/</span><span class="identifier">e</span></code> throws 667 a <code class="computeroutput"><span class="identifier">domain_error</span></code>, or returns 668 <code class="computeroutput"><span class="identifier">NaN</span></code> according to the policy. 669 </li> 670<li class="listitem"> 671 <code class="computeroutput"><span class="identifier">lambert_w0</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="identifier">T</span><span class="special">>::</span><span class="identifier">infinity</span><span class="special">())</span></code> 672 throws an <code class="computeroutput"><span class="identifier">overflow_error</span></code>. 673 </li> 674</ul></div> 675<p> 676 (An infinite argument probably indicates that something has already gone wrong, 677 but if it is desired to return infinity, this case should be handled before 678 calling <code class="computeroutput"><span class="identifier">lambert_w0</span></code>). 679 </p> 680<h6> 681<a name="math_toolkit.lambert_w.h12"></a> 682 <span class="phrase"><a name="math_toolkit.lambert_w.wm1_edges"></a></span><a class="link" href="lambert_w.html#math_toolkit.lambert_w.wm1_edges"><span class="emphasis"><em>W</em></span><sub>-1</sub> Branch</a> 683 </h6> 684<p> 685 The domain of <span class="emphasis"><em>W</em></span><sub>-1</sub> is [-<span class="emphasis"><em>e</em></span><sup>-1</sup>, 0). Numerically, 686 </p> 687<div class="itemizedlist"><ul class="itemizedlist" style="list-style-type: disc; "> 688<li class="listitem"> 689 <code class="computeroutput"><span class="identifier">lambert_wm1</span><span class="special">(-</span><span class="number">1</span><span class="special">/</span><span class="identifier">e</span><span class="special">)</span></code> is exactly -1. 690 </li> 691<li class="listitem"> 692 <code class="computeroutput"><span class="identifier">lambert_wm1</span><span class="special">(</span><span class="number">0</span><span class="special">)</span></code> returns 693 -∞ (or the nearest equivalent if <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">has_infinity</span> 694 <span class="special">==</span> <span class="keyword">false</span></code>). 695 </li> 696<li class="listitem"> 697 <code class="computeroutput"><span class="identifier">lambert_wm1</span><span class="special">(-</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="identifier">T</span><span class="special">>::</span><span class="identifier">min</span><span class="special">())</span></code> 698 returns the maximum (most negative) possible value of Lambert <span class="emphasis"><em>W</em></span> 699 for the type T. <br> For example, for <code class="computeroutput"><span class="keyword">double</span></code>: 700 lambert_wm1(-2.2250738585072014e-308) = -714.96865723796634 <br> and 701 for <code class="computeroutput"><span class="keyword">float</span></code>: lambert_wm1(-1.17549435e-38) 702 = -91.8567734 <br> 703 </li> 704<li class="listitem"> 705<p class="simpara"> 706 <code class="computeroutput"><span class="identifier">z</span> <span class="special"><</span> 707 <span class="special">-</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="identifier">T</span><span class="special">>::</span><span class="identifier">min</span><span class="special">()</span></code>, means that z is zero or denormalized 708 (if <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="identifier">T</span><span class="special">>::</span><span class="identifier">has_denorm_min</span> <span class="special">==</span> 709 <span class="keyword">true</span></code>), for example: <code class="computeroutput"><span class="identifier">r</span> <span class="special">=</span> <span class="identifier">lambert_wm1</span><span class="special">(-</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="keyword">double</span><span class="special">>::</span><span class="identifier">denorm_min</span><span class="special">());</span></code> 710 and an overflow_error exception is thrown, and will give a message like: 711 </p> 712<p class="simpara"> 713 Error in function boost::math::lambert_wm1<RealType>(<RealType>): 714 Argument z = -4.9406564584124654e-324 is too small (z < -std::numeric_limits<T>::min 715 so denormalized) for Lambert W-1 branch! 716 </p> 717</li> 718</ul></div> 719<p> 720 Denormalized values are not supported for Lambert <span class="emphasis"><em>W</em></span><sub>-1</sub> (because 721 not all floating-point types denormalize), and anyway it only covers a tiny 722 fraction of the range of possible z arguments values. 723 </p> 724<h5> 725<a name="math_toolkit.lambert_w.h13"></a> 726 <span class="phrase"><a name="math_toolkit.lambert_w.compilers"></a></span><a class="link" href="lambert_w.html#math_toolkit.lambert_w.compilers">Compilers</a> 727 </h5> 728<p> 729 The <code class="computeroutput"><span class="identifier">lambert_w</span><span class="special">.</span><span class="identifier">hpp</span></code> code has been shown to work on most C++98 730 compilers. (Apart from requiring C++11 extensions for using of <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><>::</span><span class="identifier">max_digits10</span></code> 731 in some diagnostics. Many old pre-c++11 compilers provide this extension but 732 may require enabling to use, for example using b2/bjam the lambert_w examples 733 use this command: 734 </p> 735<pre class="programlisting"><span class="special">[</span> <span class="identifier">run</span> <span class="identifier">lambert_w_basic_example</span><span class="special">.</span><span class="identifier">cpp</span> <span class="special">:</span> <span class="special">:</span> <span class="special">:</span> <span class="special">[</span> <span class="identifier">requires</span> <span class="identifier">cxx11_numeric_limits</span> <span class="special">]</span> <span class="special">]</span> 736</pre> 737<p> 738 See <a href="../../../example/Jamfile.v2" target="_top">jamfile.v2</a>.) 739 </p> 740<p> 741 For details of which compilers are expected to work see lambert_w tests and 742 examples in:<br> <a href="https://www.boost.org/development/tests/master/developer/math.html" target="_top">Boost 743 Test Summary report for master branch (used for latest release)</a><br> 744 <a href="https://www.boost.org/development/tests/develop/developer/math.html" target="_top">Boost 745 Test Summary report for latest developer branch</a>. 746 </p> 747<p> 748 As expected, debug mode is very much slower than release. 749 </p> 750<h6> 751<a name="math_toolkit.lambert_w.h14"></a> 752 <span class="phrase"><a name="math_toolkit.lambert_w.diagnostics"></a></span><a class="link" href="lambert_w.html#math_toolkit.lambert_w.diagnostics">Diagnostics 753 Macros</a> 754 </h6> 755<p> 756 Several macros are provided to output diagnostic information (potentially 757 <span class="bold"><strong>much</strong></span> output). These can be statements, for 758 example: 759 </p> 760<p> 761 <code class="computeroutput"><span class="preprocessor">#define</span> <span class="identifier">BOOST_MATH_INSTRUMENT_LAMBERT_W_TERMS</span></code> 762 </p> 763<p> 764 placed <span class="bold"><strong>before</strong></span> the <code class="computeroutput"><span class="identifier">lambert_w</span></code> 765 include statement 766 </p> 767<p> 768 <code class="computeroutput"><span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">boost</span><span class="special">/</span><span class="identifier">math</span><span class="special">/</span><span class="identifier">special_functions</span><span class="special">/</span><span class="identifier">lambert_w</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span></code>, 769 </p> 770<p> 771 or defined on the project compile command-line: <code class="computeroutput"><span class="special">/</span><span class="identifier">DBOOST_MATH_INSTRUMENT_LAMBERT_W_TERMS</span></code>, 772 </p> 773<p> 774 or defined in a jamfile.v2: <code class="computeroutput"><span class="special"><</span><span class="identifier">define</span><span class="special">></span><span class="identifier">BOOST_MATH_INSTRUMENT_LAMBERT_W_TERMS</span></code> 775 </p> 776<pre class="programlisting"><span class="comment">// #define-able macros</span> 777<span class="identifier">BOOST_MATH_INSTRUMENT_LAMBERT_W_HALLEY</span> <span class="comment">// Halley refinement diagnostics.</span> 778<span class="identifier">BOOST_MATH_INSTRUMENT_LAMBERT_W_PRECISION</span> <span class="comment">// Precision.</span> 779<span class="identifier">BOOST_MATH_INSTRUMENT_LAMBERT_WM1</span> <span class="comment">// W1 branch diagnostics.</span> 780<span class="identifier">BOOST_MATH_INSTRUMENT_LAMBERT_WM1_HALLEY</span> <span class="comment">// Halley refinement diagnostics only for W-1 branch.</span> 781<span class="identifier">BOOST_MATH_INSTRUMENT_LAMBERT_WM1_TINY</span> <span class="comment">// K > 64, z > -1.0264389699511303e-26</span> 782<span class="identifier">BOOST_MATH_INSTRUMENT_LAMBERT_WM1_LOOKUP</span> <span class="comment">// Show results from W-1 lookup table.</span> 783<span class="identifier">BOOST_MATH_INSTRUMENT_LAMBERT_W_SCHROEDER</span> <span class="comment">// Schroeder refinement diagnostics.</span> 784<span class="identifier">BOOST_MATH_INSTRUMENT_LAMBERT_W_TERMS</span> <span class="comment">// Number of terms used for near-singularity series.</span> 785<span class="identifier">BOOST_MATH_INSTRUMENT_LAMBERT_W_SINGULARITY_SERIES</span> <span class="comment">// Show evaluation of series near branch singularity.</span> 786<span class="identifier">BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES</span> 787<span class="identifier">BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES_ITERATIONS</span> <span class="comment">// Show evaluation of series for small z.</span> 788</pre> 789<h5> 790<a name="math_toolkit.lambert_w.h15"></a> 791 <span class="phrase"><a name="math_toolkit.lambert_w.implementation"></a></span><a class="link" href="lambert_w.html#math_toolkit.lambert_w.implementation">Implementation</a> 792 </h5> 793<p> 794 There are many previous implementations, each with increasing accuracy and/or 795 speed. See <a class="link" href="lambert_w.html#math_toolkit.lambert_w.references">references</a> 796 below. 797 </p> 798<p> 799 For most of the range of <span class="emphasis"><em>z</em></span> arguments, some initial approximation 800 followed by a single refinement, often using Halley or similar method, gives 801 a useful precision. For speed, several implementations avoid evaluation of 802 a iteration test using the exponential function, estimating that a single refinement 803 step will suffice, but these rarely get to the best result possible. To get 804 a better precision, additional refinements, probably iterative, are needed 805 for example, using <a class="link" href="roots_deriv.html#math_toolkit.roots_deriv.halley">Halley</a> 806 or <a class="link" href="roots_deriv.html#math_toolkit.roots_deriv.schroder">Schröder</a> methods. 807 </p> 808<p> 809 For C++, the most precise results possible, closest to the nearest <a href="http://en.wikipedia.org/wiki/Floating_point#Representable_numbers.2C_conversion_and_rounding" target="_top">representable</a> 810 for the C++ type being used, it is usually necessary to use a higher precision 811 type for intermediate computation, finally static-casting back to the smaller 812 desired result type. This strategy is used by <a href="https://www.maplesoft.com" target="_top">Maple</a> 813 and <a href="http://www.wolframalpha.com/" target="_top">Wolfram Alpha</a>, for example, 814 using arbitrary precision arithmetic, and some of their high-precision values 815 are used for testing this library. This method is also used to provide some 816 <a href="https://www.boost.org/doc/libs/release/libs/test/doc/html/index.html" target="_top">Boost.Test</a> 817 values using <a href="../../../../../libs/multiprecision/doc/html/index.html" target="_top">Boost.Multiprecision</a>, 818 typically, a 50 decimal digit type like <code class="computeroutput"><span class="identifier">cpp_bin_float_50</span></code> 819 <code class="computeroutput"><span class="keyword">static_cast</span></code> to a <code class="computeroutput"><span class="keyword">float</span></code>, <code class="computeroutput"><span class="keyword">double</span></code> 820 or <code class="computeroutput"><span class="keyword">long</span> <span class="keyword">double</span></code> 821 type. 822 </p> 823<p> 824 For <span class="emphasis"><em>z</em></span> argument values near the singularity and near zero, 825 other approximations may be used, possibly followed by refinement or increasing 826 number of series terms until a desired precision is achieved. At extreme arguments 827 near to zero or the singularity at the branch point, even this fails and the 828 only method to achieve a really close result is to cast from a higher precision 829 type. 830 </p> 831<p> 832 In practical applications, the increased computation required (often towards 833 a thousand-fold slower and requiring much additional code for <a href="../../../../../libs/multiprecision/doc/html/index.html" target="_top">Boost.Multiprecision</a>) 834 is not justified and the algorithms here do not implement this. But because 835 the Boost.Lambert_W algorithms has been tested using <a href="../../../../../libs/multiprecision/doc/html/index.html" target="_top">Boost.Multiprecision</a>, 836 users who require this can always easily achieve the nearest representation 837 for <a href="http://en.cppreference.com/w/cpp/language/types" target="_top">fundamental 838 (built-in) types</a> - if the application justifies the very large extra 839 computation cost. 840 </p> 841<h6> 842<a name="math_toolkit.lambert_w.h16"></a> 843 <span class="phrase"><a name="math_toolkit.lambert_w.evolution_of_this_implementation"></a></span><a class="link" href="lambert_w.html#math_toolkit.lambert_w.evolution_of_this_implementation">Evolution 844 of this implementation</a> 845 </h6> 846<p> 847 One compact real-only implementation was based on an algorithm by <a href="http://discovery.ucl.ac.uk/1482128/1/Luu_thesis.pdf" target="_top">Thomas 848 Luu, Thesis, University College London (2015)</a>, (see routine 11 on page 849 98 for his Lambert W algorithm) and his Halley refinement is used iteratively 850 when required. A first implementation was based on Thomas Luu's code posted 851 at <a href="https://svn.boost.org/trac/boost/ticket/11027" target="_top">Boost Trac #11027</a>. 852 It has been implemented from Luu's algorithm but templated on <code class="computeroutput"><span class="identifier">RealType</span></code> parameter and result and handles 853 both <a href="http://en.cppreference.com/w/cpp/language/types" target="_top">fundamental 854 (built-in) types</a> (<code class="computeroutput"><span class="keyword">float</span><span class="special">,</span> <span class="keyword">double</span><span class="special">,</span> 855 <span class="keyword">long</span> <span class="keyword">double</span></code>), 856 <a href="../../../../../libs/multiprecision/doc/html/index.html" target="_top">Boost.Multiprecision</a>, 857 and also has been tested successfully with a proposed fixed_point type. 858 </p> 859<p> 860 A first approximation was computed using the method of Barry et al (see references 861 5 & 6 below). This was extended to the widely used <a href="https://people.sc.fsu.edu/~jburkardt/f_src/toms443/toms443.html" target="_top">TOMS443</a> 862 FORTRAN and C++ versions by John Burkardt using Schroeder refinement(s). (For 863 users only requiring an accuracy of relative accuracy of 0.02%, Barry's function 864 alone might suffice, but a better <a href="https://en.wikipedia.org/wiki/Rational_function" target="_top">rational 865 function</a> approximation method has since been developed for this implementation). 866 </p> 867<p> 868 We also considered using <a class="link" href="roots_deriv.html#math_toolkit.roots_deriv.newton">Newton-Raphson 869 iteration</a> method. 870 </p> 871<pre class="programlisting"><span class="identifier">f</span><span class="special">(</span><span class="identifier">w</span><span class="special">)</span> <span class="special">=</span> <span class="identifier">w</span> <span class="identifier">e</span><span class="special">^</span><span class="identifier">w</span> <span class="special">-</span><span class="identifier">z</span> <span class="special">=</span> <span class="number">0</span> <span class="comment">// Luu equation 6.37</span> 872<span class="identifier">f</span><span class="char">'(w) = e^w (1 + w), Wolfram alpha (d)/(dw)(f(w) = w exp(w) - z) = e^w (w + 1) 873if (f(w) / f'</span><span class="special">(</span><span class="identifier">w</span><span class="special">)</span> <span class="special">-</span><span class="number">1</span> <span class="special"><</span> <span class="identifier">tolerance</span> 874<span class="identifier">w1</span> <span class="special">=</span> <span class="identifier">w0</span> <span class="special">-</span> <span class="special">(</span><span class="identifier">expw0</span> <span class="special">*</span> <span class="special">(</span><span class="identifier">w0</span> <span class="special">+</span> <span class="number">1</span><span class="special">));</span> <span class="comment">// Refine new Newton/Raphson estimate.</span> 875</pre> 876<p> 877 but concluded that since the Newton-Raphson method takes typically 6 iterations 878 to converge within tolerance, whereas Halley usually takes only 1 to 3 iterations 879 to achieve an result within 1 <a href="http://en.wikipedia.org/wiki/Unit_in_the_last_place" target="_top">Unit 880 in the last place (ULP)</a>, so the Newton-Raphson method is unlikely to 881 be quicker than the additional cost of computing the 2nd derivative for Halley's 882 method. 883 </p> 884<p> 885 Halley refinement uses the simplified formulae obtained from <a href="http://www.wolframalpha.com/input/?i=%5B2(z+exp(z)-w)+d%2Fdx+(z+exp(z)-w)%5D+%2F+%5B2+(d%2Fdx+(z+exp(z)-w))%5E2+-+(z+exp(z)-w)+d%5E2%2Fdx%5E2+(z+exp(z)-w)%5D" target="_top">Wolfram 886 Alpha</a> 887 </p> 888<pre class="programlisting"><span class="special">[</span><span class="number">2</span><span class="special">(</span><span class="identifier">z</span> <span class="identifier">exp</span><span class="special">(</span><span class="identifier">z</span><span class="special">)-</span><span class="identifier">w</span><span class="special">)</span> <span class="identifier">d</span><span class="special">/</span><span class="identifier">dx</span> <span class="special">(</span><span class="identifier">z</span> <span class="identifier">exp</span><span class="special">(</span><span class="identifier">z</span><span class="special">)-</span><span class="identifier">w</span><span class="special">)]</span> <span class="special">/</span> <span class="special">[</span><span class="number">2</span> <span class="special">(</span><span class="identifier">d</span><span class="special">/</span><span class="identifier">dx</span> <span class="special">(</span><span class="identifier">z</span> <span class="identifier">exp</span><span class="special">(</span><span class="identifier">z</span><span class="special">)-</span><span class="identifier">w</span><span class="special">))^</span><span class="number">2</span> <span class="special">-</span> <span class="special">(</span><span class="identifier">z</span> <span class="identifier">exp</span><span class="special">(</span><span class="identifier">z</span><span class="special">)-</span><span class="identifier">w</span><span class="special">)</span> <span class="identifier">d</span><span class="special">^</span><span class="number">2</span><span class="special">/</span><span class="identifier">dx</span><span class="special">^</span><span class="number">2</span> <span class="special">(</span><span class="identifier">z</span> <span class="identifier">exp</span><span class="special">(</span><span class="identifier">z</span><span class="special">)-</span><span class="identifier">w</span><span class="special">)]</span> 889</pre> 890<h5> 891<a name="math_toolkit.lambert_w.h17"></a> 892 <span class="phrase"><a name="math_toolkit.lambert_w.compact_implementation"></a></span><a class="link" href="lambert_w.html#math_toolkit.lambert_w.compact_implementation">Implementing 893 Compact Algorithms</a> 894 </h5> 895<p> 896 The most compact algorithm can probably be implemented using the log approximation 897 of Corless et al. followed by Halley iteration (but is also slowest and least 898 precise near zero and near the branch singularity). 899 </p> 900<h5> 901<a name="math_toolkit.lambert_w.h18"></a> 902 <span class="phrase"><a name="math_toolkit.lambert_w.faster_implementation"></a></span><a class="link" href="lambert_w.html#math_toolkit.lambert_w.faster_implementation">Implementing 903 Faster Algorithms</a> 904 </h5> 905<p> 906 More recently, the Tosio Fukushima has developed an even faster algorithm, 907 avoiding any transcendental function calls as these are necessarily expensive. 908 The current implementation of Lambert <span class="emphasis"><em>W</em></span><sub>-1</sub> is based on his 909 algorithm starting with a translation from Fukushima's FORTRAN into C++ by 910 Darko Veberic. 911 </p> 912<p> 913 Many applications of the Lambert W function make many repeated evaluations 914 for Monte Carlo methods; for these applications speed is very important. Luu, 915 and Chapeau-Blondeau and Monir provide typical usage examples. 916 </p> 917<p> 918 Fukushima improves the important observation that much of the execution time 919 of all previous iterative algorithms was spent evaluating transcendental functions, 920 usually <code class="computeroutput"><span class="identifier">exp</span></code>. He has put a lot 921 of work into avoiding any slow transcendental functions by using lookup tables 922 and bisection, finishing with a single Schroeder refinement, without any check 923 on the final precision of the result (necessarily evaluating an expensive exponential). 924 </p> 925<p> 926 Theoretical and practical tests confirm that Fukushima's algorithm gives Lambert 927 W estimates with a known small error bound (several <a href="http://en.wikipedia.org/wiki/Unit_in_the_last_place" target="_top">Unit 928 in the last place (ULP)</a>) over nearly all the range of <span class="emphasis"><em>z</em></span> 929 argument. 930 </p> 931<p> 932 A mean difference was computed to express the typical error and is often about 933 0.5 epsilon, the theoretical minimum. Using the <a href="../../../../../libs/math/doc/html/math_toolkit/next_float/float_distance.html" target="_top">Boost.Math 934 float_distance</a>, we can also express this as the number of bits that 935 are different from the nearest representable or 'exact' or 'best' value. The 936 number and distribution of these few bits differences was studied by binning, 937 including their sign. Bins for (signed) 0, 1, 2 and 3 and 4 bits proved suitable. 938 </p> 939<p> 940 However, though these give results within a few <a href="http://en.wikipedia.org/wiki/Machine_epsilon" target="_top">machine 941 epsilon</a> of the nearest representable result, they do not get as close 942 as is very often possible with further refinement, nearly always to within 943 one or two <a href="http://en.wikipedia.org/wiki/Machine_epsilon" target="_top">machine 944 epsilon</a>. 945 </p> 946<p> 947 More significantly, the evaluations of the sum of all signed differences using 948 the Fukshima algorithm show a slight bias, being more likely to be a bit or 949 few below the nearest representation than above; bias might have unwanted effects 950 on some statistical computations. 951 </p> 952<p> 953 Fukushima's method also does not cover the full range of z arguments of 'float' 954 precision and above. 955 </p> 956<p> 957 For this implementation of Lambert <span class="emphasis"><em>W</em></span><sub>0</sub>, John Maddock used 958 the Boost.Math <a href="http://en.wikipedia.org/wiki/Remez_algorithm" target="_top">Remez 959 algorithm</a> method program to devise a <a href="https://en.wikipedia.org/wiki/Rational_function" target="_top">rational 960 function</a> for several ranges of argument for the <span class="emphasis"><em>W</em></span><sub>0</sub> branch 961 of Lambert <span class="emphasis"><em>W</em></span> function. These minimax rational approximations 962 are combined for an algorithm that is both smaller and faster. 963 </p> 964<p> 965 Sadly it has not proved practical to use the same <a href="http://en.wikipedia.org/wiki/Remez_algorithm" target="_top">Remez 966 algorithm</a> method for Lambert <span class="emphasis"><em>W</em></span><sub>-1</sub> branch and so 967 the Fukushima algorithm is retained for this branch. 968 </p> 969<p> 970 An advantage of both minimax rational <a href="http://en.wikipedia.org/wiki/Remez_algorithm" target="_top">Remez 971 algorithm</a> approximations is that the <span class="bold"><strong>distribution</strong></span> 972 from the reference values is reasonably random and insignificantly biased. 973 </p> 974<p> 975 For example, table below a test of Lambert <span class="emphasis"><em>W</em></span><sub>0</sub> 10000 values 976 of argument covering the main range of possible values, 10000 comparisons from 977 z = 0.0501 to 703, in 0.001 step factor 1.05 when module 7 == 0 978 </p> 979<div class="table"> 980<a name="math_toolkit.lambert_w.lambert_w0_Fukushima"></a><p class="title"><b>Table 8.73. Fukushima Lambert <span class="emphasis"><em>W</em></span><sub>0</sub> and typical improvement from 981 a single Halley step.</b></p> 982<div class="table-contents"><table class="table" summary="Fukushima Lambert W0 and typical improvement from 983 a single Halley step."> 984<colgroup> 985<col> 986<col> 987<col> 988<col> 989<col> 990<col> 991<col> 992</colgroup> 993<thead><tr> 994<th> 995 <p> 996 Method 997 </p> 998 </th> 999<th> 1000 <p> 1001 Exact 1002 </p> 1003 </th> 1004<th> 1005 <p> 1006 One_bit 1007 </p> 1008 </th> 1009<th> 1010 <p> 1011 Two_bits 1012 </p> 1013 </th> 1014<th> 1015 <p> 1016 Few_bits 1017 </p> 1018 </th> 1019<th> 1020 <p> 1021 inexact 1022 </p> 1023 </th> 1024<th> 1025 <p> 1026 bias 1027 </p> 1028 </th> 1029</tr></thead> 1030<tbody> 1031<tr> 1032<td> 1033 <p> 1034 Schroeder <span class="emphasis"><em>W</em></span><sub>0</sub> 1035 </p> 1036 </td> 1037<td> 1038 <p> 1039 8804 1040 </p> 1041 </td> 1042<td> 1043 <p> 1044 1154 1045 </p> 1046 </td> 1047<td> 1048 <p> 1049 37 1050 </p> 1051 </td> 1052<td> 1053 <p> 1054 5 1055 </p> 1056 </td> 1057<td> 1058 <p> 1059 1243 1060 </p> 1061 </td> 1062<td> 1063 <p> 1064 -1193 1065 </p> 1066 </td> 1067</tr> 1068<tr> 1069<td> 1070 <p> 1071 after Halley step 1072 </p> 1073 </td> 1074<td> 1075 <p> 1076 9710 1077 </p> 1078 </td> 1079<td> 1080 <p> 1081 288 1082 </p> 1083 </td> 1084<td> 1085 <p> 1086 2 1087 </p> 1088 </td> 1089<td> 1090 <p> 1091 0 1092 </p> 1093 </td> 1094<td> 1095 <p> 1096 292 1097 </p> 1098 </td> 1099<td> 1100 <p> 1101 22 1102 </p> 1103 </td> 1104</tr> 1105</tbody> 1106</table></div> 1107</div> 1108<br class="table-break"><p> 1109 Lambert <span class="emphasis"><em>W</em></span><sub>0</sub> values computed using the Fukushima method with 1110 Schroeder refinement gave about 1/6 <code class="computeroutput"><span class="identifier">lambert_w0</span></code> 1111 values that are one bit different from the 'best', and < 1% that are a few 1112 bits 'wrong'. If a Halley refinement step is added, only 1 in 30 are even one 1113 bit different, and only 2 two-bits 'wrong'. 1114 </p> 1115<div class="table"> 1116<a name="math_toolkit.lambert_w.lambert_w0_plus_halley"></a><p class="title"><b>Table 8.74. Rational polynomial Lambert <span class="emphasis"><em>W</em></span><sub>0</sub> and typical improvement 1117 from a single Halley step.</b></p> 1118<div class="table-contents"><table class="table" summary="Rational polynomial Lambert W0 and typical improvement 1119 from a single Halley step."> 1120<colgroup> 1121<col> 1122<col> 1123<col> 1124<col> 1125<col> 1126<col> 1127<col> 1128</colgroup> 1129<thead><tr> 1130<th> 1131 <p> 1132 Method 1133 </p> 1134 </th> 1135<th> 1136 <p> 1137 Exact 1138 </p> 1139 </th> 1140<th> 1141 <p> 1142 One_bit 1143 </p> 1144 </th> 1145<th> 1146 <p> 1147 Two_bits 1148 </p> 1149 </th> 1150<th> 1151 <p> 1152 Few_bits 1153 </p> 1154 </th> 1155<th> 1156 <p> 1157 inexact 1158 </p> 1159 </th> 1160<th> 1161 <p> 1162 bias 1163 </p> 1164 </th> 1165</tr></thead> 1166<tbody> 1167<tr> 1168<td> 1169 <p> 1170 rational/polynomial 1171 </p> 1172 </td> 1173<td> 1174 <p> 1175 7135 1176 </p> 1177 </td> 1178<td> 1179 <p> 1180 2863 1181 </p> 1182 </td> 1183<td> 1184 <p> 1185 2 1186 </p> 1187 </td> 1188<td> 1189 <p> 1190 0 1191 </p> 1192 </td> 1193<td> 1194 <p> 1195 2867 1196 </p> 1197 </td> 1198<td> 1199 <p> 1200 -59 1201 </p> 1202 </td> 1203</tr> 1204<tr> 1205<td> 1206 <p> 1207 after Halley step 1208 </p> 1209 </td> 1210<td> 1211 <p> 1212 9724 1213 </p> 1214 </td> 1215<td> 1216 <p> 1217 273 1218 </p> 1219 </td> 1220<td> 1221 <p> 1222 3 1223 </p> 1224 </td> 1225<td> 1226 <p> 1227 0 1228 </p> 1229 </td> 1230<td> 1231 <p> 1232 279 1233 </p> 1234 </td> 1235<td> 1236 <p> 1237 5 1238 </p> 1239 </td> 1240</tr> 1241</tbody> 1242</table></div> 1243</div> 1244<br class="table-break"><p> 1245 With the rational polynomial approximation method, there are a third one-bit 1246 from the best and none more than two-bits. Adding a Halley step (or iteration) 1247 reduces the number that are one-bit different from about a third down to one 1248 in 30; this is unavoidable 'computational noise'. An extra Halley step would 1249 double the runtime for a tiny gain and so is not chosen for this implementation, 1250 but remains a option, as detailed above. 1251 </p> 1252<p> 1253 For the Lambert <span class="emphasis"><em>W</em></span><sub>-1</sub> branch, the Fukushima algorithm is 1254 used. 1255 </p> 1256<div class="table"> 1257<a name="math_toolkit.lambert_w.lambert_wm1_fukushima"></a><p class="title"><b>Table 8.75. Lambert <span class="emphasis"><em>W</em></span><sub>-1</sub> using Fukushima algorithm.</b></p> 1258<div class="table-contents"><table class="table" summary="Lambert W-1 using Fukushima algorithm."> 1259<colgroup> 1260<col> 1261<col> 1262<col> 1263<col> 1264<col> 1265<col> 1266<col> 1267</colgroup> 1268<thead><tr> 1269<th> 1270 <p> 1271 Method 1272 </p> 1273 </th> 1274<th> 1275 <p> 1276 Exact 1277 </p> 1278 </th> 1279<th> 1280 <p> 1281 One_bit 1282 </p> 1283 </th> 1284<th> 1285 <p> 1286 Two_bits 1287 </p> 1288 </th> 1289<th> 1290 <p> 1291 Few_bits 1292 </p> 1293 </th> 1294<th> 1295 <p> 1296 inexact 1297 </p> 1298 </th> 1299<th> 1300 <p> 1301 bias 1302 </p> 1303 </th> 1304</tr></thead> 1305<tbody> 1306<tr> 1307<td> 1308 <p> 1309 Fukushima <span class="emphasis"><em>W</em></span><sub>-1</sub> 1310 </p> 1311 </td> 1312<td> 1313 <p> 1314 7167 1315 </p> 1316 </td> 1317<td> 1318 <p> 1319 2704 1320 </p> 1321 </td> 1322<td> 1323 <p> 1324 129 1325 </p> 1326 </td> 1327<td> 1328 <p> 1329 0 1330 </p> 1331 </td> 1332<td> 1333 <p> 1334 2962 1335 </p> 1336 </td> 1337<td> 1338 <p> 1339 -160 1340 </p> 1341 </td> 1342</tr> 1343<tr> 1344<td> 1345 <p> 1346 plus Halley step 1347 </p> 1348 </td> 1349<td> 1350 <p> 1351 7379 1352 </p> 1353 </td> 1354<td> 1355 <p> 1356 2529 1357 </p> 1358 </td> 1359<td> 1360 <p> 1361 92 1362 </p> 1363 </td> 1364<td> 1365 <p> 1366 0 1367 </p> 1368 </td> 1369<td> 1370 <p> 1371 2713 1372 </p> 1373 </td> 1374<td> 1375 <p> 1376 549 1377 </p> 1378 </td> 1379</tr> 1380</tbody> 1381</table></div> 1382</div> 1383<br class="table-break"><h6> 1384<a name="math_toolkit.lambert_w.h19"></a> 1385 <span class="phrase"><a name="math_toolkit.lambert_w.lookup_tables"></a></span><a class="link" href="lambert_w.html#math_toolkit.lambert_w.lookup_tables">Lookup 1386 tables</a> 1387 </h6> 1388<p> 1389 For speed during the bisection, Fukushima's algorithm computes lookup tables 1390 of powers of e and z for integral Lambert W. There are 64 elements in these 1391 tables. The FORTRAN version (and the C++ translation by Veberic) computed these 1392 (once) as <code class="computeroutput"><span class="keyword">static</span></code> data. This is 1393 slower, may cause trouble with multithreading, and is slightly inaccurate because 1394 of rounding errors from repeated(64) multiplications. 1395 </p> 1396<p> 1397 In this implementation the array values have been computed using <a href="../../../../../libs/multiprecision/doc/html/index.html" target="_top">Boost.Multiprecision</a> 1398 50 decimal digit and output as C++ arrays 37 decimal digit <code class="computeroutput"><span class="keyword">long</span> 1399 <span class="keyword">double</span></code> literals using <code class="computeroutput"><span class="identifier">max_digits10</span></code> precision 1400 </p> 1401<pre class="programlisting"><span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span><span class="special">.</span><span class="identifier">precision</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="identifier">cpp_bin_float_quad</span><span class="special">>::</span><span class="identifier">max_digits10</span><span class="special">);</span> 1402</pre> 1403<p> 1404 The arrays are as <code class="computeroutput"><span class="keyword">const</span></code> and <code class="computeroutput"><span class="keyword">constexpr</span></code> and <code class="computeroutput"><span class="keyword">static</span></code> 1405 as possible (for the compiler version), using BOOST_STATIC_CONSTEXPR macro. 1406 (See <a href="../../../tools/lambert_w_lookup_table_generator.cpp" target="_top">lambert_w_lookup_table_generator.cpp</a> 1407 The precision was chosen to ensure that if used as <code class="computeroutput"><span class="keyword">long</span> 1408 <span class="keyword">double</span></code> arrays, then the values output 1409 to <a href="../../../include/boost/math/special_functions/detail/lambert_w_lookup_table.ipp" target="_top">lambert_w_lookup_table.ipp</a> 1410 will be the nearest representable value for the type chose by a <code class="computeroutput"><span class="keyword">typedef</span></code> in <a href="../../../include/boost/math/special_functions/lambert_w.hpp" target="_top">lambert_w.hpp</a>. 1411 </p> 1412<pre class="programlisting"><span class="keyword">typedef</span> <span class="keyword">double</span> <span class="identifier">lookup_t</span><span class="special">;</span> <span class="comment">// Type for lookup table (`double` or `float`, or even `long double`?)</span> 1413</pre> 1414<p> 1415 This is to allow for future use at higher precision, up to platforms that use 1416 128-bit (hardware or software) for their <code class="computeroutput"><span class="keyword">long</span> 1417 <span class="keyword">double</span></code> type. 1418 </p> 1419<p> 1420 The accuracy of the tables was confirmed using <a href="http://www.wolframalpha.com/" target="_top">Wolfram 1421 Alpha</a> and agrees at the 37th decimal place, so ensuring that the value 1422 is exactly read into even 128-bit <code class="computeroutput"><span class="keyword">long</span> 1423 <span class="keyword">double</span></code> to the nearest representation. 1424 </p> 1425<h6> 1426<a name="math_toolkit.lambert_w.h20"></a> 1427 <span class="phrase"><a name="math_toolkit.lambert_w.higher_precision"></a></span><a class="link" href="lambert_w.html#math_toolkit.lambert_w.higher_precision">Higher 1428 precision</a> 1429 </h6> 1430<p> 1431 For types more precise than <code class="computeroutput"><span class="keyword">double</span></code>, 1432 Fukushima reported that it was best to use the <code class="computeroutput"><span class="keyword">double</span></code> 1433 estimate as a starting point, followed by refinement using <a class="link" href="roots_deriv.html#math_toolkit.roots_deriv.halley">Halley</a> 1434 iterations or other methods; our experience confirms this. 1435 </p> 1436<p> 1437 Using <a href="../../../../../libs/multiprecision/doc/html/index.html" target="_top">Boost.Multiprecision</a> 1438 it is simple to compute very high precision values of Lambert W at least to 1439 thousands of decimal digits over most of the range of z arguments. 1440 </p> 1441<p> 1442 For this reason, the lookup tables and bisection are only carried out at low 1443 precision, usually <code class="computeroutput"><span class="keyword">double</span></code>, chosen 1444 by the <code class="computeroutput"><span class="keyword">typedef</span> <span class="keyword">double</span> 1445 <span class="identifier">lookup_t</span></code>. Unlike the FORTRAN version, 1446 the lookup tables of Lambert_W of integral values are precomputed as C++ static 1447 arrays of floating-point literals. The default is a <code class="computeroutput"><span class="keyword">typedef</span></code> 1448 setting the type to <code class="computeroutput"><span class="keyword">double</span></code>. To 1449 allow users to vary the precision from <code class="computeroutput"><span class="keyword">float</span></code> 1450 to <code class="computeroutput"><span class="keyword">long</span> <span class="keyword">double</span></code> 1451 these are computed to 128-bit precision to ensure that even platforms with 1452 <code class="computeroutput"><span class="keyword">long</span> <span class="keyword">double</span></code> 1453 do not lose precision. 1454 </p> 1455<p> 1456 The FORTRAN version and translation only permits the z argument to be the largest 1457 items in these lookup arrays, <code class="computeroutput"><span class="identifier">wm0s</span><span class="special">[</span><span class="number">64</span><span class="special">]</span> 1458 <span class="special">=</span> <span class="number">3.99049</span></code>, 1459 producing an error message and returning <code class="computeroutput"><span class="identifier">NaN</span></code>. 1460 So 64 is the largest possible value ever returned from the <code class="computeroutput"><span class="identifier">lambert_w0</span></code> 1461 function. This is far from the <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><>::</span><span class="identifier">max</span><span class="special">()</span></code> for even <code class="computeroutput"><span class="keyword">float</span></code>s. 1462 Therefore this implementation uses an approximation or 'guess' and Halley's 1463 method to refine the result. Logarithmic approximation is discussed at length 1464 by R.M.Corless et al. (page 349). Here we use the first two terms of equation 1465 4.19: 1466 </p> 1467<pre class="programlisting"><span class="identifier">T</span> <span class="identifier">lz</span> <span class="special">=</span> <span class="identifier">log</span><span class="special">(</span><span class="identifier">z</span><span class="special">);</span> 1468<span class="identifier">T</span> <span class="identifier">llz</span> <span class="special">=</span> <span class="identifier">log</span><span class="special">(</span><span class="identifier">lz</span><span class="special">);</span> 1469<span class="identifier">guess</span> <span class="special">=</span> <span class="identifier">lz</span> <span class="special">-</span> <span class="identifier">llz</span> <span class="special">+</span> <span class="special">(</span><span class="identifier">llz</span> <span class="special">/</span> <span class="identifier">lz</span><span class="special">);</span> 1470</pre> 1471<p> 1472 This gives a useful precision suitable for Halley refinement. 1473 </p> 1474<p> 1475 Similarly, for Lambert <span class="emphasis"><em>W</em></span><sub>-1</sub> branch, tiny values very near 1476 zero, W > 64 cannot be computed using the lookup table. For this region, 1477 an approximation followed by a few (usually 3) Halley refinements. See <a class="link" href="lambert_w.html#math_toolkit.lambert_w.wm1_near_zero">wm1_near_zero</a>. 1478 </p> 1479<p> 1480 For the less well-behaved regions for Lambert <span class="emphasis"><em>W</em></span><sub>0</sub> <span class="emphasis"><em>z</em></span> 1481 arguments near zero, and near the branch singularity at <span class="emphasis"><em>-1/e</em></span>, 1482 some series functions are used. 1483 </p> 1484<h6> 1485<a name="math_toolkit.lambert_w.h21"></a> 1486 <span class="phrase"><a name="math_toolkit.lambert_w.small_z"></a></span><a class="link" href="lambert_w.html#math_toolkit.lambert_w.small_z">Small 1487 values of argument z near zero</a> 1488 </h6> 1489<p> 1490 When argument <span class="emphasis"><em>z</em></span> is small and near zero, there is an efficient 1491 and accurate series evaluation method available (implemented in <code class="computeroutput"><span class="identifier">lambert_w0_small_z</span></code>). There is no equivalent 1492 for the <span class="emphasis"><em>W</em></span><sub>-1</sub> branch as this only covers argument <code class="computeroutput"><span class="identifier">z</span> <span class="special"><</span> <span class="special">-</span><span class="number">1</span><span class="special">/</span><span class="identifier">e</span></code>. 1493 The cutoff used <code class="computeroutput"><span class="identifier">abs</span><span class="special">(</span><span class="identifier">z</span><span class="special">)</span> <span class="special"><</span> 1494 <span class="number">0.05</span></code> is as found by trial and error by 1495 Fukushima. 1496 </p> 1497<p> 1498 Coefficients of the inverted series expansion of the Lambert W function around 1499 <code class="computeroutput"><span class="identifier">z</span> <span class="special">=</span> 1500 <span class="number">0</span></code> are computed following Fukushima using 1501 17 terms of a Taylor series computed using <a href="http://www.wolfram.com/products/mathematica/index.html" target="_top">Wolfram 1502 Mathematica</a> with 1503 </p> 1504<pre class="programlisting"><span class="identifier">InverseSeries</span><span class="special">[</span><span class="identifier">Series</span><span class="special">[</span><span class="identifier">z</span> <span class="identifier">Exp</span><span class="special">[</span><span class="identifier">z</span><span class="special">],{</span><span class="identifier">z</span><span class="special">,</span><span class="number">0</span><span class="special">,</span><span class="number">17</span><span class="special">}]]</span> 1505</pre> 1506<p> 1507 See Tosio Fukushima, Journal of Computational and Applied Mathematics 244 (2013), 1508 page 86. 1509 </p> 1510<p> 1511 To provide higher precision constants (34 decimal digits) for types larger 1512 than <code class="computeroutput"><span class="keyword">long</span> <span class="keyword">double</span></code>, 1513 </p> 1514<pre class="programlisting"><span class="identifier">InverseSeries</span><span class="special">[</span><span class="identifier">Series</span><span class="special">[</span><span class="identifier">z</span> <span class="identifier">Exp</span><span class="special">[</span><span class="identifier">z</span><span class="special">],{</span><span class="identifier">z</span><span class="special">,</span><span class="number">0</span><span class="special">,</span><span class="number">34</span><span class="special">}]]</span> 1515</pre> 1516<p> 1517 were also computed, but for current hardware it was found that evaluating a 1518 <code class="computeroutput"><span class="keyword">double</span></code> precision and then refining 1519 with Halley's method was quicker and more accurate. 1520 </p> 1521<p> 1522 Decimal values of specifications for built-in floating-point types below are 1523 21 digits precision == <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="identifier">T</span><span class="special">>::</span><span class="identifier">max_digits10</span></code> for <code class="computeroutput"><span class="keyword">long</span> 1524 <span class="keyword">double</span></code>. 1525 </p> 1526<p> 1527 Specializations for <code class="computeroutput"><span class="identifier">lambert_w0_small_z</span></code> 1528 are provided for <code class="computeroutput"><span class="keyword">float</span></code>, <code class="computeroutput"><span class="keyword">double</span></code>, <code class="computeroutput"><span class="keyword">long</span> 1529 <span class="keyword">double</span></code>, <code class="computeroutput"><span class="identifier">float128</span></code> 1530 and for <a href="../../../../../libs/multiprecision/doc/html/index.html" target="_top">Boost.Multiprecision</a> 1531 types. 1532 </p> 1533<p> 1534 The <code class="computeroutput"><span class="identifier">tag_type</span></code> selection is based 1535 on the value <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="identifier">T</span><span class="special">>::</span><span class="identifier">max_digits10</span></code> 1536 (and <span class="bold"><strong>not</strong></span> on the floating-point type T). This 1537 distinguishes between <code class="computeroutput"><span class="keyword">long</span> <span class="keyword">double</span></code> 1538 types that commonly vary between 64 and 80-bits, and also compilers that have 1539 a <code class="computeroutput"><span class="keyword">float</span></code> type using 64 bits and/or 1540 <code class="computeroutput"><span class="keyword">long</span> <span class="keyword">double</span></code> 1541 using 128-bits. 1542 </p> 1543<p> 1544 As noted in the <a class="link" href="lambert_w.html#math_toolkit.lambert_w.implementation">implementation</a> 1545 section above, it is only possible to ensure the nearest representable value 1546 by casting from a higher precision type, computed at very, very much greater 1547 cost. 1548 </p> 1549<p> 1550 For multiprecision types, first several terms of the series are tabulated and 1551 evaluated as a polynomial: (this will save us a bunch of expensive calls to 1552 <code class="computeroutput"><span class="identifier">pow</span></code>). Then our series functor 1553 is initialized "as if" it had already reached term 18, enough evaluation 1554 of built-in 64-bit double and float (and 80-bit <code class="computeroutput"><span class="keyword">long</span> 1555 <span class="keyword">double</span></code>) types. Finally the functor is 1556 called repeatedly to compute as many additional series terms as necessary to 1557 achieve the desired precision, set from <code class="computeroutput"><span class="identifier">get_epsilon</span></code> 1558 (or terminated by <code class="computeroutput"><span class="identifier">evaluation_error</span></code> 1559 on reaching the set iteration limit <code class="computeroutput"><span class="identifier">max_series_iterations</span></code>). 1560 </p> 1561<p> 1562 A little more than one decimal digit of precision is gained by each additional 1563 series term. This allows computation of Lambert W near zero to at least 1000 1564 decimal digit precision, given sufficient compute time. 1565 </p> 1566<h5> 1567<a name="math_toolkit.lambert_w.h22"></a> 1568 <span class="phrase"><a name="math_toolkit.lambert_w.near_singularity"></a></span><a class="link" href="lambert_w.html#math_toolkit.lambert_w.near_singularity">Argument 1569 z near the singularity at -1/e between branches <span class="emphasis"><em>W</em></span><sub>0</sub> and 1570 <span class="emphasis"><em>W</em></span><sub>-1</sub> </a> 1571 </h5> 1572<p> 1573 Variants of Function <code class="computeroutput"><span class="identifier">lambert_w_singularity_series</span></code> 1574 are used to handle <span class="emphasis"><em>z</em></span> arguments which are near to the singularity 1575 at <code class="computeroutput"><span class="identifier">z</span> <span class="special">=</span> 1576 <span class="special">-</span><span class="identifier">exp</span><span class="special">(-</span><span class="number">1</span><span class="special">)</span> 1577 <span class="special">=</span> <span class="special">-</span><span class="number">3.6787944</span></code> where the branches <span class="emphasis"><em>W</em></span><sub>0</sub> and 1578 <span class="emphasis"><em>W</em></span><sub>-1</sub> join. 1579 </p> 1580<p> 1581 T. Fukushima / Journal of Computational and Applied Mathematics 244 (2013) 1582 77-89 describes using <a href="http://www.wolfram.com/products/mathematica/index.html" target="_top">Wolfram 1583 Mathematica</a> 1584 </p> 1585<pre class="programlisting"><span class="identifier">InverseSeries</span><span class="special">\[</span><span class="identifier">Series</span><span class="special">\[</span><span class="identifier">sqrt</span><span class="special">\[</span><span class="number">2</span><span class="special">(</span><span class="identifier">p</span> <span class="identifier">Exp</span><span class="special">\[</span><span class="number">1</span> <span class="special">+</span> <span class="identifier">p</span><span class="special">\]</span> <span class="special">+</span> <span class="number">1</span><span class="special">)\],</span> <span class="special">{</span><span class="identifier">p</span><span class="special">,-</span><span class="number">1</span><span class="special">,</span> <span class="number">20</span><span class="special">}\]\]</span> 1586</pre> 1587<p> 1588 to provide his Table 3, page 85. 1589 </p> 1590<p> 1591 This implementation used <a href="http://www.wolfram.com/products/mathematica/index.html" target="_top">Wolfram 1592 Mathematica</a> to obtain 40 series terms at 50 decimal digit precision 1593 </p> 1594<pre class="programlisting"><span class="identifier">N</span><span class="special">\[</span><span class="identifier">InverseSeries</span><span class="special">\[</span><span class="identifier">Series</span><span class="special">\[</span><span class="identifier">Sqrt</span><span class="special">\[</span><span class="number">2</span><span class="special">(</span><span class="identifier">p</span> <span class="identifier">Exp</span><span class="special">\[</span><span class="number">1</span> <span class="special">+</span> <span class="identifier">p</span><span class="special">\]</span> <span class="special">+</span> <span class="number">1</span><span class="special">)\],</span> <span class="special">{</span> <span class="identifier">p</span><span class="special">,-</span><span class="number">1</span><span class="special">,</span><span class="number">40</span> <span class="special">}\]\],</span> <span class="number">50</span><span class="special">\]</span> 1595 1596<span class="special">-</span><span class="number">1</span><span class="special">+</span><span class="identifier">p</span><span class="special">-</span><span class="identifier">p</span><span class="special">^</span><span class="number">2</span><span class="special">/</span><span class="number">3</span><span class="special">+(</span><span class="number">11</span> <span class="identifier">p</span><span class="special">^</span><span class="number">3</span><span class="special">)/</span><span class="number">72</span><span class="special">-(</span><span class="number">43</span> <span class="identifier">p</span><span class="special">^</span><span class="number">4</span><span class="special">)/</span><span class="number">540</span><span class="special">+(</span><span class="number">769</span> <span class="identifier">p</span><span class="special">^</span><span class="number">5</span><span class="special">)/</span><span class="number">17280</span><span class="special">-(</span><span class="number">221</span> <span class="identifier">p</span><span class="special">^</span><span class="number">6</span><span class="special">)/</span><span class="number">8505</span><span class="special">+(</span><span class="number">680863</span> <span class="identifier">p</span><span class="special">^</span><span class="number">7</span><span class="special">)/</span><span class="number">43545600</span> <span class="special">...</span> 1597</pre> 1598<p> 1599 These constants are computed at compile time for the full precision for any 1600 <code class="computeroutput"><span class="identifier">RealType</span> <span class="identifier">T</span></code> 1601 using the original rationals from Fukushima Table 3. 1602 </p> 1603<p> 1604 Longer decimal digits strings are rationals pre-evaluated using <a href="http://www.wolfram.com/products/mathematica/index.html" target="_top">Wolfram 1605 Mathematica</a>. Some integer constants overflow, so largest size available 1606 is used, suffixed by <code class="computeroutput"><span class="identifier">uLL</span></code>. 1607 </p> 1608<p> 1609 Above the 14th term, the rationals exceed the range of <code class="computeroutput"><span class="keyword">unsigned</span> 1610 <span class="keyword">long</span> <span class="keyword">long</span></code> 1611 and are replaced by pre-computed decimal values at least 21 digits precision 1612 == <code class="computeroutput"><span class="identifier">max_digits10</span></code> for <code class="computeroutput"><span class="keyword">long</span> <span class="keyword">double</span></code>. 1613 </p> 1614<p> 1615 A macro <code class="computeroutput"><span class="identifier">BOOST_MATH_TEST_VALUE</span></code> 1616 (defined in <a href="../../../test/test_value.hpp" target="_top">test_value.hpp</a>) 1617 taking a decimal floating-point literal was used to allow testing with both 1618 built-in floating-point types like <code class="computeroutput"><span class="keyword">double</span></code> 1619 which have constructors taking literal decimal values like <code class="computeroutput"><span class="number">3.14</span></code>, 1620 <span class="bold"><strong>and</strong></span> also multiprecision and other User-defined 1621 Types that only provide full-precision construction from decimal digit strings 1622 like <code class="computeroutput"><span class="string">"3.14"</span></code>. (Construction 1623 of multiprecision types from built-in floating-point types only provides the 1624 precision of the built-in type, like <code class="computeroutput"><span class="keyword">double</span></code>, 1625 only 17 decimal digits). 1626 </p> 1627<div class="tip"><table border="0" summary="Tip"> 1628<tr> 1629<td rowspan="2" align="center" valign="top" width="25"><img alt="[Tip]" src="../../../../../doc/src/images/tip.png"></td> 1630<th align="left">Tip</th> 1631</tr> 1632<tr><td align="left" valign="top"><p> 1633 Be exceeding careful not to silently lose precision by constructing multiprecision 1634 types from literal decimal types, usually <code class="literal">double</code>. Use 1635 decimal digit strings like "3.1459" instead. See examples. 1636 </p></td></tr> 1637</table></div> 1638<p> 1639 Fukushima's implementation used 20 series terms; it was confirmed that using 1640 more terms does not usefully increase accuracy. 1641 </p> 1642<h6> 1643<a name="math_toolkit.lambert_w.h23"></a> 1644 <span class="phrase"><a name="math_toolkit.lambert_w.wm1_near_zero"></a></span><a class="link" href="lambert_w.html#math_toolkit.lambert_w.wm1_near_zero">Lambert 1645 <span class="emphasis"><em>W</em></span><sub>-1</sub> arguments values very near zero.</a> 1646 </h6> 1647<p> 1648 The lookup tables of Fukushima have only 64 elements, so that the z argument 1649 nearest zero is -1.0264389699511303e-26, that corresponds to a maximum Lambert 1650 <span class="emphasis"><em>W</em></span><sub>-1</sub> value of 64.0. Fukushima's implementation did not cater 1651 for z argument values that are smaller (nearer to zero), but this implementation 1652 adds code to accept smaller (but not denormalised) values of z. A crude approximation 1653 for these very small values is to take the exponent and multiply by ln[10] 1654 ~= 2.3. We also tried the approximation first proposed by Corless et al. using 1655 ln(-z), (equation 4.19 page 349) and then tried improving by a 2nd term -ln(ln(-z)), 1656 and finally the ratio term -ln(ln(-z))/ln(-z). 1657 </p> 1658<p> 1659 For a z very close to z = -1.0264389699511303e-26 when W = 64, when effect 1660 of ln(ln(-z) term, and ratio L1/L2 is greatest, the possible 'guesses' are 1661 </p> 1662<pre class="programlisting"><span class="identifier">z</span> <span class="special">=</span> <span class="special">-</span><span class="number">1.e-26</span><span class="special">,</span> <span class="identifier">w</span> <span class="special">=</span> <span class="special">-</span><span class="number">64.02</span><span class="special">,</span> <span class="identifier">guess</span> <span class="special">=</span> <span class="special">-</span><span class="number">64.0277</span><span class="special">,</span> <span class="identifier">ln</span><span class="special">(-</span><span class="identifier">z</span><span class="special">)</span> <span class="special">=</span> <span class="special">-</span><span class="number">59.8672</span><span class="special">,</span> <span class="identifier">ln</span><span class="special">(-</span><span class="identifier">ln</span><span class="special">(-</span><span class="identifier">z</span><span class="special">)</span> <span class="special">=</span> <span class="number">4.0921</span><span class="special">,</span> <span class="identifier">llz</span><span class="special">/</span><span class="identifier">lz</span> <span class="special">=</span> <span class="special">-</span><span class="number">0.0684</span> 1663</pre> 1664<p> 1665 whereas at the minimum (unnormalized) z 1666 </p> 1667<pre class="programlisting"><span class="identifier">z</span> <span class="special">=</span> <span class="special">-</span><span class="number">2.2250e-308</span><span class="special">,</span> <span class="identifier">w</span> <span class="special">=</span> <span class="special">-</span><span class="number">714.9</span><span class="special">,</span> <span class="identifier">guess</span> <span class="special">=</span> <span class="special">-</span><span class="number">714.9687</span><span class="special">,</span> <span class="identifier">ln</span><span class="special">(-</span><span class="identifier">z</span><span class="special">)</span> <span class="special">=</span> <span class="special">-</span><span class="number">708.3964</span><span class="special">,</span> <span class="identifier">ln</span><span class="special">(-</span><span class="identifier">ln</span><span class="special">(-</span><span class="identifier">z</span><span class="special">)</span> <span class="special">=</span> <span class="number">6.5630</span><span class="special">,</span> <span class="identifier">llz</span><span class="special">/</span><span class="identifier">lz</span> <span class="special">=</span> <span class="special">-</span><span class="number">0.0092</span> 1668</pre> 1669<p> 1670 Although the addition of the 3rd ratio term did not reduce the number of Halley 1671 iterations needed, it might allow return of a better low precision estimate 1672 <span class="bold"><strong>without any Halley iterations</strong></span>. For the worst 1673 case near w = 64, the error in the 'guess' is 0.008, ratio 0.0001 or 1 in 10,000 1674 digits 10 ~= 4. Two log evaluations are still needed, but is probably over 1675 an order of magnitude faster. 1676 </p> 1677<p> 1678 Halley's method was then used to refine the estimate of Lambert <span class="emphasis"><em>W</em></span><sub>-1</sub> from 1679 this guess. Experiments showed that although all approximations reached with 1680 <a href="http://en.wikipedia.org/wiki/Unit_in_the_last_place" target="_top">Unit in the 1681 last place (ULP)</a> of the closest representable value, the computational 1682 cost of the log functions was easily paid by far fewer iterations (typically 1683 from 8 down to 4 iterations for double or float). 1684 </p> 1685<h6> 1686<a name="math_toolkit.lambert_w.h24"></a> 1687 <span class="phrase"><a name="math_toolkit.lambert_w.halley"></a></span><a class="link" href="lambert_w.html#math_toolkit.lambert_w.halley">Halley 1688 refinement</a> 1689 </h6> 1690<p> 1691 After obtaining a double approximation, for <code class="computeroutput"><span class="keyword">double</span></code>, 1692 <code class="computeroutput"><span class="keyword">long</span> <span class="keyword">double</span></code> 1693 and <code class="computeroutput"><span class="identifier">quad</span></code> 128-bit precision, 1694 a single iteration should suffice because Halley iteration should triple the 1695 precision with each step (as long as the function is well behaved - and it 1696 is), and since we have at least half of the bits correct already, one Halley 1697 step is ample to get to 128-bit precision. 1698 </p> 1699<h6> 1700<a name="math_toolkit.lambert_w.h25"></a> 1701 <span class="phrase"><a name="math_toolkit.lambert_w.lambert_w_derivatives"></a></span><a class="link" href="lambert_w.html#math_toolkit.lambert_w.lambert_w_derivatives">Lambert 1702 W Derivatives</a> 1703 </h6> 1704<p> 1705 The derivatives are computed using the formulae in <a href="https://en.wikipedia.org/wiki/Lambert_W_function#Derivative" target="_top">Wikipedia</a>. 1706 </p> 1707<h5> 1708<a name="math_toolkit.lambert_w.h26"></a> 1709 <span class="phrase"><a name="math_toolkit.lambert_w.testing"></a></span><a class="link" href="lambert_w.html#math_toolkit.lambert_w.testing">Testing</a> 1710 </h5> 1711<p> 1712 Initial testing of the algorithm was done using a small number of spot tests. 1713 </p> 1714<p> 1715 After it was established that the underlying algorithm (including unlimited 1716 Halley refinements with a tight terminating criterion) was correct, some tables 1717 of Lambert W values were computed using a 100 decimal digit precision <a href="../../../../../libs/multiprecision/doc/html/index.html" target="_top">Boost.Multiprecision</a> 1718 <code class="computeroutput"><span class="identifier">cpp_dec_float_100</span></code> type and 1719 saved as a C++ program that will initialise arrays of values of z arguments 1720 and lambert_W0 (<code class="computeroutput"><span class="identifier">lambert_w_mp_high_values</span><span class="special">.</span><span class="identifier">ipp</span></code> and 1721 <code class="computeroutput"><span class="identifier">lambert_w_mp_low_values</span><span class="special">.</span><span class="identifier">ipp</span></code> ). 1722 </p> 1723<p> 1724 (A few of these pairs were checked against values computed by Wolfram Alpha 1725 to try to guard against mistakes; all those tested agreed to the penultimate 1726 decimal place, so they can be considered reliable to at least 98 decimal digits 1727 precision). 1728 </p> 1729<p> 1730 A macro <code class="computeroutput"><span class="identifier">BOOST_MATH_TEST_VALUE</span></code> 1731 was used to allow tests with any real type, both <a href="http://en.cppreference.com/w/cpp/language/types" target="_top">fundamental 1732 (built-in) types</a> and <a href="../../../../../libs/multiprecision/doc/html/index.html" target="_top">Boost.Multiprecision</a>. 1733 (This is necessary because <a href="http://en.cppreference.com/w/cpp/language/types" target="_top">fundamental 1734 (built-in) types</a> have a constructor from floating-point literals like 1735 3.1459F, 3.1459 or 3.1459L whereas <a href="../../../../../libs/multiprecision/doc/html/index.html" target="_top">Boost.Multiprecision</a> 1736 types may lose precision unless constructed from decimal digits strings like 1737 "3.1459"). 1738 </p> 1739<p> 1740 The 100-decimal digits precision pairs were then used to assess the precision 1741 of less-precise types, including <a href="../../../../../libs/multiprecision/doc/html/index.html" target="_top">Boost.Multiprecision</a> 1742 <code class="computeroutput"><span class="identifier">cpp_bin_float_quad</span></code> and <code class="computeroutput"><span class="identifier">cpp_bin_float_50</span></code>. <code class="computeroutput"><span class="keyword">static_cast</span></code>ing 1743 from the high precision types should give the closest representable value of 1744 the less-precise type; this is then be used to assess the precision of the 1745 Lambert W algorithm. 1746 </p> 1747<p> 1748 Tests using confirm that over nearly all the range of z arguments, nearly all 1749 estimates are the nearest <a href="http://en.wikipedia.org/wiki/Floating_point#Representable_numbers.2C_conversion_and_rounding" target="_top">representable</a> 1750 value, a minority are within 1 <a href="http://en.wikipedia.org/wiki/Unit_in_the_last_place" target="_top">Unit 1751 in the last place (ULP)</a> and only a very few 2 ULP. 1752 </p> 1753<div class="blockquote"><blockquote class="blockquote"><p> 1754 <span class="inlinemediaobject"><img src="../../graphs/lambert_w0_errors_graph.svg" align="middle"></span> 1755 1756 </p></blockquote></div> 1757<div class="blockquote"><blockquote class="blockquote"><p> 1758 <span class="inlinemediaobject"><img src="../../graphs/lambert_wm1_errors_graph.svg" align="middle"></span> 1759 1760 </p></blockquote></div> 1761<p> 1762 For the range of z arguments over the range -0.35 to 0.5, a different algorithm 1763 is used, but the same technique of evaluating reference values using a <a href="../../../../../libs/multiprecision/doc/html/index.html" target="_top">Boost.Multiprecision</a> 1764 <code class="computeroutput"><span class="identifier">cpp_dec_float_100</span></code> was used. 1765 For extremely small z arguments, near zero, and those extremely near the singularity 1766 at the branch point, precision can be much lower, as might be expected. 1767 </p> 1768<p> 1769 See source at: <a href="../../../example/lambert_w_simple_examples.cpp" target="_top">lambert_w_simple_examples.cpp</a> 1770 <a href="../../../test/test_lambert_w.cpp" target="_top">test_lambert_w.cpp</a> contains 1771 routine tests using <a href="https://www.boost.org/doc/libs/release/libs/test/doc/html/index.html" target="_top">Boost.Test</a>. 1772 <a href="../../../tools/lambert_w_errors_graph.cpp" target="_top">lambert_w_errors_graph.cpp</a> 1773 generating error graphs. 1774 </p> 1775<h6> 1776<a name="math_toolkit.lambert_w.h27"></a> 1777 <span class="phrase"><a name="math_toolkit.lambert_w.quadrature_testing"></a></span><a class="link" href="lambert_w.html#math_toolkit.lambert_w.quadrature_testing">Testing 1778 with quadrature</a> 1779 </h6> 1780<p> 1781 A further method of testing over a wide range of argument z values was devised 1782 by Nick Thompson (cunningly also to test the recently written quadrature routines 1783 including <a href="../../../../../libs/multiprecision/doc/html/index.html" target="_top">Boost.Multiprecision</a> 1784 !). These are definite integral formulas involving the W function that are 1785 exactly known constants, for example, LambertW0(1/(z²) == √(2π), see <a href="https://en.wikipedia.org/wiki/Lambert_W_function#Definite_integrals" target="_top">Definite 1786 Integrals</a>. Some care was needed to avoid overflow and underflow as 1787 the integral function must evaluate to a finite result over the entire range. 1788 </p> 1789<h6> 1790<a name="math_toolkit.lambert_w.h28"></a> 1791 <span class="phrase"><a name="math_toolkit.lambert_w.other_implementations"></a></span><a class="link" href="lambert_w.html#math_toolkit.lambert_w.other_implementations">Other 1792 implementations</a> 1793 </h6> 1794<p> 1795 The Lambert W has also been discussed in a <a href="http://lists.boost.org/Archives/boost/2016/09/230819.php" target="_top">Boost 1796 thread</a>. 1797 </p> 1798<p> 1799 This also gives link to a prototype version by which also gives complex results 1800 <code class="literal">(x < -exp(-1)</code>, about -0.367879). <a href="https://github.com/CzB404/lambert_w/" target="_top">Balazs 1801 Cziraki 2016</a> Physicist, PhD student at Eotvos Lorand University, ELTE 1802 TTK Institute of Physics, Budapest. has also produced a prototype C++ library 1803 that can compute the Lambert W function for floating point <span class="bold"><strong>and 1804 complex number types</strong></span>. This is not implemented here but might be 1805 completed in the future. 1806 </p> 1807<h5> 1808<a name="math_toolkit.lambert_w.h29"></a> 1809 <span class="phrase"><a name="math_toolkit.lambert_w.acknowledgements"></a></span><a class="link" href="lambert_w.html#math_toolkit.lambert_w.acknowledgements">Acknowledgements</a> 1810 </h5> 1811<div class="itemizedlist"><ul class="itemizedlist" style="list-style-type: disc; "> 1812<li class="listitem"> 1813 Thanks to Wolfram for use of their invaluable online Wolfram Alpha service. 1814 </li> 1815<li class="listitem"> 1816 Thanks for Mark Chapman for performing offline Wolfram computations. 1817 </li> 1818</ul></div> 1819<h5> 1820<a name="math_toolkit.lambert_w.h30"></a> 1821 <span class="phrase"><a name="math_toolkit.lambert_w.references"></a></span><a class="link" href="lambert_w.html#math_toolkit.lambert_w.references">References</a> 1822 </h5> 1823<div class="orderedlist"><ol class="orderedlist" type="1"> 1824<li class="listitem"> 1825 NIST Digital Library of Mathematical Functions. <a href="http://dlmf.nist.gov/4.13.F1" target="_top">http://dlmf.nist.gov/4.13.F1</a>. 1826 </li> 1827<li class="listitem"> 1828 <a href="http://www.orcca.on.ca/LambertW/" target="_top">Lambert W Poster</a>, 1829 R. M. Corless, G. H. Gonnet, D. E. G. Hare, D. J. Jeffery and D. E. Knuth, 1830 On the Lambert W function Advances in Computational Mathematics, Vol 5, 1831 (1996) pp 329-359. 1832 </li> 1833<li class="listitem"> 1834 <a href="https://people.sc.fsu.edu/~jburkardt/f_src/toms443/toms443.html" target="_top">TOMS443</a>, 1835 Andrew Barry, S. J. Barry, Patricia Culligan-Hensley, Algorithm 743: WAPR 1836 - A Fortran routine for calculating real values of the W-function,<br> 1837 ACM Transactions on Mathematical Software, Volume 21, Number 2, June 1995, 1838 pages 172-181.<br> BISECT approximates the W function using bisection 1839 (GNU licence). Original FORTRAN77 version by Andrew Barry, S. J. Barry, 1840 Patricia Culligan-Hensley, this version by C++ version by John Burkardt. 1841 </li> 1842<li class="listitem"> 1843 <a href="https://people.sc.fsu.edu/~jburkardt/f_src/toms743/toms743.html" target="_top">TOMS743</a> 1844 Fortran 90 (updated 2014). 1845 </li> 1846</ol></div> 1847<p> 1848 Initial guesses based on: 1849 </p> 1850<div class="orderedlist"><ol class="orderedlist" type="1"> 1851<li class="listitem"> 1852 R.M.Corless, G.H.Gonnet, D.E.G.Hare, D.J.Jeffrey, and D.E.Knuth, On the 1853 Lambert W function, Adv.Comput.Math., vol. 5, pp. 329 to 359, (1996). 1854 </li> 1855<li class="listitem"> 1856 D.A. Barry, J.-Y. Parlange, L. Li, H. Prommer, C.J. Cunningham, and F. 1857 Stagnitti. Analytical approximations for real values of the Lambert W-function. 1858 Mathematics and Computers in Simulation, 53(1), 95-103 (2000). 1859 </li> 1860<li class="listitem"> 1861 D.A. Barry, J.-Y. Parlange, L. Li, H. Prommer, C.J. Cunningham, and F. 1862 Stagnitti. Erratum to analytical approximations for real values of the 1863 Lambert W-function. Mathematics and Computers in Simulation, 59(6):543-543, 1864 2002. 1865 </li> 1866<li class="listitem"> 1867 C++ <a href="https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#c-cplusplus-language-support" target="_top">CUDA 1868 NVidia GPU C/C++ language support</a> version of Luu algorithm, <a href="https://github.com/thomasluu/plog/blob/master/plog.cu" target="_top">plog</a>. 1869 </li> 1870<li class="listitem"> 1871 <a href="http://discovery.ucl.ac.uk/1482128/1/Luu_thesis.pdf" target="_top">Thomas 1872 Luu, Thesis, University College London (2015)</a>, see routine 11, 1873 page 98 for Lambert W algorithm. 1874 </li> 1875<li class="listitem"> 1876 Having Fun with Lambert W(x) Function, Darko Veberic University of Nova 1877 Gorica, Slovenia IK, Forschungszentrum Karlsruhe, Germany, J. Stefan Institute, 1878 Ljubljana, Slovenia. 1879 </li> 1880<li class="listitem"> 1881 François Chapeau-Blondeau and Abdelilah Monir, Numerical Evaluation of the 1882 Lambert W Function and Application to Generation of Generalized Gaussian 1883 Noise With Exponent 1/2, IEEE Transactions on Signal Processing, 50(9) 1884 (2002) 2160 - 2165. 1885 </li> 1886<li class="listitem"> 1887 Toshio Fukushima, Precise and fast computation of Lambert W-functions without 1888 transcendental function evaluations, Journal of Computational and Applied 1889 Mathematics, 244 (2013) 77-89. 1890 </li> 1891<li class="listitem"> 1892 T.C. Banwell and A. Jayakumar, Electronic Letter, Feb 2000, 36(4), pages 1893 291-2. Exact analytical solution for current flow through diode with series 1894 resistance. <a href="https://doi.org/10.1049/el:20000301" target="_top">https://doi.org/10.1049/el:20000301</a> 1895 </li> 1896<li class="listitem"> 1897 Princeton Companion to Applied Mathematics, 'The Lambert-W function', Section 1898 1.3: Series and Generating Functions. 1899 </li> 1900<li class="listitem"> 1901 Cleve Moler, Mathworks blog <a href="https://blogs.mathworks.com/cleve/2013/09/02/the-lambert-w-function/#bfba4e2d-e049-45a6-8285-fe8b51d69ce7" target="_top">The 1902 Lambert W Function</a> 1903 </li> 1904<li class="listitem"> 1905 Digital Library of Mathematical Function, <a href="https://dlmf.nist.gov/4.13" target="_top">Lambert 1906 W function</a>. 1907 </li> 1908</ol></div> 1909</div> 1910<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr> 1911<td align="left"></td> 1912<td align="right"><div class="copyright-footer">Copyright © 2006-2019 Nikhar 1913 Agrawal, Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, 1914 Hubert Holin, Bruno Lalande, John Maddock, Jeremy Murphy, Matthew Pulver, Johan 1915 Råde, Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg, 1916 Daryle Walker and Xiaogang Zhang<p> 1917 Distributed under the Boost Software License, Version 1.0. (See accompanying 1918 file LICENSE_1_0.txt or copy at <a href="http://www.boost.org/LICENSE_1_0.txt" target="_top">http://www.boost.org/LICENSE_1_0.txt</a>) 1919 </p> 1920</div></td> 1921</tr></table> 1922<hr> 1923<div class="spirit-nav"> 1924<a accesskey="p" href="jacobi/jacobi_sn.html"><img src="../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../special.html"><img src="../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../index.html"><img src="../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="zetas.html"><img src="../../../../../doc/src/images/next.png" alt="Next"></a> 1925</div> 1926</body> 1927</html> 1928
