From de4d578fc56ff548d1222834fc0555b2ebd8122d Mon Sep 17 00:00:00 2001 From: pabristow Date: Mon, 19 Feb 2018 16:54:52 +0000 Subject: [PATCH] deconflict lambert_w.hpp --- doc/background/remez.qbk | 76 +++---- doc/html/index.html | 2 +- doc/html/indexes/s01.html | 10 +- doc/html/indexes/s02.html | 2 +- doc/html/indexes/s03.html | 2 +- doc/html/indexes/s04.html | 2 +- doc/html/indexes/s05.html | 12 +- doc/html/math_toolkit/conventions.html | 2 +- doc/html/math_toolkit/navigation.html | 2 +- .../pol_tutorial/what_is_a_policy.html | 8 +- doc/html/math_toolkit/remez.html | 2 +- doc/math.qbk | 9 +- doc/sf/lambert_w.qbk | 115 ++++++---- example/lambert_w_simple_examples.cpp | 7 - .../math/special_functions/lambert_w.hpp | 197 ++++++++++-------- .../lambert_w_lookup_table.ipp | 154 +++++++------- .../math/special_functions/lambert_w_pb.hpp | 67 +++--- test/test_lambert_w.cpp | 27 +-- 18 files changed, 384 insertions(+), 312 deletions(-) diff --git a/doc/background/remez.qbk b/doc/background/remez.qbk index 6dd2718d2..f7ebe4019 100644 --- a/doc/background/remez.qbk +++ b/doc/background/remez.qbk @@ -8,9 +8,9 @@ should consult your favorite textbook. Imagine that you want to approximate some function f(x) by way of a rational function R(x), where R(x) may be either a polynomial P(x) or a ratio of two -polynomials P(x)/Q(x) (a rational function). Initially we'll concentrate on the -polynomial case, as it's by far the easier to deal with, later we'll extend -to the full rational function case. +polynomials P(x)/Q(x) (a rational function). Initially we'll concentrate on the +polynomial case, as it's by far the easier to deal with, later we'll extend +to the full rational function case. We want to find the "best" rational approximation, where "best" is defined to be the approximation that has the least deviation @@ -61,11 +61,11 @@ the range \[-1, 1\]. Before we can begin the Remez method, we must obtain an initial value for the location of the extrema of the error function. We could "guess" -these, but a much closer first approximation can be obtained by first +these, but a much closer first approximation can be obtained by first constructing an interpolated polynomial approximation to f(x). In order to obtain the N+1 coefficients of the interpolated polynomial -we need N+1 points (x[sub 0]...x[sub N]): with our interpolated form +we need N+1 points (x[sub 0]...x[sub N]): with our interpolated form passing through each of those points that yields N+1 simultaneous equations: @@ -73,12 +73,12 @@ f(x[sub i]) = P(x[sub i]) = c[sub 0] + c[sub 1]x[sub i] ... + c[sub N]x[sub i][s Which can be solved for the coefficients c[sub 0]...c[sub N] in P(x). -Obviously this is not a minimax solution, indeed our only guarantee is that f(x) and +Obviously this is not a minimax solution, indeed our only guarantee is that f(x) and P(x) touch at N+1 locations, away from those points the error may be arbitrarily large. However, we would clearly like this initial approximation to be as close to f(x) as possible, and it turns out that using the zeros of an orthogonal polynomial -as the initial interpolation points is a good choice. In our example we'll use the -zeros of a Chebyshev polynomial as these are particularly easy to calculate, +as the initial interpolation points is a good choice. In our example we'll use the +zeros of a Chebyshev polynomial as these are particularly easy to calculate, interpolating for a polynomial of degree 4, and measuring /relative error/ we get the following error function: @@ -86,7 +86,7 @@ we get the following error function: Which has a peak relative error of 1.2x10[super -3]. -While this is a pretty good approximation already, judging by the +While this is a pretty good approximation already, judging by the shape of the error function we can clearly do better. Before starting on the Remez method propper, we have one more step to perform: locate all the extrema of the error function, and store @@ -101,7 +101,7 @@ achieve, and typically is very close to the minimax solution. However, if we want to optimise for /relative error/, or if the approximation is a rational function, then the initial Chebyshev solution can be quite far -from the ideal minimax solution. +from the ideal minimax solution. A more technical discussion of the theory involved can be found in this [@http://math.fullerton.edu/mathews/n2003/ChebyshevPolyMod.html online course].] @@ -118,7 +118,7 @@ To obtain the error term E, and the coefficients of the polynomial P(x). This gives us a new approximation to f(x) that has the same error /E/ at each of the control points, and whose error function ['alternates in sign] -at the control points. This is still not necessarily the minimax +at the control points. This is still not necessarily the minimax solution though: since the control points may not be at the extrema of the error function. After this first step here's what our approximation's error function looks like: @@ -131,10 +131,10 @@ dropped to 5.6x10[super -4]. [h4 Remez Step 2] -The second step is to locate the extrema of the new approximation, which we do +The second step is to locate the extrema of the new approximation, which we do in two stages: first, since the error function changes sign at each control point, we must have N+1 roots of the error function located between -each pair of N+2 control points. Once these roots are found by standard root finding +each pair of N+2 control points. Once these roots are found by standard root finding techniques, we know that N extrema are bracketed between each pair of roots, plus two more between the endpoints of the range and the first and last roots. The N+2 extrema can then be found using standard function minimisation techniques. @@ -195,10 +195,10 @@ polynomial alternatives. For example, if we takes our previous example of an approximation to e[super x], we obtained 5x10[super -4] accuracy with an order 4 polynomial. If we move two of the unknowns into the denominator to give a pair of order 2 polynomials, and re-minimise, then the peak relative error drops -to 8.7x10[super -5]. That's a 5 fold increase in accuracy, for the same number +to 8.7x10[super -5]. That's a 5 fold increase in accuracy, for the same number of terms overall. -[h4 Practical Considerations] +[h4:remez_practical Practical Considerations] Most treatises on approximation theory stop at this point. However, from a practical point of view, most of the work involves finding the right @@ -211,12 +211,12 @@ f(x) = R(x) But this will converge to a useful approximation only if f(x) is smooth. In addition round-off errors when evaluating the rational form mean that this -will never get closer than within a few epsilon of machine precision. +will never get closer than within a few epsilon of machine precision. Therefore this form of direct approximation is often reserved for situations where we want efficiency, rather than accuracy. The first step in improving the situation is generally to split f(x) into -a dominant part that we can compute accurately by another method, and a +a dominant part that we can compute accurately by another method, and a slowly changing remainder which can be approximated by a rational approximation. We might be tempted to write: @@ -235,7 +235,7 @@ to the constant /c/, that error will effectively get wiped out when R(x) is adde The difficult part is obviously finding the right g(x) to extract from your function: often the asymptotic behaviour of the function will give a clue, so -for example the function __erfc becomes proportional to +for example the function __erfc becomes proportional to e[super -x[super 2]]\/x as x becomes large. Therefore using: erfc(z) = (C + R(x)) e[super -x[super 2]]/x @@ -262,15 +262,15 @@ the roots, the approximation is nonetheless useless. Assuming that the approximation does not have any fatal errors, and that the only issue is converging adequately on the minimax solution, the aim is to get as close as possible to the minimax solution before beginning the Remez method. -Using the zeros of a Chebyshev polynomial for the initial interpolation is a +Using the zeros of a Chebyshev polynomial for the initial interpolation is a good start, but may not be ideal when dealing with relative errors and\/or rational (rather than polynomial) approximations. One approach is to skew the initial interpolation points to one end: for example if we raise the -roots of the Chebyshev polynomial to a positive power greater than 1 -then the roots will be skewed towards the middle of the \[-1,1\] interval, +roots of the Chebyshev polynomial to a positive power greater than 1 +then the roots will be skewed towards the middle of the \[-1,1\] interval, while a positive power less than one will skew them towards either end. More usefully, if we initially rescale the -points over \[0,1\] and then raise to a positive power, we can skew them to the left +points over \[0,1\] and then raise to a positive power, we can skew them to the left or right. Returning to our example of e[super x][space] over \[-1,1\], the initial interpolated form was some way from the minimax solution: @@ -287,7 +287,7 @@ our desired minimax solution (5x10[super -4]). [h4 Remez Method Checklist] -The following lists some of the things to check if the Remez method goes wrong, +The following lists some of the things to check if the Remez method goes wrong, it is by no means an exhaustive list, but is provided in the hopes that it will prove useful. @@ -312,7 +312,7 @@ location, try starting from a different location. * If a rational function won't converge, one can minimise a polynomial (which presents no problems), then rotate one term from the numerator to the denominator and minimise again. In theory one can continue moving -terms one at a time from numerator to denominator, and then re-minimising, +terms one at a time from numerator to denominator, and then re-minimising, retaining the last set of control points at each stage. * Try using a smaller interval. It may also be possible to optimise over one (small) interval, rescale the control points over a larger interval, @@ -325,41 +325,41 @@ over, say \[0, b\], for some smallish value /b/. The original references for the Remez Method and it's extension to rational functions are unfortunately in Russian: -Remez, E.Ya., ['Fundamentals of numerical methods for Chebyshev approximations], +Remez, E.Ya., ['Fundamentals of numerical methods for Chebyshev approximations], "Naukova Dumka", Kiev, 1969. -Remez, E.Ya., Gavrilyuk, V.T., ['Computer development of certain approaches -to the approximate construction of solutions of Chebyshev problems +Remez, E.Ya., Gavrilyuk, V.T., ['Computer development of certain approaches +to the approximate construction of solutions of Chebyshev problems nonlinearly depending on parameters], Ukr. Mat. Zh. 12 (1960), 324-338. -Gavrilyuk, V.T., ['Generalization of the first polynomial algorithm of -E.Ya.Remez for the problem of constructing rational-fractional +Gavrilyuk, V.T., ['Generalization of the first polynomial algorithm of +E.Ya.Remez for the problem of constructing rational-fractional Chebyshev approximations], Ukr. Mat. Zh. 16 (1961), 575-585. Some English language sources include: -Fraser, W., Hart, J.F., ['On the computation of rational approximations +Fraser, W., Hart, J.F., ['On the computation of rational approximations to continuous functions], Comm. of the ACM 5 (1962), 401-403, 414. -Ralston, A., ['Rational Chebyshev approximation by Remes' algorithms], +Ralston, A., ['Rational Chebyshev approximation by Remes' algorithms], Numer.Math. 7 (1965), no. 4, 322-330. -A. Ralston, ['Rational Chebyshev approximation, Mathematical -Methods for Digital Computers v. 2] (Ralston A., Wilf H., eds.), +A. Ralston, ['Rational Chebyshev approximation, Mathematical +Methods for Digital Computers v. 2] (Ralston A., Wilf H., eds.), Wiley, New York, 1967, pp. 264-284. Hart, J.F. e.a., ['Computer approximations], Wiley, New York a.o., 1968. -Cody, W.J., Fraser, W., Hart, J.F., ['Rational Chebyshev approximation +Cody, W.J., Fraser, W., Hart, J.F., ['Rational Chebyshev approximation using linear equations], Numer.Math. 12 (1968), 242-251. -Cody, W.J., ['A survey of practical rational and polynomial +Cody, W.J., ['A survey of practical rational and polynomial approximation of functions], SIAM Review 12 (1970), no. 3, 400-423. -Barrar, R.B., Loeb, H.J., ['On the Remez algorithm for non-linear +Barrar, R.B., Loeb, H.J., ['On the Remez algorithm for non-linear families], Numer.Math. 15 (1970), 382-391. -Dunham, Ch.B., ['Convergence of the Fraser-Hart algorithm for rational +Dunham, Ch.B., ['Convergence of the Fraser-Hart algorithm for rational Chebyshev approximation], Math. Comp. 29 (1975), no. 132, 1078-1082. G. L. Litvinov, ['Approximate construction of rational @@ -368,7 +368,7 @@ Russian Journal of Mathematical Physics, vol.1, No. 3, 1994. [endsect][/section:remez The Remez Method] -[/ +[/ Copyright 2006 John Maddock and Paul A. Bristow. Distributed under the Boost Software License, Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at diff --git a/doc/html/index.html b/doc/html/index.html index 9cdba317b..cf86d0e63 100644 --- a/doc/html/index.html +++ b/doc/html/index.html @@ -120,7 +120,7 @@ This manual is also available in -

Last revised: December 18, 2017 at 18:05:20 GMT

+

Last revised: February 01, 2018 at 15:50:54 GMT

diff --git a/doc/html/indexes/s01.html b/doc/html/indexes/s01.html index 8f19b405c..6de54b099 100644 --- a/doc/html/indexes/s01.html +++ b/doc/html/indexes/s01.html @@ -24,7 +24,7 @@

-Function Index

+Function Index

2 4 A B C D E F G H I J K L M N O P Q R S T U V W X Y Z

@@ -420,6 +420,10 @@
  • +

    close

    +
    +
    • +

  • complement

  • +

    step

    +
    +
    • +

  • subtraction

    • diff --git a/doc/html/indexes/s02.html b/doc/html/indexes/s02.html index 969845731..b3279a834 100644 --- a/doc/html/indexes/s02.html +++ b/doc/html/indexes/s02.html @@ -24,7 +24,7 @@

    -Class Index

    +Class Index

    A B C D E F G H I L M N O P Q R S T U W

    diff --git a/doc/html/indexes/s03.html b/doc/html/indexes/s03.html index 99dca7502..4fa43be01 100644 --- a/doc/html/indexes/s03.html +++ b/doc/html/indexes/s03.html @@ -24,7 +24,7 @@

    -Typedef Index

    +Typedef Index

    A B C D E F G H I L N O P R S T U V W

    diff --git a/doc/html/indexes/s04.html b/doc/html/indexes/s04.html index 77f3997ec..6f07c37ca 100644 --- a/doc/html/indexes/s04.html +++ b/doc/html/indexes/s04.html @@ -24,7 +24,7 @@

    -Macro Index

    +Macro Index

    B F

    diff --git a/doc/html/indexes/s05.html b/doc/html/indexes/s05.html index 2dcaf5d10..edac1cc56 100644 --- a/doc/html/indexes/s05.html +++ b/doc/html/indexes/s05.html @@ -23,7 +23,7 @@

    -Index

    +Index

    2 4 A B C D E F G H I J K L M N O P Q R S T U V W X Y Z

    @@ -1662,6 +1662,10 @@

  • chi_squared_distribution

  • +

    close

    +
    +
    • +

  • Comparing Different Compilers

    • @@ -4851,6 +4855,7 @@

  • refinement

  • required

  • small

    • +

  • step

  • W

  • zero

    • @@ -5074,6 +5079,7 @@
  • +

    step

    +
    +
    • +

  • Students t Distribution

    • accuracy

      • diff --git a/doc/html/math_toolkit/conventions.html b/doc/html/math_toolkit/conventions.html index 84747b896..8ff3da6d2 100644 --- a/doc/html/math_toolkit/conventions.html +++ b/doc/html/math_toolkit/conventions.html @@ -27,7 +27,7 @@ Document Conventions

    - +

    This documentation aims to use of the following naming and formatting conventions. diff --git a/doc/html/math_toolkit/navigation.html b/doc/html/math_toolkit/navigation.html index 955b0fcb3..2d5689a3e 100644 --- a/doc/html/math_toolkit/navigation.html +++ b/doc/html/math_toolkit/navigation.html @@ -27,7 +27,7 @@ Navigation

    - +

    Boost.Math documentation is provided in both HTML and PDF formats. diff --git a/doc/html/math_toolkit/pol_tutorial/what_is_a_policy.html b/doc/html/math_toolkit/pol_tutorial/what_is_a_policy.html index cb7e338b0..b7cf9df0c 100644 --- a/doc/html/math_toolkit/pol_tutorial/what_is_a_policy.html +++ b/doc/html/math_toolkit/pol_tutorial/what_is_a_policy.html @@ -62,13 +62,13 @@ of user-selected policies (sometimes called a type-list).

    - Many policy + Over a dozen policy defaults are provided, so most of the time you can ignore the policy - framework, but you can overwrite with your own policies to give detailed - control, for example: + framework, but you can overwrite the defaults with your own policies to give + detailed control, for example: 

    using namespace boost::math::policies;
-//
+
 // Define a policy that sets ::errno on overflow,
 // and does not promote double to long double internally,
 // and only aims for precision of only 3 decimal digits,
diff --git a/doc/html/math_toolkit/remez.html b/doc/html/math_toolkit/remez.html
index 118bee043..76a8bb907 100644
--- a/doc/html/math_toolkit/remez.html
+++ b/doc/html/math_toolkit/remez.html
@@ -297,7 +297,7 @@
     
    
 
    
 
-      Practical
+      Practical
       Considerations
     
    
 
    
diff --git a/doc/math.qbk b/doc/math.qbk
index 9165e08cb..d3288d3ea 100644
--- a/doc/math.qbk
+++ b/doc/math.qbk
@@ -113,7 +113,7 @@ and use the function's name as the link text.]
 [def __brent_minima_example [link math_toolkit.roots.brent_minima.example Brent's method example]]
 [def __newton [link math_toolkit.roots.roots_deriv.newton Newton-Raphson iteration]]
 [def __halley [link math_toolkit.roots.roots_deriv.halley Halley]]
-[def __schroder [link math_toolkit.roots.roots_deriv.schroder Schr'''ö'''der]]
+[def __schroder [link math_toolkit.roots.roots_deriv.schroder Schr'''ö'''der]] [/ link should really be __schroeder]
 [def __brent_minima_example [link math_toolkit.roots.brent_minima.example Brent minima finding example]]
 [def __bisection [link math_toolkit.roots.roots_noderiv.bisect bisection]]
 [def __bisect [link math_toolkit.roots.roots_noderiv.bisect bisect]]
@@ -222,7 +222,7 @@ and use the function's name as the link text.]
 [def __lambert_w_precision[link math_toolkit.lambert_w.precision precision]]
 [def __lambert_w_examples [link math_toolkit.lambert_w.examples examples]]
 [def __lambert_w_wm1_near_zero [link math_toolkit.lambert_w.wm1_near_zero wm1_near_zero]]
-
+[def __lambert_w_references [link math_toolkit.lambert_w.references references]]
 
 
 [/Airy Functions]
@@ -388,6 +388,8 @@ and use the function's name as the link text.]
 
 [def __real_concept [link math_toolkit.real_concepts  real concept]]
 [def __next_float [link math_toolkit.next_float Adjacent Floating-Point Values]]
+[def __remez_practical [link math_toolkit.remez.practical remez_practical]]
+
 
 [/ External links]
 [def __gsl [@http://www.gnu.org/software/gsl/ GSL-1.9]]
@@ -443,7 +445,8 @@ and use the function's name as the link text.]
 [def __Lambert_W [@http://en.wikipedia.org/wiki/Lambert_W_function Lambert W function]]
 [def __CUDA [@https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#c-cplusplus-language-support CUDA NVidia GPU C/C++ language support]]
 [def __Luu_thesis [@http://discovery.ucl.ac.uk/1482128/1/Luu_thesis.pdf Thomas Luu, Thesis, University College London (2015)]]
-
+[def __rational_function [@https://en.wikipedia.org/wiki/Rational_function rational function]]
+[def __remez [@http://en.wikipedia.org/wiki/Remez_algorithm Remez algorithm]]
 
 [/ Some composite templates]
 [template super[x]''''''[x]'''''']
diff --git a/doc/sf/lambert_w.qbk b/doc/sf/lambert_w.qbk
index 13e0a2ac5..0c07ec710 100644
--- a/doc/sf/lambert_w.qbk
+++ b/doc/sf/lambert_w.qbk
@@ -127,7 +127,7 @@ at least to thousands of decimal digits, but not shown here for brevity.
 See [@../../example/lambert_w_simple_examples.cpp lambert_w_simple_examples.cpp]
 ).
 
-Policies can be used to control precision (a tiny less precision gives a significant speedup):
+Policies can be used to control precision (a tiny bit less precision gives a significant speed-up):
 
 [lambert_w_simple_examples_precision_policies]
 
@@ -139,8 +139,8 @@ and to control what action to take on errors:
 
 [lambert_w_simple_examples_error_message_1]
 
-Show error reporting if passing a value to `lambert_w0` that is out of range,
-and probably was meant to be passed to `lambert_wm1` instead.
+Show error reporting if a value is passed to `lambert_w0` that is out of range,
+and was probably meant to be passed to `lambert_wm1` instead.
 
 [lambert_w_simple_examples_out_of_range]
 
@@ -208,7 +208,7 @@ The symbolic algebra program __Maple also computes Lambert W to an arbitrary pre
 
 [h4:precision Controlling the compromise between Precision and Speed]
 
-The compromise between precision and speed can be controlled using __precision_policy.
+This implementation allows the compromise between precision and speed to be controlled using __precision_policy.
 
 Using the default policy, all the functions usually return values within a few __ulp
 for the floating-point type, except for very small arguments very near zero,
@@ -216,10 +216,10 @@ and for arguments very close to the singularity at the branch point.
 
 By default, this implementation provides the best possible precision
 (adding a __halley refinement, but without using __multiprecision higher precision types)
-making it nearly twice as slow as Fukushima's algorithm,
+making it slower than Fukushima's or other algorithms,
 (Boost.Math generally prefers accuracy over speed).
 However, it will return at intermediate stages if the precision specified in the policy has been met,
-usually at least halving the execution time.
+usually at least halving the execution time, and so is, we believe, the fastest algorithm now available.
 
 It is convenient to define some `using` statements when using __policy_section to control precision.
 
@@ -229,7 +229,7 @@ It is convenient to define some `using` statements when using __policy_section t
     using boost::math::policies::digits10; // Precision as decimal digits.
 
 [tip Because some macros and metaprogramming are used by the implementation of policies,
-some ways of expressing these items may confuse the compiler and Visual Studio Intellisense.]
+some ways of naming these items may confuse Visual Studio Intellisense and the compiler.]
 
 Precision targets can be specified in bits using `digits2` or decimal using `digits10`.
 These are related by `digits10 = log10(2) * digits2 where log10(2) = 0.30103`.
@@ -271,10 +271,10 @@ The final Schroeder refinement is skipped by this implementation if
 `(digits2 <= (std::numeric_limits::digits - 3))`
 3 bits less than the maximum for the type `std::numeric_limits::digits`.
 
-For type `double`, the switch is at `policy `,
+For type `double`, the switch is at `policy`,
 so this is recommended for general use when speed is important.
 
-[tip Use `lambert_w0(z, policy >())` if speed is more important than the ultimate precision,
+[tip Use `lambert_w0(z, policy >())` or `lambert_w0(z, policy >())`  if speed is more important than the ultimate precision,
 for example, Monte Carlo applications.]
 
 If speed is [*very] important then `digits2<10>` only using Bisection alone might suffice:
@@ -293,7 +293,7 @@ the precision gain from the Schroeder and Halley is small.
 
 For floating-point types with precision greater than __fundamental_types,
 a `double` evaluation is used as a first approximation followed by Halley refinement,
-so that there is is little to gain trying to control precision.
+so that there is is little to gain in trying to control precision.
 Higher precision types are always going to be [*very, very much slower].
 
 For more examples of controlling precision, see
@@ -341,7 +341,7 @@ and anyway it only covers a tiny fraction of the range of possible z arguments v
 
 [h4:compilers Compilers]
 
-The code has been tested on VS 15.5.2, Clang 5.0.0 and GCC 7.1.0 but it is expected to work on all C++ compilers.
+The code has been tested on VS 15.5.5, Clang 5.0.0 and GCC 7.2.0 but it is expected to work on all C++ compilers.
 
 [h5:diagnostics Diagnostics Macros]
 
@@ -366,20 +366,21 @@ or defined in a jamfile.v2:  `BOOST_MATH_INSTRUMENT_LAMBERT_W_TERMS`
 
 ['First be right, then be fast.]
 
-There are many previous implementations with increasing accuracy and speed. See references below.
+There are many previous implementations, each with increasing accuracy and speed.
+See __lambert_w_references below.
 
 For most of the range of ['z] arguments, some initial approximation is followed by a single refinement,
-often using Halley or similar method, and gives a useful precision.
-To get a better precision, additional refinements, probably iterative, may be needed
-for example, using __halley or __schroder methods.
-For speed, several implementations avoid evaluation of a test
-using the exponential function calculating that a single refinement step will suffice,
+often using Halley or similar method, gives a useful precision.
+For speed, several implementations avoid evaluation of a test using the exponential function,
+estimating that a single refinement step will suffice,
 but these rarely get to the best result possible.
+To get a better precision, additional refinements, probably iterative, are needed
+for example, using __halley or __schroder methods.
 
-For the most precise results possible, for C++, closest to the nearest representation for the type,
+For C++, the most precise results possible, closest to the nearest __representable for the C++ type being used,
 it is usually necessary to use a higher precision type for intermediate computation,
 finally casting back to the smaller desired result type.
-This strategy is used by Maple and Wolfram, for example, using arbitrary precision arithmetic,
+This strategy is used by __Maple and WolframAlpha, for example, using arbitrary precision arithmetic,
 and some of their high-precision values are used for testing this library.
 This method is also used to provide __boost_test values using __multiprecision,
 typically, a 50 decimal digit type like `cpp_dec_float_50`
@@ -391,27 +392,26 @@ At extreme arguments near to zero or the singularity at the branch point,
 even this fails and the only method to achieve a really close result is to cast from a higher precision type.
 
 In practical applications, the increased computation required
-(often towards a thousand fold slower and requiring the additional code for multiprecision)
+(often towards a thousand-fold slower and requiring the additional code for multiprecision)
 is not justified and the algorithms here do not implement this.
 But because the Boost.Lambert_W algorithm has been tested using __multiprecision,
 users who require this can always easily achieve the nearest representation for the __fundamental_types
 if the application justifies the very large extra computation cost.
 
-One real-only implementation was based on an algorithm by __Luu_thesis,
+One compact real-only implementation was based on an algorithm by __Luu_thesis,
 (see routine 11 on page 98 for his Lambert W algorithm)
-and his Halley refinement is used iteratively in this implementation when required.
-
-This implementation is based on Thomas Luu's code posted at
+and his Halley refinement is used iteratively when required.
+A first implementation was based on Thomas Luu's code posted at
 [@https://svn.boost.org/trac/boost/ticket/11027 Boost Trac \#11027].
-
 It has been implemented from Luu's algorithm but templated on `RealType` parameter and result
 and handles both __fundamental_types (`float, double, long double`), __multiprecision,
-and also has been tested with a proposed fixed_point type.
+and also has been tested successfully with a proposed fixed_point type.
 
 A first approximation was computed using the method of Barry et al (see references 5 & 6 below).
-(For users only requiring an accuracy of relative accuracy of 0.02%, this function might suffice).
 This was extended to the widely used [@https://people.sc.fsu.edu/~jburkardt/f_src/toms443/toms443.html TOMS443]
 FORTRAN and C++ versions by John Burkardt using Schroeder refinement(s).
+(For users only requiring an accuracy of relative accuracy of 0.02%, Barry's function alone might suffice,
+but a better __rational_function approximation method has since been developed for this implementation).
 
 We also considered using __newton method.
 ``
@@ -422,11 +422,17 @@ We also considered using __newton method.
 ``
 but concluded that since the Newton-Raphson method takes typically 6 iterations to converge within tolerance,
 whereas Halley usually takes only 1 to 3 iterations to achieve an result within 1 __ulp,
+<<<<<<< HEAD
 so the Newton-Raphson method is unlikely to be quicker
 than the additional cost of computing the 2nd derivative for Halley's method,
 in this case requiring an expensive exponential function evaluation on each step.
 
 
+=======
+so Newton/Raphson's method unlikely to be quicker
+than the additional cost of computating the 2nd derivative for Halley's method,
+each iteration step requiring an expensive exponential function evaluation.
+>>>>>>> before trying to switch to develop to add test_value.hpp
 
 [h4:compact_implementation Implementing Compact Algorithms]
 
@@ -438,16 +444,16 @@ followed by Halley iteration (but is also slowest and least precise near zero an
 More recently, the Tosio Fukushima has developed an even faster algorithm,
 avoiding any transcendental function calls as these are necessarily expensive.
 The current implementation is based on his algorithm
-starting with a translation from FORTRAN into C++ by Darko Veberic.
+starting with a translation from Fukushima's FORTRAN into C++ by Darko Veberic.
 
 Many applications of the Lambert W function make many repeated evaluations for Monte Carlo methods;
 for these applications speed is very important.
-Luu and Chapeau-Blondeau and Monir provide typical usage examples.
+Luu, and Chapeau-Blondeau and Monir provide typical usage examples.
 
 Fukushima improves the important observation that much of the execution time of all
-previous iterative algorithms was spent evaluating transcendental functions, mainly `exp`.
+previous iterative algorithms was spent evaluating transcendental functions, usually `exp`.
 He has put a lot of work into avoiding any slow transcendental functions by using lookup tables and
-bisection, and finally by a single Schroeder refinement,
+bisection, finishing with a single Schroeder refinement,
 without any check on the precision of the result (necessarily evaluating an expensive exponential).
 
 Theoretical and practical tests confirm that this gives Lambert W estimates
@@ -457,22 +463,41 @@ However, though these give results within several __epsilon of the nearest repre
 they do not get as close as is very often possible with further refinement,
 usually to within one or two __epsilon.
 A mean difference was computed to express the typical error and is often about 0.5 epsilon,
-the theoretical minimum.
+the theoretical minimum.  Using the __float_distance, we can also express this as the number of
+bits that are different from the nearest representable or 'exact' value.
 
-A disadvantage is that the distribution of differences from the references value is reasonably random and unbiased,
-whereas the Schroeder stage is significantly biased.
+For this implementation of Lambert W-1, John Maddock used the Boost.Math `minimax.exe` __remez method program
+to devise a __rational_function for several ranges of argument z for the W0 branch of Lambert W0 function.
+These minimax rational approximations are combined for an algorithm that is both smaller and faster.
+Sadly it has not proved __remez_practical to use the same method for Lambert W-1 branch
+and so the Fukushima algorithm is retained for this branch.
 
+[h5:bias Bias]
 
-For speed during the bisection, Fukushima produces lookup tables of powers of e and z for integral Lambert W.
+An advantage of both minimax rational approximations and Halley iterated estimates
+is that the [*distribution] from the references value is reasonably random and unbiased,
+whereas any using the Schroeder refinement is significantly biased.
+Bias might have bad effects on statistical computations.
+For example, a test of 5000 values of argument z covering the main range of possible values,
+computed using the Fukushima method with Schroeder refinement gave about 1/3 lambert_w0 values that are
+one bit different from the best representation, and about 1% that are a few bits 'wrong'.
+If a Halley refinement step is added, only 1 in 30 are even one bit different, and none more than one bit 'wrong'.
+Starting with the rational approximation method, there is insignificant bias.
+Adding a Halley step (or iteration) reduces the number that are one bit different from about a third down to one in 30.
+This is unavoidable 'computational noise'.
+
+[h5:lookup_tables.
+
+For speed during the bisection, Fukushima computes lookup tables of powers of e and z for integral Lambert W.
 There are 64 elements in these tables. The FORTRAN version (and the C++ translation by Veberic)
 computed these (once) as `static` data. This is slower, may cause trouble with multithreading, and
 is slightly inaccurate because of rounding errors from repeated(64) multiplications. In this implementation
-the values have been computed using __multiprecision 50 decimal digit
-and output as C++ arrays 37 decimal digit long literals using
+the array values have been computed using __multiprecision 50 decimal digit
+and output as C++ arrays 37 decimal digit `long double` literals using `max_digits10` precision
 
   std::cout.precision(std::numeric_limits::max_digits10);
 
-The arrays are as const and static as possible (for the compiler version),
+The arrays are as const and constexpr and static as possible (for the compiler version),
 using BOOST_STATIC_CONSTEXPR macro.
 (See [@../../test/lambert_w_lookup_table_generator.cpp  lambert_w_lookup_table_generator.cpp]
 The precision was chosen to ensure that if used as `long double` arrays,
@@ -484,16 +509,16 @@ will be the nearest representable value for the type chose by a `typedef` in
   typedef double lookup_t; // Type for lookup table (`double` or `float`, or even `long double`?)
 
 This is to allow for future use at higher precision, up to platforms that use 128-bit
-(hardware or software) for `long double`.
+(hardware or software) for their `long double` type.
 
 The accuracy of the tables was confirmed using __WolframAlpha and agrees at the 37th decimal place,
-so ensuring that the value is read into even 128-bit `long double` to the nearest representation.
+so ensuring that the value is exactly read into even 128-bit `long double` to the nearest representation.
 
 [h5:higher_precision Higher precision]
 
 For types more precise than `double`, Fukushima showed that it is best to use the `double` estimate
-as a starting point followed by refinement using __halley iterations or other methods
-and our experience confirms this.
+as a starting point, followed by refinement using __halley iterations or other methods;
+our experience confirms this.
 
 Using __multiprecision it is simple to compute very high precision values of Lambert
 W at least to thousands of decimal digits over most of the range of z arguments.
@@ -632,7 +657,7 @@ only provides the precision of the built-in type, like `double`, only 17 decimal
 
 [tip Be exceeding careful not to silently lose precision by constructing multiprecision types from literal decimal types, usually [^double].  See examples.]
 
-Fukushima implmentation used 20 series terms and it was confirmed that using more terms does not usefully increase accuracy.
+Fukushima's implementation used 20 series terms; it was confirmed that using more terms does not usefully increase accuracy.
 
 [h5:wm1_near_zero Lambert W-1 arguments values very near zero.]
 
@@ -671,7 +696,7 @@ the computational cost of the log functions was easily paid by far fewer iterati
    double r = lambert_wm1(z, policy >());
    // lambert_wm1(-1.0000000000000000e-26) = -64.026509628385895
 
-[h5:halley]
+[h5:halley Halley refinement]
 
 After obtaining a double approximation, for `double` and quad 128-bit precision,
 a single iteration should suffice because
@@ -680,8 +705,6 @@ Halley iteration should triple the precision with each step
 and since we have at least half of the bits correct already,
 one Halley step is ample to get to 128-bit precision.
 
-
-
 [h4:testing Testing]
 
 Initial testing of the algorithm was done using a small number of spot tests.
diff --git a/example/lambert_w_simple_examples.cpp b/example/lambert_w_simple_examples.cpp
index e3e7d6cda..265aad14c 100644
--- a/example/lambert_w_simple_examples.cpp
+++ b/example/lambert_w_simple_examples.cpp
@@ -27,13 +27,6 @@
 #include  // For exp_minus_one == 3.67879441171442321595523770161460867e-01.
 #include 
 
-// Built-in/fundamental GCC float128 or Intel Quad 128-bit type, if available.
-#ifdef __GNUC__
-#include  // Not available for MSVC.
-// sets BOOST_MP_USE_FLOAT128 for GCC
-using boost::multiprecision::float128;
-#endif //# __GNUC__ and NOT _MSC_VER
-
 #include  // boost::multiprecision::cpp_dec_float_50
 using boost::multiprecision::cpp_dec_float_50; // 50 decimal digits type.
 using boost::multiprecision::cpp_dec_float_100; // 100 decimal digits type.
diff --git a/include/boost/math/special_functions/lambert_w.hpp b/include/boost/math/special_functions/lambert_w.hpp
index bdd1c029f..3f2c0e8e9 100644
--- a/include/boost/math/special_functions/lambert_w.hpp
+++ b/include/boost/math/special_functions/lambert_w.hpp
@@ -5,12 +5,8 @@
 // (See accompanying file LICENSE_1_0.txt or
 //  copy at http ://www.boost.org/LICENSE_1_0.txt).
 
-#ifndef BOOST_MATH_LAMBERTW0_HPP
-#define BOOST_MATH_LAMBERTW0_HPP
-
-// TODO change name of include file after PB/TF/DV version replaced fully.
-//#ifndef BOOST_MATH_SF_LAMBERT_W_HPP
-//#define BOOST_MATH_SF_LAMBERT_W_HPP
+#ifndef BOOST_MATH_SF_LAMBERT_W_HPP
+#define BOOST_MATH_SF_LAMBERT_W_HPP
 
 #ifdef _MSC_VER
 #  pragma warning (disable: 4127) // warning C4127: conditional expression is constant.
@@ -38,7 +34,7 @@ BOOST_MATH_INSTRUMENT_LAMBERT_W0 // W1 branch diagnostics.
 BOOST_MATH_INSTRUMENT_LAMBERT_W_HALLEY // Halley refinement diagnostics.
 BOOST_MATH_INSTRUMENT_LAMBERT_W_PRECISION // Precision.
 BOOST_MATH_INSTRUMENT_LAMBERT_W_HALLEY_W0 // Halley refinement diagnostics only for W0 branch.
-BOOST_MATH_INSTRUMENT_LAMBERT_W_HALLEY_WM1 // Halley refinement diagnostics only for W-1 branch.
+BOOST_MATH_INSTRUMENT_LAMBERT_WM1_HALLEY // Halley refinement diagnostics only for W-1 branch.
 BOOST_MATH_INSTRUMENT_LAMBERT_WM1_TINY // K > 64, z > -1.0264389699511303e-26
 BOOST_MATH_INSTRUMENT_LAMBERT_W0_HUGE // K > 64, z > 3.9904954117194348e+29
 BOOST_MATH_INSTRUMENT_LAMBERT_W_SCHROEDER // Schroeder refinement diagnostics.
@@ -1054,7 +1050,7 @@ inline T lambert_w0_imp(T z, const Policy& pol, const mpl::int_<1>&)
 template 
 inline T lambert_w0_imp(T z, const Policy& pol, const mpl::int_<2>&)
 {
-   static const char* function = "boost::math::lamber_w0<%1%>";
+   static const char* function = "boost::math::lambert_w0<%1%>";
    BOOST_MATH_STD_USING // Aid ADL of std functions.
 
    // Unusual case of 32-bit double with a wider/64-bit long double
@@ -1064,10 +1060,10 @@ inline T lambert_w0_imp(T z, const Policy& pol, const mpl::int_<2>&)
    "- or possibly edit the coefficients to have "
    "an appropriate size-suffix for 64-bit floats on your platform - L?");
 
-     if ((boost::math::isnan)(z))
-     {
-       return boost::math::policies::raise_domain_error(function, "Expected a value > -e^-1 (-0.367879...) but got %1%", z, pol);
-     }
+    if ((boost::math::isnan)(z))
+    {
+      return boost::math::policies::raise_domain_error(function, "Expected a value > -e^-1 (-0.367879...) but got %1%", z, pol);
+    }
 
    if (z >= 0.05)
    {
@@ -1207,7 +1203,7 @@ inline T lambert_w0_imp(T z, const Policy& pol, const mpl::int_<2>&)
             1.36363515125489502e-06,
             3.44200749053237945e-09,
          };
-         T log_w= log(z);
+         T log_w = log(z);
          return log_w+ Y + boost::math::tools::evaluate_rational(P, Q, log_w);
       }
       else if (z < 7.896296e+13)  // 9.2 < log(z) <= 32
@@ -1236,7 +1232,7 @@ inline T lambert_w0_imp(T z, const Policy& pol, const mpl::int_<2>&)
             4.97253225968548872e-09,
             3.39460723731970550e-12,
          };
-         T log_w= log(z);
+         T log_w = log(z);
          return log_w+ Y + boost::math::tools::evaluate_rational(P, Q, log_w);
       }
       else if (z < 2.6881171e+43) // 32 < log(z) < 100
@@ -1265,7 +1261,7 @@ inline T lambert_w0_imp(T z, const Policy& pol, const mpl::int_<2>&)
             -9.35271498075378319e-11,
             -2.60648331090076845e-14,
          };
-         T log_w= log(z);
+         T log_w = log(z);
          return log_w+ Y + boost::math::tools::evaluate_rational(P, Q, log_w);
       }
       else // 100 < log(z) < 710
@@ -1298,7 +1294,7 @@ inline T lambert_w0_imp(T z, const Policy& pol, const mpl::int_<2>&)
             1.11775518708172009e-20,
             3.78250395617836059e-25,
          };
-         T log_w= log(z);
+         T log_w = log(z);
          return log_w+ Y + boost::math::tools::evaluate_rational(P, Q, log_w);
       }
    }
@@ -1419,7 +1415,6 @@ inline T lambert_w0_imp(T z, const Policy& pol, const mpl::int_<2>&)
             -1.72845216404874299e+10,
          };
          T d = z + 0.36787944117144232159552377016146086744581113103176804;
-         // bad!
          return -d / (Y + boost::math::tools::evaluate_polynomial(P, d) / boost::math::tools::evaluate_polynomial(Q, d));
       }
       else if (z <= -0.36787944117144232159552377016146086744581113103176804) // Precision is max_digits10(cpp_bin_float_50).
@@ -1491,11 +1486,11 @@ inline T lambert_w0_imp(T z, const Policy& pol, const mpl::int_<0>&)
       return lambert_w_detail::lambert_w_singularity_series(p);
    }
 
-   // Whew!
+   // Phew!
    // If we get here we are in the normal range of the function, so
    // get a double precision approximation first, then iterate to full precision of T.
    // We define a type that is:
-   // mpl::true_  if there are so many digits required that iteration is required.
+   // mpl::true_  if there are so many digits precision wanted that iteration is necessary.
    // mpl::false_ if a single Halley step is sufficient.
    // For speed, we also cast z to type double when that is possible
    //   (boost::is_constructible() == true).
@@ -1507,12 +1502,12 @@ inline T lambert_w0_imp(T z, const Policy& pol, const mpl::int_<0>&)
       : false
    > tag_type;
 
-  // JM reused name z but Clang complains that it is hiding parameter z.
-  // T z = lambert_w0_imp(maybe_reduce_to_double(z, boost::is_constructible()), pol, mpl::int_<2>());
-  // return lambert_w_maybe_halley_iterate(z, z, tag_type());
    T w = lambert_w0_imp(maybe_reduce_to_double(z, boost::is_constructible()), pol, mpl::int_<2>());
 
    return lambert_w_maybe_halley_iterate(w, z, tag_type());
+
+ //  result = lambert_w_halley_iterate(result, z);
+
 } // T lambert_w0_imp(T z, const Policy& pol, const mpl::int_<0>&)  extended precision.
 
 
@@ -1843,20 +1838,20 @@ T lambert_wm1_imp(const T z, const Policy&  pol)
       { // Only enough digits2 required, so just return the Schroeder refined value.
         // For example, for float, if (d2 <= 22) then
         // digits-3 returns Schroeder result up to 3 bits might be wrong.
-#ifdef BOOST_MATH_INSTRUMENT_LAMBERT_WM1
+#ifdef BOOST_MATH_INSTRUMENT_LAMBERT_W_SCHROEDER
         std::cout << "Schroeder refinement estimate = " << result << std::endl;
-#endif // BOOST_MATH_INSTRUMENT_LAMBERT_WM1
+#endif // BOOST_MATH_INSTRUMENT_LAMBERT_W_SCHROEDER
         return result; // Schroeder only.
       }
       else // Perform additional Halley refinement(s) to ensure that
            // get a near as possible to correct result (usually +/- epsilon).
       {
         result = lambert_w_halley_iterate(result, z);
-#ifdef BOOST_MATH_INSTRUMENT_LAMBERT_W0_HALLEY
+#ifdef BOOST_MATH_INSTRUMENT_LAMBERT_WM1_HALLEY
         std::cout.precision(std::numeric_limits::max_digits10);
         std::cout << "Halley refinement estimate =    " << result << std::endl;
         std::cout.precision(saved_precision);
-#endif // BOOST_MATH_INSTRUMENT_LAMBERT_W0_HALLEY
+#endif // BOOST_MATH_INSTRUMENT_LAMBERT_W1_HALLEY
         return result; // Halley
       } // schroeder or Schroeder and Halley.
     } // more than 10 digits2
@@ -1867,76 +1862,106 @@ T lambert_wm1_imp(const T z, const Policy&  pol)
 /////////////////////  User Lambert w functions. //////////////////////////
 
 //! Lambert W0 using User-defined policy.
-template 
-inline
-typename tools::promote_args::type
-lambert_w0(T z, const Policy& pol)
-{
-   // Promote integer or expression template arguments to double,
-   // without doing any other internal promotion like float to double.
-   typedef typename tools::promote_args::type result_type;
-   //
-   // Work out what precision has been selected, based on the Policy and the
-   // number type,
-   typedef typename policies::precision::type precision_type;
-   // and then select the correct implementation based on that (not the type T):
-   typedef mpl::int_<
+  template 
+  inline
+    typename tools::promote_args::type
+    lambert_w0(T z, const Policy& pol)
+  {
+     // Promote integer or expression template arguments to double,
+     // without doing any other internal promotion like float to double.
+    typedef typename tools::promote_args::type result_type;
+    
+    // Work out what precision has been selected, 
+    // based on the Policy and the number type.
+    typedef typename policies::precision::type precision_type;
+    // and then select the correct implementation based on that (not the type T):
+    typedef mpl::int_<
       (precision_type::value == 0) || (precision_type::value > 53) ?
       0  // either variable precision (0), or greater than 64-bit precision.
       : (precision_type::value <= 24) ? 1 // 32-bit (probably float) precision.
       : 2  // 64-bit (probably double) precision.
-   > tag_type;
+      > tag_type;
 
-   return lambert_w_detail::lambert_w0_imp(result_type(z), pol, tag_type());
-} // lambert_w0(T z, const Policy& pol)
+    // Optionally, if policy is all T's digits (default), return after a extra Halley step.
+    // (Might be better with mpl?)
+    int p = precision_type::value;
+    int d = std::numeric_limits::digits;
 
-  //! Lambert W0 using default policy.
-template  >
-inline
-typename tools::promote_args::type
-lambert_w0(T z)
-{
-  // Promote integer or expression template arguments to double,
-  // without doing any other internal promotion like float to double.
-   typedef typename tools::promote_args::type result_type;
+    T result = lambert_w_detail::lambert_w0_imp(result_type(z), pol, tag_type()); // Pre-Halley.
 
-   // Work out what precision has been selected, based on the Policy and the number type.
-   // For the default policy version, we want the *default policy* precision for T.
-   typedef typename policies::precision::type precision_type;
-   // and then select the correct implementation based on that (not the type T):
-   typedef mpl::int_<
-     (precision_type::value == 0) || (precision_type::value > 53) ?
-     0  // either variable precision (0), or greater than 64-bit precision.
-     : (precision_type::value <= 24) ? 1 // 32-bit (probably float) precision.
-     : 2  // 64-bit (probably double) precision.
-   > tag_type;
+    bool quick = (precision_type::value < std::numeric_limits::digits);
 
-   return lambert_w_detail::lambert_w0_imp(result_type(z), policies::policy<>(), tag_type()); //
-} // lambert_w0(T z)
+    if (quick || (result == -1) || (result == 0) || (!boost::math::isnormal(result)) || (z < -0.3578))
+    {
+      return result;
+    }
+    else
+    {
+      result = lambert_w_detail::lambert_w_halley_step(result, z);
+    }
+    return result;
+    // return lambert_w_detail::lambert_w0_imp(result_type(z), pol, tag_type());
+  } // lambert_w0(T z, const Policy& pol)
 
-  //! W-1 branch (-max(z) < z <= -1/e).
+    //! Lambert W0 using default policy.
+  template  >
+  inline
+    typename tools::promote_args::type
+    lambert_w0(T z)
+  {
+    // Promote integer or expression template arguments to double,
+    // without doing any other internal promotion like float to double.
+    typedef typename tools::promote_args::type result_type;
 
-  //! Lambert W-1 using User-defined policy.
-template 
-inline
-typename tools::promote_args::type
-lambert_wm1(T z, const Policy& pol)
-{
-  // Promote integer or expression template arguments to double,
-  // without doing any other internal promotion like float to double.
-  typedef typename tools::promote_args::type result_type;
-  return lambert_w_detail::lambert_wm1_imp(result_type(z), pol); //
-}
+    // Work out what precision has been selected, based on the Policy and the number type.
+    // For the default policy version, we want the *default policy* precision for T.
+    typedef typename policies::precision::type precision_type;
+    // and then select the correct implementation based on that (not the type T):
+    typedef mpl::int_<
+      (precision_type::value == 0) || (precision_type::value > 53) ?
+      0  // either variable precision (0), or greater than 64-bit precision.
+      : (precision_type::value <= 24) ? 1 // 32-bit (probably float) precision.
+      : 2  // 64-bit (probably double) precision.
+    > tag_type;
+    // Always add Halley step for default policy.
+    T result = lambert_w_detail::lambert_w0_imp(result_type(z),  policies::policy<>(), tag_type()); // Pre-Halley.
+    if ((z < -0.3578) // Using singularity series.
+      || (!boost::math::isnormal(result))  // NaN, infinity or denormal.
+      || (result == -1) || (result == 0)) // exact
+    {
+      return result;
+    }
+    result = lambert_w_detail::lambert_w_halley_step(result, z);
+    return result;
+    // TODO this ugliness is because JM rationals return rather than result = ...
+    // and should use existing checks for -1, 0, singularity series and small_z series inside _imp
+    // to avoid repeating here.
+//    return lambert_w_detail::lambert_w0_imp(result_type(z), policies::policy<>(), tag_type()); //
+  } // lambert_w0(T z)
 
-//! Lambert W-1 using default policy.
-template 
-inline
-typename tools::promote_args::type
-lambert_wm1(T z)
-{
-  typedef typename tools::promote_args::type result_type;
-  return lambert_w_detail::lambert_wm1_imp(result_type(z), policies::policy<>());
-} // lambert_wm1(T z)
+    //! W-1 branch (-max(z) < z <= -1/e).
+
+    //! Lambert W-1 using User-defined policy.
+  template 
+  inline
+    typename tools::promote_args::type
+    lambert_wm1(T z, const Policy& pol)
+  {
+    // Promote integer or expression template arguments to double,
+    // without doing any other internal promotion like float to double.
+    typedef typename tools::promote_args::type result_type;
+    return lambert_w_detail::lambert_wm1_imp(result_type(z), pol); //
+  }
+
+  //! Lambert W-1 using default policy.
+  template 
+  inline
+    typename tools::promote_args::type
+    lambert_wm1(T z)
+  {
+    typedef typename tools::promote_args::type result_type;
+    return lambert_w_detail::lambert_wm1_imp(result_type(z), policies::policy<>());
+  } // lambert_wm1(T z)
 
 
 template 
@@ -1982,4 +2007,4 @@ lambert_wm1_prime(T z)
 
 }} //boost::math namespaces
 
-#endif // #ifdef BOOST_MATH_LAMBERTW_HPP
+#endif // #ifdef BOOST_MATH_SF_LAMBERT_W_HPP
diff --git a/include/boost/math/special_functions/lambert_w_lookup_table.ipp b/include/boost/math/special_functions/lambert_w_lookup_table.ipp
index 718c73498..4c48cdff5 100644
--- a/include/boost/math/special_functions/lambert_w_lookup_table.ipp
+++ b/include/boost/math/special_functions/lambert_w_lookup_table.ipp
@@ -14,13 +14,13 @@
 
 // Written by I:\modular-boost\libs\math\test\lambert_w_lookup_table_generator.cpp Thu Jan 25 16:52:07 2018
 
-// Sizes of arrays of z values for Lambert W[0], W[1] ... W[64]"nand W[-1], W[-2] ... W[-64].
+// Sizes of arrays of z values for Lambert W[0], W[1] ... W[64]" and W[-1], W[-2] ... W[-64].
 
 namespace boost {
 namespace math {
 namespace lambert_w_detail {
 namespace lambert_w_lookup
-{ 
+{
 BOOST_STATIC_CONSTEXPR std::size_t  noof_sqrts = 12;
 BOOST_STATIC_CONSTEXPR std::size_t  noof_halves = 12;
 BOOST_STATIC_CONSTEXPR std::size_t  noof_w0es = 65;
@@ -28,104 +28,104 @@ BOOST_STATIC_CONSTEXPR std::size_t  noof_w0zs = 65;
 BOOST_STATIC_CONSTEXPR std::size_t  noof_wm1es = 64;
 BOOST_STATIC_CONSTEXPR std::size_t  noof_wm1zs = 64;
 
-BOOST_STATIC_CONSTEXPR lookup_t halves[noof_halves] = 
+BOOST_STATIC_CONSTEXPR lookup_t halves[noof_halves] =
 { // Common to Lambert W0 and W-1 (and exactly representable).
   0.5L, 0.25L, 0.125L, 0.0625L, 0.03125L, 0.015625L, 0.0078125L, 0.00390625L, 0.001953125L, 0.0009765625L, 0.00048828125L, 0.000244140625L
 }; // halves, 0.5, 0.25, ... 0.000244140625, common to W0 and W-1.
 
-BOOST_STATIC_CONSTEXPR lookup_t sqrtw0s[noof_sqrts] = 
+BOOST_STATIC_CONSTEXPR lookup_t sqrtw0s[noof_sqrts] =
 {  // For Lambert W0 only.
   0.6065306597126334242631173765403218567L, 0.77880078307140486866846070009071995L, 0.882496902584595403104717592968701829L, 0.9394130628134757862473572557999761753L, 0.9692332344763440819139583751755278177L, 0.9844964370054084060204319075254540376L, 0.9922179382602435121227899136829802692L, 0.996101369470117490071323985506950379L, 0.9980487811074754727142805899944244847L, 0.9990239141819756622368328253791383317L, 0.9995118379398893653889967919448497792L, 0.9997558891748972165136242351259789505L
 }; // sqrtw0s
 
-BOOST_STATIC_CONSTEXPR lookup_t sqrtwm1s[noof_sqrts] = 
+BOOST_STATIC_CONSTEXPR lookup_t sqrtwm1s[noof_sqrts] =
 { // For Lambert W-1 only.
   1.648721270700128146848650787814163572L, 1.284025416687741484073420568062436458L, 1.133148453066826316829007227811793873L, 1.064494458917859429563390594642889673L, 1.031743407499102670938747815281507144L, 1.015747708586685747458535072082351749L, 1.007843097206447977693453559760123579L, 1.003913889338347573443609603903460282L, 1.001955033591002812046518898047477216L, 1.000977039492416535242845292611606506L, 1.000488400478694473126173623807163354L, 1.000244170429747854937005233924135774L
 }; // sqrtwm1s
 
-BOOST_STATIC_CONSTEXPR lookup_t w0es[noof_w0zs] = 
+BOOST_STATIC_CONSTEXPR lookup_t w0es[noof_w0zs] =
 { // Fukushima e powers array e[0] = 2.718, 1., e[2] = e^-1 = 0.135, e[3] = e^-2 = 0.133 ... e[64] = 4.3596100000630809736231248158884615452e-28.
-  2.7182818284590452353602874713526624978e+00L, 
-  1.0000000000000000000000000000000000000e+00L, 3.6787944117144232159552377016146086745e-01L, 1.3533528323661269189399949497248440341e-01L, 4.9787068367863942979342415650061776632e-02L, 
-  1.8315638888734180293718021273241242212e-02L, 6.7379469990854670966360484231484242488e-03L, 2.4787521766663584230451674308166678915e-03L, 9.1188196555451620800313608440928262647e-04L, 
-  3.3546262790251183882138912578086101931e-04L, 1.2340980408667954949763669073003382607e-04L, 4.5399929762484851535591515560550610238e-05L, 1.6701700790245659312635517360580879078e-05L, 
-  6.1442123533282097586823081788055323112e-06L, 2.2603294069810543257852772905386894694e-06L, 8.3152871910356788406398514256526229461e-07L, 3.0590232050182578837147949770228963937e-07L, 
-  1.1253517471925911451377517906012719164e-07L, 4.1399377187851666596510277189552806229e-08L, 1.5229979744712628436136629233517431862e-08L, 5.6027964375372675400129828162064630798e-09L, 
-  2.0611536224385578279659403801558209764e-09L, 7.5825604279119067279417432412681264430e-10L, 2.7894680928689248077189130306442932077e-10L, 1.0261879631701890303927527840612497760e-10L, 
-  3.7751345442790977516449695475234067792e-11L, 1.3887943864964020594661763746086856910e-11L, 5.1090890280633247198744001934792157666e-12L, 1.8795288165390832947582704184221926212e-12L, 
-  6.9144001069402030094125846587414092712e-13L, 2.5436656473769229103033856148576816666e-13L, 9.3576229688401746049158322233787067450e-14L, 3.4424771084699764583923893328515572846e-14L, 
-  1.2664165549094175723120904155965096382e-14L, 4.6588861451033973641842455436101684114e-15L, 1.7139084315420129663027203425760492412e-15L, 6.3051167601469893856390211922465427614e-16L, 
-  2.3195228302435693883122636097380800411e-16L, 8.5330476257440657942780498229412441658e-17L, 3.1391327920480296287089646522319196491e-17L, 1.1548224173015785986262442063323868655e-17L, 
-  4.2483542552915889953292347828586580179e-18L, 1.5628821893349887680908829951058341550e-18L, 5.7495222642935598066643808805734234249e-19L, 2.1151310375910804866314010070226514702e-19L, 
-  7.7811322411337965157133167292798981918e-20L, 2.8625185805493936444701216291839372068e-20L, 1.0530617357553812378763324449428108806e-20L, 3.8739976286871871129314774972691278293e-21L, 
-  1.4251640827409351062853210280340602263e-21L, 5.2428856633634639371718053028323436716e-22L, 1.9287498479639177830173428165270125748e-22L, 7.0954741622847041389832693878080734877e-23L, 
-  2.6102790696677048047026953153318648093e-23L, 9.6026800545086760302307696700074909076e-24L, 3.5326285722008070297353928101772088374e-24L, 1.2995814250075030736007134060714855303e-24L, 
-  4.7808928838854690812771770423179628939e-25L, 1.7587922024243116489558751288034363178e-25L, 6.4702349256454603261540395529264893765e-26L, 2.3802664086944006058943245888024963309e-26L, 
-  8.7565107626965203384887328007391660366e-27L, 3.2213402859925160890012477758489437534e-27L, 1.1850648642339810062850307390972809891e-27L, 4.3596100000630809736231248158884596428e-28L, 
-  
+  2.7182818284590452353602874713526624978e+00L,
+  1.0000000000000000000000000000000000000e+00L, 3.6787944117144232159552377016146086745e-01L, 1.3533528323661269189399949497248440341e-01L, 4.9787068367863942979342415650061776632e-02L,
+  1.8315638888734180293718021273241242212e-02L, 6.7379469990854670966360484231484242488e-03L, 2.4787521766663584230451674308166678915e-03L, 9.1188196555451620800313608440928262647e-04L,
+  3.3546262790251183882138912578086101931e-04L, 1.2340980408667954949763669073003382607e-04L, 4.5399929762484851535591515560550610238e-05L, 1.6701700790245659312635517360580879078e-05L,
+  6.1442123533282097586823081788055323112e-06L, 2.2603294069810543257852772905386894694e-06L, 8.3152871910356788406398514256526229461e-07L, 3.0590232050182578837147949770228963937e-07L,
+  1.1253517471925911451377517906012719164e-07L, 4.1399377187851666596510277189552806229e-08L, 1.5229979744712628436136629233517431862e-08L, 5.6027964375372675400129828162064630798e-09L,
+  2.0611536224385578279659403801558209764e-09L, 7.5825604279119067279417432412681264430e-10L, 2.7894680928689248077189130306442932077e-10L, 1.0261879631701890303927527840612497760e-10L,
+  3.7751345442790977516449695475234067792e-11L, 1.3887943864964020594661763746086856910e-11L, 5.1090890280633247198744001934792157666e-12L, 1.8795288165390832947582704184221926212e-12L,
+  6.9144001069402030094125846587414092712e-13L, 2.5436656473769229103033856148576816666e-13L, 9.3576229688401746049158322233787067450e-14L, 3.4424771084699764583923893328515572846e-14L,
+  1.2664165549094175723120904155965096382e-14L, 4.6588861451033973641842455436101684114e-15L, 1.7139084315420129663027203425760492412e-15L, 6.3051167601469893856390211922465427614e-16L,
+  2.3195228302435693883122636097380800411e-16L, 8.5330476257440657942780498229412441658e-17L, 3.1391327920480296287089646522319196491e-17L, 1.1548224173015785986262442063323868655e-17L,
+  4.2483542552915889953292347828586580179e-18L, 1.5628821893349887680908829951058341550e-18L, 5.7495222642935598066643808805734234249e-19L, 2.1151310375910804866314010070226514702e-19L,
+  7.7811322411337965157133167292798981918e-20L, 2.8625185805493936444701216291839372068e-20L, 1.0530617357553812378763324449428108806e-20L, 3.8739976286871871129314774972691278293e-21L,
+  1.4251640827409351062853210280340602263e-21L, 5.2428856633634639371718053028323436716e-22L, 1.9287498479639177830173428165270125748e-22L, 7.0954741622847041389832693878080734877e-23L,
+  2.6102790696677048047026953153318648093e-23L, 9.6026800545086760302307696700074909076e-24L, 3.5326285722008070297353928101772088374e-24L, 1.2995814250075030736007134060714855303e-24L,
+  4.7808928838854690812771770423179628939e-25L, 1.7587922024243116489558751288034363178e-25L, 6.4702349256454603261540395529264893765e-26L, 2.3802664086944006058943245888024963309e-26L,
+  8.7565107626965203384887328007391660366e-27L, 3.2213402859925160890012477758489437534e-27L, 1.1850648642339810062850307390972809891e-27L, 4.3596100000630809736231248158884596428e-28L,
+
 }; // w0es
 
-BOOST_STATIC_CONSTEXPR lookup_t w0zs[noof_w0zs] = 
-{ // z values for W[0], W[1], W[2] ... W[64] (Fukushima array Fk). 
-  0.0000000000000000000000000000000000000e+00L, 
-  2.7182818284590452353602874713526624978e+00L, 1.4778112197861300454460854921150015626e+01L, 6.0256610769563003222785588963745153691e+01L, 2.1839260013257695631244104481144351361e+02L, 
-  7.4206579551288301710557790020276139812e+02L, 2.4205727609564107356503230832603296776e+03L, 7.6764321089992101948460416680168500271e+03L, 2.3847663896333826197948736795623109390e+04L, 
-  7.2927755348178456069389970204894839685e+04L, 2.2026465794806716516957900645284244366e+05L, 6.5861555886717600300859134371483559776e+05L, 1.9530574970280470496960624587818413980e+06L, 
-  5.7513740961159665432393360873381476632e+06L, 1.6836459978306874888489314790750032292e+07L, 4.9035260587081659589527825691375819733e+07L, 1.4217776832812596218820837985250320561e+08L, 
-  4.1063419681078006965118239806655900596e+08L, 1.1818794444719492004981570586630806042e+09L, 3.3911637183005579560532906419857313738e+09L, 9.7033039081958055593821366108308111737e+09L, 
-  2.7695130424147508641409976558651358487e+10L, 7.8868082614895014356985518811525255163e+10L, 2.2413047926372475980079655175092843139e+11L, 6.3573893111624333505933989166748517618e+11L, 
-  1.8001224834346468131040337866531539479e+12L, 5.0889698451498078710141863447784789126e+12L, 1.4365302496248562650461177217211790925e+13L, 4.0495197800161304862957327843914007993e+13L, 
-  1.1400869461717722015726999684446230289e+14L, 3.2059423744573386440971405952224204950e+14L, 9.0051433962267018216365614546207459567e+14L, 2.5268147258457822451512967243234631750e+15L, 
-  7.0832381329352301326018261305316090522e+15L, 1.9837699245933465967698692976753294646e+16L, 5.5510470830970075484537561902113104381e+16L, 1.5520433569614702817608320254284931407e+17L, 
-  4.3360826779369661842459877227403719730e+17L, 1.2105254067703227363724895246493485480e+18L, 3.3771426165357561311906703760513324357e+18L, 9.4154106734807994163159964299613921804e+18L, 
-  2.6233583234732252918129199544138403574e+19L, 7.3049547543861043990576614751671879498e+19L, 2.0329709713386190214340167519800405595e+20L, 5.6547040503180956413560918381429636734e+20L, 
-  1.5720421975868292906615658755032744790e+21L, 4.3682149334771264822761478593874428627e+21L, 1.2132170565093316762294432610117848880e+22L, 3.3680332378068632345542636794533635462e+22L, 
-  9.3459982052259884835729892206738573922e+22L, 2.5923527642935362320437266614667426924e+23L, 7.1876803203773878618909930893087860822e+23L, 1.9921241603726199616378561653688236827e+24L, 
-  5.5192924995054165325072406547517121131e+24L, 1.5286067837683347062387143159276002521e+25L, 4.2321318958281094260005100745711666956e+25L, 1.1713293177672778461879598480402173158e+26L, 
-  3.2408603996214813669049988277609543829e+26L, 8.9641258264226027960478448084812796397e+26L, 2.4787141382364034104243901241243054434e+27L, 6.8520443388941057019777430988685937812e+27L, 
+BOOST_STATIC_CONSTEXPR lookup_t w0zs[noof_w0zs] =
+{ // z values for W[0], W[1], W[2] ... W[64] (Fukushima array Fk).
+  0.0000000000000000000000000000000000000e+00L,
+  2.7182818284590452353602874713526624978e+00L, 1.4778112197861300454460854921150015626e+01L, 6.0256610769563003222785588963745153691e+01L, 2.1839260013257695631244104481144351361e+02L,
+  7.4206579551288301710557790020276139812e+02L, 2.4205727609564107356503230832603296776e+03L, 7.6764321089992101948460416680168500271e+03L, 2.3847663896333826197948736795623109390e+04L,
+  7.2927755348178456069389970204894839685e+04L, 2.2026465794806716516957900645284244366e+05L, 6.5861555886717600300859134371483559776e+05L, 1.9530574970280470496960624587818413980e+06L,
+  5.7513740961159665432393360873381476632e+06L, 1.6836459978306874888489314790750032292e+07L, 4.9035260587081659589527825691375819733e+07L, 1.4217776832812596218820837985250320561e+08L,
+  4.1063419681078006965118239806655900596e+08L, 1.1818794444719492004981570586630806042e+09L, 3.3911637183005579560532906419857313738e+09L, 9.7033039081958055593821366108308111737e+09L,
+  2.7695130424147508641409976558651358487e+10L, 7.8868082614895014356985518811525255163e+10L, 2.2413047926372475980079655175092843139e+11L, 6.3573893111624333505933989166748517618e+11L,
+  1.8001224834346468131040337866531539479e+12L, 5.0889698451498078710141863447784789126e+12L, 1.4365302496248562650461177217211790925e+13L, 4.0495197800161304862957327843914007993e+13L,
+  1.1400869461717722015726999684446230289e+14L, 3.2059423744573386440971405952224204950e+14L, 9.0051433962267018216365614546207459567e+14L, 2.5268147258457822451512967243234631750e+15L,
+  7.0832381329352301326018261305316090522e+15L, 1.9837699245933465967698692976753294646e+16L, 5.5510470830970075484537561902113104381e+16L, 1.5520433569614702817608320254284931407e+17L,
+  4.3360826779369661842459877227403719730e+17L, 1.2105254067703227363724895246493485480e+18L, 3.3771426165357561311906703760513324357e+18L, 9.4154106734807994163159964299613921804e+18L,
+  2.6233583234732252918129199544138403574e+19L, 7.3049547543861043990576614751671879498e+19L, 2.0329709713386190214340167519800405595e+20L, 5.6547040503180956413560918381429636734e+20L,
+  1.5720421975868292906615658755032744790e+21L, 4.3682149334771264822761478593874428627e+21L, 1.2132170565093316762294432610117848880e+22L, 3.3680332378068632345542636794533635462e+22L,
+  9.3459982052259884835729892206738573922e+22L, 2.5923527642935362320437266614667426924e+23L, 7.1876803203773878618909930893087860822e+23L, 1.9921241603726199616378561653688236827e+24L,
+  5.5192924995054165325072406547517121131e+24L, 1.5286067837683347062387143159276002521e+25L, 4.2321318958281094260005100745711666956e+25L, 1.1713293177672778461879598480402173158e+26L,
+  3.2408603996214813669049988277609543829e+26L, 8.9641258264226027960478448084812796397e+26L, 2.4787141382364034104243901241243054434e+27L, 6.8520443388941057019777430988685937812e+27L,
   1.8936217407781711443114787060753312270e+28L, 5.2317811346197017832254642778313331353e+28L, 1.4450833904658542238325922893692265683e+29L, 3.9904954117194348050619127737142206367e+29L
-  
+
 }; // w0zs
 
-BOOST_STATIC_CONSTEXPR lookup_t wm1es[noof_wm1es] = 
+BOOST_STATIC_CONSTEXPR lookup_t wm1es[noof_wm1es] =
 { // Fukushima e array e[0] = e^1 = 2.718, e[1] = e^2 = 7.39 ... e[64] = 4.60718e+28.
-  2.7182818284590452353602874713526624978e+00L, 
-  7.3890560989306502272304274605750078132e+00L, 2.0085536923187667740928529654581717897e+01L, 5.4598150033144239078110261202860878403e+01L, 1.4841315910257660342111558004055227962e+02L, 
-  4.0342879349273512260838718054338827961e+02L, 1.0966331584284585992637202382881214324e+03L, 2.9809579870417282747435920994528886738e+03L, 8.1030839275753840077099966894327599650e+03L, 
-  2.2026465794806716516957900645284244366e+04L, 5.9874141715197818455326485792257781614e+04L, 1.6275479141900392080800520489848678317e+05L, 4.4241339200892050332610277594908828178e+05L, 
-  1.2026042841647767777492367707678594494e+06L, 3.2690173724721106393018550460917213155e+06L, 8.8861105205078726367630237407814503508e+06L, 2.4154952753575298214775435180385823880e+07L, 
-  6.5659969137330511138786503259060033569e+07L, 1.7848230096318726084491003378872270388e+08L, 4.8516519540979027796910683054154055868e+08L, 1.3188157344832146972099988837453027851e+09L, 
-  3.5849128461315915616811599459784206892e+09L, 9.7448034462489026000346326848229752776e+09L, 2.6489122129843472294139162152811882341e+10L, 7.2004899337385872524161351466126157915e+10L, 
-  1.9572960942883876426977639787609534279e+11L, 5.3204824060179861668374730434117744166e+11L, 1.4462570642914751736770474229969288569e+12L, 3.9313342971440420743886205808435276858e+12L, 
-  1.0686474581524462146990468650741401650e+13L, 2.9048849665247425231085682111679825667e+13L, 7.8962960182680695160978022635108224220e+13L, 2.1464357978591606462429776153126088037e+14L, 
-  5.8346174252745488140290273461039101900e+14L, 1.5860134523134307281296446257746601252e+15L, 4.3112315471151952271134222928569253908e+15L, 1.1719142372802611308772939791190194522e+16L, 
-  3.1855931757113756220328671701298646000e+16L, 8.6593400423993746953606932719264934250e+16L, 2.3538526683701998540789991074903480451e+17L, 6.3984349353005494922266340351557081888e+17L, 
-  1.7392749415205010473946813036112352261e+18L, 4.7278394682293465614744575627442803708e+18L, 1.2851600114359308275809299632143099258e+19L, 3.4934271057485095348034797233406099533e+19L, 
-  9.4961194206024488745133649117118323102e+19L, 2.5813128861900673962328580021527338043e+20L, 7.0167359120976317386547159988611740546e+20L, 1.9073465724950996905250998409538484474e+21L, 
-  5.1847055285870724640874533229334853848e+21L, 1.4093490824269387964492143312370168789e+22L, 3.8310080007165768493035695487861993899e+22L, 1.0413759433029087797183472933493796440e+23L, 
-  2.8307533032746939004420635480140745409e+23L, 7.6947852651420171381827455901293939921e+23L, 2.0916594960129961539070711572146737782e+24L, 5.6857199993359322226403488206332533034e+24L, 
-  1.5455389355901039303530766911174620068e+25L, 4.2012104037905142549565934307191617684e+25L, 1.1420073898156842836629571831447656302e+26L, 3.1042979357019199087073421411071003721e+26L, 
+  2.7182818284590452353602874713526624978e+00L,
+  7.3890560989306502272304274605750078132e+00L, 2.0085536923187667740928529654581717897e+01L, 5.4598150033144239078110261202860878403e+01L, 1.4841315910257660342111558004055227962e+02L,
+  4.0342879349273512260838718054338827961e+02L, 1.0966331584284585992637202382881214324e+03L, 2.9809579870417282747435920994528886738e+03L, 8.1030839275753840077099966894327599650e+03L,
+  2.2026465794806716516957900645284244366e+04L, 5.9874141715197818455326485792257781614e+04L, 1.6275479141900392080800520489848678317e+05L, 4.4241339200892050332610277594908828178e+05L,
+  1.2026042841647767777492367707678594494e+06L, 3.2690173724721106393018550460917213155e+06L, 8.8861105205078726367630237407814503508e+06L, 2.4154952753575298214775435180385823880e+07L,
+  6.5659969137330511138786503259060033569e+07L, 1.7848230096318726084491003378872270388e+08L, 4.8516519540979027796910683054154055868e+08L, 1.3188157344832146972099988837453027851e+09L,
+  3.5849128461315915616811599459784206892e+09L, 9.7448034462489026000346326848229752776e+09L, 2.6489122129843472294139162152811882341e+10L, 7.2004899337385872524161351466126157915e+10L,
+  1.9572960942883876426977639787609534279e+11L, 5.3204824060179861668374730434117744166e+11L, 1.4462570642914751736770474229969288569e+12L, 3.9313342971440420743886205808435276858e+12L,
+  1.0686474581524462146990468650741401650e+13L, 2.9048849665247425231085682111679825667e+13L, 7.8962960182680695160978022635108224220e+13L, 2.1464357978591606462429776153126088037e+14L,
+  5.8346174252745488140290273461039101900e+14L, 1.5860134523134307281296446257746601252e+15L, 4.3112315471151952271134222928569253908e+15L, 1.1719142372802611308772939791190194522e+16L,
+  3.1855931757113756220328671701298646000e+16L, 8.6593400423993746953606932719264934250e+16L, 2.3538526683701998540789991074903480451e+17L, 6.3984349353005494922266340351557081888e+17L,
+  1.7392749415205010473946813036112352261e+18L, 4.7278394682293465614744575627442803708e+18L, 1.2851600114359308275809299632143099258e+19L, 3.4934271057485095348034797233406099533e+19L,
+  9.4961194206024488745133649117118323102e+19L, 2.5813128861900673962328580021527338043e+20L, 7.0167359120976317386547159988611740546e+20L, 1.9073465724950996905250998409538484474e+21L,
+  5.1847055285870724640874533229334853848e+21L, 1.4093490824269387964492143312370168789e+22L, 3.8310080007165768493035695487861993899e+22L, 1.0413759433029087797183472933493796440e+23L,
+  2.8307533032746939004420635480140745409e+23L, 7.6947852651420171381827455901293939921e+23L, 2.0916594960129961539070711572146737782e+24L, 5.6857199993359322226403488206332533034e+24L,
+  1.5455389355901039303530766911174620068e+25L, 4.2012104037905142549565934307191617684e+25L, 1.1420073898156842836629571831447656302e+26L, 3.1042979357019199087073421411071003721e+26L,
   8.4383566687414544890733294803731179601e+26L, 2.2937831594696098790993528402686136005e+27L, 6.2351490808116168829092387089284697448e+27L
 }; // wm1es
 
-BOOST_STATIC_CONSTEXPR lookup_t wm1zs[noof_wm1zs] = 
+BOOST_STATIC_CONSTEXPR lookup_t wm1zs[noof_wm1zs] =
 { // Fukushima G array of z values for integral K, (Fukushima Gk) g[0] (k = -1) = 1 ... g[64] = -1.0264389699511303e-26.
-  -3.6787944117144232159552377016146086745e-01L, 
-  -2.7067056647322538378799898994496880682e-01L, -1.4936120510359182893802724695018532990e-01L, -7.3262555554936721174872085092964968848e-02L, -3.3689734995427335483180242115742121244e-02L, 
-  -1.4872513059998150538271004584900007349e-02L, -6.3831737588816134560219525908649783853e-03L, -2.6837010232200947105711130062468881545e-03L, -1.1106882367801159454787302165703044346e-03L, 
-  -4.5399929762484851535591515560550610238e-04L, -1.8371870869270225243899069096638966986e-04L, -7.3730548239938517104187698145666387735e-05L, -2.9384282290753706235208604777002963102e-05L, 
-  -1.1641402067449950376895791995913672125e-05L, -4.5885348075273868255721924655343445906e-06L, -1.8005627955081458322204028649620350662e-06L, -7.0378941219347833214067471222239770590e-07L, 
-  -2.7413963540482731185045932620331377352e-07L, -1.0645313231320808326024667350792279852e-07L, -4.1223072448771156559318807603116419528e-08L, -1.5923376898615004128677660806663065530e-08L, 
-  -6.1368298043116345769816086674174450569e-09L, -2.3602323152914347699033314033408744848e-09L, -9.0603229062698346039479269140561762700e-10L, -3.4719859662410051486654409365217142276e-10L, 
-  -1.3283631472964644271673440503045960993e-10L, -5.0747278046555248958473301297399200773e-11L, -1.9360320299432568426355237044475945959e-11L, -7.3766303773930764398798182830872768331e-12L, 
-  -2.8072868906520523814747496670136120235e-12L, -1.0671679036256927021016406931839827582e-12L, -4.0525329757101362313986893299088308423e-13L, -1.5374324278841211301808010293913555758e-13L, 
-  -5.8272886672428440854292491647585674200e-14L, -2.2067908660514462849736574172862899665e-14L, -8.3502821888768497979241489950570881481e-15L, -3.1572276215253043438828784344882603413e-15L, 
-  -1.1928704609782512589094065678481294667e-15L, -4.5038074274761565346423524046963087756e-16L, -1.6993417021166355981316939131434632072e-16L, -6.4078169762734539491726202799339200356e-17L, 
-  -2.4147993510032951187990399698408378385e-17L, -9.0950634616416460925150243301974013218e-18L, -3.4236981860988704669138593608831552044e-18L, -1.2881333612472271400115547331327717431e-18L, 
-  -4.8440839844747536942311292467369300505e-19L, -1.8207788854829779430777944237164900798e-19L, -6.8407875971564885101695409345634890862e-20L, -2.5690139750480973292141845983878483991e-20L, 
-  -9.6437492398195889150867140826350628738e-21L, -3.6186918227651991108814673877821174787e-21L, -1.3573451162272064984454015639725697008e-21L, -5.0894204288895982960223079251039701810e-22L, 
-  -1.9076194289884357960571121174956927722e-22L, -7.1476978375412669048039237333931704166e-23L, -2.6773000149758626855152191436980592206e-23L, -1.0025115553818576399048488234179587012e-23L, 
-  -3.7527362568743669891693429406973638384e-24L, -1.4043571811296963574776515073934728352e-24L, -5.2539064576179122030932396804434996219e-25L, -1.9650175744554348142907611432678556896e-25L, 
+  -3.6787944117144232159552377016146086745e-01L,
+  -2.7067056647322538378799898994496880682e-01L, -1.4936120510359182893802724695018532990e-01L, -7.3262555554936721174872085092964968848e-02L, -3.3689734995427335483180242115742121244e-02L,
+  -1.4872513059998150538271004584900007349e-02L, -6.3831737588816134560219525908649783853e-03L, -2.6837010232200947105711130062468881545e-03L, -1.1106882367801159454787302165703044346e-03L,
+  -4.5399929762484851535591515560550610238e-04L, -1.8371870869270225243899069096638966986e-04L, -7.3730548239938517104187698145666387735e-05L, -2.9384282290753706235208604777002963102e-05L,
+  -1.1641402067449950376895791995913672125e-05L, -4.5885348075273868255721924655343445906e-06L, -1.8005627955081458322204028649620350662e-06L, -7.0378941219347833214067471222239770590e-07L,
+  -2.7413963540482731185045932620331377352e-07L, -1.0645313231320808326024667350792279852e-07L, -4.1223072448771156559318807603116419528e-08L, -1.5923376898615004128677660806663065530e-08L,
+  -6.1368298043116345769816086674174450569e-09L, -2.3602323152914347699033314033408744848e-09L, -9.0603229062698346039479269140561762700e-10L, -3.4719859662410051486654409365217142276e-10L,
+  -1.3283631472964644271673440503045960993e-10L, -5.0747278046555248958473301297399200773e-11L, -1.9360320299432568426355237044475945959e-11L, -7.3766303773930764398798182830872768331e-12L,
+  -2.8072868906520523814747496670136120235e-12L, -1.0671679036256927021016406931839827582e-12L, -4.0525329757101362313986893299088308423e-13L, -1.5374324278841211301808010293913555758e-13L,
+  -5.8272886672428440854292491647585674200e-14L, -2.2067908660514462849736574172862899665e-14L, -8.3502821888768497979241489950570881481e-15L, -3.1572276215253043438828784344882603413e-15L,
+  -1.1928704609782512589094065678481294667e-15L, -4.5038074274761565346423524046963087756e-16L, -1.6993417021166355981316939131434632072e-16L, -6.4078169762734539491726202799339200356e-17L,
+  -2.4147993510032951187990399698408378385e-17L, -9.0950634616416460925150243301974013218e-18L, -3.4236981860988704669138593608831552044e-18L, -1.2881333612472271400115547331327717431e-18L,
+  -4.8440839844747536942311292467369300505e-19L, -1.8207788854829779430777944237164900798e-19L, -6.8407875971564885101695409345634890862e-20L, -2.5690139750480973292141845983878483991e-20L,
+  -9.6437492398195889150867140826350628738e-21L, -3.6186918227651991108814673877821174787e-21L, -1.3573451162272064984454015639725697008e-21L, -5.0894204288895982960223079251039701810e-22L,
+  -1.9076194289884357960571121174956927722e-22L, -7.1476978375412669048039237333931704166e-23L, -2.6773000149758626855152191436980592206e-23L, -1.0025115553818576399048488234179587012e-23L,
+  -3.7527362568743669891693429406973638384e-24L, -1.4043571811296963574776515073934728352e-24L, -5.2539064576179122030932396804434996219e-25L, -1.9650175744554348142907611432678556896e-25L,
   -7.3474021582506822389671905824031421322e-26L, -2.7465543000397410133825686340097295749e-26L, -1.0264389699511282259046957018510946438e-26L
 }; // wm1zs
 } // namespace lambert_w_lookup
diff --git a/include/boost/math/special_functions/lambert_w_pb.hpp b/include/boost/math/special_functions/lambert_w_pb.hpp
index 74eac4943..ed82f3259 100644
--- a/include/boost/math/special_functions/lambert_w_pb.hpp
+++ b/include/boost/math/special_functions/lambert_w_pb.hpp
@@ -8,7 +8,7 @@
 Implementation of the Fukushima algorithm for the Lambert W0 and W-1 real-only functions.
 
 This code is based on the algorithm by
-Toshio Fukushima, 
+Toshio Fukushima,
 "Precise and fast computation of Lambert W-functions without transcendental function evaluations",
 J. Comp. Appl. Math. 244 (2013) 77-89,
 and its author's FORTRAN code,
@@ -39,8 +39,8 @@ BOOST_MATH_INSTRUMENT_LAMBERT_WM1_LOOKUP // Show results from W-1 lookup table.
 //] [/boost_math_instrument_lambert_w_macros]
 */
 
-#ifndef BOOST_MATH_SF_LAMBERT_W_HPP
-#define BOOST_MATH_SF_LAMBERT_W_HPP
+#ifndef BOOST_MATH_FKDVPB_LAMBERT_W_HPP
+#define BOOST_MATH_FKDVPB_LAMBERT_W_HPP
 
 #ifdef _MSC_VER
 #  pragma warning (disable: 4127) // warning C4127: conditional expression is constant
@@ -75,8 +75,8 @@ BOOST_MATH_INSTRUMENT_LAMBERT_WM1_LOOKUP // Show results from W-1 lookup table.
 
 typedef double lookup_t; // Type for lookup table (double or float, or even long double?)
 
-//#include "J:\Cpp\Misc\lambert_w_lookup_table_generator\lambert_w_lookup_table.ipp"
- #include "lambert_w_lookup_table.ipp" // Boost.Math version.
+#include "J:\Cpp\Misc\lambert_w_lookup_table_generator\lambert_w_lookup_table.ipp" // BOTH w0 and wm1 version in namespace detail.
+// #include "lambert_w_lookup_table.ipp" // Boost.Math version.  Only Wm1 needed using namespace lambert_w_detail
 
 namespace boost
 {
@@ -85,6 +85,7 @@ namespace math
   namespace detail
   {
 
+
   //! \brief Series expansion used near the singularity/branch point z = -exp(-1) = -3.6787944.
   //! Wolfram InverseSeries[Series[sqrt[2(p Exp[1 + p] + 1)], {p,-1, 20}]]
   //! Wolfram command used to obtain 40 series terms at 50 decimal digit precision was
@@ -283,7 +284,7 @@ namespace math
   /////////////////////////////////////////////////////////////////////////////////////////////
 
   //! \brief Series expansion used near zero (abs(z) < 0.05).
-  //! \details 
+  //! \details
   //! Coefficients of the inverted series expansion of the Lambert W function around z = 0.
   //! Tosio Fukushima always uses all 17 terms of a Taylor series computed using Wolfram with
   //!   InverseSeries[Series[z Exp[z],{z,0,17}]]
@@ -295,7 +296,7 @@ namespace math
 
   //! This version is intended to allow use by user-defined types
   //! like Boost.Multiprecision quad and cpp_dec_float types.
-  //! The three specializations below for built-in float, double 
+  //! The three specializations below for built-in float, double
   //! (and perhaps long double) will be chosen in preference for these types.
 
   //! This version uses rationals computed by Wolfram as far as possible,
@@ -321,7 +322,7 @@ namespace math
   //! It assumes that max_digits10 is defined correctly or this might fail to make the correct selection.
   //! causing very small differences in computing lambert_w that would be very difficult to detect and diagnose.
   //! Cannot switch on @c std::numeric_limits<>::max() because comparison values may overflow the compiler limit.
-  //! Cannot switch on @c std::numeric_limits::max_exponent10() 
+  //! Cannot switch on @c std::numeric_limits::max_exponent10()
   //! because both 80 and 128 bit floating-point implmentations use 11 bits for the exponent.
   //! So must rely on @c std::numeric_limits::max_digits10.
 
@@ -699,7 +700,7 @@ namespace math
     using boost::math::tools::max_value;
 
     int iterations = 0;
-   // int iterations_required = 6; // 
+   // int iterations_required = 6; //
     int max_iterations = 20; // Only as a check on looping.
 
     T tolerance = boost::math::policies::get_epsilon();
@@ -730,7 +731,7 @@ namespace math
       return w0;
     }
     // relative does NOT work well, so check for improvement in diff instead.
-    //else if ((abs(diff) / w0) <= tolerance) 
+    //else if ((abs(diff) / w0) <= tolerance)
     //{
     //  return w0;
     //}
@@ -788,7 +789,7 @@ namespace math
   } // T halley_update(T w0, const T z)
 
 
-  //! Halley update version that does one update without a check, 
+  //! Halley update version that does one update without a check,
   //! and then more updates as necessary.
   // TODO Might use a template parameter to avoid duplication?
   template
@@ -881,7 +882,7 @@ namespace math
     //   result = schroeder_update(w, y);
     //where
     // w is estimate of Lambert W (from bisection or series).
-    // y is z * e^-w. 
+    // y is z * e^-w.
 
     BOOST_MATH_STD_USING // Aid argument dependent lookup of abs.
     #ifdef BOOST_MATH_INSTRUMENT_LAMBERT_W_SCHROEDER
@@ -919,7 +920,7 @@ namespace math
   } // template T schroeder_update(const T w, const T y)
 
   /////////////  namespace detail implementations of Lambert W for W0 and W-1 branches.
- 
+
   //! Compute Lambert W for W0 (or W+) branch.
   template
   T lambert_w0_imp(const T z, const Policy& pol)
@@ -931,7 +932,7 @@ namespace math
     // Integral types should been promoted to double by user Lambert w functions.
     BOOST_STATIC_ASSERT_MSG(!std::is_integral::value,
       "Must be floating-point or fixed type (not integer type), for example: lambert_w0(1.), not lambert_w0(1)!");
-     
+
     const char* function = "boost::math::lambert_w0()"; // Used for error messages.
 
 #ifdef BOOST_MATH_INSTRUMENT_LAMBERT_W0
@@ -1002,7 +1003,7 @@ namespace math
         T w_series = lambert_w_singularity_series(T(sqrt(p2)));
         // Neither Halley nor Schroeder do not improve result for builtin,
         // differences about 5e-15, so do NOT use refinement step below unless multiprecision.
-        // 
+        //
         //int d2 = policies::digits(); // Precision as digits base 2 from policy.
         //if (d2 > (std::numeric_limits::digits - 6))
         //{ // If more (fullest) accuracy required, then use refinement.
@@ -1022,11 +1023,14 @@ namespace math
       }
     } //  z < -0.35
  // z < -1/e
-    else 
+    else
     { //! Check if z so big that cannot use lookup table, so approximate and Halley iterate.
      // std::cout << " lambert_w_lookup::w0zs[63] "  << lambert_w_lookup::w0zs[63] << ' ' << lambert_w_lookup::w0zs[64]  << std::endl;
-      // lambert_w_lookup::w0zs[63] = 1.4450833904658542e+29,  lambert_w_lookup::w0zs[64]=  3.9904954117194348e+29, 
-      if (z >= lambert_w_lookup::w0zs[lambert_w_lookup::noof_w0zs - 1]) // F[64] =  3.9904954117194348e+29
+
+      using namespace boost::math::detail::lambert_w_lookup;
+
+      // lambert_w_lookup::w0zs[63] = 1.4450833904658542e+29,  lambert_w_lookup::w0zs[64]=  3.9904954117194348e+29,
+      if (z >= w0zs[noof_w0zs - 1]) // F[64] =  3.9904954117194348e+29
       { //! z so big that cannot use lookup table, so approximate and Halley iterate.
         T guess;
          // R.M.Corless, G.H.Gonnet, D.E.G.Hare, D.J.Jeffrey, and D.E.Knuth, “On the Lambert W function, ” Adv.Comput.Math., vol. 5, pp. 329–359, 1996.
@@ -1252,7 +1256,7 @@ namespace math
     } // normal range and precision.
     throw std::logic_error("No result from lambert_w0_imp!  (Please report!)");  // Should never get here.
   } //  T lambert_w0_imp(const T z, const Policy& pol)
- 
+
 
   //! Lambert W for W-1 branch, -max(z) < z <= -1/e.
   // TODO is -max(z) allowed?
@@ -1296,7 +1300,7 @@ namespace math
         return -std::numeric_limits::infinity();
       }
       else
-      { 
+      {
         return -tools::max_value();
       }
     }
@@ -1336,7 +1340,7 @@ namespace math
         z, pol);
     }
     if (z < static_cast(-0.35))
-    { // Close to singularity/branch point z = -0.3678794411714423215955237701614608727 but on W-1 branch. 
+    { // Close to singularity/branch point z = -0.3678794411714423215955237701614608727 but on W-1 branch.
       const T p2 = 2 * (boost::math::constants::e() * z + 1);
       if (p2 == 0)
       { // At the singularity at branch point.
@@ -1366,12 +1370,12 @@ namespace math
       //  z, pol);
     } // if (z < -0.35)
 
-    using lambert_w_lookup::wm1es;
-    using lambert_w_lookup::wm1zs;
-    using lambert_w_lookup::noof_wm1zs; // size == 64 
+    using boost::math::detail::lambert_w_lookup::wm1es;
+    using boost::math::detail::lambert_w_lookup::wm1zs;
+    using boost::math::detail::lambert_w_lookup::noof_wm1zs; // size == 64
 
     // std::cout <<" Wm1zs[63] (== G[64]) = " << " " << wm1zs[63] << std::endl; // Wm1zs[63] (== G[64]) =  -1.0264389699511283e-26
-    // Check that z argument value is not smaller than lookup_table G[64] 
+    // Check that z argument value is not smaller than lookup_table G[64]
     // std::cout << "(z > wm1zs[63]) = " << std::boolalpha << (z > wm1zs[63]) << std::endl;
 
     if (z >= wm1zs[63]) // wm1zs[63]  = -1.0264389699511282259046957018510946438e-26L  W = 64.00000000000000000
@@ -1452,7 +1456,7 @@ namespace math
     else // T is double or less precision.
     {
       // Use a lookup table to find the nearest integer part of Lambert W as starting point for Bisection.
-      using namespace boost::math::detail::lambert_w_lookup;
+      using boost::math::detail::lambert_w_lookup::wm1zs;
       // Bracketing sequence  n = (2, 4, 8, 16, 32, 64) for W-1 branch. (0 is -infinity)
       // Since z is probably quite small, start with lowest n (=2).
       int n = 2;
@@ -1552,8 +1556,8 @@ namespace math
       // (Avoiding using exponential function for speed).
       // Only use @c lookup_t precision, default double, for bisection (again for speed),
       // and use later Halley refinement for higher precisions.
-      using lambert_w_lookup::halves;
-      using lambert_w_lookup::sqrtwm1s;
+      using detail::lambert_w_lookup::halves;
+      using detail::lambert_w_lookup::sqrtwm1s;
       lookup_t w = -static_cast(n); // Equation 25,
       lookup_t y = static_cast(z * wm1es[n - 1]); // Equation 26,
       // Perform the bisections fractional bisections for necessary precision.
@@ -1604,6 +1608,9 @@ namespace math
   }
 } // namespace detail ////////////////////////////////////////////////////////////
 
+// Put user functions in a separate namespace.
+
+namespace FKDVPB {
   // User Lambert W0 and Lambert W-1 functions.
 
   //! W0 branch, -1/e <= z < max(z)
@@ -1656,7 +1663,9 @@ namespace math
     return detail::lambert_wm1_imp(result_type(z), policies::policy<>());
   } // lambert_wm1(T z)
 
+  } // namespace FKDVPB {
+
   } // namespace math
 } // namespace boost
 
-#endif //  #ifndef BOOST_MATH_SF_LAMBERT_W_HPP
+#endif //  #ifndef BOOST_MATH_FKDVPB_LAMBERT_W_HPP
diff --git a/test/test_lambert_w.cpp b/test/test_lambert_w.cpp
index c2798a55e..90e53fb78 100644
--- a/test/test_lambert_w.cpp
+++ b/test/test_lambert_w.cpp
@@ -40,9 +40,6 @@ using boost::math::float_prior;
 #include 
 using boost::math::policies::digits2;
 #include  // For Lambert W lambert_w function.
-
-// #include "C:\Users\Paul\Desktop\lambert_w0 - Copy before amalgamation.hpp"  // JM latest providing jm_lambert_w0 and lambert_wm1
-//#include "C:\Users\Paul\Desktop\lambert_w0.hpp"  // JM latest providing jm_lambert_w0 and lambert_wm1
 using boost::math::lambert_wm1;
 using boost::math::lambert_w0; // Use jm version instead.
 
@@ -109,7 +106,7 @@ std::cout << "MINGW64 " << __MINGW32_MAJOR_VERSION << __MINGW32_MINOR_VERSION <<
   message << "\n  STL " << BOOST_STDLIB;
   message << "\n  Boost version " << BOOST_VERSION / 100000 << "." << BOOST_VERSION / 100 % 1000 << "." << BOOST_VERSION % 100;
 
-#ifdef BOOST_HAS_FLOAT
+#ifdef BOOST_HAS_FLOAT128
   message << ",  BOOST_HAS_FLOAT128" << std::endl;
 #endif
   message << std::endl;
@@ -129,6 +126,7 @@ void test_spots(RealType)
   using boost::math::policies::policy;
 
   /*  Example of an exception-free 'ignore_all' policy (possibly ill-advised?).
+  */
   typedef policy <
     boost::math::policies::domain_error,
     boost::math::policies::overflow_error,
@@ -137,19 +135,24 @@ void test_spots(RealType)
     boost::math::policies::pole_error,
     boost::math::policies::evaluation_error
   > ignore_all_policy;
-  */
 
-//  Test some bad parameters to the function, with default policy and with ignore_all policy.
+//  Test some bad parameters to the function, with default policy and also with ignore_all policy.
 #ifndef BOOST_NO_EXCEPTIONS
   BOOST_CHECK_THROW(lambert_w0(-1.), std::domain_error);
   BOOST_CHECK_THROW(lambert_wm1(-1.), std::domain_error);
   BOOST_CHECK_THROW(lambert_w0(std::numeric_limits::quiet_NaN()), std::domain_error); // Would be NaN.
+  //BOOST_CHECK_EQUAL(lambert_w0(std::numeric_limits::quiet_NaN(), ignore_all_policy()), std::numeric_limits::quiet_NaN()); // Should be NaN.
+  // Fails as NaN != NaN by definition.
+  BOOST_CHECK(boost::math::isnan(lambert_w0(std::numeric_limits::quiet_NaN(), ignore_all_policy())));
+  //BOOST_MATH_CHECK_EQUAL(boost::math::lambert_w0(std::numeric_limits::infinity(), ignore_all_policy()), std::numeric_limits(std::numeric_limits::infinity()), std::domain_error); // Was if infinity should throw, now infinity.
   BOOST_CHECK_THROW(lambert_w0(-static_cast(0.4)), std::domain_error); // Would be complex.
-#else // No exceptions so set policy to ignore and check result is NaN
-  BOOST_MATH_CHECK_EQUAL(boost::math::lambert_w0(std::numeric_limits::quiet_NaN(), ignore_all_policy), std::numeric_limits(std::numeric_limits::infinity(), ignore_all_policy), std::numeric_limits(std::numeric_limits::infinity(), ignore_all_policy), std::numeric_limits(std::numeric_limits::quiet_NaN(), ignore_all_policy()), std::numeric_limits(std::numeric_limits::infinity(), ignore_all_policy()), std::numeric_limits(std::numeric_limits::infinity(), ignore_all_policy()), std::numeric_limits