Find all fixed points of exp[b]

[mike3]

That's strange ... can't reproduce the apparent problem. Anyways, this code:

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 = (0.639 + 0.152*I)*log(z + oneovere) + (0.390 - 0.343*I);
     );
   /* Large asymptotic for large values of z */
   if(abs(z + oneovere) >= 10,
      r = log(z) - log(log(z)) + (log(log(z))/log(z));
     );

   /* ensure purely real result for positive and not-too-negative real inputs */
   if((imag(z) == 0) && (real(z) >= -oneovere), r = real(r));

   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 + 1) - ((r + 2)*(r*expr - z))/(2*r + 2));
        );

   /* 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);

fix(n, b) = Wbr(n, -log(b))/(-log(b));

has an improved initial guess for the principal branch. Use "fix(n, b)" to find the nth fixed point of base b.

This method is highly efficient, being able to compute the fixed points to thousands of digits in seconds on my 2.4GHz Core 2 machine. E.g. fix(4, 2) gives

Code:

5.2738486587969421374330315218571433423534146415852330401066294382423556546501676001880571055973904853663679469274813764620000693663488878468878992575133848048931917567327477719613549216850000740670739245408605739020265724724505221922149165357780497139655091703495699587502138683716726663081682755519618348414500732021219035238961163860776131773637452101798940688378851967487213699777902485555983113779509422469571944880078361570690436015178078493425261482043433724684114279337464744511646424440703191308122783279 -
38.327787228780575754183706624252395882735744965637005425264867291957073919505982312610787328357650749246081503608230774425033102729695628681018727886229348965817703601243873934949849312466874132625859591878503716732102631968981649986915683789960123279085687250569540219407639546682348231751332415437790604182399983123646641624765667027493935467428011331004110282142632651966995879343798360074353334590171964725501188490595138754927696831217983525656552968160304620509175861943833510599624314954439880911317940197*I

512 digits of precision almost instantly.