Hmm, coming a bit late... However - this fiddling reminds me nicely of my own struggles, when I started to look at tetration. So I'll do the proposal to introduce matrix-notation, and especially that of Carleman-matrices, into that study. It has also the advantage, that the formulae can then immediately be (numerically approximate) checked by use of the matrix-features in Pari/GP.
What I introduced first was the notation of a "vandermondevector" V(x) of ideally infinite, practically truncated size to size =n (and using that Vector- and matrixfunctions implicitely by the global variable "n") . In Pari/GP
V(x,dim=n) = vector(dim,r,x^(r-1))
Then with a vector A_1, which contains the coefficients a0,a1,... a_{n-1} of some powerseries one can write down easily the approximate formula
f(x) = V(x) * A_1 ~ \\ approximate due to truncation to "n" terms
This shall already indicate your initial definition of the A-coefficients, and of course, A shall be so that
f(x) = b^^x = V(x) * A_1~
If one needs also the powers of f(x) and thus some manipulations on A_1 (as in your post) one does with advantage using a Carlemanmatrix A, which is the concatenation of columnvectors, such that for some function f(x)
V(x) * A = [1,f(x),f(x)^2, f(x)^3,... ] = V(f(x))
and again this is ideally of infinite columns but practically (in Pari/GP) truncated to n columns.
With your idea of some coefficients a_k and now the matrix A instead, A is searched such that
V(x) * A = [1, b^^x , (b^^x)^2 , (b^^ x)^3, ...] = V(b^^x)
(I took the liberty to rename the base variable for the exponentiation to "b" (from "base") and to avoid confusion with the capital letter A and the small letters "a" for its entries) .
Also I have a standardmatrix "Bb" which is a Carlemanmatrix for the exponentiation
V(x) * Bb = V(b^x)
beta = log(b)
Bb = matrix(n,n,r,c, (beta*(c-1))^(r-1)/(r-1)!)
With size n=32 and b=sqrt(2) it is easy to see the map x -> b^x with the visible digits (I do by default 200 dec digits internal float precision and 12 digits to show)
Your convolution-formula is then simply written
V(x) * A * Bb = V(b^^x) * Bb = V(b^(b^^x))
and the initial idea to be used is much concisely expressed in my notation here:
V(x+1)*A = V(x)*A*Bb \\ searched for A, such that this equation holds
You've done even more: you've introduced the binomial-rules, and -in principle- discovered, that the matrix of binomial-coefficients is of importance here. I call this matrix usually P, and here it should be used as upper triangular Pascalmatrix, such that (with its property being a Carlemanmatrix as well!) (in Pari/GP: P=matpascal(n-1) ~ )
V(x) * P = V(x+1)
V(x) * P^h = V(x+h) \\ where this must be implemented programmatically in Pari/GP for fractional h
We arrive now at the conditions of your first post:
Thus your goal is: let us find a carlemanmatrix A such that
V(x)*P*A = V(x)*A*Bb
As long as the involved dot-products are convergent, there is a theorem, that if for a continuous variable x there is a nonzero interval, where V(x)*A = V(x)*B (for infinite size, thus for a power series), then A=B. The existence of this condition is important here, because (only) if A is of the required form, then also without the leading V(x)-vector it must be equal:
P*A = A*Bb
And now we can go on using knowledge from matrix-operations to try to solve this euqation for A. It reminds of an eigenvalue-system, however neither P nor Bb are diagonal here.
I'm not yet stepping further here, I've limited time this weeks, but I'll come back later to this. Just for a short approximation of Sheldon's coefficients:
Fractional iteration can be approximated by fractional powers of Bb:
V(x) * Bb^h = V(exp_b^h(x))
V(1) * Bb^h = V(exp_b^h(1)) = V(b^^h)
and fractional powers of Bb can be generated if Bb is diagonalizable. Now, for the approximation by the finite truncation, say n=16 and size = nxn Pari/GP can do this in a blink:
tmpM=mateigen(Bb)
tmpW=tmpM^-1
tmpD= diag(tmpW * Bb * tmpM )
and check, that indeed
tmpM*matdiagonal(tmpD)*tmpW - Bb \\ approximately zero
Fractional powers of Bb can now be constructed by fractional powers of the entries of tmpD only! For this I have also a standard-procedure:
dpow(D,ex=1)= matdiagonal(vector(#D,r,D[r]^ex))
so
tmpM*dpow(tmpD,h)*tmpW \\ is equal to Bb^h
Because we want
V(1)* B^x \\ as a powerseries in x
I just did with the symbolic variable x and powerseries expansion to 32 terms (by setting \ps 32)
V(1) * tmpM * dpow(tmpD,x) * tmpW[,2]
to get the powerseries, with coefficients very near to Sheldon's coefficients.
Unfortunately, for sizes nxn=32x32 the mateigen-procedure needs high precision (perhaps 800 digits or more) and the 64x64 version I could only compute with precision of 3200 internal dec digits, so this method is only for toying around and for the intuition. But one can see a comparision to Sheldons solution to base b=4 in my "comparision of methods" -(e-) paper which I've linked to in older postings here.
V(1) ~ * (tmpM*(dpow(tmpD,x)*tmpW[,2] ))
%910 = 1.00000000000 + 1.09176762501*x + 0.271483635676*x^2 + 0.212452965751*x^3 + 0.0695402029399*x^4
+ 0.0442918971124*x^5 + 0.0147365956688*x^6 + 0.00866880451890*x^7 + 0.00279641322771*x^8 + 0.00161065399837*x^9
+ O(x^10)
(I'll continue this another day... and wish so far much joy and success on your detective-journey for the composition of power series!)
Gottfried
Gottfried
What I introduced first was the notation of a "vandermondevector" V(x) of ideally infinite, practically truncated size to size =n (and using that Vector- and matrixfunctions implicitely by the global variable "n") . In Pari/GP
V(x,dim=n) = vector(dim,r,x^(r-1))
Then with a vector A_1, which contains the coefficients a0,a1,... a_{n-1} of some powerseries one can write down easily the approximate formula
f(x) = V(x) * A_1 ~ \\ approximate due to truncation to "n" terms
This shall already indicate your initial definition of the A-coefficients, and of course, A shall be so that
f(x) = b^^x = V(x) * A_1~
If one needs also the powers of f(x) and thus some manipulations on A_1 (as in your post) one does with advantage using a Carlemanmatrix A, which is the concatenation of columnvectors, such that for some function f(x)
V(x) * A = [1,f(x),f(x)^2, f(x)^3,... ] = V(f(x))
and again this is ideally of infinite columns but practically (in Pari/GP) truncated to n columns.
With your idea of some coefficients a_k and now the matrix A instead, A is searched such that
V(x) * A = [1, b^^x , (b^^x)^2 , (b^^ x)^3, ...] = V(b^^x)
(I took the liberty to rename the base variable for the exponentiation to "b" (from "base") and to avoid confusion with the capital letter A and the small letters "a" for its entries) .
Also I have a standardmatrix "Bb" which is a Carlemanmatrix for the exponentiation
V(x) * Bb = V(b^x)
beta = log(b)
Bb = matrix(n,n,r,c, (beta*(c-1))^(r-1)/(r-1)!)
With size n=32 and b=sqrt(2) it is easy to see the map x -> b^x with the visible digits (I do by default 200 dec digits internal float precision and 12 digits to show)
Your convolution-formula is then simply written
V(x) * A * Bb = V(b^^x) * Bb = V(b^(b^^x))
and the initial idea to be used is much concisely expressed in my notation here:
V(x+1)*A = V(x)*A*Bb \\ searched for A, such that this equation holds
You've done even more: you've introduced the binomial-rules, and -in principle- discovered, that the matrix of binomial-coefficients is of importance here. I call this matrix usually P, and here it should be used as upper triangular Pascalmatrix, such that (with its property being a Carlemanmatrix as well!) (in Pari/GP: P=matpascal(n-1) ~ )
V(x) * P = V(x+1)
V(x) * P^h = V(x+h) \\ where this must be implemented programmatically in Pari/GP for fractional h
We arrive now at the conditions of your first post:
Thus your goal is: let us find a carlemanmatrix A such that
V(x)*P*A = V(x)*A*Bb
As long as the involved dot-products are convergent, there is a theorem, that if for a continuous variable x there is a nonzero interval, where V(x)*A = V(x)*B (for infinite size, thus for a power series), then A=B. The existence of this condition is important here, because (only) if A is of the required form, then also without the leading V(x)-vector it must be equal:
P*A = A*Bb
And now we can go on using knowledge from matrix-operations to try to solve this euqation for A. It reminds of an eigenvalue-system, however neither P nor Bb are diagonal here.
I'm not yet stepping further here, I've limited time this weeks, but I'll come back later to this. Just for a short approximation of Sheldon's coefficients:
Fractional iteration can be approximated by fractional powers of Bb:
V(x) * Bb^h = V(exp_b^h(x))
V(1) * Bb^h = V(exp_b^h(1)) = V(b^^h)
and fractional powers of Bb can be generated if Bb is diagonalizable. Now, for the approximation by the finite truncation, say n=16 and size = nxn Pari/GP can do this in a blink:
tmpM=mateigen(Bb)
tmpW=tmpM^-1
tmpD= diag(tmpW * Bb * tmpM )
and check, that indeed
tmpM*matdiagonal(tmpD)*tmpW - Bb \\ approximately zero
Fractional powers of Bb can now be constructed by fractional powers of the entries of tmpD only! For this I have also a standard-procedure:
dpow(D,ex=1)= matdiagonal(vector(#D,r,D[r]^ex))
so
tmpM*dpow(tmpD,h)*tmpW \\ is equal to Bb^h
Because we want
V(1)* B^x \\ as a powerseries in x
I just did with the symbolic variable x and powerseries expansion to 32 terms (by setting \ps 32)
V(1) * tmpM * dpow(tmpD,x) * tmpW[,2]
to get the powerseries, with coefficients very near to Sheldon's coefficients.
Unfortunately, for sizes nxn=32x32 the mateigen-procedure needs high precision (perhaps 800 digits or more) and the 64x64 version I could only compute with precision of 3200 internal dec digits, so this method is only for toying around and for the intuition. But one can see a comparision to Sheldons solution to base b=4 in my "comparision of methods" -(e-) paper which I've linked to in older postings here.
V(1) ~ * (tmpM*(dpow(tmpD,x)*tmpW[,2] ))
%910 = 1.00000000000 + 1.09176762501*x + 0.271483635676*x^2 + 0.212452965751*x^3 + 0.0695402029399*x^4
+ 0.0442918971124*x^5 + 0.0147365956688*x^6 + 0.00866880451890*x^7 + 0.00279641322771*x^8 + 0.00161065399837*x^9
+ O(x^10)
(I'll continue this another day... and wish so far much joy and success on your detective-journey for the composition of power series!)
Gottfried
Gottfried
Gottfried Helms, Kassel
