// Copyright Paul A. Bristow 2016, 2017.

// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or
//  copy at http ://www.boost.org/LICENSE_1_0.txt).

// Build and run a simple examples of Lambert W function.

// Some macros that will show some(or much) diagnostic values if #defined.
//#define-able macros
//#define BOOST_MATH_INSTRUMENT_LAMBERT_W0 // W0 branch diagnostics.
//#define BOOST_MATH_INSTRUMENT_LAMBERT_Wm1 // W1 branch diagnostics.
//#define BOOST_MATH_INSTRUMENT_LAMBERT_W_HALLEY // Halley refinement diagnostics.
//#define BOOST_MATH_INSTRUMENT_LAMBERT_W_SCHROEDER // Schroeder refinement diagnostics.
//#define BOOST_MATH_INSTRUMENT_LAMBERT_W_TERMS // Number of terms used for near-singularity series.
//#define BOOST_MATH_INSTRUMENT_LAMBERT_W0_NOT_BUILTIN // higher than built-in precision types approximation and refinement.
//#define BOOST_MATH_INSTRUMENT_LAMBERT_W_SINGULARITY_SERIES // Show evaluation of series near branch singularity.
//#define BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES_ITERATIONS  // Show evaluation of series for small z.
//#define BOOST_MATH_INSTRUMENT_LAMBERT_W0_LOOKUP // Show results from lookup table.

#include <boost/config.hpp> // for BOOST_PLATFORM, BOOST_COMPILER,  BOOST_STDLIB ...
#include <boost/version.hpp>   // for BOOST_MSVC versions.
#include <boost/cstdint.hpp>
#include <boost/exception/exception.hpp>  // boost::exception
#include <boost/math/constants/constants.hpp> // For exp_minus_one == 3.67879441171442321595523770161460867e-01.
#include <boost/math/policies/policy.hpp>

#include <boost/multiprecision/cpp_dec_float.hpp> // boost::multiprecision::cpp_dec_float_50
using boost::multiprecision::cpp_dec_float_50; // 50 decimal digits type.
using boost::multiprecision::cpp_dec_float_100; // 100 decimal digits type.
using boost::multiprecision::backends::cpp_dec_float;
using boost::multiprecision::number;
typedef number<cpp_dec_float<1000> > cpp_dec_float_1000; // 1000 decimal digit types

#include <boost/multiprecision/cpp_bin_float.hpp>
using boost::multiprecision::cpp_bin_float_double; // == double
using boost::multiprecision::cpp_bin_float_double_extended; // 80-bit long double emulation.
using boost::multiprecision::cpp_bin_float_quad; // 128-bit quad precision.

//[lambert_w_simple_examples_includes
#include <boost/math/special_functions/lambert_w.hpp> // For lambert_w function.

using boost::math::lambert_w0;
using boost::math::lambert_wm1;
//] //[/lambert_w_simple_examples_includes]

#include <iostream>
// using std::cout;
// using std::endl;
#include <exception>
#include <stdexcept>
#include <string>
#include <limits>  // For std::numeric_limits.

//! Show value of z to the full possibly-significant max_digits10 precision of type T.
template<typename T>
void show_value(T z)
{
  std::streamsize precision = std::cout.precision(std::numeric_limits<T>::max_digits10);  // Save.
  std::cout.precision(std::numeric_limits<T>::max_digits10); // Show all posssibly significant digits.
  std::ios::fmtflags flags(std::cout.flags());
  std::cout.setf(std::ios_base::showpoint); // Include any trailing zeros.
  std::cout << z;
  // Restore:
  std::cout.precision(precision);
  std::cout.flags(flags);
} // template<typename T> void show_value(T z)

