Constructing Taylor Series from a Recursive function in Pari-GP

This is a continuation of my questions:

Declaring a functional recursive sequence in Matlab

Is there a more efficient way of nesting logarithms?

But I'll keep this question self contained. I have made a coding project for myself; which is to program a working simple calculator for a tetration function I've constructed. This tetration function is holomorphic, and stated not to be Kneser's solution (as to all the jargon, ignore); long story short, I need to run the numbers; to win over the nay-sayers.

As to this, I have to use Pari-GP; as this is a fantastic language for handling large numbers and algebraic expressions. As we are dealing with tetration (think numbers of the order e^e^e^e^e^e); this language is, of the few that exist, the best for such affairs. It is the favourite when doing iterated exponential computations.

Now, the trouble I am facing is odd. It is not so much that my code doesn't work; it's that it's overflowing because it should over flow (think, we're getting inputs like e^e^e^e^e^e; and no computer can handle it properly). I'll post the first batch of code, before I dive deeper.

The following code works perfectly; and does everything I want. The trouble is with the next batch of code. This produces all the numbers I want.

\\This is the asymptotic solution to tetration. z is the variable, l is the multiplier, and n is the depth of recursion
\\Warning: z with large real part looks like tetration; and therefore overflows very fast. Additionally there are singularities which occur where l*(z-j) = (2k+1)*Pi*I.
\\j,k are integers

beta_function(z,l,n) =
{
    my(out = 0);
    for(i=0,n-1,
        out = exp(out)/(exp(l*(n-i-z)) +1));
    out;
}

\\This is the error between the asymptotic tetration and the tetration. This is pretty much good for 200 digit accuracy if you need.
\\modify the 0.000000001 to a bigger number to make this go faster and receive less precision. When graphing 0.0001 is enough
\\Warning: This will blow up at some points. This is part of the math; these functions have singularities/branch cuts.

tau(z,l,n)={
    if(1/real(beta_function(z,l,n)) <= 0.000000001, //this is where we'll have problems; if I try to grab a taylor series with this condition we error out
        -log(1+exp(-l*z)),
        log(1 + tau(z+1,l,n)/beta_function(z+1,l,n)) - log(1+exp(-l*z))
    )
}

\\This is the sum function. I occasionally modify it; to make better graphs, but the basis is this.

Abl(z,l,n) = {
    beta_function(z,l,n) + tau(z,l,n)
}

Plugging this in, you get the following expressions:

Abl(1,log(2),100)
   realprecision = 28 significant digits (20 digits displayed)
%109 = 0.15201551563214167060
exp(Abl(0,log(2),100))
%110 = 0.15201551563214167060
Abl(1+I,2+0.5*I,100)
%111 = 0.28416643148885326261 + 0.80115283113944703984*I
exp(Abl(0+I,2+0.5*I,100))
%112 = 0.28416643148885326261 + 0.80115283113944703984*I

And so on and so forth; where Abl(z,l,n) = exp(Abl(z-1,l,n)). There's no problem with this code. Absolutely none at all; we can set this to 200 precision and it'll still produce correct results. The graphs behave exactly as the math says they should behave. The problem is, in my construction of tetration (the one we actually want); we have to sort of paste together the solutions of Abl(z,l,n) across the value l. Now, you don't have to worry about any of that at all; but, mathematically, this is what we're doing.

This is the second batch of code; which is designed to "paste together" all these Abl(z,l,n) into one function.

//This is the modified asymptotic solution to the Tetration equation.
beta(z,n) = {
    beta_function(z,1/sqrt(1+z),n);
}

//This is the Tetration function.
Tet(z,n) ={
    if(1/abs(beta_function(z,1/sqrt(1+z),n)) <= 0.00000001,//Again, we see here this if statement; and we can't have this.
        beta_function(z,1/sqrt(1+z),n),
        log(Tet(z+1,n))
    )
}

