libs/math/tools/lambert_w_lookup_table_generator.cpp

   1 // Lambert W lookup table generator lambert_w_lookup_table_generator.cpp
   2
   3 //! \file
   4 //! Output a table of precomputed array values for Lambert W0 and W-1,
   5 //! and square roots and halves, and powers of e,
   6 // as a .ipp file for use by lambert_w.hpp.
   7
   8 //! \details Output as long double precision (suffix L) using Boost.Multiprecision
   9 //! to 34 decimal digits precision to cater for platforms that have 128-bit long double.
  10 //! The function bisection can then use any built-in floating-point type,
  11 //! which may have different precision and speed on different platforms.
  12 //! The actual builtin floating-point type of the arrays is chosen by a
  13 //! typedef in \modular-boost\libs\math\include\boost\math\special_functions\lambert_w.hpp
  14 //! by default, for example: typedef double lookup_t;
  15
  16
  17 // This includes lookup tables for both branches W0 and W-1.
  18 // Only W-1 is needed by current code that uses JM rational Polynomials,
  19 // but W0 is kept (for now) to allow comparison with the previous FKDVPB version
  20 // that uses lookup for W0 branch as well as W-1.
  21
  22 #include <boost/config.hpp>
  23 #include <boost/math/constants/constants.hpp> // For exp_minus_one == 3.67879441171442321595523770161460867e-01.
  24 using boost::math::constants::exp_minus_one; // 0.36787944
  25 using boost::math::constants::root_e; // 1.64872
  26 #include <boost/multiprecision/cpp_bin_float.hpp>
  27 using boost::multiprecision::cpp_bin_float_quad;
  28 using boost::multiprecision::cpp_bin_float_50;
  29
  30 #include <iostream>
  31 #include <fstream>
  32 #include <typeinfo>
  33
  34 /*
  35 typedef double lookup_t; // Type for lookup table (double or float?)
  36
  37 BOOST_STATIC_CONSTEXPR std::size_t noof_sqrts = 12;
  38 BOOST_STATIC_CONSTEXPR lookup_t a[noof_sqrts] = // 0.6065 0.7788, 0.8824 ... 0.9997, sqrt of previous elements.
  39 {
  40 0.60653065971263342, 0.77880078307140487, 0.8824969025845954, 0.93941306281347579, 0.96923323447634408, 0.98449643700540841,
  41 0.99221793826024351, 0.99610136947011749, 0.99804878110747547, 0.99902391418197566, 0.99951183793988937, 0.99975588917489722
  42 };
  43 BOOST_STATIC_CONSTEXPR  lookup_t b[noof_sqrts] = // 0.5 0.25 0.125, 0.0625 ...  0.0002441, halves of previous elements.
  44 { 0.5, 0.25, 0.125, 0.0625, 0.03125, 0.015625, 0.0078125, 0.00390625, 0.001953125, 0.0009765625, 0.00048828125, 0.000244140625 };
  45
  46 BOOST_STATIC_CONSTEXPR size_t noof_w0zs = 65;
  47 BOOST_STATIC_CONSTEXPR lookup_t g[noof_w0zs] = // lambert_w[k] for W0 branch. 0 2.7182 14.77 60.2566 ... 1.445e+29 3.990e+29.
  48 { 0., 2.7182818284590452, 14.7781121978613, 60.256610769563003, 218.39260013257696, 742.06579551288302, 2420.5727609564107, 7676.4321089992102,
  49 23847.663896333826, 72927.755348178456, 220264.65794806717, 658615.558867176, 1953057.4970280471, 5751374.0961159665, 16836459.978306875, 49035260.58708166,
  50 142177768.32812596, 410634196.81078007, 1181879444.4719492, 3391163718.300558, 9703303908.1958056, 27695130424.147509, 78868082614.895014, 224130479263.72476,
  51 635738931116.24334, 1800122483434.6468, 5088969845149.8079, 14365302496248.563, 40495197800161.305, 114008694617177.22, 320594237445733.86,
  52 900514339622670.18, 2526814725845782.2, 7083238132935230.1, 19837699245933466., 55510470830970076., 1.5520433569614703e+17, 4.3360826779369662e+17,
  53 1.2105254067703227e+18, 3.3771426165357561e+18, 9.4154106734807994e+18, 2.6233583234732253e+19, 7.3049547543861044e+19, 2.032970971338619e+20,
  54 5.6547040503180956e+20, 1.5720421975868293e+21, 4.3682149334771265e+21, 1.2132170565093317e+22, 3.3680332378068632e+22, 9.3459982052259885e+22,
  55 2.5923527642935362e+23, 7.1876803203773879e+23, 1.99212416037262e+24, 5.5192924995054165e+24, 1.5286067837683347e+25, 4.2321318958281094e+25,
  56 1.1713293177672778e+26, 3.2408603996214814e+26, 8.9641258264226028e+26, 2.4787141382364034e+27, 6.8520443388941057e+27, 1.8936217407781711e+28,
  57 5.2317811346197018e+28, 1.4450833904658542e+29, 3.9904954117194348e+29
  58 };
  59
  60 BOOST_STATIC_CONSTEXPR std::size_t noof_wm1zs = 66;
  61 BOOST_STATIC_CONSTEXPR lookup_t e[noof_wm1zs] = // lambert_w[k] for W-1 branch. 2.7182 1. 0.3678 0.1353 0.04978 ... 4.359e-28 1.603e-28
  62 {
  63 2.7182818284590452, 1., 0.36787944117144232, 0.13533528323661269, 0.049787068367863943, 0.01831563888873418, 0.0067379469990854671,
  64 0.0024787521766663584, 0.00091188196555451621, 0.00033546262790251184, 0.00012340980408667955, 4.5399929762484852e-05, 1.6701700790245659e-05,
  65 6.1442123533282098e-06, 2.2603294069810543e-06, 8.3152871910356788e-07, 3.0590232050182579e-07, 1.1253517471925911e-07, 4.1399377187851667e-08,
  66 1.5229979744712628e-08, 5.6027964375372675e-09, 2.0611536224385578e-09, 7.5825604279119067e-10, 2.7894680928689248e-10, 1.026187963170189e-10,
  67 3.7751345442790977e-11, 1.3887943864964021e-11, 5.1090890280633247e-12, 1.8795288165390833e-12, 6.914400106940203e-13, 2.5436656473769229e-13,
  68 9.3576229688401746e-14, 3.4424771084699765e-14, 1.2664165549094176e-14, 4.6588861451033974e-15, 1.713908431542013e-15, 6.3051167601469894e-16,
  69 2.3195228302435694e-16, 8.5330476257440658e-17, 3.1391327920480296e-17, 1.1548224173015786e-17, 4.248354255291589e-18, 1.5628821893349888e-18,
  70 5.7495222642935598e-19, 2.1151310375910805e-19, 7.7811322411337965e-20, 2.8625185805493936e-20, 1.0530617357553812e-20, 3.8739976286871871e-21,
  71 1.4251640827409351e-21, 5.2428856633634639e-22, 1.9287498479639178e-22, 7.0954741622847041e-23, 2.6102790696677048e-23, 9.602680054508676e-24,
  72 3.532628572200807e-24, 1.2995814250075031e-24, 4.7808928838854691e-25, 1.7587922024243116e-25, 6.4702349256454603e-26, 2.3802664086944006e-26,
  73 8.7565107626965203e-27, 3.2213402859925161e-27, 1.185064864233981e-27, 4.359610000063081e-28, 1.6038108905486378e-28
  74 };
  75
  76 lambert_w0 version of array from Fukushima
  77
  78 // lambert_w[k] for W-1 branch. 2.7182 1. 0.3678 0.1353 0.04978 ... 4.359e-28 1.603e-28
  79 e: 2.7182818284590452, 1., 0.36787944117144232, 0.13533528323661269, 0.049787068367863943, 0.01831563888873418, 0.0067379469990854671,
  80 0.0024787521766663584, 0.00091188196555451621, 0.00033546262790251184, 0.00012340980408667955, 4.5399929762484852e-05, 1.6701700790245659e-05,
  81 6.1442123533282098e-06, 2.2603294069810543e-06, 8.3152871910356788e-07, 3.0590232050182579e-07, 1.1253517471925911e-07, 4.1399377187851667e-08,
  82 1.5229979744712628e-08, 5.6027964375372675e-09, 2.0611536224385578e-09, 7.5825604279119067e-10, 2.7894680928689248e-10, 1.026187963170189e-10,
  83 3.7751345442790977e-11, 1.3887943864964021e-11, 5.1090890280633247e-12, 1.8795288165390833e-12, 6.914400106940203e-13, 2.5436656473769229e-13,
  84 9.3576229688401746e-14, 3.4424771084699765e-14, 1.2664165549094176e-14, 4.6588861451033974e-15, 1.713908431542013e-15, 6.3051167601469894e-16,
  85 2.3195228302435694e-16, 8.5330476257440658e-17, 3.1391327920480296e-17, 1.1548224173015786e-17, 4.248354255291589e-18, 1.5628821893349888e-18,
  86 5.7495222642935598e-19, 2.1151310375910805e-19, 7.7811322411337965e-20, 2.8625185805493936e-20, 1.0530617357553812e-20, 3.8739976286871871e-21,
  87 1.4251640827409351e-21, 5.2428856633634639e-22, 1.9287498479639178e-22, 7.0954741622847041e-23, 2.6102790696677048e-23, 9.602680054508676e-24,
  88 3.532628572200807e-24, 1.2995814250075031e-24, 4.7808928838854691e-25, 1.7587922024243116e-25, 6.4702349256454603e-26, 2.3802664086944006e-26,
  89 8.7565107626965203e-27, 3.2213402859925161e-27, 1.185064864233981e-27, 4.359610000063081e-28, 1.6038108905486378e-28
  90
  91 // lambert_w[k] for W0 branch. 0 2.7182 14.77 60.2566 ... 1.445e+29 3.990e+29.
  92
  93 g: 0, 2.7182818284590452, 14.7781121978613, 60.256610769563003, 218.39260013257696, 742.06579551288302, 2420.5727609564107, 7676.4321089992102,
  94 23847.663896333826, 72927.755348178456, 220264.65794806717, 658615.558867176, 1953057.4970280471, 5751374.0961159665, 16836459.978306875, 49035260.58708166,
  95 142177768.32812596, 410634196.81078007, 1181879444.4719492, 3391163718.300558, 9703303908.1958056, 27695130424.147509, 78868082614.895014, 224130479263.72476,
  96 635738931116.24334, 1800122483434.6468, 5088969845149.8079, 14365302496248.563, 40495197800161.305, 114008694617177.22, 320594237445733.86,
  97 900514339622670.18, 2526814725845782.2, 7083238132935230.1, 19837699245933466, 55510470830970076, 1.5520433569614703e+17, 4.3360826779369662e+17,
  98 1.2105254067703227e+18, 3.3771426165357561e+18, 9.4154106734807994e+18, 2.6233583234732253e+19, 7.3049547543861044e+19, 2.032970971338619e+20,
  99 5.6547040503180956e+20, 1.5720421975868293e+21, 4.3682149334771265e+21, 1.2132170565093317e+22, 3.3680332378068632e+22, 9.3459982052259885e+22,
 100 2.5923527642935362e+23, 7.1876803203773879e+23, 1.99212416037262e+24, 5.5192924995054165e+24, 1.5286067837683347e+25, 4.2321318958281094e+25,
 101 1.1713293177672778e+26, 3.2408603996214814e+26, 8.9641258264226028e+26, 2.4787141382364034e+27, 6.8520443388941057e+27, 1.8936217407781711e+28,
 102 5.2317811346197018e+28, 1.4450833904658542e+29, 3.9904954117194348e+29
 103
 104
 105 lambert_wm1 version of arrays from Fukushima
 106 e: 2.7182817459106445 7.3890557289123535 20.085535049438477 54.59814453125 148.41314697265625 403.42874145507813 1096.6329345703125 2980.957275390625 8103.08154296875 22026.458984375 59874.12109375 162754.734375 442413.21875 1202603.75 3269015.75 8886106 24154940 65659932 178482192 485164896 1318814848 3584910336 9744796672 26489102336 72004845568 195729457152 532047822848 1446255919104 3931331100672 10686465835008 29048824659968 78962889850880 214643389759488 583461240832000 1586012102852608 4311227773747200 11719131799748608 31855901283450880 86593318145753088 2.3538502982225101e+17 6.398428560008151e+17 1.7392731886358364e+18 4.7278345784949473e+18 1.2851586685678387e+19 3.493423319351296e+19 9.4961089747571704e+19 2.581309902546461e+20 7.0167278463083348e+20 1.9073443887231177e+21 5.1846992652160593e+21 1.4093473476000776e+22 3.831003235981371e+22 1.0413746376682761e+23 2.8307496154307266e+23 7.6947746628514896e+23 2.0916565667371597e+24 5.6857119515524837e+24 1.5455367020327599e+25 4.2012039964445827e+25 1.1420056438012293e+26 3.1042929865047826e+26 8.4383428037470738e+26 2.2937792813113457e+27 6.2351382164292627e+27
 107
 108 g: -0.36787945032119751 -0.27067059278488159 -0.14936122298240662 -0.073262564837932587 -0.033689741045236588 -0.014872515574097633 -0.0063831745646893978 -0.0026837014593183994 -0.0011106884339824319 -0.00045399941154755652 -0.00018371877376921475 -7.3730567237362266e-05 -2.9384291337919421e-05 -1.1641405762929935e-05 -4.5885362851549871e-06 -1.8005634956352878e-06 -7.0378973759943619e-07 -2.7413975089984888e-07 -1.0645318582191976e-07 -4.122309249510181e-08 -1.5923385277005764e-08 -6.1368328196920174e-09 -2.3602335641470518e-09 -9.0603280433754207e-10 -3.471987974901225e-10 -1.3283640853956058e-10 -5.0747316071575455e-11 -1.9360334516105304e-11 -7.3766357605586919e-12 -2.8072891233854591e-12 -1.0671687058344537e-12 -4.0525363013271809e-13 -1.5374336461045079e-13 -5.8272932648966574e-14 -2.206792725173521e-14 -8.3502896573240185e-15 -3.1572303958374423e-15 -1.192871523299666e-15 -4.5038112940094517e-16 -1.699343306816689e-16 -6.4078234365689933e-17 -2.4148019279880996e-17 -9.095073346605316e-18 -3.4237017961279004e-18 -1.2881348671140216e-18 -4.8440896082993503e-19 -1.8207810463107454e-19 -6.8407959442757565e-20 -2.569017156788846e-20 -9.6437611040447661e-21 -3.6186962678628536e-21 -1.357346940624028e-21 -5.0894276378983633e-22 -1.9076220526102576e-22 -7.1477077345829229e-23 -2.6773039821769189e-23 -1.0025130740057213e-23 -3.7527418826161672e-24 -1.4043593713279384e-24 -5.2539147015754201e-25 -1.9650207139502987e-25 -7.3474141096711539e-26 -2.7465588329293218e-26 -1.0264406957471058e-26
 109
 110
 111 a: 1.6487212181091309 1.2840254306793213 1.1331484317779541 1.0644944906234741 1.0317434072494507 1.0157476663589478 1.007843017578125 1.0039138793945313 1.0019550323486328 1.0009770393371582 1.0004884004592896 1.000244140625
 112
 113 // These are common to both W0 and W-1
 114 b: 0.5 0.25 0.125 0.0625 0.03125 0.015625 0.0078125 0.00390625 0.001953125 0.0009765625 0.00048828125 0.000244140625
 115
 116 */
 117
 118 // Creates if no file exists, & uses default overwrite/ ios::replace.
 119 //const char filename[] = // "lambert_w_lookup_table.ipp";  // Write to same folder as generator:
 120 //"I:/modular-boost/libs/math/include/boost/math/special_functions/lambert_w_lookup_table.ipp";
 121 const char filename[] = "lambert_w_lookup_table.ipp";
 122
 123 std::ofstream fout(filename, std::ios::out); // File output stream.
 124
 125 // 128-bit precision type (so that full precision if long double type uses 128-bit).
 126 // typedef cpp_bin_float_quad table_lookup_t; // Output using max_digits10 for 37 decimal digit precision.
 127 // (This is the precision for the tables output as a C++ program,
 128 // not the precision used by the lambert_w.hpp, that defines another typedef lookup_t, default double.
 129
 130 typedef cpp_bin_float_50 table_lookup_t; // Compute tables to 50 decimla digit precision to avoid slight inaccuracy from repeated multiply.
 131
 132 // But Output using max_digits10 for 37 decimal digit precision.
 133
 134 int main()
 135 {
 136   std::cout << "Lambert W table lookup values." << std::endl;
 137   if (!fout.is_open())
 138   {  // File failed to open OK.
 139     std::cerr << "Open file " << filename << " failed!" << std::endl;
 140     std::cerr << "errno " << errno << std::endl;
 141     return -1;
 142   }
 143   try
 144   {
 145     std::cout << "Lambert W test values writing to file " << filename << std::endl;
 146     int output_precision = std::numeric_limits<cpp_bin_float_quad>::max_digits10; // 37 decimal digits.
 147     fout.precision(output_precision);
 148     fout <<
 149       "// Copyright Paul A. Bristow 2017."   "\n"
 150       "// Distributed under the Boost Software License, Version 1.0." "\n"
 151       "// (See accompanying file LICENSE_1_0.txt" "\n"
 152       "// or copy at http://www.boost.org/LICENSE_1_0.txt)" "\n"
 153       "\n"
 154       "// " << filename << "\n\n"
 155       "// A collection of 128-bit precision integral z argument Lambert W values computed using "
 156       << output_precision << " decimal digits precision.\n"
 157       "// C++ floating-point precision is 128-bit long double.\n"
 158       "// Output as "
 159       << std::numeric_limits<table_lookup_t>::max_digits10
 160       << " decimal digits, suffixed L.\n"
 161       "\n"
 162       "// C++ floating-point type is provided by lambert_w.hpp typedef."  "\n"
 163       "// For example: typedef lookup_t double; (or float or long double)"  "\n"
 164
 165       "\n"
 166       "// Written by " << __FILE__ << " " << __TIMESTAMP__ << "\n"
 167       << std::endl;
 168
 169     fout << "// Sizes of arrays of z values for Lambert W[0], W[1] ... W[64]"
 170       "\"n""and W[-1], W[-2] ... W[-64]." << std::endl;
 171
 172     fout << "\nnamespace boost {\nnamespace math {\nnamespace lambert_w_detail {\nnamespace lambert_w_lookup\n{ \n";
 173
 174     BOOST_STATIC_CONSTEXPR std::size_t noof_sqrts = 12;
 175     BOOST_STATIC_CONSTEXPR std::size_t noof_halves = 12;
 176     fout << "BOOST_STATIC_CONSTEXPR std::size_t  noof_sqrts = " << noof_sqrts << ";" << std::endl;
 177     fout << "BOOST_STATIC_CONSTEXPR std::size_t  noof_halves = " << noof_halves << ";" << std::endl; // Common to both branches.
 178
 179     BOOST_STATIC_CONSTEXPR std::size_t noof_w0zs = 65; // F[k] 0 <= k <= 64. f[0] = F[0], f[64] = F[64]
 180     BOOST_STATIC_CONSTEXPR std::size_t noof_w0es = 66; //  noof_w0zs +1 for gratuitous extra power.
 181     BOOST_STATIC_CONSTEXPR std::size_t noof_wm1zs = 64; // G[k] 1 <= k <= 64. (W-1 = 0 would be z = -infinity, so not stored in array) g[0] == G[1], g[63] = G[64]
 182     BOOST_STATIC_CONSTEXPR std::size_t noof_wm1es = 64; //
 183
 184     fout << "BOOST_STATIC_CONSTEXPR std::size_t  noof_w0es = " << noof_w0zs << ";" << std::endl;
 185     fout << "BOOST_STATIC_CONSTEXPR std::size_t  noof_w0zs = " << noof_w0zs << ";" << std::endl;
 186     fout << "BOOST_STATIC_CONSTEXPR std::size_t  noof_wm1es = " << noof_wm1zs << ";" << std::endl;
 187     fout << "BOOST_STATIC_CONSTEXPR std::size_t  noof_wm1zs = " << noof_wm1zs << ";" << std::endl;
 188
 189     // Defining actual lookup table sqrts of e^k, e^-k = 1/e, etc.
 190     table_lookup_t halves[noof_halves]; // 0.5 0.25 0.125, 0.0625 ...  0.0002441, halves of previous elements.
 191     table_lookup_t sqrtw0s[noof_sqrts]; // 0.6065 0.7788, 0.8824 ... 0.9997, sqrt of previous elements.
 192     table_lookup_t sqrtwm1s[noof_sqrts]; // 1.6487, 1.2840 1.1331  ... 1.00024 , sqrt of previous elements.
 193     table_lookup_t w0es[noof_w0es]; // lambert_w[k] for W0 branch. 2.7182, 1, 0.3678, 0.1353, ... 1.6038e-28
 194     table_lookup_t w0zs[noof_w0zs]; // lambert_w[k] for W0 branch. 0. , 2.7182, 14.77, 60.2566 ... 1.445e+29, 3.990e+29.
 195     table_lookup_t wm1es[noof_wm1es];  // lambert_w[k] for W-1 branch. 2.7182 7.38905 20.085 ... 6.235e+27
 196     table_lookup_t wm1zs[noof_wm1zs];  // lambert_w[k] for W-1 branch. -0.3678 ... -1.0264e-26
 197
 198     // e values lambert_w[k] for W-1 branch. 2.7182 1. 0.3678 0.1353 0.04978 ... 4.359e-28 1.603e-28
 199
 200     using boost::math::constants::e;
 201     using boost::math::constants::exp_minus_one;
 202
 203     {  // z values for integral W F[k] and powers for W0 branch.
 204       table_lookup_t ej = 1; //
 205       w0es[0] = e<table_lookup_t>(); // e = 2.7182 exp(-1) - 1/e exp_minus_one = 0.36787944.
 206       w0es[1] = 1; // e^0
 207       w0zs[0] = 0; // F[0] = 0 or W0 branch.
 208       for (int j = 1, jj = 2; jj != noof_w0es; ++jj)
 209       {
 210         w0es[jj] = w0es[j] * exp_minus_one<table_lookup_t>(); // previous * 0.36787944.
 211         ej *= e<table_lookup_t>(); // * 2.7182
 212         w0zs[j] = j * ej; // For W0 branch.
 213         j = jj; // Previous.
 214       } // for
 215     }
 216     // Checks on accuracy of W0 exponents.
 217
 218     // Checks on e power w0es
 219
 220     // w0es[64] =       4.3596100000630809736231248158884615452e-28
 221     // N[e ^ -63, 37] = 4.359610000063080973623124815888459643*10^-28
 222     // So slight loss at last decimal place.
 223
 224     // Checks on accuracy of z for integral W0 w0zs
 225     // w0zs[0] = 0,  z = -infinity expected? but = zero
 226     // w0zs[1] = 2.7182818284590452353602874713526623144
 227     // w0[2] z = 14.778112197861300454460854921150012956
 228     // w0zs[64] = 3.9904954117194348050619127737142022705e+29
 229     // N[productlog(0, 3.9904954117194348050619127737142022705 10^+29), 37]
 230     //   =  63.99999999999999999999999999999999547
 231     //   = 64.0 to 34 decimal digits, so exact. :-)
 232
 233     { // z values for integral powers G[k] and e^-k for W-1 branch.
 234       // Fukushima indexing of G (k-1) differs by 1 from(k).
 235       // G[0] = -infinity, so his first item in array g[0] is -0.3678 which is G[1]
 236       // and last is g[63] = G[64] = 1.026e-26
 237       table_lookup_t e1 = 1. / e<table_lookup_t>(); // 1/e = 0.36787944117144233
 238       table_lookup_t ej = e1;
 239       wm1es[0] = e<table_lookup_t>(); // e = 2.7182
 240       wm1zs[0] = -e1; // -1/e = - 0.3678
 241       for (int j = 0, jj = 1; jj != noof_wm1zs; ++jj)
 242       {
 243         ej *= e1; // * 0.3678..
 244         wm1es[jj] = wm1es[j] * e<table_lookup_t>();
 245         wm1zs[jj] = -(jj + 1) * ej;
 246         j = jj; // Previous.
 247       } // for
 248     }
 249
 250     // Checks on W-1 branch accuracy wm1es by comparing with Wolfram.
 251     // exp powers:
 252     // N[e ^ 1, 37]  2.718281828459045235360287471352662498
 253     // wm1es[0] =   2.7182818284590452353602874713526623144  - close enough.
 254     // N[e ^ 3, 37]         20.08553692318766774092852965458171790
 255     // computed wm1es[2]    2.0085536923187667740928529654581712847e+01L  OK
 256     // e ^ 66 =        4.6071866343312915426773184428060086893349003037096040 * 10^28
 257     // N[e ^ 66, 34] = 4.607186634331291542677318442806009 10^28
 258     // computed        4.6071866343312915426773184428059867859e+28L
 259     // N[e ^ 66, 37] = 4.607186634331291542677318442806008689*10^28
 260     // so suffering some loss of precision by repeated multiplication computation.
 261     // :-(
 262
 263     // Repeat with cpp_bin_float_50 and correct to 37th decimal digit.
 264     //                 4.60718663433129154267731844280600868933490030370929982
 265     // output std::cout.precision(std::numeric_limits<cpp_bin_float_quad>::max_digits10) as 37 decimal digits.
 266     //                 4.6071866343312915426773184428060086893e+28L
 267     // N[e ^ 66, 37] = 4.607186634331291542677318442806008689*10^28
 268     // Agrees exactly for 37th place, so should be read in to nearest representable value.
 269
 270
 271     // Checks W-1 branch z values wm1zs
 272     // W-1[0] = -2.7067056647322538378799898994496883858e-01
 273     // w-1[1] = -1.4936120510359182893802724695018536337e-01
 274
 275     // wm1zs[65] -1.4325445274604020119111357113179868158e-27
 276
 277     // N[productlog(-1, -1.4325445274604020119111357113179868158* 10^-27), 37]
 278     // = -65.99999999999999999999999999999999955
 279     // = -66 accurately, so this is OK.
 280     // z = 66 * e^66 =
 281     // =N[-66*e ^ -66, 37]
 282     //           -1.432544527460402011911135711317986177*10^-27
 283     // wm1zs[65] -1.4325445274604020119111357113179868158e-27
 284     // which agrees well enough to 34 decimal digits.
 285     // last wm1zs[65] = 0 is unused.
 286
 287     // Halves, common to both W0 and W-1.
 288     halves[0] = static_cast<table_lookup_t>(0.5); // Exactly representable.
 289     for (int j = 0; j != noof_sqrts -1; ++j)
 290     {
 291       halves[j+1] = halves[j] / 2;  // Half previous element (/2 will be optimised better?).
 292     } // for j
 293
 294     // W0 sqrts
 295     sqrtw0s[0] = static_cast<table_lookup_t>(0.606530659712633423603799534991180453441918135487186955682L);
 296     for (int j = 0; j != noof_sqrts -1; ++j)
 297     {
 298       sqrtw0s[j+1] = sqrt(sqrtw0s[j]); //  Sqrt of previous element. sqrt(1/e), sqrt(sqrt(1/e)) ...
 299     } // for j
 300
 301     // W-1 sqrts
 302     sqrtwm1s[0] = root_e<table_lookup_t>();
 303     for (int j = 0; j != noof_sqrts -1; ++j)
 304     {
 305       sqrtwm1s[j+1] = sqrt(sqrtwm1s[j]); //  Sqrt of previous element. sqrt(1/e), sqrt(sqrt(1/e)) ...
 306     } // for j
 307
 308     // Output values as C++ arrays,
 309     // using BOOST_STATIC_CONSTEXPR as static and constexpr as possible for platform.
 310     fout << std::noshowpoint; // Do show NOT trailing zeros for halves and sqrts values.
 311
 312     fout <<
 313       "\n" "BOOST_STATIC_CONSTEXPR lookup_t halves[noof_halves] = " //
 314       "\n" "{ // Common to Lambert W0 and W-1 (and exactly representable)." << "\n  ";
 315       for (int i = 0; i != noof_halves; i++)
 316       {
 317         fout << halves[i] << 'L';
 318         if (i != noof_halves - 1)
 319         { // Omit trailing comma on last element.
 320           fout << ", ";
 321         }
 322         else
 323         {
 324           fout << std::endl;
 325         }
 326       }
 327       fout << "}; // halves, 0.5, 0.25, ... 0.000244140625, common to W0 and W-1." << std::endl;
 328
 329       fout <<
 330         "\n" "BOOST_STATIC_CONSTEXPR lookup_t sqrtw0s[noof_sqrts] = " //
 331         "\n" "{  // For Lambert W0 only." << "\n  ";
 332       for (int i = 0; i != noof_sqrts; i++)
 333       {
 334         fout << sqrtw0s[i] << 'L';
 335         if (i != noof_sqrts - 1)
 336         { // Omit trailing comma on last element.
 337           fout << ", ";
 338         }
 339         else
 340         {
 341           fout << std::endl;
 342         }
 343       }
 344       fout << "}; // sqrtw0s" << std::endl;
 345
 346       fout <<
 347         "\n" "BOOST_STATIC_CONSTEXPR lookup_t sqrtwm1s[noof_sqrts] = " //
 348         "\n" "{ // For Lambert W-1 only." << "\n  ";
 349       for (int i = 0; i != noof_sqrts; i++)
 350       {
 351         fout << sqrtwm1s[i] << 'L';
 352         if (i != noof_sqrts - 1)
 353         { // Omit trailing comma on last element.
 354           fout << ", ";
 355         }
 356         else
 357         {
 358           fout << std::endl;
 359         }
 360       }
 361       fout << "}; // sqrtwm1s" << std::endl;
 362
 363       fout << std::scientific // May be needed to avoid very large dddddddddddddddd.ddddddddddddddd output?
 364         << std::showpoint; // Do show trailing zeros for sqrts and halves.
 365
 366       // Two W0 arrays
 367
 368       fout << // W0 e values.
 369         // Fukushima code generates an extra unused power, but it is not output by using noof_w0zs instead of noof_w0es.
 370         "\n" "BOOST_STATIC_CONSTEXPR lookup_t w0es[noof_w0zs] = " //
 371         "\n" "{ // Fukushima e powers array e[0] = 2.718, 1., e[2] = e^-1 = 0.135, e[3] = e^-2 = 0.133 ... e[64] = 4.3596100000630809736231248158884615452e-28." << "\n  ";
 372       for (int i = 0; i != noof_w0zs; i++)
 373       {
 374         fout << w0es[i] << 'L';
 375         if (i != noof_w0es - 1)
 376         { // Omit trailing comma on last element.
 377           fout << ", ";
 378         }
 379         if (i % 4 == 0)
 380         {
 381           fout << "\n  ";
 382         }
 383       }
 384       fout << "\n}; // w0es" << std::endl;
 385
 386       fout << // W0 z values for W[1], .
 387         "\n" "BOOST_STATIC_CONSTEXPR lookup_t w0zs[noof_w0zs] = " //
 388         "\n" "{ // z values for W[0], W[1], W[2] ... W[64] (Fukushima array Fk). " << "\n  ";
 389       for (int i = 0; i != noof_w0zs; i++)
 390       {
 391         fout << w0zs[i] << 'L';
 392         if (i != noof_w0zs - 1)
 393         { // Omit trailing comma on last element.
 394           fout << ", ";
 395         }
 396         if (i % 4 == 0)
 397         {
 398           fout << "\n  ";
 399         }
 400       }
 401       fout << "\n}; // w0zs" << std::endl;
 402
 403       // Two arrays for w-1
 404
 405       fout << // W-1 e values.
 406         "\n" "BOOST_STATIC_CONSTEXPR lookup_t wm1es[noof_wm1es] = " //
 407         "\n" "{ // Fukushima e array e[0] = e^1 = 2.718, e[1] = e^2 = 7.39 ... e[64] = 4.60718e+28." << "\n  ";
 408       for (int i = 0; i != noof_wm1es; i++)
 409       {
 410         fout << wm1es[i] << 'L';
 411         if (i != noof_wm1es - 1)
 412         { // Omit trailing comma on last element.
 413           fout << ", ";
 414         }
 415         if (i % 4 == 0)
 416         {
 417           fout << "\n  ";
 418         }
 419       }
 420       fout << "\n}; // wm1es" << std::endl;
 421
 422       fout << // Wm1 z values for integral K.
 423         "\n" "BOOST_STATIC_CONSTEXPR lookup_t wm1zs[noof_wm1zs] = " //
 424         "\n" "{ // Fukushima G array of z values for integral K, (Fukushima Gk) g[0] (k = -1) = 1 ... g[64] = -1.0264389699511303e-26." << "\n  ";
 425       for (int i = 0; i != noof_wm1zs; i++)
 426       {
 427         fout << wm1zs[i] << 'L';
 428         if (i != noof_wm1zs - 1)
 429         { // Omit trailing comma on last element.
 430           fout << ", ";
 431         }
 432         if (i % 4 == 0)
 433         { // 4 values per line.
 434           fout << "\n  ";
 435         }
 436       }
 437       fout << "\n}; // wm1zs" << std::endl;
 438
 439       fout << "} // namespace lambert_w_lookup\n} // namespace lambert_w_detail\n} // namespace math\n} // namespace boost" << std::endl;
 440     }
 441     catch (std::exception& ex)
 442     {
 443       std::cout << "Exception " << ex.what() << std::endl;
 444     }
 445     fout.close();
 446     return 0;
 447
 448 } // int main()
 449
 450 /*
 451
 452 Original arrays as output by Veberic/Fukushima code:
 453
 454 w0 branch
 455
 456 W-1 branch
 457
 458 e: 2.7182818284590451 7.3890560989306495 20.085536923187664 54.598150033144229 148.41315910257657 403.42879349273500 1096.6331584284583 2980.9579870417274 8103.0839275753815 22026.465794806707 59874.141715197788 162754.79141900383 442413.39200892020 1202604.2841647759 3269017.3724721079 8886110.5205078647 24154952.753575277 65659969.137330450 178482300.96318710 485165195.40978980 1318815734.4832134 3584912846.1315880 9744803446.2488918 26489122129.843441 72004899337.385788 195729609428.83853 532048240601.79797 1446257064291.4734 3931334297144.0371 10686474581524.447 29048849665247.383 78962960182680.578 214643579785915.75 583461742527454.00 1586013452313428.3 4311231547115188.5 11719142372802592. 31855931757113704. 86593400423993600. 2.3538526683701958e+17 6.3984349353005389e+17 1.7392749415204982e+18 4.7278394682293381e+18 1.2851600114359284e+19 3.4934271057485025e+19 9.4961194206024286e+19 2.5813128861900616e+20 7.0167359120976157e+20 1.9073465724950953e+21 5.1847055285870605e+21 1.4093490824269355e+22 3.8310080007165677e+22 1.0413759433029062e+23 2.8307533032746866e+23 7.6947852651419974e+23 2.0916594960129907e+24 5.6857199993359170e+24 1.5455389355900996e+25 4.2012104037905024e+25 1.1420073898156810e+26 3.1042979357019109e+26 8.4383566687414291e+26 2.2937831594696028e+27 6.2351490808115970e+27
 459
 460 g: -0.36787944117144233 -0.27067056647322540 -0.14936120510359185 -0.073262555554936742 -0.033689734995427351 -0.014872513059998156 -0.0063831737588816162 -0.0026837010232200957 -0.0011106882367801162 -0.00045399929762484866 -0.00018371870869270232 -7.3730548239938541e-05 -2.9384282290753722e-05 -1.1641402067449956e-05 -4.5885348075273889e-06 -1.8005627955081467e-06 -7.0378941219347870e-07 -2.7413963540482742e-07 -1.0645313231320814e-07 -4.1223072448771179e-08 -1.5923376898615014e-08 -6.1368298043116385e-09 -2.3602323152914367e-09 -9.0603229062698418e-10 -3.4719859662410078e-10 -1.3283631472964657e-10 -5.0747278046555293e-11 -1.9360320299432585e-11 -7.3766303773930841e-12 -2.8072868906520550e-12 -1.0671679036256938e-12 -4.0525329757101402e-13 -1.5374324278841227e-13 -5.8272886672428505e-14 -2.2067908660514491e-14 -8.3502821888768594e-15 -3.1572276215253082e-15 -1.1928704609782527e-15 -4.5038074274761624e-16 -1.6993417021166378e-16 -6.4078169762734621e-17 -2.4147993510032983e-17 -9.0950634616416589e-18 -3.4236981860988753e-18 -1.2881333612472291e-18 -4.8440839844747606e-19 -1.8207788854829806e-19 -6.8407875971564987e-20 -2.5690139750481013e-20 -9.6437492398196038e-21 -3.6186918227652047e-21 -1.3573451162272088e-21 -5.0894204288896066e-22 -1.9076194289884390e-22 -7.1476978375412793e-23 -2.6773000149758669e-23 -1.0025115553818592e-23 -3.7527362568743735e-24 -1.4043571811296988e-24 -5.2539064576179218e-25 -1.9650175744554385e-25 -7.3474021582506962e-26 -2.7465543000397468e-26 -1.0264389699511303e-26
 461
 462 a: 1.6487212707001282 1.2840254166877414 1.1331484530668263 1.0644944589178595 1.0317434074991028 1.0157477085866857 1.0078430972064480 1.0039138893383477 1.0019550335910028 1.0009770394924165 1.0004884004786945 1.0002441704297478
 463
 464 b: 0.50000000000000000 0.25000000000000000 0.12500000000000000 0.062500000000000000 0.031250000000000000 0.015625000000000000 0.0078125000000000000 0.0039062500000000000 0.0019531250000000000 0.00097656250000000000 0.00048828125000000000 0.00024414062500000000
 465
 466
 467 */
 468
 469