int main()
{
  try
  {
    std::cout << "Lambert W simple examples." << std::endl;

    using boost::math::constants::exp_minus_one; //-1/e, the branch point, a singularity ~= -0.367879.

    // using statements needed for changing error handling policy.
    using boost::math::policies::policy;
    using boost::math::policies::make_policy;
    using boost::math::policies::evaluation_error;
    using boost::math::policies::domain_error;
    using boost::math::policies::overflow_error;
    using boost::math::policies::ignore_error;
    using boost::math::policies::throw_on_error;

  {
//[lambert_w_simple_examples_0
    std::cout.precision(std::numeric_limits<double>::max_digits10);
    // Show all potentially significant decimal digits,
    std::cout << std::showpoint << std::endl;
    // and show significant trailing zeros too.

    double z = 10.;
    double r = lambert_w0(z); // Default policy for double.
    std::cout << "lambert_w0(z) = " << r << std::endl;
    // lambert_w0(z) = 1.7455280027406994
//] [/lambert_w_simple_examples_0]
  }
  {
    // Other floating-point types can be used too, here float.
    // It is convenient to use a function like `show_value`
    // to display all potentially significant decimal digits
    // for the type, including any significant trailing zeros.
    //[lambert_w_simple_examples_1
    float z = 10.F;
    float r;
    r = lambert_w0(z);        // Default policy digits10 = 7, digits2 = 24
    std::cout << "lambert_w0(";
    show_value(z);
    std::cout << ") = ";
    show_value(r);
    std::cout << std::endl;   // lambert_w0(10.0000000) = 1.74552798
   //] //[/lambert_w_simple_examples_1]
  }
   {
     // Example of an integer argument to lambert_w,
     // showing that an integer is correctly promoted to a double.
//[lambert_w_simple_examples_2
     std::cout.precision(std::numeric_limits<double>::max_digits10);
     double r = lambert_w0(10);                           // Pass an int argument "10" that should be promoted to double argument.
     std::cout << "lambert_w0(10) = " << r << std::endl;  // lambert_w0(10) = 1.7455280027406994
     double rp = lambert_w0(10);
     std::cout << "lambert_w0(10) = " << rp << std::endl;
     // lambert_w0(10) = 1.7455280027406994
     auto rr = lambert_w0(10);                            // C++11 needed.
     std::cout << "lambert_w0(10) = " << rr << std::endl;
     // lambert_w0(10) = 1.7455280027406994 too, showing that rr has been promoted to double.
//] //[/lambert_w_simple_examples_2]
   }
   {
     // Using multiprecision types to get much higher precision is painless.
     //[lambert_w_simple_examples_3
     cpp_dec_float_50 z("10");
     // Note construction using a decimal digit string "10",
     // NOT a floating-point double literal 10.
     cpp_dec_float_50 r;
     r = lambert_w0(z);
     std::cout << "lambert_w0("; show_value(z); std::cout << ") = ";
     show_value(r);
     std::cout << std::endl;
     // lambert_w0(10.000000000000000000000000000000000000000000000000000000000000000000000000000000) =
     //   1.7455280027406993830743012648753899115352881290809413313533156980404446940000000
     //] //[/lambert_w_simple_examples_3]
   }
   // Using multiprecision types to get multiprecision precision wrong!
   {
     //[lambert_w_simple_examples_4
     cpp_dec_float_50 z(0.7777777777777777777777777777777777777777777777777777777777777777777777777);
     // Compiler evaluates the nearest double-precision binary representation,
     // from the max_digits10 of the floating_point literal double 0.7777777777777777777777777777...,
     // so any extra digits in the multiprecision type
     // beyond max_digits10 (usually 17) are random and meaningless.
     cpp_dec_float_50 r;
     r = lambert_w0(z);
     std::cout << "lambert_w0(";
     show_value(z);
     std::cout << ") = "; show_value(r);
     std::cout << std::endl;
     // lambert_w0(0.77777777777777779011358916250173933804035186767578125000000000000000000000000000)
     //   = 0.48086152073210493501934682309060873341910109230469724725005039758139532631901386
     //] //[/lambert_w_simple_examples_4]
   }
   {
     //[lambert_w_simple_examples_4a
     cpp_dec_float_50 z(0.9); // Construct from floating_point literal double 0.9.
     cpp_dec_float_50 r;
     r = lambert_w0(0.9);
     std::cout << "lambert_w0(";
     show_value(z);
     std::cout << ") = "; show_value(r);
     std::cout << std::endl;
     // lambert_w0(0.90000000000000002220446049250313080847263336181640625000000000000000000000000000)
     //   = 0.52983296563343440510607251781038939952850341796875000000000000000000000000000000
     std::cout << "lambert_w0(0.9) = " << lambert_w0(static_cast<double>(0.9))
     // lambert_w0(0.9)
     //   = 0.52983296563343441
       << std::endl;
     //] //[/lambert_w_simple_examples_4a]
   }
   {
     // Using multiprecision types to get multiprecision precision right!
     //[lambert_w_simple_examples_4b
     cpp_dec_float_50 z("0.9");     // Construct from decimal digit string.
     cpp_dec_float_50 r;
     r = lambert_w0(z);
     std::cout << "lambert_w0(";
     show_value(z);
     std::cout << ") = "; show_value(r);
     std::cout << std::endl;
     // 0.90000000000000000000000000000000000000000000000000000000000000000000000000000000)
     // = 0.52983296563343441213336643954546304857788132269804249284012528304239956413801252
     //] //[/lambert_w_simple_examples_4b]
   }
   // Getting extreme precision (1000 decimal digits) Lambert W values.
   {
     std::cout.precision(std::numeric_limits<cpp_dec_float_1000>::digits10);
     cpp_dec_float_1000 z("2.0");
     cpp_dec_float_1000 r;
     r = lambert_w0(z);
     std::cout << "lambert_w0(z) = " << r << std::endl;
     // 0.8526055020137254913464724146953174668984533001514035087721073946525150656742630448965773783502494847334503972691804119834761668851953598826198984364998343940330324849743119327028383008883133161249045727544669202220292076639777316648311871183719040610274221013237163543451621208284315007250267190731048119566857455987975973474411544571619699938899354169616378479326962044241495398851839432070255805880208619490399218130868317114428351234208216131218024303904457925834743326836272959669122797896855064630871955955318383064292191644322931561534814178034773896739684452724587331245831001449498844495771266728242975586931792421997636537572767708722190588748148949667744956650966402600446780664924889043543203483210769017254907808218556111831854276511280553252641907484685164978750601216344998778097446525021666473925144772131644151718261199915247932015387685261438125313159125475113124470774926288823525823567568542843625471594347837868505309329628014463491611881381186810879712667681285740515197493390563
     // Wolfram alpha command N[productlog[0, 2.0],1000] gives the identical result:
     // 0.8526055020137254913464724146953174668984533001514035087721073946525150656742630448965773783502494847334503972691804119834761668851953598826198984364998343940330324849743119327028383008883133161249045727544669202220292076639777316648311871183719040610274221013237163543451621208284315007250267190731048119566857455987975973474411544571619699938899354169616378479326962044241495398851839432070255805880208619490399218130868317114428351234208216131218024303904457925834743326836272959669122797896855064630871955955318383064292191644322931561534814178034773896739684452724587331245831001449498844495771266728242975586931792421997636537572767708722190588748148949667744956650966402600446780664924889043543203483210769017254907808218556111831854276511280553252641907484685164978750601216344998778097446525021666473925144772131644151718261199915247932015387685261438125313159125475113124470774926288823525823567568542843625471594347837868505309329628014463491611881381186810879712667681285740515197493390563
   }
   {
//[lambert_w_simple_examples_error_policies
      // Define an error handling policy:
      typedef policy<
        domain_error<throw_on_error>,
        overflow_error<ignore_error> // possibly unwise?
      > my_throw_policy;

      std::cout.precision(std::numeric_limits<double>::max_digits10);
      // Show all potentially significant decimal digits,
      std::cout << std::showpoint << std::endl;
      // and show significant trailing zeros too.
      double z = +1;
      std::cout << "Lambert W (" << z << ") = " << lambert_w0(z) << std::endl;
      // Lambert W (1.0000000000000000) = 0.56714329040978384
      std::cout << "\nLambert W (" << z << ", my_throw_policy()) = "
        << lambert_w0(z, my_throw_policy()) << std::endl;
      // Lambert W (1.0000000000000000, my_throw_policy()) = 0.56714329040978384
    //] //[/lambert_w_simple_example_error_policies]
    }
    {
      // Show error reporting from passing a value to lambert_wm1 that is out of range,
      // (and probably was meant to be passed to lambert_0 instead).
//[lambert_w_simple_examples_out_of_range
      double z = +1.;
      double r = lambert_wm1(z);
      std::cout << "lambert_wm1(+1.) = " << r << std::endl;
 //] [/lambert_w_simple_examples_out_of_range]
     // Error in function boost::math::lambert_wm1<RealType>(<RealType>):
      // Argument z = 1 is out of range (z <= 0) for Lambert W-1 branch! (Try Lambert W0 branch?)
    }
  }
  catch (std::exception& ex)
  {
    std::cout << ex.what() << std::endl;
  }
}  // int main()

   /*

   Output:
//[lambert_w_simple_examples_error_message_1
Error in function boost::math::lambert_wm1<RealType>(<RealType>):
Argument z = 1 is out of range (z <= 0) for Lambert W-1 branch! (Try Lambert W0 branch?)
//] [/lambert_w_simple_examples_error_message_1]

   */