gfunc(z)={ local(n); n=real(120*log((z-gtcenter)/gtradius)/(Pi*I)); if (n>0, return(vv[floor(0.51+n)]); , n=n+240; return(vv[floor(0.51+n)]); ); } /* for height=2; nearest singularity for migfunc is at (1/e)^(1/e)~=0.6922 where z^z has f'(z)=0; */ gen_sroot(z,r,height) = { local(n,vb); quietmode=1; vv=vector(240); default(format,"g0.13"); gtcenter=z; gtradius=r; gtheight=height; for (n=1,120, vb=z+r*exp(Pi*I*(n-0.5)/120); sexpinit(vb); vv[n]=sexp(height); vv[241-n]=conj(sexp(conj(height))); print(n" "vb" "vv[n]); ); gt=gtaylor(gtcenter,gtradius); if (imag(height)==0, gt=real(gt)); gm=O(x^200); for (n=1,199,gm=gm+polcoeff(gt,n)*x^n); gm=serreverse(gm); gm=Pol(gm)+gtcenter; gmcenter=polcoeff(gt,0); default(format,"g0.32"); return(0); } ztth(z)=subst(gt,x,(z-gtcenter)); zzztth(z)=z^z^subst(gt,x,(z-gtcenter)); sroot(z)=subst(gm,x,(z-gmcenter)); bfromp(z) = { /* bfromp(z) = { /* ploth(t=0,2,z=bfromp(t);z2=3.3+1.8*exp(t*Pi*I);[real(z),imag(z),real(z2),imag(z2)],1); */ if (z<1,z=2-z;return(conj(bfromp(z)))); exp(exp(kfromp(z)-1)); } helps() = { print ("helps(); superroot help menu"); print ("\\r superroot.gp requires \"\\r fatou.gp\" to be loaded first"); print ("newtonsroot(goal,est,height,nlim); uses Newton's method to interate "); print ("gen_sroot(center,radius,height); samples 240 points; center must be real"); print (" gen_sroot(4,2,1.5);" ); print ("sroot(z); gm(x); gmcenter; suchthat z=sroot(z)^^height"); print ("ztth(z); gt(x); gtcenter; z^^height "); print ("ztth(z); gt(x); gtcenter; z^^height "); print ("zzztth(z); gt(x); gtcenter; z^z^z^^height "); print ("invzzztth(z); Newton's method, and then sexpinit "); } helps(); newtonsroot(goal,est,height,nlim) = { local(d,z,n,ert,y1,y2); n=0; if (nlim==0, nlim=20); ert=10; quietmode=1; default(format,"g0.13"); while ((n1E-20), n++; sexpinit(est+1E-12); y2=sexp(height); sexpinit(est); y1=sexp(height); print(n" "est" "y1); d=(y2-y1)/1E-12; ert=(y1-goal); est=est-ert/d; ert=abs(ert); ); return(est); } invzzztth(goal,est,nlim,nprt) = { local(d,z,n,ert,y,y1,y2); n=0; if (nlim==0, nlim=20); if (est==0, est=3; if (real(goal)<0, est=est+3*I)); ert=10; quietmode=1; default(format,"g0.13"); while ((n1E-20), n++; if (abs(est-gtcenter)>gtradius, est=gtcenter+gtradius*exp(I*arg(est-gtcenter))); /* cya; don't let the function get out of bounds */ z = est+x+O(x^2); y1=z^z^subst(gt,x,(z-gtcenter)); d=polcoeff(y1,1); y=polcoeff(y1,0); ert=(y-goal); est=est-ert/d; ert=abs(ert); if (nprt==0, print(n" "est" "ert)); ); /*n=0; */ /*ert=10; */ /*while ((n1E-20), */ /* n++; */ /* sexpinit(est+1E-12); */ /* y2=sexp(gtheight+2); */ /* sexpinit(est); */ /* y1=sexp(gtheight+2); */ /* d=(y2-y1)/1E-12; */ /* ert=(y1-goal); */ /* est=est-ert/d; */ /* ert=abs(ert); */ /* print(n" "est" "ert); */ /*); */ return(est); }