Tetration series for integer exponent. Can you find the pattern?

[sheldonison]

Despite Sheldonison kind attempts to explain fatou.gp algorithm, it goes over my head, so I cannot tweak fatou for more decimals, and must treat it like a black box.

^^I found my problem with fatou.gp

because I cranked up precision, and calculation of tetrations took hours, I was saving the results on a file using the function write(path,sexp(z)), so I could reutilize his results multiple times.

But fatou's sexp() function is formatted to display only 15 decimals. Saving to disk with write() only saves 15 decimals, so it destroys all the hours of work.

I was believing that fatou was fundamentally limited to 15 decimals and the algorithm only showed meaningful digits...

Well, I don't really know, so I will ask on fatou's thread.

EDIT: oh, it actually have 15 digits of precision.

A 512x512 matrix gives results accurate to about 15-16 decimal digits for base(e), \( z\mapsto \exp(z) \). The fatou.gp pari-program uses theta mapping to iteratively solve p(z), to dramatically improve convergence.

fatou currently defaults to disabling the theta mapping when you're anywhere near the attracting/repelling boundary, the "shell thron" boundary, where the function is approximately indifferent. Without a theta mapping means precision is limited to about 15-16 decimal digits. Allowing a theta mapping actual works fine unless you're almost exactly on the indifferent boundary, which is rarely the case. So you need some patches/updates, for fatou.gp

I'm waiting a little while more before I update the fatou.gp in the computation thread, before deciding what defaults make the most sense, but here it is for now. fatou.gp update:
  fatou.gp (Size: 46.67 KB / Downloads: 938)

Also, you need something to make a taylor series. Here's some additional code that samples 100 points centered on the real axis, and uses the conjugate of those 100 points in the negative half of the real axis, to accurately interpolate the half(b) function centered at b=3, with radius 1.5 accurate to about 32 decimal digits. It takes about 6 mintues to run. Most of the bases can be calculated in a few seconds, but the one closet to the Shel Thron boundary needs a couple of minutes. The theta mapping allows you to have exponential convergence, even if it does get slow down really badly near the Shel Thron boundary. Without the theta mapping, you only get x^4.5 or so convergence, so you need a lot more iterations and a lot more sample points, and then you still wind up with only 16 decimal digits unless you have an infinite amount of time.

Code:

/* in fatou.gp p 38 gives you about 33-35 decimal digits of usable precision */
\p 38
default(format,"g0.15");
/* ffunc returns the sexp(-0.5) for base z */
ffunc(z) = {local(y,y2,t2);
  quietmode=1;
  gettime();
  setdefaults();
  print1(z" ");
  y2=loop(log(log(z))+1); /* loop is like sexpinit(z) using loop(k) */
  y=sexp(-0.5);
  t2=gettime();
  print(y2[1]" "t2);
  return(y);
}
/* ftaylor generates a taylor series centered at w, radius r, samples/2 */
ftaylor( w,r,samples) = {
  local(rinv,s,t,x1,y,z,tot,t_est,tcrc,halfsamples,wtaylor,terms,s2);
  if (samples==0, samples=200);
  halfsamples=samples/2;
  terms = floor(samples*200/240);
  t_est    = vector (samples,i,0);
  tcrc     = vector (samples,i,0);
  if (r==0,r=1);
  rinv = 1/r;
  wtaylor=0;
  s2=samples/2;
  for(s=1, samples, x1=-1/(samples)+(s/halfsamples); tcrc[s]=exp(Pi*I*x1); );
  for (t=1,s2, t_est[t] = ffunc(w+r*tcrc[t]); );
  for (t=s2+1,samples, t_est[t] = conj(t_est[samples+1-t]));

  for (s=0,terms-1,
    tot=0;
    for (t=1,samples,
      tot=tot+t_est[t];
      t_est[t]=t_est[t]*conj(tcrc[t]);
    );
    tot=tot/samples;
    if (s>=1, tot=tot*(rinv)^s);
    wtaylor=wtaylor+tot*x^s;
  );
  wtaylor=precision(wtaylor,precis);
  return(real(wtaylor)); /* real valued; ignore imaginary part errors */
}
ft=ftaylor(3,1.45,200); /* centered at 3, radius=1.45, 100 sample points */
default(format,"g0.32");
prtpoly(ft,150,"xhalf");  /* print out accurate to 32 decimal digits */