Thread Rating:
  • 0 Vote(s) - 0 Average
  • 1
  • 2
  • 3
  • 4
  • 5
Find all fixed points of exp[b]
#9
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.
Reply


Messages In This Thread
Find all fixed points of exp[b] - by MorgothV8 - 08/30/2014, 01:41 PM
RE: Find all fixed points of exp[b] - by jaydfox - 09/17/2014, 11:10 PM
RE: Find all fixed points of exp[b] - by jaydfox - 09/26/2014, 04:47 PM
RE: Find all fixed points of exp[b] - by mike3 - 09/29/2014, 01:49 AM
RE: Find all fixed points of exp[b] - by mike3 - 09/29/2014, 02:41 PM
RE: Find all fixed points of exp[b] - by mike3 - 10/07/2014, 01:32 AM
RE: Find all fixed points of exp[b] - by mike3 - 10/07/2014, 01:50 AM

Possibly Related Threads...
Thread Author Replies Views Last Post
  fixed point formula sheldonison 6 9,912 05/23/2015, 04:32 AM
Last Post: mike3
  Attempt to find a limit point but each step needs doubling the precision... Gottfried 15 16,860 11/09/2014, 10:25 PM
Last Post: tommy1729
  sexp(strip) is winding around the fixed points Kouznetsov 8 12,182 06/29/2009, 10:05 AM
Last Post: bo198214
  An error estimate for fixed point computation of b^x bo198214 0 2,446 05/31/2008, 04:11 PM
Last Post: bo198214



Users browsing this thread: 1 Guest(s)