libs/math/tools/lambert_w_high_reference_values.cpp

   1 // \modular-boost\libs\math\test\lambert_w_high_reference_values.cpp
   2
   3 // Copyright Paul A. Bristow 2017.
   4 // Distributed under the Boost Software License, Version 1.0.
   5 // (See accompanying file LICENSE_1_0.txt
   6 // or copy at http://www.boost.org/LICENSE_1_0.txt)
   7
   8 // Write a C++ file \lambert_w_mp_hi_values.ipp
   9 // containing arrays of z arguments and 100 decimal digit precision lambert_w0(z) reference values.
  10 // These can be used in tests of precision of less-precise types like
  11 // built-in float, double, long double and quad and cpp_dec_float_50.
  12
  13 // These cover the range from 0.5 to (std::numeric_limits<>::max)();
  14 // The Fukushima algorithm changes from a series function for all z > 0.5.
  15
  16 // Multiprecision types:
  17 //#include <boost/multiprecision/cpp_bin_float.hpp>
  18 #include <boost/multiprecision/cpp_dec_float.hpp> // boost::multiprecision::cpp_dec_float_100
  19 using  boost::multiprecision::cpp_dec_float_100;
  20
  21 #include <boost/math/special_functions/lambert_w.hpp> //
  22 using boost::math::lambert_w0;
  23 using boost::math::lambert_wm1;
  24
  25 #include <iostream>
  26 // using std::cout; using std::std::endl; using std::ios; using std::std::cerr;
  27 #include <iomanip>
  28 using std::setprecision;
  29 using std::showpoint;
  30 #include <fstream>
  31 using std::ofstream;
  32 #include <cassert>
  33 #include <cfloat> // for DBL_EPS etc
  34 #include <limits> // for numeric limits.
  35 //#include <ctype>
  36 #include <string>
  37
  38 const long double eps = std::numeric_limits<long double>::epsilon();
  39
  40 // Creates if no file exists, & uses default overwrite/ ios::replace.
  41 const char filename[] = "lambert_w_high_reference_values.ipp"; //
  42 std::ofstream fout(filename, std::ios::out);
  43
  44 typedef cpp_dec_float_100 RealType;  // 100 decimal digits for the value fed to macro BOOST_MATH_TEST_VALUE.
  45 // Could also use cpp_dec_float_50 or cpp_bin_float types.
  46
  47 const int no_of_tests = 450; // 500 overflows float.
  48
  49 static const float min_z = 0.5F; // for element[0]
  50
  51 int main()
  52 {  // Make C++ file containing Lambert W test values.
  53   std::cout << filename << " ";
  54   std::cout << std::endl;
  55   std::cout << "Lambert W0 decimal digit precision values for high z argument values." << std::endl;
  56
  57   if (!fout.is_open())
  58   {  // File failed to open OK.
  59     std::cerr << "Open file " << filename << " failed!" << std::endl;
  60     std::cerr << "errno " << errno << std::endl;
  61     return -1;
  62   }
  63   try
  64   {
  65     int output_precision = std::numeric_limits<RealType>::digits10;
  66     // cpp_dec_float_100 is ample precision and
  67     // has several extra bits internally so max_digits10 are not needed.
  68     fout.precision(output_precision);
  69     fout << std::showpoint << std::endl; // Do show trailing zeros.
  70
  71     // Intro for RealType values.
  72     std::cout << "Lambert W test values written to file " << filename << std::endl;
  73     fout <<
  74       "\n"
  75       "// A collection of big Lambert W test values computed using "
  76       << output_precision << " decimal digits precision.\n"
  77       "// C++ floating-point type is " << "RealType."  "\n"
  78       "\n"
  79       "// Written by " << __FILE__ << " " << __TIMESTAMP__ << "\n"
  80
  81       "\n"
  82       "// Copyright Paul A. Bristow 2017."   "\n"
  83       "// Distributed under the Boost Software License, Version 1.0." "\n"
  84       "// (See accompanying file LICENSE_1_0.txt" "\n"
  85       "// or copy at http://www.boost.org/LICENSE_1_0.txt)" "\n"
  86       << std::endl;
  87
  88     fout << "// Size of arrays of arguments z and Lambert W" << std::endl;
  89     fout << "static const unsigned int noof_tests = " << no_of_tests << ";" << std::endl;
  90
  91     // Declare arrays of z and Lambert W.
  92     fout << "\n// Declare arrays of arguments z and Lambert W(z)" << std::endl;
  93     fout <<
  94       "\n"
  95       "template <typename RealType>""\n"
  96       "static RealType zs[" << no_of_tests << "];"
  97       << std::endl;
  98
  99     fout <<
 100       "\n"
 101       "template <typename RealType>""\n"
 102       "static RealType ws[" << no_of_tests << "];"
 103       << std::endl;
 104
 105     fout << "// The values are defined using the macro BOOST_MATH_TEST_VALUE to ensure\n"
 106       "// that both built-in and multiprecision types are correctly initialiased with full precision.\n"
 107       "// built-in types like float, double require a floating-point literal like 3.14,\n"
 108       "// but multiprecision types require a decimal digit string like \"3.14\".\n"
 109       "// Numerical values are chosen to avoid exactly representable values."
 110       << std::endl;
 111
 112     static const RealType min_z = 0.6; // for element[0]
 113
 114     const RealType max_z = (std::numeric_limits<float>::max)() / 10; // (std::numeric_limits<float>::max)() to make sure is OK for all floating-point types.
 115     // Less a bit as lambert_w0(max) may be inaccurate.
 116     const RealType step_size = 0.5F; // Increment step size.
 117     const RealType step_factor = 2.f; // Multiple factor, typically 2, 5 or 10.
 118     const int step_modulo = 5;
 119
 120     RealType z = min_z;
 121
 122     // Output function to initialize array of arguments z and Lambert W.
 123     fout <<
 124       "\n"
 125       << "template <typename RealType>\n"
 126       "void init_zws()\n"
 127       "{\n";
 128
 129     for (size_t index = 0; (index != no_of_tests); index++)
 130     {
 131       fout
 132         << "  zs<RealType>[" << index << "] = BOOST_MATH_TEST_VALUE(RealType, "
 133         << z // Since start with converting a float may get lots of usefully random digits.
 134         << ");"
 135         << std::endl;
 136
 137       fout
 138         << "  ws<RealType>[" << index << "] = BOOST_MATH_TEST_VALUE(RealType, "
 139         << lambert_w0(z)
 140         << ");"
 141         << std::endl;
 142
 143       if ((index % step_modulo) == 0)
 144       {
 145         z *= step_factor; //
 146       }
 147       z += step_size;
 148       if (z >= max_z)
 149       { // Don't go over max for float.
 150         std::cout << "too big z" << std::endl;
 151         break;
 152       }
 153     } // for index
 154     fout << "};" << std::endl;
 155
 156     fout << "// End of lambert_w_mp_high_values.ipp " << std::endl;
 157   }
 158   catch (std::exception& ex)
 159   {
 160     std::cout << "Exception " << ex.what() << std::endl;
 161   }
 162
 163   fout.close();
 164
 165   std::cout << no_of_tests << " Lambert_w0 values written to files " << __TIMESTAMP__ << std::endl;
 166   return 0;
 167 }  // main
 168
 169
 170 /*
 171 A few spot checks again Wolfram:
 172
 173   zs<RealType>[1] = BOOST_MATH_TEST_VALUE(RealType, 1.6999999999999999555910790149937383830547332763671875);
 174   ws<RealType>[1] = BOOST_MATH_TEST_VALUE(RealType, 0.7796011225311008662356536916883580556792500749037209859530390902424444585607630246126725241921761054);
 175   Wolfram                                            0.7796011225311008662356536916883580556792500749037209859530390902424444585607630246126725241921761054
 176
 177   zs<RealType>[99] = BOOST_MATH_TEST_VALUE(RealType, 3250582.599999999976716935634613037109375);
 178   ws<RealType>[99] = BOOST_MATH_TEST_VALUE(RealType, 12.47094339016839065212822905567651460418204106065566910956134121802725695306834966790193342511971825);
 179   Wolfram                                            12.47094339016839065212822905567651460418204106065566910956134121802725695306834966790193342511971825
 180
 181 */
 182