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 main()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 183