# include <cmath>
# include <ctime>
# include <iomanip>
# include <iostream>

using namespace std;

# include "lambert_w.hpp"

int main ( );
void lambert_w_test01 ( );
void lambert_w_values ( int &n_data, double &x, double &fx, int &b );
void timestamp ( );

//****************************************************************************80

int main ( )

//****************************************************************************80
//
//  Purpose:
//
//    lambert_w_test() tests lambert_w().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    21 June 2023
//
//  Author:
//
//    John Burkardt.
//
{
  timestamp ( );
  cout << "\n";
  cout << "lambert_w_test():\n";
  cout << "  C++ version\n";
  cout << "  Test lambert_w().\n";

  lambert_w_test01 ( );
/*
  Terminate.
*/
  cout << "\n";
  cout << "lambert_w_test():\n";
  cout << "  Normal end of execution.\n";
  cout << "\n";
  timestamp ( );

  return 0;
}
//****************************************************************************80

void lambert_w_test01 ( )

//****************************************************************************80
//
//  Purpose:
//
//    lambert_w_test01() compares lambert_w() to stored values, and a built in function
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    21 June 2023
//
//  Author:
//
//    John Burkardt.
//
{
  int b;
  double fx1;
  double fx2;
  int n_data;
  double x;

  cout << "\n";
  cout << "lambert_w_values_test():\n";
  cout << "  Compare a stored value to the lambert_w(X) function.\n";
  cout << "\n";
  cout << "                X            Branch           W(X)    lambert_w(X)\n";
  cout << "\n";

  n_data = 0;

  for ( ; ; )
  {
    lambert_w_values ( n_data, x, fx1, b );

    if ( n_data == 0 )
    {
      break;
    }

    fx2 = lambert_w ( x, b, 0 );

    cout << "  " << setw(24) << setprecision ( 16 ) << x
         << "  " << setw(2) << b
         << "  " << setw(24) << setprecision ( 16 ) << fx1
         << "  " << setw(24) << setprecision ( 16 ) << fx2 << "\n";
  }
  return;
}
//****************************************************************************80

void lambert_w_values ( int &n_data, double &x, double &fx, int &b )