This code works perfectly for real-values; and for complex values. Some sample values,

Tet(1+I,100)
%113 = 0.12572857262453957030 - 0.96147559586703141524*I
exp(Tet(0+I,100))
%114 = 0.12572857262453957030 - 0.96147559586703141524*I
Tet(0.5,100)
%115 = -0.64593666417664607364
exp(Tet(0.5,100))
%116 = 0.52417133958039107545
Tet(1.5,100)
%117 = 0.52417133958039107545

We can also effectively graph this object on the real-line. Which just looks like the following,

ploth(X=0,4,Tet(X,100))

Now, you may be asking; What's the problem then?

If you try and plot this function in the complex plane, it's doomed to fail. The nested logarithms produce too many singularities near the real line. For imaginary arguments away from the real-line, there's no problem. And I've produced some nice graphs; but the closer you get to the real line; the more it misbehaves and just short circuits. You may be thinking; well then, the math is wrong! But, no, the reason this is happening is because Kneser's tetration is the only tetration that is stable about the principal branch of the logarithm. Since this tetration IS NOT Kneser's tetration, it's inherently unstable about the principal branch of the logarithm. Of course, Pari just chooses the principal branch. So when I do log(log(log(log(log(beta(z+5,100)))))); the math already says this will diverge. But on the real line; it's perfectly adequate. And for values of z with an imaginary argument away from zero, we're fine too.

So, how I want to solve this, is to grab the Taylor series at Tet(1+z,100); which Pari-GP is perfect for. The trouble?

Tet(1+z,100)
  ***   at top-level: Tet(1+z,100)
  ***                 ^------------
  ***   in function Tet: ...unction(z,1/sqrt(1+z),n))<=0.00000001,beta_fun
  ***                                                ^---------------------
  *** _<=_: forbidden comparison t_SER , t_REAL.

The numerical comparison I've done doesn't translate to a comparison between t_SER and t_REAL.

So, my question, at long last: what is an effective strategy at getting the Taylor series of Tet(1+z,100) using only real inputs. The complex inputs near z=0 are erroneous; the real values are not. And if my math is right; we can take the derivatives along the real-line and get the right result. Then, we can construct a Tet_taylor(z,n) which is just the Taylor Series expansion. Which; will most definitely have no errors when trying to graph.

Any help, questions, comments, suggestions--anything, is greatly appreciated! I really need some outside eyes on this.

Thanks so much if you got to the bottom of this post. This one is bugging me.

Regards, James

EDIT:

I should add that a Tet(z+c,100) for some number c is the actual tetration function we want. There is a shifting constant I haven't talked about yet. Nonetheless; this is spurious to the question, and is more a mathematical point.

Solution

So I've successfully answered my question. I haven't programmed in so long; I'm kind of shoddy. But I figured it out after enough coffee. I created 3 new functions, which allow me to grab the Taylor series.

\\This function attempts to find the number of iterations we need.

Tet_GRAB_k(A,n) ={
    my(k=0);
    while( 1/real(beta(A+k,n)) >= 0.0001, k++);
    return(k);
}


\\This function will run and produce the same results as Tet; but it's slower; but it let's us estimate Taylor coefficients. 
\\You have to guess which k to use for whatever accuracy before overflowing; which is what the last function is good for.

Tet_taylor(z,n,k) = {
    my(val = beta(z+k,n));
    for(i=1,k,val = log(val));
    return(val);
}

\\This function produces an array of all the coefficients about a value A.

TAYLOR_SERIES(A,n) = {
    my(ser = vector(40,i,0));
    for(i=1,40, ser[i] = polcoeff(Tet_taylor(A+z,n,Tet_GRAB_k(A,n)),i-1,z));
    return(ser);
}

After running the numbers, I'm confident this works. The Taylor series is converging; albeit rather slowly and slightly less accurately than desired; but this will have to do.

Thanks to anyone who read this. I'm just answering this question for completeness.