Find all fixed points of exp[b]

[mike3]

This problem -- how to compute the various branches of \( W(z) \) -- is already a solved problem. See here:

http://www.aqit.uwo.ca/~rcorless/frames/...index.html
http://www.aqit.uwo.ca/~rcorless/frames/...ambertW.ps

The following is a simple Pari/GP code snippet I came up with that will compute any branch of \( W(z) \). It is based on the asymptotic from the above paper and others found online and guessed, as well as some old Lambert W code of mine. Hopefully it should work for all values of \( z \):

Code:

/* Computes all branches of the Lambert W function. W(z) is the principal branch W_0(z). Wbr(k, z) is any branch W_k(z). */
stirling1(n, k) = ((2*n-k)!)/((k-1)!) * sum(kk=0,n-k,1/((n+kk)*(n-k-kk)!*(n-k+kk)!) * sum(j=0,kk,((((-1)^j)*(j^(n-k+kk)))/(j!*(kk-j)!))))
c(k, m) = 1/m! * ((-1)^k) * stirling1(k+m, k+1)
brlog(k, z) = log(z) + 2*k*Pi*I;
Wapprox(k, z) = brlog(k, z) - log(brlog(k, z)) + sum(kk=0,5,sum(m=1,5,c(kk, m)*log(brlog(k, z))^m * brlog(k, z)^(-kk-m)))
W0approx(z) = {
   local(r);
   local(e = exp(1));
   local(oneovere = exp(-1));

   /* Asymptotic expansion at z = -1/e */
   if(abs(z + oneovere) < 0.6,
      r = -1 + sqrt(2*e)*sqrt(z + oneovere) - (2*e)/3*(z + oneovere);
     );
   /* Crude fit for modest values of z */
   if((abs(z + oneovere) >= 0.6) && (abs(z + oneovere) < 10),
      r = log(z + oneovere);
     );
   /* Large asymptotic for large values of z */
   if(abs(z + oneovere) >= 10,
      r = log(z) - log(log(z)) + (log(log(z))/log(z));
     );

   return(r);
}

Wbr(k, z) = {
   local(r);
   local(expr);

   if(x == 0, return(0));
   if((k == 0) && (z == -1/exp(1)), return(-1.0));

   /* Create initial guess. */
   if(k == 0,
      r = W0approx(z),
      r = Wapprox(k, z)
     );

   /* Use Newton's method to refine to full precision */
   rold = r + 1;
   while(abs(rold - r) > 10^-(precision(0.01)-1),
         expr = exp(r);
         rold = r;
         r = r - (r*expr - z)/(expr + r*expr);
        );

   /* Do one last iteration for insurance */
   expr = exp(r);
   r = r - (r*expr - z)/(expr + r*expr);

   /* Done! */
   return(r);
}

W(z) = Wbr(0, z);