//****************************************************************************80
//
//  Purpose:
//
//    lambert_w_values() returns some values of the Lambert W function.
//
//  Discussion:
//
//    The function W(X) is defined implicitly by:
//
//      W(X) * e^W(X) = X
//
//    The function is defined for -1/e <= x.
//
//    There are two branches, joining at -1/e = x.
//    The lower branch extends from -1/e <= x < 0
//    The upper branch extends from -1/e <= x
//
//    The function is also known as the "Omega" function.
//
//    In Mathematica, the function can be evaluated by:
//      W = ProductLog [ X ]
//
//    In MATLAB,
//      W = lambertw ( b, x )
//    where b = -1 for lower branch, 0 for upper branch.
//
//    In Python,
//      W = scipy.special.lambertw ( x, b )
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    20 June 2023
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Brian Hayes,
//    "Why W?",
//    The American Scientist,
//    Volume 93, March-April 2005, pages 104-108.
//
//    Eric Weisstein,
//    "Lambert's W-Function",
//    CRC Concise Encyclopedia of Mathematics,
//    CRC Press, 1998.
//
//    Stephen Wolfram,
//    The Mathematica Book,
//    Fourth Edition,
//    Cambridge University Press, 1999,
//    ISBN: 0-521-64314-7,
//    LC: QA76.95.W65.
//
//  Input:
//
//    int &n_data: the user sets N_DATA to 0 before the first call.
//
//  Output:
//
//    int &n_data: on each call, the routine increments N_DATA by 1, and
//    returns the corresponding data; when there is no more data, the
//    output value of N_DATA will be 0 again.
//
//    double &x: the argument of the function.
//
//    double &fx: the value of the function.
//
//    int b: -1 (lower branch) or 0 (upper branch).
//
{
# define N_MAX 41

  static int branch_vec[N_MAX] = {
    -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 
     0,  0,  0,  0,  0,  0,  0,  0,  0,  0, 
     0,  0,  0,  0,  0,  0,  0,  0,  0,  0, 
     0,  0,  0,  0,  0,  0,  0,  0,  0,  0, 
     0 };

  static double fx_vec[N_MAX] = {
   -4.889720169867429,
   -3.994308347002122,
   -3.439216483280204,
   -3.022313245324657,
   -2.678346990016661,
   -2.376421342062887,
   -2.097349210703492,
   -1.824388309032984,
   -1.531811608389612,
   -1.000000000000000,
   -0.608341284733432,
   -0.471671909743522,
   -0.374493134019498,
   -0.297083462446424,
   -0.231960952986534,
   -0.175356500529299,
   -0.125066982982524,
   -0.079678160511477,
   -0.038221241746799,
    0.0000000000000000E+00,
    0.3517337112491958E+00,
    0.5671432904097839E+00,
    0.7258613577662263E+00,
    0.8526055020137255E+00,
    0.9585863567287029E+00,
    0.1000000000000000E+01,
    0.1049908894964040E+01,
    0.1130289326974136E+01,
    0.1202167873197043E+01,
    0.1267237814307435E+01,
    0.1326724665242200E+01,
    0.1381545379445041E+01,
    0.1432404775898300E+01,
    0.1479856830173851E+01,
    0.1524345204984144E+01,
    0.1566230953782388E+01,
    0.1605811996320178E+01,
    0.1745528002740699E+01,
    0.3385630140290050E+01,
    0.5249602852401596E+01,
    0.1138335808614005E+02 };

  static double x_vec[N_MAX] = {
   -0.036787944117144,
   -0.073575888234288,
   -0.110363832351433,
   -0.147151776468577,
   -0.183939720585721,
   -0.220727664702865,
   -0.257515608820010,
   -0.294303552937154,
   -0.331091497054298,
   -0.367879441171442,
   -0.331091497054298,
   -0.294303552937154,
   -0.257515608820010,
   -0.220727664702865,
   -0.183939720585721,
   -0.147151776468577,
   -0.110363832351433,
   -0.073575888234288,
   -0.036787944117144,
    0.0000000000000000E+00,
    0.5000000000000000E+00,
    0.1000000000000000E+01,
    0.1500000000000000E+01,
    0.2000000000000000E+01,
    0.2500000000000000E+01,
    0.2718281828459045E+01,
    0.3000000000000000E+01,
    0.3500000000000000E+01,
    0.4000000000000000E+01,
    0.4500000000000000E+01,
    0.5000000000000000E+01,
    0.5500000000000000E+01,
    0.6000000000000000E+01,
    0.6500000000000000E+01,
    0.7000000000000000E+01,
    0.7500000000000000E+01,
    0.8000000000000000E+01,
    0.1000000000000000E+02,
    0.1000000000000000E+03,
    0.1000000000000000E+04,
    0.1000000000000000E+07 };

  if ( n_data < 0 )
  {
    n_data = 0;
  }

  n_data = n_data + 1;

  if ( N_MAX < n_data )
  {
    n_data = 0;
    x = 0.0;
    fx = 0.0;
    b = 0;
  }
  else
  {
    x  = x_vec[n_data-1];
    fx = fx_vec[n_data-1];
    b = branch_vec[n_data-1];
  }

  return;
# undef N_MAX
}
//****************************************************************************80

void timestamp ( )

//****************************************************************************80
//
//  Purpose:
//
//    timestamp() prints the current YMDHMS date as a time stamp.
//
//  Example:
//
//    31 May 2001 09:45:54 AM
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    19 March 2018
//
//  Author:
//
//    John Burkardt
//
{
# define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct std::tm *tm_ptr;
  std::time_t now;

  now = std::time ( NULL );
  tm_ptr = std::localtime ( &now );

  std::strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm_ptr );

  std::cout << time_buffer << "\n";

  return;
# undef TIME_SIZE
}