\\ Lambert's W function.
\\ Method: solve the equation w*exp(w)=x by a handcrafted iterative method.
\\ lambertw(x,w): initial value supplied.
\\ lambertw(x): initial value is ln(1+x).  This is OK for x>=0 and may well
\\         work elsewhere.  

\\ We actually have three functions.  The basic engine is lambertw1,
\\   which does a single iteration.  A really sophisticated user would
\\   monitor his own calls to lambertw1, but we provide an automatic
\\   driver lambertw2, which generates an initial value if none is
\\   supplied and iterates to convergence, usually making an unnecessary
\\   evaluation of exp(w) as a check on the final value.  The recommended
\\   highest level routine lambertw calls lambertw2 for a low-precision
\\   initial value that nevertheless is accurate enough so that one call
\\   to lambertw1 should suffice.

{lambertw(x,w)=local(low,high);
  high=precision(0.); low=max(9,ceil(high/7.5)); 
  default(realprecision,low); w=lambertw2(precision(x,low),precision(w,low)); 
  print(precision(w));
  default(realprecision,high); w=lambertw1(x,precision(w,high));
  w}

global LAMBERTW_CONVERGE;
{lambertw2(x,w)=local(err,oldw,olderr);
    if (x==0,return(x));
    if (w==0, w=lplog1(x)); err=100;
    while (1,
      oldw=w; w=lambertw1(x,w); olderr=err; err=abs(w-oldw);
      if (LAMBERTW_CONVERGE | err==0 | err>=olderr, break);
     );
  w
}

\\ improve approximate value w using Taylor series up to N-th derivative; 
\\ N=0 means stop after smallest term, but use at most 7 terms.
{lambertw1(x,w,N)=local(s,a,h,hn,p,d,dn,c,term,oldterm,olds,vary);
   if (x==0, return(x)); s=w; a=w*exp(w); h=x-a; if (h==0, return(w));
   hn=1; p=1; d=w/(1+w)/a; dn=d; c=d; term=1e100;
   vary=N==0; if (vary, N=7);
   for (n=1,N,
     hn=hn*h/n; oldterm=term; term=hn*dn; olds=s; s=s+term;
     LAMBERTW_CONVERGE=s==olds;
     if (LAMBERTW_CONVERGE | (vary & abs(term)>=abs(oldterm)), s=olds; break);
     p=(1-3*n-n*q_)*p+(1+q_)*deriv(p,q_);
     c=c*w/(1+w)^2/a; dn=c*subst(p,q_,w));
   s
}

\\ low precision log(1+x)
{lplog1(x)=if (abs(x)<0.5, x-x^2/2, 
/*else*/ precision(log(precision(1+x,9)),precision(0.)))}