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