123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581158215831584158515861587158815891590159115921593159415951596159715981599160016011602160316041605160616071608160916101611161216131614161516161617161816191620162116221623162416251626162716281629163016311632163316341635163616371638163916401641164216431644164516461647164816491650165116521653165416551656165716581659166016611662166316641665166616671668166916701671167216731674167516761677167816791680168116821683168416851686168716881689169016911692169316941695169616971698169917001701170217031704170517061707170817091710171117121713171417151716171717181719172017211722172317241725172617271728172917301731173217331734173517361737173817391740174117421743174417451746174717481749175017511752175317541755175617571758175917601761176217631764176517661767176817691770177117721773177417751776177717781779178017811782178317841785178617871788178917901791179217931794179517961797179817991800180118021803180418051806180718081809181018111812181318141815181618171818181918201821182218231824182518261827182818291830183118321833183418351836183718381839184018411842184318441845184618471848184918501851185218531854185518561857185818591860186118621863186418651866186718681869187018711872187318741875187618771878187918801881188218831884188518861887188818891890189118921893189418951896189718981899190019011902190319041905190619071908190919101911191219131914191519161917191819191920192119221923192419251926192719281929193019311932193319341935193619371938193919401941194219431944194519461947194819491950195119521953195419551956195719581959196019611962196319641965196619671968196919701971197219731974197519761977197819791980198119821983198419851986198719881989199019911992199319941995199619971998199920002001200220032004200520062007200820092010201120122013201420152016201720182019202020212022202320242025202620272028202920302031203220332034203520362037203820392040204120422043204420452046204720482049205020512052205320542055205620572058205920602061206220632064206520662067206820692070207120722073207420752076207720782079208020812082208320842085208620872088208920902091209220932094209520962097209820992100210121022103210421052106210721082109211021112112211321142115211621172118211921202121212221232124212521262127212821292130213121322133213421352136213721382139214021412142214321442145214621472148214921502151215221532154215521562157215821592160216121622163216421652166216721682169217021712172217321742175217621772178217921802181218221832184218521862187218821892190219121922193219421952196 |
- #ifndef BOOST_MATH_SF_LAMBERT_W_HPP
- #define BOOST_MATH_SF_LAMBERT_W_HPP
- #ifdef _MSC_VER
- #pragma warning(disable : 4127)
- #endif
- #include <boost/math/tools/config.hpp>
- #include <boost/math/policies/error_handling.hpp>
- #include <boost/math/policies/policy.hpp>
- #include <boost/math/tools/promotion.hpp>
- #include <boost/math/special_functions/fpclassify.hpp>
- #include <boost/math/special_functions/log1p.hpp>
- #include <boost/math/constants/constants.hpp>
- #include <boost/math/special_functions/next.hpp>
- #include <boost/math/special_functions/pow.hpp>
- #include <boost/math/tools/series.hpp>
- #include <boost/math/tools/rational.hpp>
- #include <boost/math/tools/precision.hpp>
- #include <boost/math/tools/big_constant.hpp>
- #include <boost/math/tools/cxx03_warn.hpp>
- #ifndef BOOST_MATH_STANDALONE
- #include <boost/lexical_cast.hpp>
- #endif
- #include <limits>
- #include <cmath>
- #include <limits>
- #include <exception>
- #include <type_traits>
- #include <cstdint>
- #include <iostream>
- #include <typeinfo>
- #include <boost/math/special_functions/next.hpp>
- using lookup_t = double;
- #include <boost/math/special_functions/detail/lambert_w_lookup_table.ipp>
- #if defined(__GNUC__) && defined(BOOST_MATH_USE_FLOAT128)
- #pragma GCC system_header
- #endif
- namespace boost {
- namespace math {
- namespace lambert_w_detail {
- template <typename T>
- inline T lambert_w_halley_step(T w_est, const T z)
- {
- BOOST_MATH_STD_USING
- T e = exp(w_est);
- w_est -= 2 * (w_est + 1) * (e * w_est - z) / (z * (w_est + 2) + e * (w_est * (w_est + 2) + 2));
- return w_est;
- }
- template <typename T>
- inline T lambert_w_halley_iterate(T w_est, const T z)
- {
- BOOST_MATH_STD_USING
- static const T max_diff = boost::math::tools::root_epsilon<T>() * fabs(w_est);
- T w_new = lambert_w_halley_step(w_est, z);
- T diff = fabs(w_est - w_new);
- while (diff > max_diff)
- {
- w_est = w_new;
- w_new = lambert_w_halley_step(w_est, z);
- diff = fabs(w_est - w_new);
- }
- return w_new;
- }
- template <typename T>
- inline T lambert_w_maybe_halley_iterate(T z, T w, std::false_type const&)
- {
- return lambert_w_halley_step(z, w);
- }
- template <typename T>
- inline T lambert_w_maybe_halley_iterate(T z, T w, std::true_type const&)
- {
- return lambert_w_halley_iterate(z, w);
- }
- template <typename T>
- inline double maybe_reduce_to_double(const T& z, const std::true_type&)
- {
- return static_cast<double>(z);
- }
- template <typename T>
- inline T maybe_reduce_to_double(const T& z, const std::false_type&)
- {
- return z;
- }
- template <typename T>
- inline double must_reduce_to_double(const T& z, const std::true_type&)
- {
- return static_cast<double>(z);
- }
- template <typename T>
- inline double must_reduce_to_double(const T& z, const std::false_type&)
- {
- #ifndef BOOST_MATH_STANDALONE
- #ifdef BOOST_MATH_USE_CHARCONV_FOR_CONVERSION
-
- if constexpr (std::is_arithmetic_v<T>)
- {
- return static_cast<double>(z);
- }
- else
- {
- return boost::lexical_cast<double>(z);
- }
- #else
-
- return boost::lexical_cast<double>(z);
-
- #endif
- #else
- static_assert(sizeof(T) == 0, "Unsupported in standalone mode: don't know how to cast your number type to a double.");
- return 0.0;
- #endif
- }
- template<typename T>
- inline T schroeder_update(const T w, const T y)
- {
-
-
-
-
-
-
-
- BOOST_MATH_STD_USING
- #ifdef BOOST_MATH_INSTRUMENT_LAMBERT_W_SCHROEDER
- std::streamsize saved_precision = std::cout.precision(std::numeric_limits<T>::max_digits10);
- using boost::math::float_distance;
- T fd = float_distance<T>(w, y);
- std::cout << "Schroder ";
- if (abs(fd) < 214748000.)
- {
- std::cout << " Distance = "<< static_cast<int>(fd);
- }
- else
- {
- std::cout << "Difference w - y = " << (w - y) << ".";
- }
- std::cout << std::endl;
- #endif
-
- const T f0 = w - y;
- const T f1 = 1 + y;
- const T f00 = f0 * f0;
- const T f11 = f1 * f1;
- const T f0y = f0 * y;
- const T result =
- w - 4 * f0 * (6 * f1 * (f11 + f0y) + f00 * y) /
- (f11 * (24 * f11 + 36 * f0y) +
- f00 * (6 * y * y + 8 * f1 * y + f0y));
- #ifdef BOOST_MATH_INSTRUMENT_LAMBERT_W_SCHROEDER
- std::cout << "Schroeder refined " << w << " " << y << ", difference " << w-y << ", change " << w - result << ", to result " << result << std::endl;
- std::cout.precision(saved_precision);
- #endif
- return result;
- }
-
-
-
-
-
-
-
-
- template<typename T>
- T lambert_w_singularity_series(const T p)
- {
- #ifdef BOOST_MATH_INSTRUMENT_LAMBERT_W_SINGULARITY_SERIES
- std::size_t saved_precision = std::cout.precision(3);
- std::cout << "Singularity_series Lambert_w p argument = " << p << std::endl;
- std::cout
-
-
-
- << std::endl;
- std::cout.precision(saved_precision);
- #endif
- static const T q[] =
- {
- -static_cast<T>(1),
- +T(1),
- -T(1) / 3,
- +T(11) / 72,
- -T(43) / 540,
- +T(769) / 17280,
- -T(221) / 8505,
-
-
- +T(680863uLL) / 43545600uLL,
-
- -T(1963uLL) / 204120uLL,
-
- +T(226287557uLL) / 37623398400uLL,
- -T(5776369uLL) / 1515591000uLL,
-
- +T(169709463197uLL) / 69528040243200uLL,
-
- -T(1118511313uLL) / 709296588000uLL,
- +T(667874164916771uLL) / 650782456676352000uLL,
-
- -T(500525573uLL) / 744761417400uLL,
-
-
-
- BOOST_MATH_BIG_CONSTANT(T, 64, +0.000442473061814620910),
-
- BOOST_MATH_BIG_CONSTANT(T, 64, -0.000292677224729627445),
-
- BOOST_MATH_BIG_CONSTANT(T, 64, 0.000194387276054539318),
-
- BOOST_MATH_BIG_CONSTANT(T, 64, -0.000129574266852748819),
-
- BOOST_MATH_BIG_CONSTANT(T, 64, +0.0000866503580520812717),
-
-
-
-
- BOOST_MATH_BIG_CONSTANT(T, 113, -0.000058113607504413816772205464778828177256611844221913)
-
-
-
-
-
-
-
-
-
- };
-
-
-
-
-
-
-
- BOOST_MATH_STD_USING
- const T absp = abs(p);
- #ifdef BOOST_MATH_INSTRUMENT_LAMBERT_W_TERMS
- {
- int terms = 20;
- if (absp < 0.01159)
- {
- terms = 6;
- }
- else if (absp < 0.0766)
- {
- terms = 10;
- }
- std::streamsize saved_precision = std::cout.precision(3);
- std::cout << "abs(p) = " << absp << ", terms = " << terms << std::endl;
- std::cout.precision(saved_precision);
- }
- #endif
- if (absp < T(0.01159))
- {
- return
- -1 +
- p * (1 +
- p * (q[2] +
- p * (q[3] +
- p * (q[4] +
- p * (q[5] +
- p * q[6]
- )))));
- }
- else if (absp < T(0.0766))
- {
- return
- -1 +
- p * (1 +
- p * (q[2] +
- p * (q[3] +
- p * (q[4] +
- p * (q[5] +
- p * (q[6] +
- p * (q[7] +
- p * (q[8] +
- p * (q[9] +
- p * q[10]
- )))))))));
- }
-
- return
- -1 +
- p * (1 +
- p * (q[2] +
- p * (q[3] +
- p * (q[4] +
- p * (q[5] +
- p * (q[6] +
- p * (q[7] +
- p * (q[8] +
- p * (q[9] +
- p * (q[10] +
- p * (q[11] +
- p * (q[12] +
- p * (q[13] +
- p * (q[14] +
- p * (q[15] +
- p * (q[16] +
- p * (q[17] +
- p * (q[18] +
- p * (q[19] +
- p * q[20]
- )))))))))))))))))));
-
-
-
-
-
-
-
-
- }
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- template <typename T, typename Policy>
- T lambert_w0_small_z(T x, const Policy&, std::integral_constant<int, 0> const&);
- template <typename T, typename Policy>
- T lambert_w0_small_z(T x, const Policy&, std::integral_constant<int, 1> const&);
- template <typename T, typename Policy>
- T lambert_w0_small_z(T x, const Policy&, std::integral_constant<int, 2> const&);
- template <typename T, typename Policy>
- T lambert_w0_small_z(T x, const Policy&, std::integral_constant<int, 3> const&);
- template <typename T, typename Policy>
- T lambert_w0_small_z(T x, const Policy&, std::integral_constant<int, 4> const&);
- template <typename T, typename Policy>
- T lambert_w0_small_z(T x, const Policy&, std::integral_constant<int, 5> const&);
-
- template <typename T, typename Policy>
- T lambert_w0_small_z(T x, const Policy& pol)
- {
- using tag_type = std::integral_constant<int,
- std::numeric_limits<T>::is_specialized == 0 ? 5 :
- #ifndef BOOST_NO_CXX11_NUMERIC_LIMITS
- std::numeric_limits<T>::max_digits10 <= 9 ? 0 :
- std::numeric_limits<T>::max_digits10 <= 17 ? 1 :
- std::numeric_limits<T>::max_digits10 <= 22 ? 2 :
- std::numeric_limits<T>::max_digits10 < 37 ? 4
- #else
- std::numeric_limits<T>::radix != 2 ? 5 :
- std::numeric_limits<T>::digits <= 24 ? 0 :
- std::numeric_limits<T>::digits <= 53 ? 1 :
- std::numeric_limits<T>::digits <= 64 ? 2 :
- std::numeric_limits<T>::digits <= 113 ? 4
- #endif
- : 5>;
-
- return lambert_w0_small_z(x, pol, tag_type());
- }
-
-
-
-
-
- template <typename T, typename Policy>
- T lambert_w0_small_z(T z, const Policy&, std::integral_constant<int, 0> const&)
- {
- #ifdef BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES
- std::streamsize prec = std::cout.precision(std::numeric_limits<T>::max_digits10);
- std::cout << "\ntag_type 0 float lambert_w0_small_z called with z = " << z << " using " << 9 << " terms of precision "
- << std::numeric_limits<float>::max_digits10 << " decimal digits. " << std::endl;
- #endif
- T result =
- z * (1 -
- z * (1 -
- z * (static_cast<float>(3uLL) / 2uLL -
- z * (2.6666666666666666667F -
- z * (5.2083333333333333333F -
- z * (10.8F -
- z * (23.343055555555555556F -
- z * (52.012698412698412698F -
- z * 118.62522321428571429F))))))));
- #ifdef BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES
- std::cout << "return w = " << result << std::endl;
- std::cout.precision(prec);
- #endif
- return result;
- }
-
-
-
-
- template <typename T, typename Policy>
- T lambert_w0_small_z(const T z, const Policy&, std::integral_constant<int, 1> const&)
- {
- #ifdef BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES
- std::streamsize prec = std::cout.precision(std::numeric_limits<T>::max_digits10);
- std::cout << "\ntag_type 1 double lambert_w0_small_z called with z = " << z << " using " << 17 << " terms of precision, "
- << std::numeric_limits<double>::max_digits10 << " decimal digits. " << std::endl;
- #endif
- T result =
- z * (1. -
- z * (1. -
- z * (1.5 -
- z * (2.6666666666666666667 -
- z * (5.2083333333333333333 -
- z * (10.8 -
- z * (23.343055555555555556 -
- z * (52.012698412698412698 -
- z * (118.62522321428571429 -
- z * (275.57319223985890653 -
- z * (649.78717234347442681 -
- z * (1551.1605194805194805 -
- z * (3741.4497029592385495 -
- z * (9104.5002411580189358 -
- z * (22324.308512706601434 -
- z * (55103.621972903835338 -
- z * 136808.86090394293563))))))))))))))));
- #ifdef BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES
- std::cout << "return w = " << result << std::endl;
- std::cout.precision(prec);
- #endif
- return result;
- }
-
-
-
-
-
- template <typename T, typename Policy>
- T lambert_w0_small_z(const T z, const Policy&, std::integral_constant<int, 2> const&)
- {
- #ifdef BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES
- std::streamsize precision = std::cout.precision(std::numeric_limits<T>::max_digits10);
- std::cout << "\ntag_type 2 long double (80-bit double extended) lambert_w0_small_z called with z = " << z << " using " << 21 << " terms of precision, "
- << std::numeric_limits<long double>::max_digits10 << " decimal digits. " << std::endl;
- #endif
- T result =
- z * (1.L -
- z * (1.L -
- z * (1.500000000000000000000000000000000L -
- z * (2.666666666666666666666666666666666L -
- z * (5.208333333333333333333333333333333L -
- z * (10.80000000000000000000000000000000L -
- z * (23.34305555555555555555555555555555L -
- z * (52.01269841269841269841269841269841L -
- z * (118.6252232142857142857142857142857L -
- z * (275.5731922398589065255731922398589L -
- z * (649.7871723434744268077601410934744L -
- z * (1551.160519480519480519480519480519L -
- z * (3741.449702959238549516327294105071L -
- z * (9104.500241158018935796713574491352L -
- z * (22324.308512706601434280005708577137L -
- z * (55103.621972903835337697771560205422L -
- z * (136808.86090394293563342215789305736L -
- z * (341422.05066583836331735491399356945L -
- z * (855992.9659966075514633630250633224L -
- z * (2.154990206091088289321708745358647e6L
- ))))))))))))))))))));
- #ifdef BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES
- std::cout << "return w = " << result << std::endl;
- std::cout.precision(precision);
- #endif
- return result;
- }
- template <typename T, typename Policy>
- T lambert_w0_small_z(const T z, const Policy&, std::integral_constant<int, 3> const&)
- {
- #ifdef BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES
- std::streamsize precision = std::cout.precision(std::numeric_limits<T>::max_digits10);
- std::cout << "\ntag_type 3 long double (128-bit) lambert_w0_small_z called with z = " << z << " using " << 17 << " terms of precision, "
- << std::numeric_limits<double>::max_digits10 << " decimal digits. " << std::endl;
- #endif
- T result =
- z * (1.L -
- z * (1.L -
- z * (1.5L -
- z * (2.6666666666666666666666666666666666L -
- z * (5.2052083333333333333333333333333333L -
- z * (10.800000000000000000000000000000000L -
- z * (23.343055555555555555555555555555555L -
- z * (52.0126984126984126984126984126984126L -
- z * (118.625223214285714285714285714285714L -
- z * (275.57319223985890652557319223985890L -
- z * (649.78717234347442680776014109347442680776014109347L -
- z * (1551.1605194805194805194805194805194805194805194805L -
- z * (3741.4497029592385495163272941050718828496606274384L -
- z * (9104.5002411580189357967135744913522691300469078247L -
- z * (22324.308512706601434280005708577137148565719994291L -
- z * (55103.621972903835337697771560205422639285073147507L -
- z * 136808.86090394293563342215789305736395683485630576L
- ))))))))))))))));
- #ifdef BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES
- std::cout << "return w = " << result << std::endl;
- std::cout.precision(precision);
- #endif
- return result;
- }
- #ifdef BOOST_HAS_FLOAT128
- template <typename T, typename Policy>
- T lambert_w0_small_z(const T z, const Policy&, std::integral_constant<int, 4> const&)
- {
- #ifdef BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES
- std::streamsize precision = std::cout.precision(std::numeric_limits<T>::max_digits10);
- std::cout << "\ntag_type 4 128-bit quad float128 lambert_w0_small_z called with z = " << z << " using " << 34 << " terms of precision, "
- << std::numeric_limits<float128>::max_digits10 << " max decimal digits." << std::endl;
- #endif
- T result =
- z * (1.Q -
- z * (1.Q -
- z * (1.500000000000000000000000000000000Q -
- z * (2.666666666666666666666666666666666Q -
- z * (5.208333333333333333333333333333333Q -
- z * (10.80000000000000000000000000000000Q -
- z * (23.34305555555555555555555555555555Q -
- z * (52.01269841269841269841269841269841Q -
- z * (118.6252232142857142857142857142857Q -
- z * (275.5731922398589065255731922398589Q -
- z * (649.7871723434744268077601410934744Q -
- z * (1551.160519480519480519480519480519Q -
- z * (3741.449702959238549516327294105071Q -
- z * (9104.500241158018935796713574491352Q -
- z * (22324.308512706601434280005708577137Q -
- z * (55103.621972903835337697771560205422Q -
- z * (136808.86090394293563342215789305736Q -
- z * (341422.05066583836331735491399356945Q -
- z * (855992.9659966075514633630250633224Q -
- z * (2.154990206091088289321708745358647e6Q -
- z * (5.445552922314462431642316420035073e6Q -
- z * (1.380733000216662949061923813184508e7Q -
- z * (3.511704498513923292853869855945334e7Q -
- z * (8.956800256102797693072819557780090e7Q -
- z * (2.290416846187949813964782641734774e8Q -
- z * (5.871035041171798492020292225245235e8Q -
- z * (1.508256053857792919641317138812957e9Q -
- z * (3.882630161293188940385873468413841e9Q -
- z * (1.001394313665482968013913601565723e10Q -
- z * (2.587356736265760638992878359024929e10Q -
- z * (6.696209709358073856946120522333454e10Q -
- z * (1.735711659599198077777078238043644e11Q -
- z * (4.505680465642353886756098108484670e11Q -
- z * (1.171223178256487391904047636564823e12Q
- ))))))))))))))))))))))))))))))))));
- #ifdef BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES
- std::cout << "return w = " << result << std::endl;
- std::cout.precision(precision);
- #endif
- return result;
- }
- #else
- template <typename T, typename Policy>
- inline T lambert_w0_small_z(const T z, const Policy& pol, std::integral_constant<int, 4> const&)
- {
- return lambert_w0_small_z(z, pol, std::integral_constant<int, 5>());
- }
- #endif
- template <typename T>
- struct lambert_w0_small_z_series_term
- {
- using result_type = T;
-
-
-
-
- lambert_w0_small_z_series_term(T _z, T _term, int _k)
- : k(_k), z(_z), term(_term) { }
- T operator()()
- {
- using std::pow;
- ++k;
- term *= -z / k;
-
- T result = term * pow(T(k), T(-1 + k));
-
- return result;
- }
- private:
- int k;
- T z;
- T term;
- };
-
- template <typename T, typename Policy>
- inline T lambert_w0_small_z(T z, const Policy& pol, std::integral_constant<int, 5> const&)
- {
- #ifdef BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES
- std::streamsize precision = std::cout.precision(std::numeric_limits<T>::max_digits10);
- std::cout << "Generic lambert_w0_small_z called with z = " << z << " using as many terms needed for precision." << std::endl;
- std::cout << "Argument z is of type " << typeid(T).name() << std::endl;
- #endif
-
-
-
-
-
- static const T coeff[] =
- {
- 0,
- 1,
- -1,
- static_cast<T>(3uLL) / 2uLL,
- -static_cast<T>(8uLL) / 3uLL,
- static_cast<T>(125uLL) / 24uLL,
- -static_cast<T>(54uLL) / 5uLL,
- static_cast<T>(16807uLL) / 720uLL,
- -static_cast<T>(16384uLL) / 315uLL,
- static_cast<T>(531441uLL) / 4480uLL,
- -static_cast<T>(156250uLL) / 567uLL,
- static_cast<T>(2357947691uLL) / 3628800uLL,
- -static_cast<T>(2985984uLL) / 1925uLL,
- static_cast<T>(1792160394037uLL) / 479001600uLL,
- -static_cast<T>(7909306972uLL) / 868725uLL,
- static_cast<T>(320361328125uLL) / 14350336uLL,
- -static_cast<T>(35184372088832uLL) / 638512875uLL,
- static_cast<T>(2862423051509815793uLL) / 20922789888000uLL,
- -static_cast<T>(5083731656658uLL) / 14889875uLL,
-
-
-
-
-
-
-
-
-
-
-
- };
-
- using boost::math::policies::get_epsilon;
- using boost::math::tools::sum_series;
- using boost::math::tools::evaluate_polynomial;
-
-
- T result = evaluate_polynomial(coeff, z);
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- lambert_w0_small_z_series_term<T> s(z, -pow<18>(z) / 6402373705728000uLL, 18);
-
-
-
-
-
-
-
-
- std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
- #ifdef BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES
- std::cout << "max iter from policy = " << max_iter << std::endl;
-
- #endif
- result = sum_series(s, get_epsilon<T, Policy>(), max_iter, result);
-
-
-
-
-
-
-
-
- policies::check_series_iterations<T>("boost::math::lambert_w0_small_z<%1%>(%1%)", max_iter, pol);
- #ifdef BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES_ITERATIONS
- std::cout << "z = " << z << " needed " << max_iter << " iterations." << std::endl;
- std::cout.precision(prec);
- #endif
- return result;
- }
- template <typename T>
- inline T lambert_w0_approx(T z)
- {
- BOOST_MATH_STD_USING
- T lz = log(z);
- T llz = log(lz);
- T w = lz - llz + (llz / lz);
- return w;
-
- }
-
- template <typename T>
- inline T do_get_near_singularity_param(T z)
- {
- BOOST_MATH_STD_USING
- const T p2 = 2 * (boost::math::constants::e<T>() * z + 1);
- const T p = sqrt(p2);
- return p;
- }
- template <typename T, typename Policy>
- inline T get_near_singularity_param(T z, const Policy)
- {
- using value_type = typename policies::evaluation<T, Policy>::type;
- return static_cast<T>(do_get_near_singularity_param(static_cast<value_type>(z)));
- }
- template <typename T>
- T lambert_w_positive_rational_float(T z)
- {
- BOOST_MATH_STD_USING
- if (z < 2)
- {
- if (z < T(0.5))
- {
-
-
-
- static const T Y = 8.196592331e-01f;
- static const T P[] = {
- 1.803388345e-01f,
- -4.820256838e-01f,
- -1.068349741e+00f,
- -3.506624319e-02f,
- };
- static const T Q[] = {
- 1.000000000e+00f,
- 2.871703469e+00f,
- 1.690949264e+00f,
- };
- return z * (Y + boost::math::tools::evaluate_polynomial(P, z) / boost::math::tools::evaluate_polynomial(Q, z));
- }
- else
- {
-
- static const T Y = 5.503368378e-01f;
- static const T P[] = {
- 4.493332766e-01f,
- 2.543432707e-01f,
- -4.808788799e-01f,
- -1.244425316e-01f,
- };
- static const T Q[] = {
- 1.000000000e+00f,
- 2.780661241e+00f,
- 1.830840318e+00f,
- 2.407221031e-01f,
- };
- return z * (Y + boost::math::tools::evaluate_rational(P, Q, z));
- }
- }
- else if (z < 6)
- {
-
-
- static const T Y = 1.162393570e+00f;
- static const T P[] = {
- -1.144183394e+00f,
- -4.712732855e-01f,
- 1.563162512e-01f,
- 1.434010911e-02f,
- };
- static const T Q[] = {
- 1.000000000e+00f,
- 1.192626340e+00f,
- 2.295580708e-01f,
- 5.477869455e-03f,
- };
- return Y + boost::math::tools::evaluate_rational(P, Q, z);
- }
- else if (z < 18)
- {
-
-
- static const T Y = 1.809371948e+00f;
- static const T P[] = {
- -1.689291769e+00f,
- -3.337812742e-01f,
- 3.151434873e-02f,
- 1.134178734e-03f,
- };
- static const T Q[] = {
- 1.000000000e+00f,
- 5.716915685e-01f,
- 4.489521292e-02f,
- 4.076716763e-04f,
- };
- return Y + boost::math::tools::evaluate_rational(P, Q, z);
- }
- else if (z < T(9897.12905874))
- {
-
- static const T Y = -1.402973175e+00f;
- static const T P[] = {
- 1.966174312e+00f,
- 2.350864728e-01f,
- -5.098074353e-02f,
- -1.054818339e-02f,
- };
- static const T Q[] = {
- 1.000000000e+00f,
- 4.388208264e-01f,
- 8.316639634e-02f,
- 3.397187918e-03f,
- -1.321489743e-05f,
- };
- T log_w = log(z);
- return log_w + Y + boost::math::tools::evaluate_polynomial(P, log_w) / boost::math::tools::evaluate_polynomial(Q, log_w);
- }
- else if (z < T(7.896296e+13))
- {
-
- static const T Y = -2.735729218e+00f;
- static const T P[] = {
- 3.424903470e+00f,
- 7.525631787e-02f,
- -1.427309584e-02f,
- -1.435974178e-05f,
- };
- static const T Q[] = {
- 1.000000000e+00f,
- 2.514005579e-01f,
- 6.118994652e-03f,
- -1.357889535e-05f,
- 7.312865624e-08f,
- };
- T log_w = log(z);
- return log_w + Y + boost::math::tools::evaluate_polynomial(P, log_w) / boost::math::tools::evaluate_polynomial(Q, log_w);
- }
-
- static const T Y = -4.012863159e+00f;
- static const T P[] = {
- 4.431629226e+00f,
- 2.756690487e-01f,
- -2.992956930e-03f,
- -4.912259384e-05f,
- };
- static const T Q[] = {
- 1.000000000e+00f,
- 2.015434591e-01f,
- 4.949426142e-03f,
- 1.609659944e-05f,
- -5.111523436e-09f,
- };
- T log_w = log(z);
- return log_w + Y + boost::math::tools::evaluate_polynomial(P, log_w) / boost::math::tools::evaluate_polynomial(Q, log_w);
- }
- template <typename T, typename Policy>
- T lambert_w_negative_rational_float(T z, const Policy& pol)
- {
- BOOST_MATH_STD_USING
- if (z > T(-0.27))
- {
- if (z < T(-0.051))
- {
-
-
- static const T Y = 1.255809784e+00f;
- static const T P[] = {
- -2.558083412e-01f,
- -2.306524098e+00f,
- -5.630887033e+00f,
- -3.803974556e+00f,
- };
- static const T Q[] = {
- 1.000000000e+00f,
- 5.107680783e+00f,
- 7.914062868e+00f,
- 3.501498501e+00f,
- };
- return z * (Y + boost::math::tools::evaluate_rational(P, Q, z));
- }
- else
- {
-
- return lambert_w0_small_z(z, pol);
- }
- }
- else if (z > T(-0.3578794411714423215955237701))
- {
-
- static const T Y = 1.220928431e-01f;
- static const T P[] = {
- -1.221787446e-01f,
- -6.816155875e+00f,
- 7.144582035e+01f,
- 1.128444390e+03f,
- };
- static const T Q[] = {
- 1.000000000e+00f,
- 6.480326790e+01f,
- 1.869145243e+02f,
- -1.361804274e+03f,
- 1.117826726e+03f,
- };
- T d = z + 0.367879441171442321595523770161460867445811f;
- return -d / (Y + boost::math::tools::evaluate_polynomial(P, d) / boost::math::tools::evaluate_polynomial(Q, d));
- }
- return lambert_w_singularity_series(get_near_singularity_param(z, pol));
- }
- template <typename T, typename Policy>
- inline T lambert_w0_imp(T z, const Policy& pol, const std::integral_constant<int, 1>&)
- {
- static const char* function = "boost::math::lambert_w0<%1%>";
- BOOST_MATH_STD_USING
- if ((boost::math::isnan)(z))
- {
- return boost::math::policies::raise_domain_error<T>(function, "Expected a value > -e^-1 (-0.367879...) but got %1%.", z, pol);
- }
- if ((boost::math::isinf)(z))
- {
- return boost::math::policies::raise_overflow_error<T>(function, "Expected a finite value but got %1%.", z, pol);
- }
- if (z >= T(0.05))
-
- {
- return lambert_w_positive_rational_float(z);
- }
- else if (z <= -0.3678794411714423215955237701614608674458111310f)
- {
- if (z < -0.3678794411714423215955237701614608674458111310f)
- return boost::math::policies::raise_domain_error<T>(function, "Expected z >= -e^-1 (-0.367879...) but got %1%.", z, pol);
- return -1;
- }
- return lambert_w_negative_rational_float(z, pol);
- }
- template <typename T>
- T lambert_w_positive_rational_double(T z)
- {
- BOOST_MATH_STD_USING
- if (z < 2)
- {
- if (z < 0.5)
- {
-
- static const T offset = 8.19659233093261719e-01;
- static const T P[] = {
- 1.80340766906685177e-01,
- 3.28178241493119307e-01,
- -2.19153620687139706e+00,
- -7.24750929074563990e+00,
- -7.28395876262524204e+00,
- -2.57417169492512916e+00,
- -2.31606948888704503e-01
- };
- static const T Q[] = {
- 1.00000000000000000e+00,
- 7.36482529307436604e+00,
- 2.03686007856430677e+01,
- 2.62864592096657307e+01,
- 1.59742041380858333e+01,
- 4.03760534788374589e+00,
- 2.91327346750475362e-01
- };
- return z * (offset + boost::math::tools::evaluate_polynomial(P, z) / boost::math::tools::evaluate_polynomial(Q, z));
- }
- else
- {
-
- static const T offset = 5.50335884094238281e-01;
- static const T P[] = {
- 4.49664083944098322e-01,
- 1.90417666196776909e+00,
- 1.99951368798255994e+00,
- -6.91217310299270265e-01,
- -1.88533935998617058e+00,
- -7.96743968047750836e-01,
- -1.02891726031055254e-01,
- -3.09156013592636568e-03
- };
- static const T Q[] = {
- 1.00000000000000000e+00,
- 6.45854489419584014e+00,
- 1.54739232422116048e+01,
- 1.72606164253337843e+01,
- 9.29427055609544096e+00,
- 2.29040824649748117e+00,
- 2.21610620995418981e-01,
- 5.70597669908194213e-03
- };
- return z * (offset + boost::math::tools::evaluate_rational(P, Q, z));
- }
- }
- else if (z < 6)
- {
-
-
- static const T Y = 1.16239356994628906e+00;
- static const T P[] = {
- -1.16230494982099475e+00,
- -3.38528144432561136e+00,
- -2.55653717293161565e+00,
- -3.06755172989214189e-01,
- 1.73149743765268289e-01,
- 3.76906042860014206e-02,
- 1.84552217624706666e-03,
- 1.69434126904822116e-05,
- };
- static const T Q[] = {
- 1.00000000000000000e+00,
- 3.77187616711220819e+00,
- 4.58799960260143701e+00,
- 2.24101228462292447e+00,
- 4.54794195426212385e-01,
- 3.60761772095963982e-02,
- 9.25176499518388571e-04,
- 4.43611344705509378e-06,
- };
- return Y + boost::math::tools::evaluate_rational(P, Q, z);
- }
- else if (z < 18)
- {
-
-
- static const T offset = 1.80937194824218750e+00;
- static const T P[] =
- {
- -1.80690935424793635e+00,
- -3.66995929380314602e+00,
- -1.93842957940149781e+00,
- -2.94269984375794040e-01,
- 1.81224710627677778e-03,
- 2.48166798603547447e-03,
- 1.15806592415397245e-04,
- 1.43105573216815533e-06,
- 3.47281483428369604e-09
- };
- static const T Q[] = {
- 1.00000000000000000e+00,
- 2.57319080723908597e+00,
- 1.96724528442680658e+00,
- 5.84501352882650722e-01,
- 7.37152837939206240e-02,
- 3.97368430940416778e-03,
- 8.54941838187085088e-05,
- 6.05713225608426678e-07,
- 8.17517283816615732e-10
- };
- return offset + boost::math::tools::evaluate_rational(P, Q, z);
- }
- else if (z < 9897.12905874)
- {
-
- static const T Y = -1.40297317504882812e+00;
- static const T P[] = {
- 1.97011826279311924e+00,
- 1.05639945701546704e+00,
- 3.33434529073196304e-01,
- 3.34619153200386816e-02,
- -5.36238353781326675e-03,
- -2.43901294871308604e-03,
- -2.13762095619085404e-04,
- -4.85531936495542274e-06,
- -2.02473518491905386e-08,
- };
- static const T Q[] = {
- 1.00000000000000000e+00,
- 8.60107275833921618e-01,
- 4.10420467985504373e-01,
- 1.18444884081994841e-01,
- 2.16966505556021046e-02,
- 2.24529766630769097e-03,
- 9.82045090226437614e-05,
- 1.36363515125489502e-06,
- 3.44200749053237945e-09,
- };
- T log_w = log(z);
- return log_w + Y + boost::math::tools::evaluate_rational(P, Q, log_w);
- }
- else if (z < 7.896296e+13)
- {
-
- static const T Y = -2.73572921752929688e+00;
- static const T P[] = {
- 3.30547638424076217e+00,
- 1.64050071277550167e+00,
- 4.57149576470736039e-01,
- 4.03821227745424840e-02,
- -4.99664976882514362e-04,
- -1.28527893803052956e-04,
- -2.95470325373338738e-06,
- -1.76662025550202762e-08,
- -1.98721972463709290e-11,
- };
- static const T Q[] = {
- 1.00000000000000000e+00,
- 6.91472559412458759e-01,
- 2.48154578891676774e-01,
- 4.60893578284335263e-02,
- 3.60207838982301946e-03,
- 1.13001153242430471e-04,
- 1.33690948263488455e-06,
- 4.97253225968548872e-09,
- 3.39460723731970550e-12,
- };
- T log_w = log(z);
- return log_w + Y + boost::math::tools::evaluate_rational(P, Q, log_w);
- }
- else if (z < 2.6881171e+43)
- {
-
- static const T Y = -4.01286315917968750e+00;
- static const T P[] = {
- 5.07714858354309672e+00,
- -3.32994414518701458e+00,
- -8.61170416909864451e-01,
- -4.01139705309486142e-02,
- -1.85374201771834585e-04,
- 1.08824145844270666e-05,
- 1.17216905810452396e-07,
- 2.97998248101385990e-10,
- 1.42294856434176682e-13,
- };
- static const T Q[] = {
- 1.00000000000000000e+00,
- -4.85840770639861485e-01,
- -3.18714850604827580e-01,
- -3.20966129264610534e-02,
- -1.06276178044267895e-03,
- -1.33597828642644955e-05,
- -6.27900905346219472e-08,
- -9.35271498075378319e-11,
- -2.60648331090076845e-14,
- };
- T log_w = log(z);
- return log_w + Y + boost::math::tools::evaluate_rational(P, Q, log_w);
- }
- else
- {
-
- static const T Y = -5.70115661621093750e+00;
- static const T P[] = {
- 6.42275660145116698e+00,
- 1.33047964073367945e+00,
- 6.72008923401652816e-02,
- 1.16444069958125895e-03,
- 7.06966760237470501e-06,
- 5.48974896149039165e-09,
- -7.00379652018853621e-11,
- -1.89247635913659556e-13,
- -1.55898770790170598e-16,
- -4.06109208815303157e-20,
- -2.21552699006496737e-24,
- };
- static const T Q[] = {
- 1.00000000000000000e+00,
- 3.34498588416632854e-01,
- 2.51519862456384983e-02,
- 6.81223810622416254e-04,
- 7.94450897106903537e-06,
- 4.30675039872881342e-08,
- 1.10667669458467617e-10,
- 1.31012240694192289e-13,
- 6.53282047177727125e-17,
- 1.11775518708172009e-20,
- 3.78250395617836059e-25,
- };
- T log_w = log(z);
- return log_w + Y + boost::math::tools::evaluate_rational(P, Q, log_w);
- }
- }
- template <typename T, typename Policy>
- T lambert_w_negative_rational_double(T z, const Policy& pol)
- {
- BOOST_MATH_STD_USING
- if (z > -0.1)
- {
- if (z < -0.051)
- {
-
-
-
-
- static const T Y = 1.08633995056152344e+00;
- static const T P[] = {
- -8.63399505615014331e-02,
- -1.64303871814816464e+00,
- -7.71247913918273738e+00,
- -1.41014495545382454e+01,
- -1.02269079949257616e+01,
- -2.17236002836306691e+00,
- };
- static const T Q[] = {
- 1.00000000000000000e+00,
- 7.44775406945739243e+00,
- 2.04392643087266541e+01,
- 2.51001961077774193e+01,
- 1.31256080849023319e+01,
- 2.11640324843601588e+00,
- };
- return z * (Y + boost::math::tools::evaluate_rational(P, Q, z));
- }
- else
- {
-
- return lambert_w0_small_z(z, pol);
- }
- }
- else if (z > -0.2)
- {
-
-
-
-
- static const T Y = 1.20359611511230469e+00;
- static const T P[] = {
- -2.03596115108465635e-01,
- -2.95029082937201859e+00,
- -1.54287922188671648e+01,
- -3.81185809571116965e+01,
- -4.66384358235575985e+01,
- -2.59282069989642468e+01,
- -4.70140451266553279e+00,
- };
- static const T Q[] = {
- 1.00000000000000000e+00,
- 9.57921436074599929e+00,
- 3.60988119290234377e+01,
- 6.73977699505546007e+01,
- 6.41104992068148823e+01,
- 2.82060127225153607e+01,
- 4.10677610657724330e+00,
- };
- return z * (Y + boost::math::tools::evaluate_rational(P, Q, z));
- }
- else if (z > -0.3178794411714423215955237)
- {
-
- static const T Y = 3.49680423736572266e-01;
- static const T P[] = {
- -3.49729841718749014e-01,
- -6.28207407760709028e+01,
- -2.57226178029669171e+03,
- -2.50271008623093747e+04,
- 1.11949239154711388e+05,
- 1.85684566607844318e+06,
- 4.80802490427638643e+06,
- 2.76624752134636406e+06,
- };
- static const T Q[] = {
- 1.00000000000000000e+00,
- 1.82717661215113000e+02,
- 8.00121119810280100e+03,
- 1.06073266717010129e+05,
- 3.22848993926057721e+05,
- -8.05684814514171256e+05,
- -2.59223192927265737e+06,
- -5.61719645211570871e+05,
- 6.27765369292636844e+04,
- };
- T d = z + 0.367879441171442321595523770161460867445811;
- return -d / (Y + boost::math::tools::evaluate_polynomial(P, d) / boost::math::tools::evaluate_polynomial(Q, d));
- }
- else if (z > -0.3578794411714423215955237701)
- {
-
- static const T Y = 5.00126481056213379e-02;
- static const T P[] = {
- -5.00173570682372162e-02,
- -4.44242461870072044e+01,
- -9.51185533619946042e+03,
- -5.88605699015429386e+05,
- -1.90760843597427751e+06,
- 5.79797663818311404e+08,
- 1.11383352508459134e+10,
- 5.67791253678716467e+10,
- 6.32694500716584572e+10,
- };
- static const T Q[] = {
- 1.00000000000000000e+00,
- 9.08910517489981551e+02,
- 2.10170163753340133e+05,
- 1.67858612416470327e+07,
- 4.90435561733227953e+08,
- 4.54978142622939917e+09,
- 2.87716585708739168e+09,
- -4.59414247951143131e+10,
- -1.72845216404874299e+10,
- };
- T d = z + 0.36787944117144232159552377016146086744581113103176804;
- return -d / (Y + boost::math::tools::evaluate_polynomial(P, d) / boost::math::tools::evaluate_polynomial(Q, d));
- }
- else
- {
-
- const T p2 = 2 * (boost::math::constants::e<T>() * z + 1);
- const T p = sqrt(p2);
- return lambert_w_detail::lambert_w_singularity_series(p);
- }
- }
- template <typename T, typename Policy>
- inline T lambert_w0_imp(T z, const Policy& pol, const std::integral_constant<int, 2>&)
- {
- static const char* function = "boost::math::lambert_w0<%1%>";
- BOOST_MATH_STD_USING
-
- static_assert(std::numeric_limits<double>::digits >= 53,
- "Our double precision coefficients will be truncated, "
- "please file a bug report with details of your platform's floating point types "
- "- 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<T>(function, "Expected a value > -e^-1 (-0.367879...) but got %1%.", z, pol);
- }
- if ((boost::math::isinf)(z))
- {
- return boost::math::policies::raise_overflow_error<T>(function, "Expected a finite value but got %1%.", z, pol);
- }
- if (z >= 0.05)
- {
- return lambert_w_positive_rational_double(z);
- }
- else if (z <= -0.36787944117144232159552377016146086744581113103176804)
- {
- if (z < -0.36787944117144232159552377016146086744581113103176804)
- {
- return boost::math::policies::raise_domain_error<T>(function, "Expected z >= -e^-1 (-0.367879...) but got %1%.", z, pol);
- }
- return -1;
- }
- else
- {
- return lambert_w_negative_rational_double(z, pol);
- }
- }
- template <typename T, typename Policy>
- inline T lambert_w0_imp(T z, const Policy& pol, const std::integral_constant<int, 0>&)
- {
- static const char* function = "boost::math::lambert_w0<%1%>";
- BOOST_MATH_STD_USING
-
- if ((boost::math::isnan)(z))
- {
- return boost::math::policies::raise_domain_error<T>(function, "Expected z >= -e^-1 (-0.367879...) but got %1%.", z, pol);
- }
- if (fabs(z) <= 0.05f)
- {
-
- return lambert_w0_small_z(z, pol);
- }
- if (z > (std::numeric_limits<double>::max)())
- {
- if ((boost::math::isinf)(z))
- {
- return policies::raise_overflow_error<T>(function, nullptr, pol);
-
-
- }
-
-
- T w = lambert_w0_approx(z);
-
-
-
- return lambert_w_halley_iterate(w, z);
- }
- if (z < -0.3578794411714423215955237701)
- {
- if (z <= -boost::math::constants::exp_minus_one<T>())
- {
- if (z == -boost::math::constants::exp_minus_one<T>())
- {
- return -1;
- }
- return boost::math::policies::raise_domain_error<T>(function, "Expected z >= -e^-1 (-0.367879...) but got %1%.", z, pol);
- }
-
-
- const T p2 = 2 * (boost::math::constants::e<T>() * z + 1);
- const T p = sqrt(p2);
- T w = lambert_w_detail::lambert_w_singularity_series(p);
- return lambert_w_halley_iterate(w, z);
- }
-
-
-
-
-
- using precision_type = typename policies::precision<T, Policy>::type;
- using tag_type = std::integral_constant<bool,
- (precision_type::value == 0) || (precision_type::value > 113) ?
- true
- : false
- >;
-
-
- T w = lambert_w0_imp(maybe_reduce_to_double(z, std::is_constructible<double, T>()), pol, std::integral_constant<int, 2>());
- return lambert_w_maybe_halley_iterate(w, z, tag_type());
- }
-
-
-
- template<typename T, typename Policy>
- T lambert_wm1_imp(const T z, const Policy& pol)
- {
-
-
-
-
-
-
-
- static_assert(!std::is_integral<T>::value,
- "Must be floating-point or fixed type (not integer type), for example: lambert_wm1(1.), not lambert_wm1(1)!");
- BOOST_MATH_STD_USING
- const char* function = "boost::math::lambert_wm1<RealType>(<RealType>)";
-
- if ((boost::math::isnan)(z))
- {
- return policies::raise_domain_error(function,
- "Argument z is NaN!",
- z, pol);
- }
- if ((boost::math::isinf)(z))
- {
- return policies::raise_domain_error(function,
- "Argument z is infinite!",
- z, pol);
- }
- if (z == static_cast<T>(0))
- {
- if (std::numeric_limits<T>::has_infinity)
- {
- return -std::numeric_limits<T>::infinity();
- }
- else
- {
- return -tools::max_value<T>();
- }
- }
- if (boost::math::detail::has_denorm_now<T>())
- {
- if (!(boost::math::isnormal)(z))
- {
- return policies::raise_overflow_error(function,
- "Argument z = %1% is denormalized! (must be z > (std::numeric_limits<RealType>::min)() or z == 0)",
- z, pol);
- }
- }
- if (z > static_cast<T>(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 > -boost::math::tools::min_value<T>())
- {
-
-
- return policies::raise_overflow_error(function,
- "Argument z = %1% is too small (z < -std::numeric_limits<T>::min so denormalized) for Lambert W-1 branch!",
- z, pol);
- }
- if (z == -boost::math::constants::exp_minus_one<T>())
- {
- return -static_cast<T>(1);
- }
-
- if (z < -boost::math::constants::exp_minus_one<T>())
- {
- 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 < static_cast<T>(-0.35))
- {
- const T p2 = 2 * (boost::math::constants::e<T>() * z + 1);
- if (p2 == 0)
- {
- return -1;
- }
- if (p2 > 0)
- {
- T w_series = lambert_w_singularity_series(T(-sqrt(p2)));
- if (boost::math::tools::digits<T>() > 53)
- {
- w_series = lambert_w_detail::lambert_w_halley_iterate(w_series, z);
- #ifdef BOOST_MATH_INSTRUMENT_LAMBERT_WM1_NOT_BUILTIN
- std::streamsize saved_precision = std::cout.precision(std::numeric_limits<T>::max_digits10);
- std::cout << "Lambert W-1 Halley updated to " << w_series << std::endl;
- std::cout.precision(saved_precision);
- #endif
- }
- return w_series;
- }
-
- return policies::raise_domain_error(function,
- "Argument z = %1% is out of range for Lambert W-1 branch. (Should not get here - please report!)",
- z, pol);
- }
- using lambert_w_lookup::wm1es;
- using lambert_w_lookup::wm1zs;
- using lambert_w_lookup::noof_wm1zs;
-
-
-
- if (z >= T(wm1zs[63]))
- {
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- T guess;
- T lz = log(-z);
- T llz = log(-lz);
- guess = lz - llz + (llz / lz);
- #ifdef BOOST_MATH_INSTRUMENT_LAMBERT_WM1_TINY
- std::streamsize saved_precision = std::cout.precision(std::numeric_limits<T>::max_digits10);
- std::cout << "z = " << z << ", guess = " << guess << ", ln(-z) = " << lz << ", ln(-ln(-z) = " << llz << ", llz/lz = " << (llz / lz) << std::endl;
-
-
-
- int d10 = policies::digits_base10<T, Policy>();
- int d2 = policies::digits<T, Policy>();
- std::cout << "digits10 = " << d10 << ", digits2 = " << d2
- << std::endl;
- std::cout.precision(saved_precision);
- #endif
- if (policies::digits<T, Policy>() < 12)
- {
- return guess;
- }
- T result = lambert_w_detail::lambert_w_halley_iterate(guess, z);
- return result;
-
-
-
-
-
- }
-
- if (boost::math::tools::digits<T>() > 53)
- {
-
-
-
- using boost::math::policies::precision;
- using boost::math::policies::digits10;
- using boost::math::policies::digits2;
- using boost::math::policies::policy;
-
- T double_approx(static_cast<T>(lambert_wm1_imp(must_reduce_to_double(z, std::is_constructible<double, T>()), policy<digits2<50>>())));
- #ifdef BOOST_MATH_INSTRUMENT_LAMBERT_WM1_NOT_BUILTIN
- std::streamsize saved_precision = std::cout.precision(std::numeric_limits<T>::max_digits10);
- std::cout << "Lambert_wm1 Argument Type " << typeid(T).name() << " approximation double = " << double_approx << std::endl;
- std::cout.precision(saved_precision);
- #endif
-
-
- T result = lambert_w_halley_iterate(double_approx, z);
- #ifdef BOOST_MATH_INSTRUMENT_LAMBERT_WM1
- std::streamsize saved_precision = std::cout.precision(std::numeric_limits<T>::max_digits10);
- std::cout << "Result " << typeid(T).name() << " precision Halley refinement = " << result << std::endl;
- std::cout.precision(saved_precision);
- #endif
- return result;
- }
- else
- {
- using namespace boost::math::lambert_w_detail::lambert_w_lookup;
-
-
- int n = 2;
- if (T(wm1zs[n - 1]) > z)
- {
- goto bisect;
- }
- for (int j = 1; j <= 5; ++j)
- {
- n *= 2;
- if (T(wm1zs[n - 1]) > z)
- {
- goto overshot;
- }
- }
-
-
- return policies::raise_domain_error(function,
- "Argument z = %1% is too small (< -1.026439e-26) (logic error - please report!)",
- z, pol);
- overshot:
- {
- int nh = n / 2;
- for (int j = 1; j <= 5; ++j)
- {
- nh /= 2;
- if (nh <= 0)
- {
- break;
- }
- if (T(wm1zs[n - nh - 1]) > z)
- {
- n -= nh;
- }
- }
- }
- bisect:
- --n;
-
-
- #ifdef BOOST_MATH_INSTRUMENT_LAMBERT_WM1_LOOKUP
- std::streamsize saved_precision = std::cout.precision(std::numeric_limits<T>::max_digits10);
- std::cout << "Result lookup W-1(" << z << ") bisection between wm1zs[" << n - 1 << "] = " << wm1zs[n - 1] << " and wm1zs[" << n << "] = " << wm1zs[n]
- << ", bisect mean = " << (wm1zs[n - 1] + wm1zs[n]) / 2 << std::endl;
- std::cout.precision(saved_precision);
- #endif
-
-
-
-
- int bisections = 11;
- if (n >= 8)
- {
- bisections = 8;
- }
- else if (n >= 3)
- {
- bisections = 9;
- }
- else if (n >= 2)
- {
- bisections = 10;
- }
-
-
-
-
- using lambert_w_lookup::halves;
- using lambert_w_lookup::sqrtwm1s;
- using calc_type = typename std::conditional<std::is_constructible<lookup_t, T>::value, lookup_t, T>::type;
- calc_type w = -static_cast<calc_type>(n);
- calc_type y = static_cast<calc_type>(z * T(wm1es[n - 1]));
-
- for (int j = 0; j < bisections; ++j)
- {
- calc_type wj = w - halves[j];
- calc_type yj = y * sqrtwm1s[j];
- if (wj < yj)
- {
- w = wj;
- y = yj;
- }
- }
- return static_cast<T>(schroeder_update(w, y));
- }
- }
- }
- template <typename T, typename Policy>
- inline
- typename boost::math::tools::promote_args<T>::type
- lambert_w0(T z, const Policy& pol)
- {
-
-
- using result_type = typename tools::promote_args<T>::type;
-
-
- using precision_type = typename policies::precision<result_type, Policy>::type;
-
- using tag_type = std::integral_constant<int,
- (precision_type::value == 0) || (precision_type::value > 53) ?
- 0
- : (precision_type::value <= 24) ? 1
- : 2
- >;
- return lambert_w_detail::lambert_w0_imp(result_type(z), pol, tag_type());
- }
-
- template <typename T>
- inline
- typename tools::promote_args<T>::type
- lambert_w0(T z)
- {
-
-
- using result_type = typename tools::promote_args<T>::type;
-
-
- using precision_type = typename policies::precision<result_type, policies::policy<>>::type;
-
- using tag_type = std::integral_constant<int,
- (precision_type::value == 0) || (precision_type::value > 53) ?
- 0
- : (precision_type::value <= 24) ? 1
- : 2
- >;
- return lambert_w_detail::lambert_w0_imp(result_type(z), policies::policy<>(), tag_type());
- }
-
-
- template <typename T, typename Policy>
- inline
- typename tools::promote_args<T>::type
- lambert_wm1(T z, const Policy& pol)
- {
-
-
- using result_type = typename tools::promote_args<T>::type;
- return lambert_w_detail::lambert_wm1_imp(result_type(z), pol);
- }
-
- template <typename T>
- inline
- typename tools::promote_args<T>::type
- lambert_wm1(T z)
- {
- using result_type = typename tools::promote_args<T>::type;
- return lambert_w_detail::lambert_wm1_imp(result_type(z), policies::policy<>());
- }
-
- template <typename T, typename Policy>
- inline typename tools::promote_args<T>::type
- lambert_w0_prime(T z, const Policy& pol)
- {
- using result_type = typename tools::promote_args<T>::type;
- using std::numeric_limits;
- if (z == 0)
- {
- return static_cast<result_type>(1);
- }
-
-
- if (z == - boost::math::constants::exp_minus_one<result_type>())
- {
- return numeric_limits<result_type>::has_infinity ? numeric_limits<result_type>::infinity() : boost::math::tools::max_value<result_type>();
- }
-
- result_type w = lambert_w0(result_type(z), pol);
-
-
-
-
-
-
-
- return w / (z * (1 + w));
- }
- template <typename T>
- inline typename tools::promote_args<T>::type
- lambert_w0_prime(T z)
- {
- return lambert_w0_prime(z, policies::policy<>());
- }
- template <typename T, typename Policy>
- inline typename tools::promote_args<T>::type
- lambert_wm1_prime(T z, const Policy& pol)
- {
- using std::numeric_limits;
- using result_type = typename tools::promote_args<T>::type;
-
-
-
-
-
- if (z == 0 || z == - boost::math::constants::exp_minus_one<result_type>())
- {
- return numeric_limits<result_type>::has_infinity ? -numeric_limits<result_type>::infinity() : -boost::math::tools::max_value<result_type>();
- }
- result_type w = lambert_wm1(z, pol);
- return w/(z*(1+w));
- }
- template <typename T>
- inline typename tools::promote_args<T>::type
- lambert_wm1_prime(T z)
- {
- return lambert_wm1_prime(z, policies::policy<>());
- }
- }}
- #endif
|