Road testing Ansus' continuum product formula

[mike3]

Well, I've put together something that uses the same Bernoulli formula but with a divergent summation technique thrown in to see if it could be used to make it converge on something.

I had a thread on the sci.math newsgroup talking about continuum sums (as it is about more general math questions than just "Tetration") and if divergent summation could be used to extend the range of the Bernoulli formula:

http://groups.google.com/group/sci.math/...ee12?hl=en

Using the methods I was told about, I put together this PARI/GP code. For base e and parameters k = 1.4, h = 0.3, I got \( ^{1/2} e \) as approximately 1.646354234 (trunc from series evaluation - residual = 1.6463542340349813076692026966465493046 or series evaluation + residual = 1.6463542341294601977391298318188518176). How does that agree with other methods of tetration? Choosing all the parameters (incl. number of terms in the series) to get good accuracy seems real tricky. Especially for bases less than 1 where complex numbers are used (haven't yet figured out the right settings for the fabled \( 0 < b < e^{-e} \) range, namely b = 0.04). Any ideas on that? Can you squeeze more accuracy out (how about trying b = sqrt(2) and trying to see if it can be gotten to agree w/ the regular iteration to, say, 24 places or not)?

Code:

base = exp(1); /* base of tetration */
pdegree = 65; /* number of terms to use in the series minus 1 */
\ps 65 /* should equal pdegree */
dim = pdegree+2;
DR = matrix(dim, dim, r, c, if(c <= r, 1, 0));
PkPow(k,h) = matrix(dim,dim,r,c,\
                     if(r>=c,\
                          binomial(r-1,c-1)*h^(r-c) * ((r-1)!/(c-1)!)^(k-1)\
                    )  )

rownorm(A) = {
       local(R=A);
       for(r=1,matsize(A)[1],
           R[r,] = A[r,]/sum(c=1,matsize(A)[2],A[r,c]));
       return(R);
}

PkPowSum(k, h) = rownorm(PkPow(k, h)) * DR;

/* Bernoulli number B_n */
bernoulli(n) = if(n==0, 1, (-1)^(n+1) * n * zeta(1 - n));

/* Fill array with sum terms for a coefficient b_k. */
MakeBernTermArray(k, p) = {
                 local(nmax = dim);
                 local(v=vector(nmax));
                 v[1] = 0;
                 for(n=1, nmax-1, v[n+1] = polcoeff(p, n-1)/n*binomial(n,k)*bernoulli(n-k));
                 return(v);
}

MakeBernTermColVector(k, p) = mattranspose(MakeBernTermArray(k, p));

/* Divergent sum based continuum sum */
ContSumPoly(p, k, h) = {
            local(p2 = 0);
            local(PPS = PkPowSum(k, h));
            for(k=0,dim-1,p2=p2+((PPS*MakeBernTermColVector(k, p))[dim]*x^k));
            return(p2 - subst(p2, x, 0));
}

/* Compute exp of a power series */
exppoly(p) = truncate(base^p);

/* Multiply by log(base)^x */
mullogb(p) = { local(multpoly = truncate(log(base)^x)*p, outp = 0);
               for(i=0,pdegree,outp=outp+polcoeff(multpoly,i,x)*x^i);
               return(outp); }

/* Compute definite integral from -1 to x of a power series */
defint(p) = { local(integ=intformal(p)); return(integ - subst(integ, x, -1)); }

/* Do one iteration of the formula */
iterate(p, k, h) = defint(mullogb(exppoly(ContSumPoly(p, k, h))));

/* Normalized iteration */
iternorm(p, k, h) = { local(pi=iterate(p, k, h)); return(pi/polcoeff(pi, 0)); }

/* N iterations on an initial guess */
niters(n, p, k, h) = { local(result=p); for(i=1,n,result=iterate(result,k,h)); return(result); }

/* N normalized iterations on an initial guess */
nnormiters(n, p, k, h) = { local(result=p); for(i=1,n,result=iternorm(result,k,h)); return(result); }

/* Use polynomial */
sumpoly(p, evalpoint) = subst(p, x, evalpoint);

/* Compute residual of tetration */
residual(p) = log(sumpoly(p, 0.5))/log(base) - sumpoly(p, -0.5);

/* Real tetration */
tetcrit(X) = sumpoly(tetpoly, X)
realtet(X) = {
       if(X < -0.5, return(log(realtet(X+1))/log(base)));
       if(X > 0.5, return(base^realtet(X-1)));
       return(tetcrit(X));
}

/* To use: run tetpoly = nnormiters(<number of iterations>, <initial guesstimate for series>, <k>, <h>). For
* the test run I used k = 1.4, h = 0.3, 32 iterations, initial guess = 1. residual(tetpoly) gives
* the residual (how well the function solves F(x+1) = exp(F(x))). Then realtet(X) will give the approximation
* of the tetration ^X b using the computed series.
*/