Hi Sheldon,
first to your second question, number 1): yes, the function is real on the real axis because we use just the eigendecomposition of a real matrix - no complex coefficients anywhere.
About number 2) I've no idea yet, perhaps its too late for me tonight...
I'll come back to this tomorrow.
We do not have a true Schröder- and inverse Schröder-function because the diagonal after the diagonalization (the set of eigenvalues) does not contain the consecutive powers of the log of the fixpoint but are more or less arbitrary - depending on the size of the matrix. So I don't see at the moment how to model the sexp() and inverse sexp()-functionality.
Gottfried
Here is code to reproduce the behaviour - no optimization is made, just the basic functionality. If you want to exceed the range of convergence of the "power series" you should add explicite iteration over the integer-part of the heights.
first to your second question, number 1): yes, the function is real on the real axis because we use just the eigendecomposition of a real matrix - no complex coefficients anywhere.
About number 2) I've no idea yet, perhaps its too late for me tonight...
I'll come back to this tomorrow.
We do not have a true Schröder- and inverse Schröder-function because the diagonal after the diagonalization (the set of eigenvalues) does not contain the consecutive powers of the log of the fixpoint but are more or less arbitrary - depending on the size of the matrix. So I don't see at the moment how to model the sexp() and inverse sexp()-functionality.
Gottfried
Here is code to reproduce the behaviour - no optimization is made, just the basic functionality. If you want to exceed the range of convergence of the "power series" you should add explicite iteration over the integer-part of the heights.
Code:
n=24 \\ set size for matrices, note this needs ~ 200 digits real precision,
\\ n=64 needs about 4000 digits precision and 2 hours to be computed!
default(precision,400) \\ when n~ 24
\\ procedures ===========================================================
\\ several global variables/matrices will be created/changed
{init(base,matsize,prec) = local(a);
n=matsize; \\ set size for matrices, note size=24x24 needs ~ 200 digits real precision,
\\ n=64 needs about 4000 digits precision and 2 hours to be computed!
default(precision,prec,1); \\ when n~ 24 then prec ~ 200
b = base;bl=log(b); \\ sets the base for tetration x-> 4^x
\\ bl must have the required precision depending on matrix-size n
\\ thus must be recomputed if size and precision is increased!
B = matrix(n,n,r,c,((c-1)*bl)^(r-1)/(r-1)!) ;
\\ the Carlemanmatrix for x-> b^x
\\ the resulting coefficients for the power-series
\\ are in the second column
\\ Diagonalization such that B = M * D * W (where W=M^-1)
M = mateigen(B);
W = M^-1 ;
D = diag(W * B * M);
POLY= M * matdiagonal(W[,2]) \\ this matrix is a "kernel" which needs only the parameter
\\ for the height to give the associated power series
print("matrices are created");
return(1);}
f(x) = sum(k=0,n-1,x^k * B[1+k,2]) \\ just to show how this works for the
\\ function f(x) itself
\\ then we can the h'th fractional power of Bb approximate by M * D^h * W (= B^h)
make_powerseries(h) = return( POLY * vectorv(n,r, D[r]^h) );
{ tet_polynomial(x,h) = local(pc);
pc = POLY*vectorv(n, r , D[r]^h ) ; \\ creates coefficients for power series
result = sum(k=0,n-1, x^k * pc[1+k]); \\ evaluates the powerseries
return(result);}
\\ note that M and W do not really implement the Schröder-functions because
\\ the entries in D are not consecutive powers of one argument
\\ then you start:=================================================================
init(4,24,200) \\ for matrixsize 24x24 you'll need at least 200 digits precision
x0 = 0.1*I
x05= tet_polynomial( x0 , 0.5 ) \\ halfiterate from x0=0.1*I
x1 = tet_polynomial( x05, 0.5 ) \\ halfiterate from x05 should equal one whole iteration:
x1 - b^x0 \\ check the difference
\\ bigger matrixsize
init(4,32,600) \\ for matrixsize 32x32 you'll need at least 600 digits precision
x0 = 0.1*I
x05= tet_polynomial( x0 , 0.5 ) \\ halfiterate from x0=0.1*I
x1 = tet_polynomial( x05, 0.5 ) \\ halfiterate from x05 should equal one whole iteration:
x1 - b^x0 \\ check the difference
\\ bigger matrixsize
init(4,64,4000) \\ for matrixsize 64x64 you'll need at least 4000 digits precision
\\ this needs 2 hours for computation!
x0 = 0.1*I
x05= tet_polynomial( x0 , 0.5 ) \\ halfiterate from x0=0.1*I
x1 = tet_polynomial( x05, 0.5 ) \\ halfiterate from x05 should equal one whole iteration:
x1 - b^x0 \\ check the difference
Gottfried Helms, Kassel
