#include <math.h>

// to be used for -0.2 <= x <= infinity
double LambertW (double z)
{
  double x;
  double x_new;
  int i;
  double NewtonStep;
  if (z < -exp(-1.0))
    return(1.0 / sin(0.0));
  i = 0;
  if (z <= 2.0)
    x = 0.5 * z;
  else
    x = log(z);
  for (i = 0; i <= 8; i++)
  {
    if (z <= 2.0)
      NewtonStep = x - (x * exp(x) - z) / (1.0 + x) / exp(x);
      // NewtonStep= x-(xexpx-z)/(expx+xexpx);
    else
      NewtonStep = x * (1.0 + log(z / x)) / (1.0 + x);
    x_new = (double)(NewtonStep);
    if ((double)(x) == (double)(x_new))
      break;
    x = x_new;
    i = i + 1;
  }
  return(x);
}