libs/math/doc/html/math_toolkit/lambert_w.html

   1 <html>
   2 <head>
   3 <meta http-equiv="Content-Type" content="text/html; charset=US-ASCII">
   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.11.0">
   8 <link rel="up" href="../special.html" title="Chapter&#160;8.&#160;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">&lt;</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">&gt;</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">&lt;</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">&gt;</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">&lt;</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">&gt;</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">&lt;</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">&gt;</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">&lt;</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">&gt;</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">&lt;</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&#160;20.&#160;Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">&gt;</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&#160;20.&#160;Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">&amp;);</span>         <span class="comment">// W0 with policy.</span>
  49   <span class="keyword">template</span> <span class="special">&lt;</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&#160;20.&#160;Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">&gt;</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&#160;20.&#160;Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">&amp;);</span>        <span class="comment">// W-1 with policy.</span>
  51   <span class="keyword">template</span> <span class="special">&lt;</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&#160;20.&#160;Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">&gt;</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&#160;20.&#160;Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">&amp;);</span>   <span class="comment">// W0 derivative with policy.</span>
  53   <span class="keyword">template</span> <span class="special">&lt;</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&#160;20.&#160;Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">&gt;</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&#160;20.&#160;Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">&amp;);</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, &#8734;), 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> &#8773; <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>
 103 with 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&#160;20.&#160;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&#160;20.&#160;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">&lt;</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">&gt;</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">&lt;</span><span class="keyword">double</span><span class="special">&gt;::</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">&lt;&lt;</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">showpoint</span> <span class="special">&lt;&lt;</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">&lt;&lt;</span> <span class="string">"lambert_w0(z) = "</span> <span class="special">&lt;&lt;</span> <span class="identifier">r</span> <span class="special">&lt;&lt;</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">&lt;</span><span class="identifier">T</span><span class="special">&gt;::</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">&lt;&lt;</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">&lt;&lt;</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">&lt;&lt;</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">&lt;</span><span class="keyword">double</span><span class="special">&gt;::</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">&lt;&lt;</span> <span class="string">"lambert_w0(10) = "</span> <span class="special">&lt;&lt;</span> <span class="identifier">r</span> <span class="special">&lt;&lt;</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">&lt;&lt;</span> <span class="string">"lambert_w0(10) = "</span> <span class="special">&lt;&lt;</span> <span class="identifier">rp</span> <span class="special">&lt;&lt;</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">&lt;&lt;</span> <span class="string">"lambert_w0(10) = "</span> <span class="special">&lt;&lt;</span> <span class="identifier">rr</span> <span class="special">&lt;&lt;</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">&lt;&lt;</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">&lt;&lt;</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">&lt;&lt;</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">&lt;&lt;</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">&lt;&lt;</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">&lt;&lt;</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">&lt;&lt;</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">&lt;&lt;</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">&lt;&lt;</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">&lt;&lt;</span> <span class="string">"lambert_w0(0.9) = "</span> <span class="special">&lt;&lt;</span> <span class="identifier">lambert_w0</span><span class="special">(</span><span class="keyword">static_cast</span><span class="special">&lt;</span><span class="keyword">double</span><span class="special">&gt;(</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">&lt;&lt;</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">&lt;&lt;</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">&lt;&lt;</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">&lt;&lt;</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">&lt;</span>
 291   <span class="identifier">domain_error</span><span class="special">&lt;</span><span class="identifier">throw_on_error</span><span class="special">&gt;,</span>
 292   <span class="identifier">overflow_error</span><span class="special">&lt;</span><span class="identifier">ignore_error</span><span class="special">&gt;</span> <span class="comment">// possibly unwise?</span>
 293 <span class="special">&gt;</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">&lt;</span><span class="keyword">double</span><span class="special">&gt;::</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">&lt;&lt;</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">showpoint</span> <span class="special">&lt;&lt;</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">&lt;&lt;</span> <span class="string">"Lambert W ("</span> <span class="special">&lt;&lt;</span> <span class="identifier">z</span> <span class="special">&lt;&lt;</span> <span class="string">") = "</span> <span class="special">&lt;&lt;</span> <span class="identifier">lambert_w0</span><span class="special">(</span><span class="identifier">z</span><span class="special">)</span> <span class="special">&lt;&lt;</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">&lt;&lt;</span> <span class="string">"\nLambert W ("</span> <span class="special">&lt;&lt;</span> <span class="identifier">z</span> <span class="special">&lt;&lt;</span> <span class="string">", my_throw_policy()) = "</span>
 303   <span class="special">&lt;&lt;</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">&lt;&lt;</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">&lt;</span><span class="identifier">RealType</span><span class="special">&gt;(&lt;</span><span class="identifier">RealType</span><span class="special">&gt;):</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">&lt;=</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">&lt;&lt;</span> <span class="string">"lambert_wm1(+1.) = "</span> <span class="special">&lt;&lt;</span> <span class="identifier">r</span> <span class="special">&lt;&lt;</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) = &#956; V<sub>T</sub>/ R <sub>S</sub> &#8228; W<sub>0</sub>(I<sub>0</sub> R<sub>S</sub> / (&#956; 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">&lt;&lt;</span> <span class="string">"V thermal "</span> <span class="special">&lt;&lt;</span> <span class="identifier">vt</span> <span class="special">&lt;&lt;</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">&lt;&lt;</span> <span class="string">"Isat = "</span> <span class="special">&lt;&lt;</span> <span class="identifier">isat</span> <span class="special">&lt;&lt;</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">&lt;&lt;</span> <span class="string">"voltage = "</span> <span class="special">&lt;&lt;</span> <span class="identifier">v</span> <span class="special">&lt;&lt;</span> <span class="string">", current = "</span> <span class="special">&lt;&lt;</span> <span class="identifier">icalc</span> <span class="special">&lt;&lt;</span> <span class="string">", "</span> <span class="special">&lt;&lt;</span> <span class="identifier">log</span><span class="special">(</span><span class="identifier">icalc</span><span class="special">)</span> <span class="special">&lt;&lt;</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 &lt;&lt; "nu * vt / rsat = " &lt;&lt; i &lt;&lt; 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 &lt;&lt; "isat * rsat / (nu * vt) = " &lt;&lt; x &lt;&lt; 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 &lt;&lt; "(v + isat * rsat) / (nu * vt) = " &lt;&lt; eterm &lt;&lt; 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 &lt;&lt; "exp(eterm) = " &lt;&lt; e &lt;&lt; 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 &lt;&lt; "w0 = " &lt;&lt; w0 &lt;&lt; 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">&lt;</span><span class="identifier">cpp_bin_float_50</span><span class="special">&gt;::</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">&lt;&lt;</span> <span class="string">"lambert_w0(z) cpp_bin_float_50  = "</span> <span class="special">&lt;&lt;</span> <span class="identifier">r</span> <span class="special">&lt;&lt;</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">&lt;</span><span class="keyword">double</span><span class="special">&gt;::</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">&lt;&lt;</span> <span class="string">"lambert_w0(z) static_cast from cpp_bin_float_50  = "</span>
 506    <span class="special">&lt;&lt;</span> <span class="keyword">static_cast</span><span class="special">&lt;</span><span class="keyword">double</span><span class="special">&gt;(</span><span class="identifier">r</span><span class="special">)</span> <span class="special">&lt;&lt;</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">&lt;&lt;</span> <span class="string">"bits different from Wolfram = "</span>
 510   <span class="special">&lt;&lt;</span> <span class="keyword">static_cast</span><span class="special">&lt;</span><span class="keyword">int</span><span class="special">&gt;(</span><span class="identifier">float_distance</span><span class="special">(</span><span class="keyword">static_cast</span><span class="special">&lt;</span><span class="keyword">double</span><span class="special">&gt;(</span><span class="identifier">r</span><span class="special">),</span> <span class="number">1.7455280027406994</span><span class="special">))</span>
 511   <span class="special">&lt;&lt;</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">&lt;&lt;</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">showpoint</span> <span class="special">&lt;&lt;</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">&lt;</span><span class="keyword">double</span><span class="special">&gt;::</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">&lt;</span><span class="keyword">double</span><span class="special">&gt;(</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">&lt;</span><span class="identifier">cpp_bin_float_50</span><span class="special">&gt;::</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">&lt;&lt;</span> <span class="string">"Reference Lambert W ("</span> <span class="special">&lt;&lt;</span> <span class="identifier">z</span> <span class="special">&lt;&lt;</span> <span class="string">") =\n                                              "</span>
 539   <span class="special">&lt;&lt;</span> <span class="identifier">w50</span> <span class="special">&lt;&lt;</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">&lt;</span><span class="keyword">double</span><span class="special">&gt;::</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">&lt;</span><span class="keyword">double</span><span class="special">&gt;(</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">&lt;&lt;</span> <span class="string">"Reference Lambert W ("</span> <span class="special">&lt;&lt;</span> <span class="identifier">z</span> <span class="special">&lt;&lt;</span> <span class="string">") =    "</span> <span class="special">&lt;&lt;</span> <span class="identifier">wr</span> <span class="special">&lt;&lt;</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">&lt;&lt;</span> <span class="string">"Rat/poly Lambert W  ("</span> <span class="special">&lt;&lt;</span> <span class="identifier">z</span> <span class="special">&lt;&lt;</span> <span class="string">")  =   "</span> <span class="special">&lt;&lt;</span> <span class="identifier">lambert_w0</span><span class="special">(</span><span class="identifier">z</span><span class="special">)</span> <span class="special">&lt;&lt;</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">&lt;&lt;</span> <span class="string">"Halley Step Lambert W ("</span> <span class="special">&lt;&lt;</span> <span class="identifier">z</span> <span class="special">&lt;&lt;</span> <span class="string">") =  "</span> <span class="special">&lt;&lt;</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">&lt;&lt;</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">&lt;&lt;</span> <span class="string">"absolute difference from Halley step = "</span> <span class="special">&lt;&lt;</span> <span class="identifier">w</span> <span class="special">-</span> <span class="identifier">ww</span> <span class="special">&lt;&lt;</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">&lt;&lt;</span> <span class="string">"relative difference from Halley step = "</span> <span class="special">&lt;&lt;</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">&lt;&lt;</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">&lt;&lt;</span> <span class="string">"epsilon difference from Halley step  = "</span> <span class="special">&lt;&lt;</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">&lt;&lt;</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">&lt;&lt;</span> <span class="string">"epsilon for float =                    "</span> <span class="special">&lt;&lt;</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special">&lt;</span><span class="keyword">double</span><span class="special">&gt;::</span><span class="identifier">epsilon</span><span class="special">()</span> <span class="special">&lt;&lt;</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">&lt;&lt;</span> <span class="string">"bits different from Halley step  =     "</span> <span class="special">&lt;&lt;</span> <span class="keyword">static_cast</span><span class="special">&lt;</span><span class="keyword">int</span><span class="special">&gt;(</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">&lt;&lt;</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">&lt;&lt;</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">showpoint</span> <span class="special">&lt;&lt;</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">&lt;</span><span class="keyword">double</span><span class="special">&gt;::</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">&lt;</span><span class="keyword">double</span><span class="special">&gt;(</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">&lt;</span><span class="identifier">cpp_bin_float_50</span><span class="special">&gt;::</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">&lt;&lt;</span> <span class="string">"Reference Lambert W-1 ("</span> <span class="special">&lt;&lt;</span> <span class="identifier">z</span> <span class="special">&lt;&lt;</span> <span class="string">") =\n                                                  "</span>
 586   <span class="special">&lt;&lt;</span> <span class="identifier">wm1_50</span> <span class="special">&lt;&lt;</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">&lt;</span><span class="keyword">double</span><span class="special">&gt;::</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">&lt;</span><span class="keyword">double</span><span class="special">&gt;(</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">&lt;&lt;</span> <span class="string">"Reference Lambert W-1 ("</span> <span class="special">&lt;&lt;</span> <span class="identifier">z</span> <span class="special">&lt;&lt;</span> <span class="string">") =    "</span> <span class="special">&lt;&lt;</span> <span class="identifier">wr</span> <span class="special">&lt;&lt;</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">&lt;&lt;</span> <span class="string">"Rat/poly Lambert W-1 ("</span> <span class="special">&lt;&lt;</span> <span class="identifier">z</span> <span class="special">&lt;&lt;</span> <span class="string">")  =    "</span> <span class="special">&lt;&lt;</span> <span class="identifier">lambert_wm1</span><span class="special">(</span><span class="identifier">z</span><span class="special">)</span> <span class="special">&lt;&lt;</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">&lt;&lt;</span> <span class="string">"Halley Step Lambert W ("</span> <span class="special">&lt;&lt;</span> <span class="identifier">z</span> <span class="special">&lt;&lt;</span> <span class="string">") =    "</span> <span class="special">&lt;&lt;</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">&lt;&lt;</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">&lt;&lt;</span> <span class="string">"absolute difference from Halley step = "</span> <span class="special">&lt;&lt;</span> <span class="identifier">w</span> <span class="special">-</span> <span class="identifier">ww</span> <span class="special">&lt;&lt;</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">&lt;&lt;</span> <span class="string">"relative difference from Halley step = "</span> <span class="special">&lt;&lt;</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">&lt;&lt;</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">&lt;&lt;</span> <span class="string">"epsilon difference from Halley step  = "</span> <span class="special">&lt;&lt;</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">&lt;&lt;</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">&lt;&lt;</span> <span class="string">"epsilon for float =                    "</span> <span class="special">&lt;&lt;</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special">&lt;</span><span class="keyword">double</span><span class="special">&gt;::</span><span class="identifier">epsilon</span><span class="special">()</span> <span class="special">&lt;&lt;</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">&lt;&lt;</span> <span class="string">"bits different from Halley step  =     "</span> <span class="special">&lt;&lt;</span> <span class="keyword">static_cast</span><span class="special">&lt;</span><span class="keyword">int</span><span class="special">&gt;(</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">&lt;&lt;</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">&lt;</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>, &#8734;). 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">&lt;</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">&lt;</span><span class="identifier">T</span><span class="special">&gt;::</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           -&#8734; (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">&lt;</span><span class="identifier">T</span><span class="special">&gt;::</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">&lt;</span>
 707           <span class="special">-</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special">&lt;</span><span class="identifier">T</span><span class="special">&gt;::</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">&lt;</span><span class="identifier">T</span><span class="special">&gt;::</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">&lt;</span><span class="keyword">double</span><span class="special">&gt;::</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&lt;RealType&gt;(&lt;RealType&gt;):
 714           Argument z = -4.9406564584124654e-324 is too small (z &lt; -std::numeric_limits&lt;T&gt;::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">&lt;&gt;::</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">&lt;</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">&gt;</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">&lt;</span><span class="identifier">define</span><span class="special">&gt;</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 &gt; 64, z &gt; -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&#246;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 &amp; 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)
 873 if (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">&lt;</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, nrealy 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&#160;8.73.&#160;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 &lt; 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&#160;8.74.&#160;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&#160;8.75.&#160;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">&lt;</span><span class="identifier">cpp_bin_float_quad</span><span class="special">&gt;::</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">&lt;&gt;::</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 &gt; 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">&lt;</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">&lt;</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">&lt;</span><span class="identifier">T</span><span class="special">&gt;::</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">&lt;</span><span class="identifier">T</span><span class="special">&gt;::</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       achive 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 contructors 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 evalutations 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&#178;) == &#8730;(2&#960;), 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 &lt; -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&#231;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 &#169; 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&#229;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>