• 0 Vote(s) - 0 Average
• 1
• 2
• 3
• 4
• 5
 Find all fixed points of exp[b] mike3 Long Time Fellow Posts: 368 Threads: 44 Joined: Sep 2009 10/07/2014, 01:32 AM (This post was last modified: 10/07/2014, 02:13 AM by 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. « Next Oldest | Next Newest »

 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 MorgothV8 - 09/17/2014, 04:17 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 MorgothV8 - 09/18/2014, 04:31 AM RE: Find all fixed points of exp[b] - by MorgothV8 - 09/24/2014, 03:53 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 RE: Find all fixed points of exp[b] - by Gottfried - 10/07/2014, 11:00 AM

 Possibly Related Threads... Thread Author Replies Views Last Post Kneser-iteration on n-periodic-points (base say \sqrt(2)) Gottfried 11 5,006 05/05/2021, 04:53 AM Last Post: Gottfried On n-periodic points of the exp() - A discussion with pictures and methods Gottfried 1 3,013 06/10/2020, 09:34 AM Last Post: Gottfried fixed point formula sheldonison 6 19,097 05/23/2015, 04:32 AM Last Post: mike3 Attempt to find a limit point but each step needs doubling the precision... Gottfried 15 34,324 11/09/2014, 10:25 PM Last Post: tommy1729 sexp(strip) is winding around the fixed points Kouznetsov 8 21,164 06/29/2009, 10:05 AM Last Post: bo198214 An error estimate for fixed point computation of b^x bo198214 0 4,188 05/31/2008, 04:11 PM Last Post: bo198214

Users browsing this thread: 1 Guest(s)