From f51d987acd21f4635bf8d75373b6bb4c14280577 Mon Sep 17 00:00:00 2001 From: pabristow Date: Fri, 27 Oct 2017 18:18:06 +0100 Subject: [PATCH] added much on W-1 branch (handling tiny z), and more docs. --- doc/distributions/non_members.qbk | 75 ++-- doc/html/index.html | 2 +- doc/html/indexes/s01.html | 19 +- doc/html/indexes/s02.html | 2 +- doc/html/indexes/s03.html | 2 +- doc/html/indexes/s04.html | 6 +- doc/html/indexes/s05.html | 37 +- doc/html/math_toolkit/conventions.html | 2 +- doc/html/math_toolkit/navigation.html | 2 +- .../changing_policy_defaults.html | 6 +- doc/html/math_toolkit/tuning.html | 12 +- doc/math.qbk | 9 +- doc/sf/lambert_w.qbk | 140 ++++++- .../math/special_functions/lambert_w.hpp | 393 +++++++++++------- test/lambert_w_high_reference_values.cpp | 16 +- test/test_lambert_w.cpp | 186 ++++++++- test/test_lambert_w0_precision_high.cpp | 8 +- test/test_lambert_w0_precision_low.cpp | 5 +- 18 files changed, 677 insertions(+), 245 deletions(-) diff --git a/doc/distributions/non_members.qbk b/doc/distributions/non_members.qbk index 653d24979..8937b3141 100644 --- a/doc/distributions/non_members.qbk +++ b/doc/distributions/non_members.qbk @@ -1,11 +1,11 @@ [section:nmp Non-Member Properties] -Properties that are common to all distributions are accessed via non-member +Properties that are common to all distributions are accessed via non-member getter functions: non-membership allows more of these functions to be added over time, as the need arises. Unfortunately the literature uses many different and confusing names to refer to a rather small number of actual concepts; refer -to the [link math_toolkit.dist_ref.nmp.concept_index concept index] to find the property you -want by the name you are most familiar with. +to the [link math_toolkit.dist_ref.nmp.concept_index concept index] to find the property you +want by the name you are most familiar with. Or use the [link math_toolkit.dist_ref.nmp.function_index function index] to go straight to the function you want if you already know its name. @@ -62,8 +62,8 @@ to go straight to the function you want if you already know its name. template RealType cdf(const ``['Distribution-Type]``& dist, const RealType& x); - -The __cdf is the probability that + +The __cdf is the probability that the variable takes a value less than or equal to x. It is equivalent to the integral from -infinity to x of the __pdf. @@ -79,11 +79,11 @@ normal distribution: template RealType cdf(const ``['Unspecified-Complement-Type]``& comp); - -The complement of the __cdf -is the probability that + +The complement of the __cdf +is the probability that the variable takes a value greater than x. It is equivalent -to the integral from x to infinity of the __pdf, or 1 minus the __cdf of x. +to the integral from x to infinity of the __pdf, or 1 minus the __cdf of x. This is also known as the survival function. @@ -118,7 +118,7 @@ the defined range for the distribution. [equation hazard] [caution -Some authors refer to this as the conditional failure +Some authors refer to this as the conditional failure density function rather than the hazard function.] [h4:chf Cumulative Hazard Function] @@ -133,14 +133,14 @@ the defined range for the distribution. [equation chf] -[caution +[caution Some authors refer to this as simply the "Hazard Function".] [h4:mean mean] template RealType mean(const ``['Distribution-Type]``& dist); - + Returns the mean of the distribution /dist/. This function may return a __domain_error if the distribution does not have @@ -150,14 +150,14 @@ a defined mean (for example the Cauchy distribution). template RealType median(const ``['Distribution-Type]``& dist); - + Returns the median of the distribution /dist/. [h4:mode mode] template RealType mode(const ``['Distribution-Type]``& dist); - + Returns the mode of the distribution /dist/. This function may return a __domain_error if the distribution does not have @@ -167,14 +167,14 @@ a defined mode. template RealType pdf(const ``['Distribution-Type]``& dist, const RealType& x); - -For a continuous function, the probability density function (pdf) returns -the probability that the variate has the value x. -Since for continuous distributions the probability at a single point is actually zero, + +For a continuous function, the probability density function (pdf) returns +the probability that the variate has the value x. +Since for continuous distributions the probability at a single point is actually zero, the probability is better expressed as the integral of the pdf between two points: see the __cdf. -For a discrete distribution, the pdf is the probability that the +For a discrete distribution, the pdf is the probability that the variate takes the value x. This function may return a __domain_error if the random variable is outside @@ -188,14 +188,14 @@ For example, for a standard normal distribution the pdf looks like this: template std::pair range(const ``['Distribution-Type]``& dist); - + Returns the valid range of the random variable over distribution /dist/. [h4:quantile Quantile] template RealType quantile(const ``['Distribution-Type]``& dist, const RealType& p); - + The quantile is best viewed as the inverse of the __cdf, it returns a value /x/ such that `cdf(dist, x) == p`. @@ -217,7 +217,7 @@ See also [link math_toolkit.stat_tut.overview.complements complements]. template RealType quantile(const ``['Unspecified-Complement-Type]``& comp); - + This is the inverse of the __ccdf. It is calculated by wrapping the arguments in a call to the quantile function in a call to /complement/. For example: @@ -250,8 +250,8 @@ distribution: template RealType standard_deviation(const ``['Distribution-Type]``& dist); - -Returns the standard deviation of distribution /dist/. + +Returns the standard deviation of distribution /dist/. This function may return a __domain_error if the distribution does not have a defined standard deviation. @@ -260,7 +260,7 @@ a defined standard deviation. template std::pair support(const ``['Distribution-Type]``& dist); - + Returns the supported range of random variable over the distribution /dist/. The distribution is said to be 'supported' over a range that is @@ -274,7 +274,7 @@ Outside are uninteresting zones where the pdf is zero, and the cdf zero or unity template RealType variance(const ``['Distribution-Type]``& dist); - + Returns the variance of the distribution /dist/. This function may return a __domain_error if the distribution does not have @@ -284,7 +284,7 @@ a defined variance. template RealType skewness(const ``['Distribution-Type]``& dist); - + Returns the skewness of the distribution /dist/. This function may return a __domain_error if the distribution does not have @@ -294,7 +294,7 @@ a defined skewness. template RealType kurtosis(const ``['Distribution-Type]``& dist); - + Returns the 'proper' kurtosis (normalized fourth moment) of the distribution /dist/. kertosis = [beta][sub 2][space]= [mu][sub 4][space] / [mu][sub 2][super 2] @@ -326,7 +326,7 @@ a defined kurtosis. template RealType kurtosis_excess(const ``['Distribution-Type]``& dist); - + Returns the kurtosis excess of the distribution /dist/. kurtosis excess = [gamma][sub 2][space]= [mu][sub 4][space] / [mu][sub 2][super 2][space]- 3 = kurtosis - 3 @@ -334,7 +334,7 @@ kurtosis excess = [gamma][sub 2][space]= [mu][sub 4][space] / [mu][sub 2][super Where [mu][sub i][space] is the i'th central moment of the distribution, and in particular [mu][sub 2][space] is the variance of the distribution. -The kurtosis excess is a measure of the "peakedness" of a distribution, and +The kurtosis excess is a measure of the "peakedness" of a distribution, and is more widely used than the "kurtosis proper". It is defined so that the kurtosis excess of a normal distribution is zero. @@ -344,7 +344,7 @@ a defined kurtosis excess. Kurtosis excess can have a value from -2 to + infinity. kurtosis = kurtosis_excess +3; - + The kurtosis excess of a normal distribution is zero. [h4:cdfPQ P and Q] @@ -361,12 +361,12 @@ the __quantile. [h4:cdf_inv Inverse CDF Function.] -The inverse of the cumulative distribution function, is the same as the +The inverse of the cumulative distribution function, is the same as the __quantile. [h4:survival_inv Inverse Survival Function.] -The inverse of the survival function, is the same as computing the +The inverse of the survival function, is the same as computing the [link math_toolkit.dist_ref.nmp.quantile_c quantile from the complement of the probability]. @@ -380,13 +380,13 @@ while the term __pdf applies to continuous distributions. [h4:lower_critical Lower Critical Value.] The lower critical value calculates the value of the random variable -given the area under the left tail of the distribution. +given the area under the left tail of the distribution. It is equivalent to calculating the __quantile. -[h4: upper_critical Upper Critical Value.] +[h4:upper_critical Upper Critical Value.] The upper critical value calculates the value of the random variable -given the area under the right tail of the distribution. It is equivalent to +given the area under the right tail of the distribution. It is equivalent to calculating the [link math_toolkit.dist_ref.nmp.quantile_c quantile from the complement of the probability]. @@ -394,8 +394,7 @@ probability]. Refer to the __ccdf. -[endsect][/section:nmp Non-Member Properties] - +[endsect] [/section:nmp Non-Member Properties] [/ non_members.qbk Copyright 2006 John Maddock and Paul A. Bristow. diff --git a/doc/html/index.html b/doc/html/index.html index ab4cb0b7d..a02e24638 100644 --- a/doc/html/index.html +++ b/doc/html/index.html @@ -120,7 +120,7 @@ This manual is also available in -

Last revised: October 10, 2017 at 14:01:12 GMT

+

Last revised: October 27, 2017 at 17:06:21 GMT

diff --git a/doc/html/indexes/s01.html b/doc/html/indexes/s01.html index bb0d7a41a..e3e6e8fa5 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

@@ -307,6 +307,10 @@
  • +

    bisection

    +
    +
    • +

  • bool

    • @@ -1556,6 +1560,7 @@

  • Additional Implementation Notes

  • Error Handling Example

  • Examples Where Root Finding Goes Wrong

    • +

  • Lambert W function

  • Students t Distribution

    • @@ -1888,6 +1893,13 @@
  • +

    ln

    +
    +
    • +

  • location

  • +

    point

    +
    +
    • +

  • polar

    • @@ -2631,6 +2647,7 @@

  • Bessel Functions of the First and Second Kinds

  • expm1

  • Inverse Chi-Squared Distribution Bayes Example

    • +

  • Lambert W function

  • Modified Bessel Functions of the First and Second Kinds

  • The Incomplete Beta Function Inverses

    • diff --git a/doc/html/indexes/s02.html b/doc/html/indexes/s02.html index d3fa7c973..8c8d8c8b4 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 90073b6f9..6840371fc 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 1b76863d8..b09c7a093 100644 --- a/doc/html/indexes/s04.html +++ b/doc/html/indexes/s04.html @@ -24,7 +24,7 @@

    -Macro Index

    +Macro Index

    B F

    @@ -298,6 +298,10 @@
  • +

    BOOST_MATH_TEST_VALUE

    +
    +
    • +

  • BOOST_MATH_UNDERFLOW_ERROR_POLICY

    • diff --git a/doc/html/indexes/s05.html b/doc/html/indexes/s05.html index b84430d40..3704c87cf 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

    @@ -669,6 +669,10 @@
  • +

    bisection

    +
    +
    • +

  • bool

    • @@ -971,6 +975,10 @@
  • +

    BOOST_MATH_TEST_VALUE

    +
    +
    • +

  • BOOST_MATH_UNDERFLOW_ERROR_POLICY

    • @@ -1444,6 +1452,10 @@
  • +

    Calculation of the Type of the Result

    +
    +
    • +

  • called

    • @@ -4432,6 +4444,7 @@

  • Additional Implementation Notes

  • Error Handling Example

  • Examples Where Root Finding Goes Wrong

    • +

  • Lambert W function

  • Students t Distribution

    • @@ -4814,15 +4827,20 @@

  • accuracy

  • al

  • alpha

    • +

  • bisection

    • +

  • BOOST_MATH_TEST_VALUE

  • branch

  • constants

  • expression

  • float

    • +

  • infinity

  • interval

  • lambert_w0

  • lambert_wm1

    • +

  • ln

  • range

  • refinement

    • +

  • small

  • W

    • @@ -5034,6 +5052,13 @@
  • +

    ln

    +
    +
    • +

  • Locating Function Minima using Brent's algorithm

    • accuracy

      • @@ -5249,7 +5274,10 @@

    • Mathematical Constants

      -
      +

    • Mathematically Undefined Function Policies

      @@ -5901,6 +5929,10 @@
  • +

    point

    +
    +
    • +

  • poisson

    • @@ -6698,6 +6730,7 @@

  • Bessel Functions of the First and Second Kinds

  • expm1

  • Inverse Chi-Squared Distribution Bayes Example

    • +

  • Lambert W function

  • Modified Bessel Functions of the First and Second Kinds

  • The Incomplete Beta Function Inverses

    • diff --git a/doc/html/math_toolkit/conventions.html b/doc/html/math_toolkit/conventions.html index b63b41239..da31c6f67 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 fa6234915..9ecf45760 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/changing_policy_defaults.html b/doc/html/math_toolkit/pol_tutorial/changing_policy_defaults.html index 5c290f413..508491501 100644 --- a/doc/html/math_toolkit/pol_tutorial/changing_policy_defaults.html +++ b/doc/html/math_toolkit/pol_tutorial/changing_policy_defaults.html @@ -55,9 +55,9 @@ 

    #define BOOST_MATH_ASSERT_UNDEFINED_POLICY false
 #define BOOST_MATH_OVERFLOW_ERROR_POLICY errno_on_error
    -

    - may be ignored. -

    +
    • + may be ignored*. +

    The compiler command line shows:

    diff --git a/doc/html/math_toolkit/tuning.html b/doc/html/math_toolkit/tuning.html index 51f8ab00f..36c3c1aae 100644 --- a/doc/html/math_toolkit/tuning.html +++ b/doc/html/math_toolkit/tuning.html @@ -109,15 +109,15 @@ be stored as tables of floating point values instead. If mixed arithmetic is slow then add:

    -

    - #define BOOST_MATH_INT_TABLE_TYPE(RT, IT) RT -

    +
    1. + define BOOST_MATH_INT_TABLE_TYPE(RT, IT) RT +

    to boost/math/tools/user.hpp, otherwise the default of:

    -

    - #define BOOST_MATH_INT_TABLE_TYPE(RT, IT) IT -

    +
    1. + define BOOST_MATH_INT_TABLE_TYPE(RT, IT) IT +

    Set in boost/math/config.hpp is fine, and may well result in smaller code. diff --git a/doc/math.qbk b/doc/math.qbk index 4db9e2c34..9165e08cb 100644 --- a/doc/math.qbk +++ b/doc/math.qbk @@ -1,5 +1,5 @@ [book Math Toolkit - [quickbook 1.6] + [quickbook 1.7] [copyright 2006, 2007, 2008, 2009, 2010, 2012, 2013, 2014, 2017 Nikhar Agrawal, Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, Hubert Holin, Bruno Lalande, John Maddock, Jeremy Murphy, Johan RÃ¥de, Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg, Daryle Walker and Xiaogang Zhang] [/purpose ISBN 0-9504833-2-X 978-0-9504833-2-0, Classification 519.2-dc22] [license @@ -219,8 +219,11 @@ and use the function's name as the link text.] [/Lambert W function] [def __lambert_w [link math_toolkit.sf_lambert_w.lambert_w_function lambert_w]] [def __lambert_w_implementation [link math_toolkit.lambert_w.implementation implementation]] -[def __lambert_w_implementation [link math_toolkit.lambert_w.precision precision]] -[def __lambert_w_implementation [link math_toolkit.lambert_w.examples examples]] +[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]] + + [/Airy Functions] [def __airy_ai [link math_toolkit.airy.ai airy_ai]] diff --git a/doc/sf/lambert_w.qbk b/doc/sf/lambert_w.qbk index 4d262d04f..306cba4a4 100644 --- a/doc/sf/lambert_w.qbk +++ b/doc/sf/lambert_w.qbk @@ -85,7 +85,6 @@ and to control what action to take on errors: [tip Always use [^try'n'catch] blocks to ensure you get an error message like this:] - [lambert_w_simple_examples_error_message_1] [lambert_w_simple_examples_out_of_range] @@ -189,6 +188,37 @@ or defined in a jamfile.v2: `BOOST_MATH_INSTRUMENT_LAMBERT_W_TERMS` BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES_ITERATIONS // Show result of estimates where small z is near zero. `` +[h4:edge_cases Edge and Corner cases] + +[h5:w0_edges w0 Branch] + +[h5:wm1_edges w-1 Branch] + +The range mathematically is -exp(-1) <= z <=0, but numerically, + +* `z == -exp(-1)` (for the type T) returns exactly -1. +* `z == 0` returns infinity (or the nearest equivalent if `std::has_infinity == false`). +* `z == std::numeric_limits::min()` returns the maximum possible value of Lambert W for the type T.[br] +For example, for `double`: lambert_wm1(-2.2250738585072014e-308) = -714.96865723796634[br] +and for `float`: lambert_wm1(-1.17549435e-38) = -91.8567734 + +* `z < std::numeric_limits::min()`, means that z is denormalized (if `std::numeric_limits::has_denorm_min == true`), +for example: + + r = lambert_wm1(std::numeric_limits::denorm_min()); + +and will give a message like: + + Error in function boost::math::lambert_wm1(): + Argument z = -4.9406564584124654e-324 is too small + (z < -std::numeric_limits::min so denormalized) for Lambert W-1 branch! + +Denormalized values of z are not supported (because not all floating-point types denormalize), +and anyway it only covers a tiny fraction of the range of possible z arguments values. + + + + [h4:implemention Implementation] There are many previous implementations with increasing accuracy and speed. See references below. @@ -253,8 +283,8 @@ 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. -Many applications of the Lambert W function make many repeated evaluations for Monte Carlo methods, -for which applications speed is important. +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. Fukushima makes the important observation that much of the execution time of all @@ -274,15 +304,27 @@ as a starting point followed by refinement using __halley iterations or other me For this reason, the lookup tables and `bisection` are only carried out at low precision, usually `double`, chosen by the `typedef double lookup_t`. Unlike the FORTRAN version, -the lookup tables of Lambert_W of integral values are precomputed as C++ static arrays. +the lookup tables of Lambert_W of integral values are precomputed as C++ static arrays of floating-point literals. +The default is a `typedef` setting the type to `double`. +To allow users to vary the precision from `float` to `long double` these are computed to +128-bit precision to ensure that even platforms with `long double` do not lose precision. -For the less well-behaved regions for ['z] arguments near zero and near the branch singularity at -1/e, -some series functions are provided. +The FORTRAN version and translation only permits the z argument to be the largest +items in these lookup arrays, `wm0s[64] = 3.99049`, producing an error message and returning `NaN`. +So 64 is the largest possible value returned from the `lambert_w0` function. +This is far from the `std::numeric_limits<>::max()` for even `float`s. + +For the less well-behaved regions for Lambert W0 ['z] arguments near zero, +and near the branch singularity at -1/e, some series functions are provided. + +For Lambert W-1 branch, tiny values very near zero where W > 64 cannot be computed using the lookup table. +For this region, an approximation followed by a few (usually 3) Halley refinements. +See __lambert_w_wm1_near_zero. So, as usual, there are compromises to consider between execution speed and accuracy. This implementation allows users to get closer than Fukushima, at some small loss of speed (but does not guarantee the nearest representable result) -or to get faster but less precise evalations controlled by __policies. +or to get faster but less precise evalutions controlled by __precision_policy. [h4:small_z Small values of argument z near zero] @@ -297,7 +339,7 @@ computed using Wolfram with InverseSeries[Series[z Exp[z],{z,0,17}]] -See Tosio Fukushima, Journal of Computational and Applied Mathematics 244 (2013) page 86. +See Tosio Fukushima, Journal of Computational and Applied Mathematics 244 (2013), page 86. To provide higher precision constants (34 decimal digits) for types larger than `long double`, @@ -350,18 +392,18 @@ which are near to the singularity at `z = -exp(-1) = -3.6787944` where the branc T. Fukushima / Journal of Computational and Applied Mathematics 244 (2013) Page 85, Table 3 described using Wolfram -InverseSeries\[Series\[sqrt\[2(p Exp\[1 + p\] + 1)\], {p,-1, 20}\]\] + InverseSeries\[Series\[sqrt\[2(p Exp\[1 + p\] + 1)\], {p,-1, 20}\]\] to provide Table 3. This implementation used Wolfram to obtain 40 series terms at 50 decimal digit precision - +`` N\[InverseSeries\[Series\[Sqrt\[2(p Exp\[1 + p\] + 1)\], { p,-1,40 }\]\], 50\] -1+p-p^2/3+(11 p^3)/72-(43 p^4)/540+(769 p^5)/17280-(221 p^6)/8505+(680863 p^7)/43545600 ... - -These constants are computed at compile time for the full precision for any RealType T -using original rationals from Table 3. +`` +These constants are computed at compile time for the full precision for any `RealType T` +using the original rationals from Fukushima Table 3. Longer decimal digits strings are rationals pre-evaluated using Wolfram. Some integer constants overflow, so use largest size available, suffixed by `uLL`. @@ -379,7 +421,74 @@ 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].] -Fukushima used 20 series terms and it was confirmed that using more terms does not usefully increase accuracy. +Fukushima implmentation used 20 series terms and it was confirmed that using more terms does not usefully increase accuracy. + +[h5:wm1_near_zero Lambert W-1 arguments values very near zero.] + +The lookup tables of Fukushima have only 64 elements, so that the z argument nearest zero is -1.0264389699511303e-26, +that corresponds to a maximum Lambert W-1 value of 64.0. +His implementation did not cater for z argument values that are smaller (nearer to zero), +but this implementation adds code to accept smaller (but not denormalised) values of z. +A crude approximation for these very small values is to take the exponent and multiply by ln[10] ~= 2.3. +We also tried the approximation first proposed by Corless et al. using ln(-z), and then tried improving by a 2nd term -ln(ln(-z)), +and finally the ratio term -ln(ln(-z))/ln(-z). + +For a z very close to z = -1.0264389699511303e-26 when W = 64, when effect of ln(ln(-z) term, and ratio L1/L2 is greatest, +the possible 'guessed' are + + z = -1.e-26, w = -64.02, guess = -64.0277, ln(-z) = -59.8672, ln(-ln(-z) = 4.0921, llz/lz = -0.0684 + +whereas at the minimum (unnormalized) z + +z = -2.2250e-308, w = -714.9, guess = -714.9687, ln(-z) = -708.3964, ln(-ln(-z) = 6.5630, llz/lz = -0.0092 + +Although the addition of the 3rd ratio term did not reduce the number of Halley iterations needed, +it allows return of a better low precision estimate [*without any Halley iterations] using the __precision_policy. +For the worst case near w = 64, the error in the 'guess' is 0.008, ratio 0.0001 or 1 in 10,000 digits 10 ~= 4, +so we opt to skip Halley iterations and return the 'guess' if digits2 < 12. +Two log evalutations are still needed, but is probably over an order of magnitude faster. + +Halley's method was then used to refine the estimate of Lambert W-1 from this guess. +Experiments showed that although all approximations reached with __ulp of the closest representable value, +the computational cost of the log functions was easily paid by far fewer iterations +(typically from 8 down to 4 iterations for double or float). For example: + + using boost::math::policies::digits2; + // Very close to z = -1.0264389699511303e-26 when W = 64, when effect of ln(ln(-z) term and ratio L1/L2 is greatest. + double z = z = -1e-26; + double r = lambert_wm1(z, policy >()); + // lambert_wm1(-1.0000000000000000e-26) = -64.026509628385895 + +[h4:edge_cases Edge and Corner cases] + +[h5:w0_edges w0 Branch] + +[h5:wm1_edges w-1 Branch] + +The range mathematically is -exp(-1) <= z <=0, but numerically, + +* `z == -exp(-1)` (for the type T) returns exactly -1. +* `z == 0` returns infinity (or the nearest equivalent if `std::has_infinity == false`). +* `z == std::numeric_limits::min()` returns the maximum possible value of Lambert W for the type T. [br] +For example, for `double`: lambert_wm1(-2.2250738585072014e-308) = -714.96865723796634 [br] +and for `float`: lambert_wm1(-1.17549435e-38) = -91.8567734 [br] + +* `z < std::numeric_limits::min()`, means that z is denormalized (if `std::numeric_limits::has_denorm_min == true`), +for example: + + r = lambert_wm1(std::numeric_limits::denorm_min()); + +and will give a message like: + + Error in function boost::math::lambert_wm1(): + Argument z = -4.9406564584124654e-324 is too small + (z < -std::numeric_limits::min so denormalized) for Lambert W-1 branch! + +Denormalized values of z are not supported (because not all floating-point types denormalize), +and anyway it only covers a tiny fraction of the range of possible z arguments values. + + + [h4:testing Testing] @@ -447,6 +556,9 @@ this version by C++ version by John Burkardt. Initial guesses based on: +# 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 to 359, (1996). + # D.A. Barry, J.-Y. Parlange, L. Li, H. Prommer, C.J. Cunningham, and F. Stagnitti. Analytical approximations for real values of the Lambert W-function. Mathematics and Computers in Simulation, 53(1), 95-103 (2000). diff --git a/include/boost/math/special_functions/lambert_w.hpp b/include/boost/math/special_functions/lambert_w.hpp index cdfbdde70..b250d4411 100644 --- a/include/boost/math/special_functions/lambert_w.hpp +++ b/include/boost/math/special_functions/lambert_w.hpp @@ -21,6 +21,9 @@ BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES BOOST_MATH_INSTRUMENT_LAMBERT_W0 // W0 branch diagnostics. BOOST_MATH_INSTRUMENT_LAMBERT_W0 // W1 branch diagnostics. BOOST_MATH_INSTRUMENT_LAMBERT_W_HALLEY // Halley refinement diagnostics. +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_TINY // K > 64, z > -1.0264389699511303e-26 BOOST_MATH_INSTRUMENT_LAMBERT_W_SCHROEDER // Schroeder refinement diagnostics. BOOST_MATH_INSTRUMENT_LAMBERT_W_TERMS // Number of terms used for near-singularity series. BOOST_MATH_INSTRUMENT_LAMBERT_W0_NOT_BUILTIN // higher than built-in precision types approximation and refinement. @@ -39,7 +42,6 @@ BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES_ITERATIONS // Show evaluation of #include #include -#include // For exp_minus_one == 3.67879441171442321595523770161460867e-01. #include #include #include // for log (1 + x) @@ -51,6 +53,8 @@ BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES_ITERATIONS // Show evaluation of #include // evaluate_polynomial. #include #include +#include // boost::math::tools::max_value(). +#include // Macro BOOST_MATH_TEST_VALUE. #include #include @@ -63,60 +67,21 @@ BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES_ITERATIONS // Show evaluation of #include // For float_distance. //#include // For create_test_value and macro BOOST_MATH_TEST_VALUE. +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" // TODO copy final version. namespace boost { namespace math -{ +{ + namespace detail { - namespace lambert_w_lookup - { - - typedef double lookup_t; // Type for lookup table (double or float?) - - BOOST_STATIC_CONSTEXPR std::size_t noof_wm1s = 66; - BOOST_STATIC_CONSTEXPR lookup_t e[noof_wm1s] = // lambert_w[k] for W-1 branch. 2.7182 1. 0.3678 0.1353 0.04978 ... 4.359e-28 1.603e-28 - { - 2.7182818284590452, 1., 0.36787944117144232, 0.13533528323661269, 0.049787068367863943, 0.01831563888873418, 0.0067379469990854671, - 0.0024787521766663584, 0.00091188196555451621, 0.00033546262790251184, 0.00012340980408667955, 4.5399929762484852e-05, 1.6701700790245659e-05, - 6.1442123533282098e-06, 2.2603294069810543e-06, 8.3152871910356788e-07, 3.0590232050182579e-07, 1.1253517471925911e-07, 4.1399377187851667e-08, - 1.5229979744712628e-08, 5.6027964375372675e-09, 2.0611536224385578e-09, 7.5825604279119067e-10, 2.7894680928689248e-10, 1.026187963170189e-10, - 3.7751345442790977e-11, 1.3887943864964021e-11, 5.1090890280633247e-12, 1.8795288165390833e-12, 6.914400106940203e-13, 2.5436656473769229e-13, - 9.3576229688401746e-14, 3.4424771084699765e-14, 1.2664165549094176e-14, 4.6588861451033974e-15, 1.713908431542013e-15, 6.3051167601469894e-16, - 2.3195228302435694e-16, 8.5330476257440658e-17, 3.1391327920480296e-17, 1.1548224173015786e-17, 4.248354255291589e-18, 1.5628821893349888e-18, - 5.7495222642935598e-19, 2.1151310375910805e-19, 7.7811322411337965e-20, 2.8625185805493936e-20, 1.0530617357553812e-20, 3.8739976286871871e-21, - 1.4251640827409351e-21, 5.2428856633634639e-22, 1.9287498479639178e-22, 7.0954741622847041e-23, 2.6102790696677048e-23, 9.602680054508676e-24, - 3.532628572200807e-24, 1.2995814250075031e-24, 4.7808928838854691e-25, 1.7587922024243116e-25, 6.4702349256454603e-26, 2.3802664086944006e-26, - 8.7565107626965203e-27, 3.2213402859925161e-27, 1.185064864233981e-27, 4.359610000063081e-28, 1.6038108905486378e-28 - }; - - BOOST_STATIC_CONSTEXPR size_t noof_w0s = 65; - BOOST_STATIC_CONSTEXPR lookup_t g[noof_w0s] = // lambert_w[k] for W0 branch. 0 2.7182 14.77 60.2566 ... 1.445e+29 3.990e+29. - { 0., 2.7182818284590452, 14.7781121978613, 60.256610769563003, 218.39260013257696, 742.06579551288302, 2420.5727609564107, 7676.4321089992102, - 23847.663896333826, 72927.755348178456, 220264.65794806717, 658615.558867176, 1953057.4970280471, 5751374.0961159665, 16836459.978306875, 49035260.58708166, - 142177768.32812596, 410634196.81078007, 1181879444.4719492, 3391163718.300558, 9703303908.1958056, 27695130424.147509, 78868082614.895014, 224130479263.72476, - 635738931116.24334, 1800122483434.6468, 5088969845149.8079, 14365302496248.563, 40495197800161.305, 114008694617177.22, 320594237445733.86, - 900514339622670.18, 2526814725845782.2, 7083238132935230.1, 19837699245933466., 55510470830970076., 1.5520433569614703e+17, 4.3360826779369662e+17, - 1.2105254067703227e+18, 3.3771426165357561e+18, 9.4154106734807994e+18, 2.6233583234732253e+19, 7.3049547543861044e+19, 2.032970971338619e+20, - 5.6547040503180956e+20, 1.5720421975868293e+21, 4.3682149334771265e+21, 1.2132170565093317e+22, 3.3680332378068632e+22, 9.3459982052259885e+22, - 2.5923527642935362e+23, 7.1876803203773879e+23, 1.99212416037262e+24, 5.5192924995054165e+24, 1.5286067837683347e+25, 4.2321318958281094e+25, - 1.1713293177672778e+26, 3.2408603996214814e+26, 8.9641258264226028e+26, 2.4787141382364034e+27, 6.8520443388941057e+27, 1.8936217407781711e+28, - 5.2317811346197018e+28, 1.4450833904658542e+29, 3.9904954117194348e+29 - }; - - BOOST_STATIC_CONSTEXPR std::size_t noof_sqrts = 12; - BOOST_STATIC_CONSTEXPR lookup_t a[noof_sqrts] = // 0.6065 0.7788, 0.8824 ... 0.9997, sqrt of previous elements. - { - 0.60653065971263342, 0.77880078307140487, 0.8824969025845954, 0.93941306281347579, 0.96923323447634408, 0.98449643700540841, - 0.99221793826024351, 0.99610136947011749, 0.99804878110747547, 0.99902391418197566, 0.99951183793988937, 0.99975588917489722 - }; - BOOST_STATIC_CONSTEXPR lookup_t b[noof_sqrts] = // 0.5 0.25 0.125, 0.0625 ... 0.0002441, halves of previous elements. - { 0.5, 0.25, 0.125, 0.0625, 0.03125, 0.015625, 0.0078125, 0.00390625, 0.001953125, 0.0009765625, 0.00048828125, 0.000244140625 }; - } // namespace Lambert_w_lookup //! Series expansion used near the singularity/branch point z = -exp(-1) = -3.6787944. - // Some integer constants overflow so use largest size available. + // Some integer constants overflow so use largest integer size available (LL). // Wolfram InverseSeries[Series[sqrt[2(p Exp[1 + p] + 1)], {p,-1, 20}]] // T. Fukushima / Journal of Computational and Applied Mathematics 244 (2013) Page 85, Table 3. // Wolfram command used to obtain 40 series terms at 50 decimal digit precision was @@ -724,19 +689,19 @@ namespace math // since we can simplify the expressions algebraically, // and don't need most of the error checking of the boilerplate version // as we know in advance that the function is reasonably well-behaved, - // and in any case the derivatives require evaluation of Lambert W! + // (and in any case the derivatives require evaluation of Lambert W!) BOOST_MATH_STD_USING // Aid argument dependent (ADL) lookup of abs. using boost::math::constants::exp_minus_one; // 0.36787944 using boost::math::tools::max_value; - T tolerance = boost::math::policies::get_epsilon(); + T tolerance = boost::math::policies::get_epsilon(); // TODO should get_precision here. int iterations = 0; int iterations_required = 6; - int max_iterations = 10; + int max_iterations = 20; T w1 = w0; // Refined estimate. T previous_diff = boost::math::tools::max_value(); @@ -751,7 +716,7 @@ namespace math std::cout.precision(std::numeric_limits::max_digits10); // Show all posssibly significant digits. std::ios::fmtflags flags(std::cout.flags()); // Save. std::cout.setf(std::ios_base::showpoint); // Include any trailing zeros. - std::cout << "w = " << w0 << ", z = " << z << ", exp(w) = " << expw0 << ", diff = " << diff << std::endl; + std::cout << "w = " << w0 << ", z = " << z << ", w * exp(w) = " << w0 * expw0 << ", diff = " << diff << std::endl; std::cout.precision(precision); // Restore. std::cout.flags(flags); #endif // BOOST_MATH_INSTRUMENT_LAMBERT_W_HALLEY @@ -779,13 +744,16 @@ if (diff == 0) // Exact result - common. w1 = w0 // Refine a new estimate using Halley's method using Luu equation 6.39. - diff / ((expw0 * (w0 + 1) - (w0 + 2) * diff / (w0 + w0 + 2))); - diff = (w1 * exp(w1)) - z; #ifdef BOOST_MATH_INSTRUMENT_LAMBERT_W_HALLEY T dis = boost::math::float_distance(w0, w1); + std::cout << "w = " << w0 << ", z = " << z << ", w * exp(w) = " << w0 * expw0 << ", diff = " << diff << std::endl; int d = static_cast(dis); - std::cout << "float_distance = " << d << std::endl; + if (abs(dis) < 2147483647) + { + std::cout << "float_distance = " << d << std::endl; + } #endif // BOOST_MATH_INSTRUMENT_LAMBERT_W_HALLEY if (diff == 0) // Exact. @@ -803,7 +771,7 @@ if (diff == 0) // Exact result - common. w0 = w1; expw0 = exp(w0); } - while ((iterations < iterations_required) || (iterations <= max_iterations)); // Absolute limit during testing - looping if need this. + while ((iterations < iterations_required) || (iterations <= max_iterations)); // Absolute limit during testing - is looping if need this. return w1; } // T halley_update(T w0, const T z) @@ -951,7 +919,7 @@ if (diff == 0) // Exact result - common. // or static_casted integer, for example: static_cast(1) or static_cast(1). // Want to allow fixed_point types too, so do not just test for floating-point. // 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: W(1.), not W(1)!"); + 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)!"); #ifdef BOOST_MATH_INSTRUMENT_LAMBERT_W0 { @@ -964,25 +932,32 @@ if (diff == 0) // Exact result - common. } #endif // BOOST_MATH_INSTRUMENT_LAMBERT_W0 - // Test for edge and corner cases first, + const char* function = "boost::math::lambert_w0()"; // Used for error messages. + // Test for edge and corner cases first: if (boost::math::isnan(z)) { - const char* function = "boost::math::lambert_w0()"; return policies::raise_domain_error(function, "Argument z is NaN!", z, pol); } if (boost::math::isinf(z)) { - const char* function = "boost::math::lambert_w0()"; - return policies::raise_domain_error(function, - "Argument z is infinite!", - z, pol); + if (std::numeric_limits::has_infinity) + { + return std::numeric_limits::infinity(); + } + else + { // Arbitrary precision type. + return tools::max_value(); + } + // Or might opt to throw an exception? + //return policies::raise_domain_error(function, + // "Argument z is infinite!", + // z, pol); } if (z > (std::numeric_limits::max)()/4) - { // TODO not sure this is right - what is largest possible for which can compute Lambert_w? - const char* function = "boost::math::lambert_w0()"; + { // TODO not sure this is right - what is largest possible for which can compute Lambert_w0? return policies::raise_domain_error(function, "Argument z %1% is too big!", z, pol); @@ -1008,13 +983,12 @@ if (diff == 0) // Exact result - common. } else if (z < -boost::math::constants::exp_minus_one()) // z < -1/e so out of range of W0 branch (should using W1 branch instead). { - const char* function = "boost::math::lambert_w0()"; return policies::raise_domain_error(function, "Argument z = %1% out of range (-1/e <= z < (std::numeric_limits::max)()) for Lambert W0 branch (use W-1 branch?).", z, pol); } else if (z < static_cast(-0.35)) - { // Near singularity/branch point at z = -0.36787944... + { // Near singularity/branch point at z = -0.3678794411714423215955237701614608727 const T p2 = 2 * (boost::math::constants::e() * z + 1); if (p2 > 0) { // Use near-singularity series expansion. @@ -1132,22 +1106,23 @@ if (diff == 0) // Exact result - common. // Perform additional Halley refinement(s) to ensure that // get a near as possible to correct result (usually +/- one epsilon). T result = halley_update(double_approx, z, pol); -#ifdef BOOST_MATH_INSTRUMENT_LAMBERT_W0_NOT_BUILTIN - std::cout.precision(std::numeric_limits::max_digits10); +#ifdef BOOST_MATH_INSTRUMENT_LAMBERT_W0 + std::streamsize saved_precision = std::cout.precision(std::numeric_limits::max_digits10); std::cout << "Result " << typeid(T).name() << " precision Halley refinement = " << result << std::endl; + std::cout.precision(saved_precision); #endif // BOOST_MATH_INSTRUMENT_LAMBERT_W0 return result; - } // digits > 53 + } // digits > 53 - higher precision than double. 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; // Test sequence is n is (0, 1, 2, 4, 8, 16, 32, 64) for W0 branch. // Since z is probably quite small, start with lowest (n = 0). - int n; // Indexing W0 lookup table g[n] of z. + int n; // Indexing W0 lookup table w0s[n] of z. for (n = 0; n <= 2; ++n) { // Try 1st few. - if (g[n] > z) + if (w0s[n] > z) { // goto bisect; } @@ -1156,13 +1131,13 @@ if (diff == 0) // Exact result - common. for (int j = 1; j <= 5; ++j) { n *= 2; // Try bigger steps. - if (g[n] > z) // z <= g[n] + if (w0s[n] > z) // z <= w0s[n] { goto overshot; // } } // for // get here if - std::cout << "Too big " << z << std::endl; + //std::cout << "Too big " << z << std::endl; //// else z too large :-( //const char* function = "boost::math::lambert_w0()"; return policies::raise_domain_error("boost::math::lambert_w0()", @@ -1179,7 +1154,7 @@ if (diff == 0) // Exact result - common. { break; } - if (g[n - nh] > z) + if (w0s[n - nh] > z) { n -= nh; } @@ -1187,8 +1162,8 @@ if (diff == 0) // Exact result - common. } bisect: - --n; // g[n] is nearest below, so n is integer part of W, - // and g[n+1] is next integral part of W. + --n; // w0s[n] is nearest below, so n is integer part of W, + // and w0s[n+1] is next integral part of W. // These are used as initial values for bisection. @@ -1207,15 +1182,15 @@ if (diff == 0) // Exact result - common. #endif // BOOST_MATH_INSTRUMENT_LAMBERT_W0 #ifdef BOOST_MATH_INSTRUMENT_LAMBERT_W0_LOOKUP - std::cout << "Result lookup W(" << z << ") bisection between g[" << n - 1 << "] = " << g[n - 1] << " and g[" << n << "] = " << g[n] - << ", bisect mean = " << (g[n - 1] + g[n]) / 2 << std::endl; + std::cout << "Result lookup W(" << z << ") bisection between w0s[" << n - 1 << "] = " << w0s[n - 1] << " and w0s[" << n << "] = " << w0s[n] + << ", bisect mean = " << (w0s[n - 1] + w0s[n]) / 2 << std::endl; #endif // BOOST_MATH_INSTRUMENT_LAMBERT_W0_LOOKUP if (d2 <= 7) - { // Only 7 binary digits precision required to hold integer size of g[65], + { // Only 7 binary digits precision required to hold integer size of w0s[65], // so just return the mean of two nearby integral values. // This is a *very* approximate estimate, and perhaps not very useful? - return static_cast((g[n - 1] + g[n]) / 2); + return static_cast((w0s[n - 1] + w0s[n]) / 2); } // Compute the number of bisections (jmax) that ensure that result is close enough that // a single 5th order Schroeder update is sufficient to achieve near double (53-bit) accuracy. @@ -1236,15 +1211,18 @@ if (diff == 0) // Exact result - common. { // Small z, so need just 1 extra bisection. jmax = 9; } - lookup_t dz = static_cast(z); // Only use double precision for bisection. - + lookup_t dz = static_cast(z); + // Only use lookup_t precision, default double, for bisection. + // Use Halley refinement for higher precisions. // Perform the jmax fractional bisections for necessary precision. - lookup_t y = dz * e[n + 1]; // - lookup_t w = static_cast(n); // Integral Lambert W estimate. + // Avoid using exponential function for speed. + lookup_t w = static_cast(n); // Equation 25, Integral Lambert W. + // lookup_t y = dz * e[n + 1]; // Integral +1 Lambert W. + lookup_t y = dz * static_cast(wm1s[n + 1]); // Equation 26, Integral +1 Lambert W. for (int j = 0; j < jmax; ++j) - { - const lookup_t wj = w + b[j]; // Add 1/2, 1/4, 1/8 ... - const lookup_t yj = y * a[j]; // Multiply by sqrt(1/e), ... + { // Equation 27. + const lookup_t wj = w + static_cast(halves[j]); // Add 1/2, 1/4, 1/8 ... + const lookup_t yj = y * static_cast(sqrts[j]); // Multiply by sqrt(1/e), ... if (wj < yj) { w = wj; @@ -1257,11 +1235,11 @@ if (diff == 0) // Exact result - common. // so just return the nearest bisection. // (Might make this test depend on size of T?) #ifdef BOOST_MATH_INSTRUMENT_LAMBERT_W0_BISECTION - std::cout << "Bisection estimate = " << w << std::endl; + std::cout << "Bisection estimate = " << w << " " << y << std::endl; #endif // BOOST_MATH_INSTRUMENT_LAMBERT_W0_BISECTION return static_cast(w); // Bisection only. } - else // More than 10 digits2 wanted so continue with Fukushima's Schroeder refinement. + else // More than 10 digits2 wanted, so continue with Fukushima's Schroeder refinement. { T result = static_cast(schroeder_update(w, y)); if (d2 <= (std::numeric_limits::digits - 3)) @@ -1291,78 +1269,198 @@ if (diff == 0) // Exact result - common. } // 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? // TODO check changes to W0 branch have also been made here. template - T lambert_wm1_imp(const T z, const Policy& /* pol */) + T lambert_wm1_imp(const T z, const Policy& pol) { // Catch providing an integer value as parameter x to lambert_w, for example, lambert_w(1). // Need to ensure it is a floating-point type (of the desired type, float 1.F, double 1., or long double 1.L), // or static_casted integer, for example: static_cast(1) or static_cast(1). // Want to allow fixed_point types too, so do not just test for floating-point. // Integral types should be promoted to double by user Lambert w functions. - // If integral type provided to user function lambert_w0 or _wm1, + // If integral type provided to user function lambert_w0 or lambert_wm1, // then should already have been promoted to double. - - BOOST_STATIC_ASSERT_MSG(!std::is_integral::value, "Must be floating-point or fixed type (not integer type), for example: W(1.), not W(1)!"); + BOOST_STATIC_ASSERT_MSG(!std::is_integral::value, "Must be floating-point or fixed type (not integer type), for example: lambert_wm1(1.), not lambert_wm1(1)!"); BOOST_MATH_STD_USING // Aid argument dependent lookup (ADL) of abs. - // else z is too large for the w-1 branch. - // TODO should this be the singularity value? And others are wrong too. - if (z >= 0) + using boost::math::tools::max_value; + + const char* function = "boost::math::lambert_wm1()"; // Used for error messages. + + if (z == static_cast(0)) + { // z is exactly zero.return -std::numeric_limits::infinity(); + return -std::numeric_limits::infinity(); + } + if (z > static_cast(0)) + { // + return policies::raise_domain_error(function, + "Argument z = %1% is out of range (z > 0) for Lambert W-1 branch! (Try Lambert W0 branch?)", + z, pol); + } + if(z > -(std::numeric_limits::min)()) + { // z is denormalized, so cannot be computed. + // -std::numeric_limits::min() is smallest for type T, + // for example, for double: lambert_wm1(-2.2250738585072014e-308) = -714.96865723796634 + return policies::raise_domain_error(function, + "Argument z = %1% is too small (z < -std::numeric_limits::min so denormalized) for Lambert W-1 branch!", + z, pol); + } + if (z == -boost::math::constants::exp_minus_one()) // == singularity/branch point z = -exp(-1) = -3.6787944. + { // At singularity, so return exactly -1. + return -static_cast(1); + } + // z is too negative for the W-1 (or W0) branch. + if (z < -boost::math::constants::exp_minus_one()) // > singularity/branch point z = -exp(-1) = -3.6787944. { - std::cerr << "lambert_wm1 argument out of range, z = " << z << std::endl; - return std::numeric_limits::quiet_NaN(); + return policies::raise_domain_error(function, + "Argument z = %1% is out of range (z < -exp(-1) = -3.6787944... <= 0) for Lambert W-1 (or W0) branch!", + z, pol); } if (z < -0.35) - { // Close to singularity/branch point. + { // Close to singularity/branch point but on w-1 branch. TODO. const T p2 = 2 * (boost::math::constants::e() * z + 1); if (p2 > 0) { return lambert_w_singularity_series(-sqrt(p2)); - // TODO Halley options here. + // TODO Halley refinement options here. } if (p2 == 0) - { + { // At the singularity at branch point. return -1; } - std::cerr << "(lambert_wm1) Argument out of range, z = " << z << std::endl; - return std::numeric_limits::quiet_NaN(); - } + return policies::raise_domain_error(function, + "Argument z = %1% is out of range!", + z, pol); + } // if (z < -0.35) + + // Create static arrays of Lambert Wm1 function + // TODO these are NOT the same size or values as the w0 constant array values????? + + // double LambertW0(const double z) + // static double e[66]; + // static double g[65]; + // static double a[12]; + // static double b[12]; // Create static arrays of Lambert W function for Wm1 // with lambert_w values of integers. // Used for bisection of W-1 branch. - static const size_t noof_sqrts = 12; + static const size_t noof_sqrts = 12; // and halves + // F[k] is 0 <= k <= 64 so size is 65 (but is 64 below????). + // G[k] is 1 <= k <= 64 so size is 64 + // So e and F are not the same! - static T e[64]; - static T g[64]; - static T a[noof_sqrts]; - static T b[noof_sqrts]; + static const size_t noof_wm1s = 64; + + static lookup_t e[noof_wm1s]; // e[0] to e[63] + static lookup_t g[noof_wm1s]; // g[0] == Gk when k = 1, to g[63] == Gk=64 + static lookup_t a[noof_sqrts]; + static lookup_t b[noof_sqrts]; if (!e[0]) { - const T e1 = 1 / boost::math::constants::e(); - T ej = e1; - e[0] = boost::math::constants::e(); + const lookup_t e1 = 1 / boost::math::constants::e(); + lookup_t ej = e1; // 1/e= 0.36787945 + e[0] = boost::math::constants::e(); g[0] = -e1; for (int j = 0, jj = 1; jj < 64; ++jj) - { + { // jj is j+1 ej *= e1; - e[jj] = e[j] * boost::math::constants::e(); + e[jj] = e[j] * boost::math::constants::e(); g[jj] = -(jj + 1) * ej; j = jj; } - a[0] = sqrt(boost::math::constants::e()); - b[0] = static_cast(0.5); + a[0] = sqrt(boost::math::constants::e()); + b[0] = static_cast(0.5); for (int j = 0, jj = 1; jj < noof_sqrts; ++jj) { - a[jj] = sqrt(a[j]); - b[jj] = b[j] * static_cast(0.5); + a[jj] = sqrt(a[j]); + // b[jj] = b[j] * static_cast(0.5); + b[jj] = b[j] / 2; // More efficient for multiprecision? j = jj; - } - } // static arrays filled. + } // for j + +//std::cout << " W-1 e[0] = " << e[0] << ", e[1] = " << e[1] << "... e[63] = " << e[63] << std::endl; +// W-1 e[0] = 2.71828175, e[1] = 7.38905573... e[63] = 6.23513822e+27 +// W-1 e[0] = 2.7182818284590451, e[1] = 7.3890560989306495... e[63] = 6.235149080811597e+27, e[64] = 58265.164084420219 +//std::cout << " W-1 g[0] = " << g[0] << ", g[1] = " << g[1] << "... g[63] = " << g[63] << std::endl; +// W-1 g[0] = -0.36787945 k=1, -1e^-1 = -0.3678794 , +// g[1] = -0.270670593 ... g[63] = -1.0264407e-26 +// W-1 g[0] = -0.36787944117144233, g[1] = -0.2706705664732254... g[63] = -1.0264389699511303e-26 + // TODO temporary output. + //std::cout << "lambert_wm1 version of arrays - NOT same as w0 version." << std::endl; + // print_collection(e, "e", " "); // e[0] = 2.71828175, e[1] = 1, e[2] = 0.36787945, 0.135335296 0.0497870743 0.0183156412 0.00673794793 ... e[63] = 2.293783159, e[64]= 6.235149080e+27 + // print_collection(g, "g", " "); // g[0] = -0.3678794, g[1] = -0.2706705... , g[63] = -1.026438e-26 + //print_collection(a, "a", " "); + //print_collection(b, "b", " "); + } // if (!e[0]) +// static arrays filled. + + // TODO move up when using precomputed array. + //if (abs(z) < abs(g[63])) + +if (z > g[63]) + { // z > -1.0264389699511303e-26 (but != 0 and >= std::numeric_limits::min() and so NOT denormalized). + std::streamsize saved_precision = std::cout.precision(std::numeric_limits::max_digits10); + + // std::cout << "abs(z) " << abs(z) << ", abs(g[63] = " << abs(g[63]) << std::endl; + + // std::cout << "-std::numeric_limits::min() = " << -(std::numeric_limits::min)() << std::endl; + // std::cout << "-std::numeric_limits::min() = " << -(std::numeric_limits::min)() << std::endl; + // -std::numeric_limits::min() = -1.1754943508222875e-38 + // -std::numeric_limits::min() = -2.2250738585072014e-308 + + // N[productlog(-1, -1.1754943508222875 * 10^-38 ), 50] = -91.856775324595479509567756730093823993834155027858 + + // N[productlog(-1, -2.2250738585072014e-308 * 10^-308 ), 50] = -1424.8544521230553853558132180518404363617968042942 + + T guess ; // bisect lowest possible Gk[=64] (for lookup_t type?) + + //int exponent = ilogb(z); + //std::cout << exponent << std::endl; + //guess = static_cast(exponent); + + // 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. + // François Chapeau-Blondeau and Abdelilah Monir + // Numerical Evaluation of the Lambert W Function + // IEEE Transactions On Signal Processing, VOL. 50, NO. 9, Sep 2002 + // https://pdfs.semanticscholar.org/7a5a/76a9369586dd0dd34dda156d8f2779d1fd59.pdf + // Estimate Lambert W using ln(-z) ... +// This is roughly the power of ten * ln(10) ~= 2.3. n ~= 10^n +// and improve by adding a second term -ln(ln(-z)) + T lz = log(-z); + T llz = log(-lz); + guess = lz - llz + (llz/lz); // Chapeau-Blondeau equation 20, page 2162. +#ifdef BOOST_MATH_INSTRUMENT_LAMBERT_WM1_TINY + std::cout << "z = " << z << ", guess = " << guess << ", ln(-z) = " << lz << ", ln(-ln(-z) = " << llz << ", llz/lz = " << (llz / lz) << std::endl; + // z = -1.0000000000000001e-30, guess = -73.312782616731482, ln(-z) = -69.077552789821368, ln(-ln(-z) = 4.2352298269101114, llz/lz = -0.061311231447304194 + // z = -9.9999999999999999e-91, guess = -212.56650048504233, ln(-z) = -207.23265836946410, ln(-ln(-z) = 5.3338421155782205, llz/lz = -0.025738424423764311 + // >z = -2.2250738585072014e-308, guess = -714.95942238244606, ln(-z) = -708.39641853226408, ln(-ln(-z) = 6.5630038501819854, llz/lz = -0.0092645920821846622 + int d10 = policies::digits_base10(); // policy template parameter digits10 + int d2 = policies::digits(); // digits base 2 from policy. + std::cout << "digits10 = " << d10 << ", digits2 = " << d2 // For example: digits10 = 1, digits2 = 5 + << std::endl; +#endif // BOOST_MATH_INSTRUMENT_LAMBERT_WM1_TINY + if (policies::digits() < 12) + { // For the worst case near w = 64, the error in the 'guess' is ~0.008, ratio ~ 0.0001 or 1 in 10,000 digits 10 ~= 4, or digits2 ~= 12. + return guess; + } + T result = halley_update(guess, z, pol); + std::cout.precision(saved_precision); + return result; + + //// G[k=64] == g[63] == -1.02643897e-26 + //return policies::raise_domain_error(function, + // "Argument z = %1% is too small (< -1.02643897e-26) !", + // z, pol); + } + + + + // Use a lookup table to find the nearest integer part of Lambert W as starting point for Bisection. // Bracketing sequence n = (2, 4, 8, 16, 32, 64) for W-1 branch. // Since z is probably quite small, start with lowest n (=2). int n = 2; @@ -1378,25 +1476,33 @@ if (diff == 0) // Exact result - common. goto overshot; } } - //else - // TODO use policy here. - std::cerr << "(lambert_wm1) Argument too small, z = " << z << std::endl; - return std::numeric_limits::quiet_NaN(); + // else z < g[63] == -1.0264389699511303e-26, so Lambert W integer part > 64. + + // Might assign the lower bracket to -1.0264389699511303e-26 and upper to (std::numeric_limits::min)() +// and jump to schroeder to halley? But that would only allow use with lookup_t precision. + // So try using Halley direct? + return policies::raise_domain_error(function, + "Argument z = %1% is too small (< -1.026439e-26) !", + z, pol); overshot: { - int nh = n / 2; + int nh = n / 2; for (int j = 1; j <= 5; ++j) { - nh /= 2; + nh /= 2; // halve step size. if (nh <= 0) - break; + { + break; // goto bisect; + } if (g[n - nh - 1] > z) + { n -= nh; + } } } bisect: - --n; // g[n] now holds lambert W of floor integer n. + --n; // g[n] now holds lambert W of floor integer n and g[n+1] the ceil part. // Check precision specified by policy. int d2 = policies::digits(); @@ -1411,9 +1517,10 @@ if (diff == 0) // Exact result - common. return static_cast(n); } - // jmax is the number of bisections such that a single application of the fifth-order update formula - // after the bisections is enough to evaluate W-1 with 53-bit accuracy. - int jmax = 11; // + // jmax is the number of bisections computed from n, + // such that a single application of the fifth-order Schroeder update formula + // after the bisections is enough to evaluate Lambert W-1 with 53-bit accuracy. + int jmax = 11; // Assume maximum number of bisections will be needed (most common case). if (n >= 8) { jmax = 8; @@ -1426,28 +1533,34 @@ if (diff == 0) // Exact result - common. { jmax = 10; } - T w = - static_cast(n); - T y = z * e[n - 1]; + // Bracketing, Fukushima section 2.3, page 82: + // (Avoid using exponential function for speed). + // Only use lookup_t precision, default double, for bisection (for speed). + lookup_t w = -static_cast(n); // Equation 25, + lookup_t y = static_cast(z * e[n - 1]); // Equation 26, for (int j = 0; j < jmax; ++j) - { - const T wj = w - b[j]; - const T yj = y * a[j]; + { // Equation 27. + lookup_t wj = w - b[j]; // Subtract 1/2, 1/4, 1/8 ... + lookup_t yj = y * a[j]; // Multiply by sqrt(1/e), ... if (wj < yj) { w = wj; y = yj; } - } - T result = schroeder_update(w, y); // Schroeder 5th order method refinement. - result = halley_update(result, z); + } // for j + + T result = static_cast(schroeder_update(w, y)); // Schroeder 5th order method refinement. + // Only use Halley refinement (using exp) for higher precisions. + result = halley_update(result, z, pol); return result; } // template T lambert_wm1_imp(const T z) } // namespace detail //////////////////////////////////////////////////////////// - // User Lambert W functions. + // User Lambert W0 and Lambert W-1 functions. - //! W0 branch, -1/e < z < max(z) + //! W0 branch, -1/e <= z < max(z) + //! W-1 branch -1/e >= z > min(z) //! Lambert W0 using User-defined policy. template diff --git a/test/lambert_w_high_reference_values.cpp b/test/lambert_w_high_reference_values.cpp index 5085a2dac..3c908cc6e 100644 --- a/test/lambert_w_high_reference_values.cpp +++ b/test/lambert_w_high_reference_values.cpp @@ -7,18 +7,18 @@ // Write a C++ file \lambert_w_mp_hi_values.ipp // containing arrays of z arguments and 100 decimal digit precision lambert_w0(z) reference values. -// These can be used in tests of precision of less-precise types like +// These can be used in tests of precision of less-precise types like // built-in float, double, long double and quad and cpp_dec_float_50. // These cover the range from 0.5 to (std::numeric_limits<>::max)(); -// The Fukushima algorithm changes from a series function for all z > 0.5. +// The Fukushima algorithm changes from a series function for all z > 0.5. // Multiprecision types: //#include #include // boost::multiprecision::cpp_dec_float_100 using boost::multiprecision::cpp_dec_float_100; -#include // +#include // using boost::math::lambert_w0; using boost::math::lambert_wm1; @@ -38,9 +38,9 @@ using std::ofstream; const long double eps = std::numeric_limits::epsilon(); // Creates if no file exists, & uses default overwrite/ ios::replace. -const char filename[] = "lambert_w_high_reference_values.ipp"; // +const char filename[] = "lambert_w_high_reference_values.ipp"; // std::ofstream fout(filename, std::ios::out); - + typedef cpp_dec_float_100 RealType; // 100 decimal digits for the value fed to macro BOOST_MATH_TEST_VALUE. // Could also use cpp_dec_float_50 or cpp_bin_float types. @@ -63,10 +63,10 @@ int main() try { int output_precision = std::numeric_limits::digits10; - // cpp_dec_float_100 is ample precision and + // cpp_dec_float_100 is ample precision and // has several extra bits internally so max_digits10 are not needed. fout.precision(output_precision); - fout << std::showpoint << std::endl; // Don't show trailing zeros. + fout << std::showpoint << std::endl; // Do show trailing zeros. // Intro for RealType values. std::cout << "Lambert W test values written to file " << filename << std::endl; @@ -130,7 +130,7 @@ int main() { fout << " zs[" << index << "] = BOOST_MATH_TEST_VALUE(RealType, " - << z // Since start with converting a float may get lots of usefully random digits. + << z // Since start with converting a float may get lots of usefully random digits. << ");" << std::endl; diff --git a/test/test_lambert_w.cpp b/test/test_lambert_w.cpp index a2612b3d2..8ca0deb5d 100644 --- a/test/test_lambert_w.cpp +++ b/test/test_lambert_w.cpp @@ -8,7 +8,7 @@ // test_lambertw.cpp //! \brief Basic sanity tests for Lambert W function plog -//! or lambert_w using algorithm informed by Thomas Luu, Veberic and Fukushima. +//! or lambert_w using algorithm informed by Thomas Luu, Veberic and Tosio Fukushima. #include // for real_concept #define BOOST_TEST_MAIN @@ -84,6 +84,14 @@ std::string show_versions() //PRINT_MACRO(LONG_MAX); #endif // __GNUC__ +#ifdef __MINGW64__ +std::cout << "MINGW64 " << __MINGW64_MAJOR_VERSION__ << __MINGW64_MINOR_VERSION << std::endl; +#endif // __MINGW64__ + +#ifdef __MINGW32__ +std::cout << 2MINGW64 " << __MINGW32_MAJOR_VERSION__ << __MINGW32_MINOR_VERSION << std::endl; +#endif // __MINGW32__ + message << "\n STL " << BOOST_STDLIB; message << "\n Boost version " << BOOST_VERSION / 100000 << "." << BOOST_VERSION / 100 % 1000 << "." << BOOST_VERSION % 100; @@ -121,7 +129,7 @@ void test_spots(RealType) #ifndef BOOST_NO_EXCEPTIONS BOOST_CHECK_THROW(boost::math::lambert_w0(-1.), std::domain_error); BOOST_CHECK_THROW(lambert_w0(std::numeric_limits::quiet_NaN()), std::domain_error); // Would be NaN. - BOOST_CHECK_THROW(lambert_w0(std::numeric_limits::infinity()), std::domain_error); // Would be infinite. + // BOOST_CHECK_THROW(lambert_w0(std::numeric_limits::infinity()), std::domain_error); // If infinity should throw. BOOST_CHECK_THROW(lambert_w0(-static_cast(0.4)), std::domain_error); // Would be complex. #else @@ -165,9 +173,14 @@ void test_spots(RealType) // Tests with some spot values computed using // https://www.wolframalpha.com/input - // For example: N[lambert_w0[1], 50] outputs: + // For example: N[lambert_w[1], 50] outputs: // 0.56714329040978387299996866221035554975381578718651 + BOOST_CHECK_CLOSE_FRACTION( // Check -exp(-1) ~= -0.367879450 == -1 + lambert_w0(BOOST_MATH_TEST_VALUE(RealType, -0.36787944117144232159552377016146086744581113103176783450783680169746149574489980335714727434591964374662732527)), + BOOST_MATH_TEST_VALUE(RealType, -1.), + tolerance); + BOOST_CHECK_CLOSE_FRACTION(lambert_w0(BOOST_MATH_TEST_VALUE(RealType, 0.1)), BOOST_MATH_TEST_VALUE(RealType, 0.091276527160862264299895721423179568653119224051472), // Output from https://www.wolframalpha.com/input/?i=lambert_w0(0.2) @@ -214,6 +227,145 @@ void test_spots(RealType) // Output from https://www.wolframalpha.com/input/?i=lambert_w0(100) tolerance); + if (std::numeric_limits::has_infinity) + { + // BOOST_CHECK_THROW(lambert_w0(std::numeric_limits::infinity()), std::domain_error); // If should throw exception. + //BOOST_CHECK_EQUAL(lambert_w0(std::numeric_limits::infinity()), +std::numeric_limits::infinity()); // message is: + // Error in "test_types": class boost::exception_detail::clone_impl > : + // Error in function boost::math::lambert_w0() : Argument z is infinite! + BOOST_CHECK_EQUAL(lambert_w0(std::numeric_limits::infinity()), +std::numeric_limits::infinity()); // Infinity allowed. + } + if (std::numeric_limits::has_quiet_NaN) + { // Argument Z == NaN is always an throwable error. + // BOOST_CHECK_EQUAL(lambert_w0(std::numeric_limits::quiet_NaN()), +std::numeric_limits::infinity()); // message is: + // Error in function boost::math::lambert_w0(): Argument z is NaN! + BOOST_CHECK_THROW(lambert_w0(std::numeric_limits::quiet_NaN()), std::domain_error); + BOOST_CHECK_THROW(lambert_wm1(std::numeric_limits::quiet_NaN()), std::domain_error); + } + + // Some tests of Lambert W-1 branch. + + BOOST_CHECK_CLOSE_FRACTION( // Check -exp(-1) ~= -0.367879450 == -1 at the singularity branch point. + lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -0.36787944117144232159552377016146086744581113103176783450783680169746149574489980335714727434591964374662732527)), + BOOST_MATH_TEST_VALUE(RealType, -1.), + tolerance); + + // Near singularity and using series approximation. + //N[productlog(-1, -0.36), 50] = -1.2227701339785059531429380734238623131735264411311 + BOOST_CHECK_CLOSE_FRACTION(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -0.36)), + BOOST_MATH_TEST_VALUE(RealType, -1.2227701339785059531429380734238623131735264411311), + tolerance); + + // Near singularity and NOT using series approximation (switch at -0.35) + // N[productlog(-1, -0.34), 50] + BOOST_CHECK_CLOSE_FRACTION(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -0.34)), + BOOST_MATH_TEST_VALUE(RealType, -1.4512014851325470735077533710339268100722032730024), + tolerance); + + // Decreasing z until near zero. + //N[productlog(-1, -0.3), 50] = -1.7813370234216276119741702815127452608215583564545 + BOOST_CHECK_CLOSE_FRACTION(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -0.3)), + BOOST_MATH_TEST_VALUE(RealType, -1.7813370234216276119741702815127452608215583564545), + tolerance); + + //N[productlog(-1, -0.2), 50] = -2.5426413577735264242938061566618482901614749075294 + BOOST_CHECK_CLOSE_FRACTION(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -0.2)), + BOOST_MATH_TEST_VALUE(RealType, -2.5426413577735264242938061566618482901614749075294), + tolerance); + + BOOST_CHECK_CLOSE_FRACTION(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -0.1)), + BOOST_MATH_TEST_VALUE(RealType, -3.577152063957297218409391963511994880401796257793), + tolerance); + + //N[productlog(-1, -0.01), 50] = -6.4727751243940046947410578927244880371043455902257 + BOOST_CHECK_CLOSE_FRACTION(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -0.01)), + BOOST_MATH_TEST_VALUE(RealType, -6.4727751243940046947410578927244880371043455902257), + tolerance); + + // N[productlog(-1, -0.001), 50] = -9.1180064704027401212583371820468142742704349737639 + BOOST_CHECK_CLOSE_FRACTION(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -0.001)), + BOOST_MATH_TEST_VALUE(RealType, -9.1180064704027401212583371820468142742704349737639), + tolerance); + + // N[productlog(-1, -0.000001), 50] = -16.626508901372473387706432163984684996461726803805 + BOOST_CHECK_CLOSE_FRACTION(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -0.000001)), + BOOST_MATH_TEST_VALUE(RealType, -16.626508901372473387706432163984684996461726803805), + tolerance); + + BOOST_CHECK_CLOSE_FRACTION(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -1e-12)), + BOOST_MATH_TEST_VALUE(RealType, -31.067172842017230842039496250208586707880448763222), + tolerance); + + BOOST_CHECK_CLOSE_FRACTION(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -1e-25)), + BOOST_MATH_TEST_VALUE(RealType, -61.686695602074505366866968627049381352503620377944), + tolerance); + + // z nearly too small + BOOST_CHECK_CLOSE_FRACTION(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -2e-26)), + BOOST_MATH_TEST_VALUE(RealType, -63.322302839923597803393585145387854867226970485197), + tolerance* 2); + + // z very nearly too small. G(k=64) g[63] = -1.0264389699511303e-26 to using 1.027e-26 + BOOST_CHECK_CLOSE_FRACTION(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -1.027e-26)), + BOOST_MATH_TEST_VALUE(RealType, -63.999444896732265186957073549916026532499356695343), + tolerance); + // So -64 is the most negative value that can be determined using lookup. + // N[productlog(-1, -1.0264389699511303 * 10^-26 ), 50] -63.999999999999997947255011093606206983577811736472 == -64 + // G[k=64] = g[63] = -1.0264389699511303e-26 + + // z too small for G(k=64) g[63] = -1.0264389699511303e-26 to using 1.027e-26 + // N[productlog(-1, -10 ^ -26), 50] = -31.067172842017230842039496250208586707880448763222 + BOOST_CHECK_CLOSE_FRACTION(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -1e-26)), + BOOST_MATH_TEST_VALUE(RealType, -64.026509628385889681156090340691637712441162092868), + tolerance); + + + + // z == zero. + + if (std::numeric_limits::has_infinity) + { + BOOST_CHECK_EQUAL(lambert_wm1(0.L), -std::numeric_limits::infinity()); + + } + if (std::numeric_limits::has_quiet_NaN) + { + // BOOST_CHECK_EQUAL(lambert_w0(std::numeric_limits::quiet_NaN()), +std::numeric_limits::infinity()); // message is: + // Error in function boost::math::lambert_w0(): Argument z is NaN! + BOOST_CHECK_THROW(lambert_wm1(std::numeric_limits::quiet_NaN()), std::domain_error); + } + + + // N[productlog(-1, -10 ^ -26), 50] = -31.067172842017230842039496250208586707880448763222 + //BOOST_CHECK_CLOSE_FRACTION(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -1e-26)), + // BOOST_MATH_TEST_VALUE(RealType, -64.026509628385889681156090340691637712441162092868L), + // tolerance); + + // 1e-28 is too small + // N[productlog(-1, -10 ^ -28), 50] = -31.067172842017230842039496250208586707880448763222 + //BOOST_CHECK_CLOSE_FRACTION(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -1e-28)), + // BOOST_MATH_TEST_VALUE(RealType, -68.702163291525429160769761667024460023336801014578L), + // tolerance); + + + + // Check for overflow when using a double (including when using for approximate value for refinement for higher precision). + + // N[productlog(-1, -10 ^ -30), 50] = -73.373110313822976797067478758120874529181611813766 + //BOOST_CHECK_CLOSE_FRACTION(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -1e30)), + // BOOST_MATH_TEST_VALUE(RealType, -73.373110313822976797067478758120874529181611813766), + // tolerance); + //unknown location : fatal error : in "test_types" : + //class boost::exception_detail::clone_impl > + // : Error in function boost::math::lambert_wm1() : + // Argument z = -1.00000002e+30 out of range(z < -exp(-1) = -3.6787944) for Lambert W - 1 branch! + + BOOST_CHECK_THROW(lambert_wm1(-0.5), std::domain_error); + + + + + // This fails for fixed_point type used for other tests because out of range? //BOOST_CHECK_CLOSE_FRACTION(lambert_w0(BOOST_MATH_TEST_VALUE(RealType, 1.0e6)), //BOOST_MATH_TEST_VALUE(RealType, 11.383358086140052622000156781585004289033774706019), @@ -275,12 +427,17 @@ void test_spots(RealType) // Checks on input that should throw. - /* This should throw when implemented. - BOOST_CHECK_CLOSE_FRACTION(lambert_w0(-0.5), - BOOST_MATH_TEST_VALUE(RealType, ), - // Output from https://www.wolframalpha.com/input/?i=lambert_w0(-0.5) - tolerance); - */ + // This should throw when implemented. + + //BOOST_CHECK_CLOSE_FRACTION(lambert_w0(-0.5), + //BOOST_MATH_TEST_VALUE(RealType, 0.0), tolerance); + // unknown location : fatal error : in "test_types": + //class boost::exception_detail::clone_impl >: + // Error in function boost::math::lambert_w0(): + // Argument z = -0.5 out of range (-1/e <= z < (std::numeric_limits::max)()) for Lambert W0 branch (use W-1 branch?). + BOOST_CHECK_THROW(lambert_w0(-0.5), std::domain_error); + + } //template void test_spots(RealType) BOOST_AUTO_TEST_CASE(test_types) @@ -299,15 +456,12 @@ BOOST_AUTO_TEST_CASE(test_types) //{ // Avoid pointless re-testing if double and long double are identical (for example, MSVC). // test_spots(0.0L); // long double //} - // Built-in/fundamental Quad 128-bit type, if available. -#ifdef BOOST_HAS_FLOAT128 - test_spots(static_cast(0)); -#endif + // TODO - renable these tests when lambert_wm1 lookup table complete // Multiprecision types: - test_spots(static_cast(0)); - test_spots(static_cast(0)); - test_spots(static_cast(0)); + //test_spots(static_cast(0)); + //test_spots(static_cast(0)); + //test_spots(static_cast(0)); // Fixed-point types: diff --git a/test/test_lambert_w0_precision_high.cpp b/test/test_lambert_w0_precision_high.cpp index 31221461d..0d1d19193 100644 --- a/test/test_lambert_w0_precision_high.cpp +++ b/test/test_lambert_w0_precision_high.cpp @@ -15,7 +15,7 @@ #include // boost::exception #include // For exp_minus_one == 3.67879441171442321595523770161460867e-01. #include -//#include "test_value.hpp" // for create_test_value and macro BOOST_MATH_TEST_VALUE. +#include "test_value.hpp" // for create_test_value and macro BOOST_MATH_TEST_VALUE. #include // boost::multiprecision::cpp_dec_float_50 using boost::multiprecision::cpp_dec_float_50; // 50 decimal digits type. @@ -51,9 +51,7 @@ using boost::lexical_cast; // Include 'Known good' Lambert w values using N[productlog(-0.3), 100] // evaluated to full precision of RealType (up to 100 decimal digits). // Checked for a few z values using WolframAlpha command: N[productlog(-0.3), 100] -//#include "J:\Cpp\Misc\lambert_w_1000_test_values\lambert_w_mp_values.ipp" -#include "J:\Cpp\Misc\lambert_w_high_reference_values\lambert_w_mp_high_values.ipp" -// #include "lambert_w_mp_values.ipp" +#include "lambert_w_high_reference_values.ipp" // Optional functions to list z arguments and reference lambert w values. template @@ -229,7 +227,7 @@ int main() { std::cout << "Lambert W tests,\n" << noof_tests << - " single spot tests comparing with 100 decimal digits precision reference values." << std::endl; + " tests comparing with 100 decimal digits precision reference values." << std::endl; std::cout.precision(std::numeric_limits::max_digits10); std::cout << std::showpoint << std::endl; // Trailing zeros. diff --git a/test/test_lambert_w0_precision_low.cpp b/test/test_lambert_w0_precision_low.cpp index 9a8a39479..0ee840c42 100644 --- a/test/test_lambert_w0_precision_low.cpp +++ b/test/test_lambert_w0_precision_low.cpp @@ -15,7 +15,7 @@ #include // boost::exception #include // For exp_minus_one == 3.67879441171442321595523770161460867e-01. #include -//#include "test_value.hpp" // for create_test_value and macro BOOST_MATH_TEST_VALUE. +#include "test_value.hpp" // for create_test_value and macro BOOST_MATH_TEST_VALUE. // Built-in/fundamental Quad 128-bit type, if available. #ifdef BOOST_HAS_FLOAT128 @@ -60,8 +60,7 @@ using boost::lexical_cast; // Include 'Known good' Lambert w values using N[productlog(-0.3), 100] // evaluated to full precision of RealType (up to 100 decimal digits). // Checked for a few z values using WolframAlpha command: N[productlog(-0.3), 100] -#include "J:\Cpp\Misc\lambert_w_1000_test_values\lambert_w_mp_values.ipp" -// #include "lambert_w_mp_values.ipp" +#include "lambert_w_low_reference_values.ipp" // Optional functions to list z arguments and reference lambert w values. template