Continuous summation question
mike3
It seems to be possible, via "Faulhaber's formula", to get a
"continuous sum" (generalization of summation to cover real and
complex numbers of terms) of an analytic function Namely, if f(z) =
sum_{n=0...inf} a_n z^n then
sum_{k=0...z-1} f(k) = sum_{n=0...inf} a_n (B_(n+1)(z) - B_(n+1)(0))/(n
+1)
where B_n are Bernoulli polynomials.
However this doesn't seem to converge for functions like f(z) = log(1
+ z) or something (try putting in Mercator's series and see what
happens.). if it did, it should give a power series for log(z!) = log
(gamma(z+1)) expanded about 0 (as sum_{k=0...z-1} log(1+k) = log(z!) =
log(gamma(z+1))), but it doesn't and diverges instead. Is there a way
to extend it to cover such series as well?
Gottfried Helms
simple convergence should only be possible, when the coefficients
of the powerseries decrease with the same characteristic, say some
exponential-series, where the denominators contain factorials.
If this is not the case, as in your example of the marcator-series or
even in powerseries with coefficients in geometric progressions,
then you need certainly a divergent summation method.
In this case Cesaro- or even the ordinary Euler-summation should not be
sufficient, because their range of summability covers only series
with geometric growthrate (which means also powerseries with coefficents
in at most geometric progression, Cesaro is even weaker).
But Borel-summation or the idea of N�rlund means gives stronger
procedures.
Examples:
With some order of N�rlund-means (PkPowSum(order1,order2) in my Pari/GP)
I seem to be able to sum the Mercator-series this way; with the Gp-matrix
(=original Faulhaber/Bernoulli-matrix for summing of like powers;
an adaption of the matrix of bernoulli-polynomials which includes
the above correction for the denominators n in your above formula)
I get with direct summing of the logarithms compared to that
of N�rlund-means using the Bernoulli-numbers-formula
the following results:
a)
sum k=2,3, log(1+k)) (exact: 2.484906649788000...)
N�rlund means (orders (1.5,1.2)): ~ 2.484906649788000
the last few of 64 partial sums
...
2.484906649787997
2.484906649787999
2.484906649787999
2.484906649788000
2.484906649788000
2.484906649788000
2.484906649788000
2.484906649788000
2.484906649788000
b)
sum k=2,13, log(1+k)) (exact: 24.49807400217874...)
N�rlund means (orders (1.6,0.84)): ~ 24.498074002
the last few of 64 partial sums
...
24.49807400203999
24.49807400208613
24.49807400209763
24.49807400210822
24.49807400212399
24.49807400213819
24.49807400214765
c)and for fractional bounds
c1) symbolic: sum k=2.5,3.5 log(1+k)
practically: sum k=2,3,log(1+k+0.5) (exact: 2.756840365271642...)
N�rlund means (orders (1.32,1.1)): ~ 2.75684036527
the last few of 64 partial sums
...
2.756840365266885
2.756840365268246
2.756840365269214
2.756840365269904
2.756840365270396
2.756840365270747
2.756840365270998
2.756840365271179
c2) symbolic: sum k=1,0.5 log(1+k)
practically (???)
N�rlund means (orders (1.32,1.4985)): ~ 0.2846828704
the last few of 64 partial sums
...
0.2846828703763550
0.2846828703996567
0.2846828704145107
0.2846828704268095
0.2846828704413400
0.2846828704532148
0.2846828704535546
Gottfried Helms
mike3
> procedures.
>
> Examples:
>
> With some order of Nörlund-means (PkPowSum(order1,order2) in my Pari/GP)
> I seem to be able to sum the Mercator-series this way; with the Gp-matrix
> (=original Faulhaber/Bernoulli-matrix for summing of like powers;
> an adaption of the matrix of bernoulli-polynomials which includes
> the above correction for the denominators n in your above formula)
>
How does this Norlund method work, exactly? How does one choose the
proper
"orders" to get the "right" answer? And what about series with terms
that
are complex numbers?
Gottfried Helms
>
> How does this Norlund method work, exactly? How does one choose the
> proper
> "orders" to get the "right" answer? And what about series with terms
> that
> are complex numbers?
>
even if the partial sums oscisllate and diverge, they can
a) be "averaged" or "smoothed" by some compositional weighting and
b) still converge to a result, which is meaningful, aka is consistent
with other relations of the series, which are possibly known.
It is a method which can be expressed by matrix-transformations, that's
why I like it (actually I came to it without knowing its name).
1) Assume a sequence to be summed, in a vector format, say
A = columnvector([a,b,c,d,...])
Then premultiplied by a triangular matrix DR
DR = 1 0 0 0 ...
1 1 0 0 ...
1 1 1 0 ...
...
where DR and A are thought as of infinite size (but practically truncated)
S = DR * A
whe have in the rows of S the partial sums of A, and hopefully the partial
sums converge, so all entries of S from a certain row-index r on are near the
limit below a certain epsilon of difference.
2) Now, if A is a sequence with periodically alternating signs and also diverging
in absolute values, the also in S we will have alternating signs and diverging entries.
This is a classical problem ("what is 1-2+4-8+..."), and led to the concept of
averaging the partial sums.
How to average the entries in S? This is possible, if we build the partial sums
of S, but divide the result by the number or involved terms.
Let's define a column-vector Z as column(1,2,3,4,...) up to the truncation size n
of all our matrices and vectors, and dZ its version as diagonal-matrix.
Then
S1 = (Z^-1 * DR) * S
gives the averages of the partials sums (or "partial means") in S, each
comprising the partial sums in S from row index 1 up to row-index r.
It is accepted, that, if this converges, we can take this as value for the
divergent summation-process for A. This is somehow the basic-agreement for
the whole concept of divergent summation.
Actually, this can be repeated to some "order", so the means of the means of
the means ... of the partial sums are also ok, if that converges.
This is the most simple idea of H�lder-means. The number of repetition of computing
the means can be called the order; and by analysis of the involved matrices we can
even define fractional orders by fractional powers of DR and such.
3) But this can also be generalized. Not only Z^-1*DR can be applied to S to get
converging partial means, but nearly any matrix - well: provided some conditions.
For instance the binomial-matrix P and its arbitrary powers can be used, if its
rows are also normalized to sum up to 1. This gives one form of Euler-summation,
where I also call its order the power to which the binomial-matrix was taken.
So if the power of the matrix P used is 1, then the Euler-order is 1.
With this matrix we can sum any sequence with alternating sign in A with at most
geometric growth-rate.
So if A contains [1,-2,4,-8,16,...] we can sum it using
S1 = rownorm(P) * S
= rownorm(P) * DR * A
and get in the ~ 30'th row of S1 8 good decimals, let's say we used "Euler-order 1".
If we take
S1 = rownorm(P^2) * S
= rownorm(P^2) * DR * A
we have the correct result already in row 2 of S1 and we used "Euler-order 2"
In priciple, this is also using H�lder-means because H�lder defined it in some
general manner.
4) H�lder shows, that not only DR or P, but mearly any matrix can be used as
weighting/averaging matrix, given some conditions (I might miss some here,
you find a good overview in the springer online math-cyclopedia)
One condition is:
4.1) all row sums of the matrix Q, where
S1 = Q * S
is 1.
Another condition is: (I discuss lower triangular matrices here)
4.2) Q_{n,n} /rowsum_n) -> 0 for n-> inf
Another
4,3) Q_{r,c}>=0 (I'm not sure about this at the moment, but all my matrices
which I use have only positive entries)
4.4) ... (some more conditions)
Condition 4.2) restricts the form of matrices. Especially, some matrices Q can be
"too strong", Q_{n,n} cannot be too big compared to the complete rowsum.
For instance, Euler-summation of order 24 to a sequence in A with q=-2 is not
meaningful because the entries of a row in Q decrease too fast.
There is a lot more to say, but I think, this is the most important overview.
5) Since Euler-summation (using matrix P) can only sum divergent series with at most
geometric growth (and alternating signs!) it is not sufficient to sum, for instance,
alternating factorials, bernoulli-numbers, or even the powerseries for fractional
iterates in tetration. For the latter even Borel-summation is not sufficient.
For such cases I configured Q-matrices, which are derived from the form of the
binomial P-matrix.
Remember, that the log L of the P-matrix is [1,2,3,4,...] in the first subdiagonal.
If you insert squares, cubes or general powers, and take the matrix-exponential
you get P-matrices of -as I call it-: k-orders.Say L[k] for the k'th powers of these
entries. Actually the matrix-exponentials are factorially similarity scalings of
the P-matrix. Call F= diag([0!,1!,2!,3!,,,,,)
then
P_k = Exp(L(k)) = F^(k-1) * P * F^(-(k-1))
For k=2 we find this matrix for instance known as matrix of Laguerre-polynomials.
For k>1 we see also, that the entries of one row decrease sharply.
Thus this factorial scaling modifies the entries in one row of P in a way, that they
compensate the strong growth even of partial sums of factorials with alternating
signs (or of bernoulli-numbers) which we may have in S. (But: a selected k may be
too strong for slower growing series and the sum will then not be correct)
6) Now we have the order k for the P-matrices and as well its ordinary powers;
this gives a set of two order-parameters
PkPow(k,e) = rownorm( P_k^e )
where the function rownorm takes a matrix and norms it rows such that the sum
of its absolute values is 1 (requirement 4.1)
and usually, for sake of convenience I also include the DR matrix into it
and compute the entries of that matrix directly by the entries-definition
PkPowSum(k,e) = PkPow(k,e)*DR
and can then directly apply to a given set of coefficients in a columnvector A
PkPowSum(k,e) * A
and watch the partial sums, whether k and e are too strong or not sufficient
to produce a converging sequence of partial sums.
For instance, with B containing the first 64 Bernoulli-numbers we can use
PkPowSum(1.4,1.18) * B
to get the partial sums from some 58 up to 64 terms
...
1.64493406395
1.64493406440
1.64493406453
1.64493406509
1.64493406615
1.64493406666
nicely converging to zeta(2) (=1.64493406685...) as expected, although the
Bernoulli-numbers have hypergeometric growthrate and cannot be -for instance-
summed by Cesaro-sums or Euler-sums (note: don't confuse this term "Euler-sums"
with the more complicated form of Euler-Maclaurin-summation).
7) Earlier I called this method "modified Euler-summation" because I derived
it this way via Euler-summation, but later found, that this method of averaging
was in the general application settled by H�lder (and independently by Voronoi,
as I learnt in the Springer-online-encyclopedia)
PkPow can -by its construction- assume not only integer, but any fractional
or even complex order in both parameters. If A contains complex values,
likely PkPow must assume complex orders as well to be able to approximate the
divergent "sum" of A with finitely many terms (and especially in the range of
truncation!) - but this is a tedious task of manually finetuning so far...(perhaps
one can find a simple rule using some growthrate of the angle in the complex plane
of the partial sums, or the like...)
Gottfried
mike3
<snip>
Thanks for the post. However, I'm curious: how do you compute
fractional, real, complex matrix powers as are needed in this
method when the matrix has lots and lots of entries, without needing
extreme amounts of numerical precision?
Also, do you have a code for this method I could see? Could you post
it here?
Gottfried Helms
Well I think, integer powers is just a Pari/gp-command, but the fractional
powers. This wil be a longer post, but... ok.
Let's begin with the DR-matrix.
It is good to think in terms of matrix-Log and exponential, such that
DR^h = DRH = MExp( h* MLog(DR))
which can be done symbolically for the symbol h, if you use a small
dimension, from which you get a discernible pattern.
The symbolic h'th power of DR is for size 5
1 . . . .
h 1 . . .
1/2*h^2+1/2*h h 1 . .
1/6*h^3+1/2*h^2+1/3*h 1/2*h^2+1/2*h h 1 .
1/24*h^4+1/4*h^3+11/24*h^2+1/4*h 1/6*h^3+1/2*h^2+1/3*h 1/2*h^2+1/2*h h 1
What we can see (especially if that matrix has more rows/columns) is that along
the diagonals the entries are equal. So we need only compute the first column.
We can smooth the coefficients if the rows are scaled (multiplied) by the factorials
and the columns by their reciprocals (F = matdiagonal(vector(n,r,(r-1)!))
F * DRH * F^-1 gives for the only relevant first column
1
h
h^2+h
h^3+3*h^2+2*h
h^4+6*h^3+11*h^2+6*h
and the coefficents of h are the coefficients of the stirling matrix 1'st
kind. So they can easily be computed when expressed as factors:
1
h
(h+1)h
(h^2+1)(h+1)h
(h^3+1)(h^2+1)(h+1)h
and also the introduced factorials have to be removed.
This makes an easy procedure:
{ DrPow(h, dim=6) = local(res, f);
res=matid(dim); \\ set ones in the diagonal of the result-matrix
f=h; \\ compute leading entry in first column(=f), second row
for(r=2,dim,
for(c=1,r, res[r+c-1,c]=f);\\ copy it along the whole diagonal
f=f*(r-1+h)/r; \\ compute leading entry in first column(=f) for next row
);
return(res);}
which Pari/GP can numerically compute "on the fly" up to, say, dim=128 which is enough
for all my tests
------------------------------------------------------
For the Euler-summation we need powers of P, which is especially simple.
Again the reference is
PPow(h) = P^h = PH = MExp( h * MLog(P))
where MLog(P) may be stored as constant, since it is simply
MLog(P) = L = subdiagonal(1,[1,2,3,4,...])
The MExp is then simply the exponential-series with h*L as parameter.
But here we can do even easier. Assume we have P as lower triangular
binomial-matrix stored as matrix-constant of a certain size. Also assume the
vandermonde-vector V(x) = column(1,x,x^2,x^3,...) and dV(x) as its
version as diagonalmatrix. Then
PPow(h) = PH = dV(h) * P * dV(1/h)
and thus the entries in PH are just the binomials of P multiplied by
powers of h:
PPow(h,dim=6) = matrix(dim,dim,r,c, P[r,c] * h^(r-c) )
or more general
PPow(h,dim=6) = matrix(dim,dim,r,c, if(r>=c,binomial(r-1,c-1) * h^(r-c)))
which I'd rewrite into recursion to avoid the explicite computation of binomials.
The Pk-Powers are furtherly similarity-scalings of P (or PPow) by powers of
the factorials. Assume
dF(e) = matdiagonal(vector(dim,r,(r-1)!^e)) \\ e'th powers of factorials
then
PkPow(k,h) = dF(k-1) * PPow(h) * dF(-(k-1))
Which can be seen from the exponentiation of the accordingly modified matrix-logarithm
of P
L(k) = subdiagonal(1, [1^k ,2^k, 3^k,...])
= matrix(dim,dim,r,c,if(r-1==c, c^k ))
and
PkPow(k,h) = MExp( h * L(k) )
The entries in PkPow are then simply the entries of PPow scaled by (r-1)!/(c-1)!)^(k-1)
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)
) )
which can also be optimized when computed recursively instead of explicite
computation of binomials, factorials and h'th, k'th powers. Then Pari/Gp
gives this in a blink and one can experiment with some fractional power of
the Bell matrix of tetration, or the vector of Bernoulli-numbers or ...
print( PkPow(1, 1) *DR* B[,1] )
print( PkPow(1, 2.2) *DR* B[,1] )
print( PkPow(1.4 , 1.3)*DR* B[,1] )
(or put the results into a sandbox, using PariTTY)
%box >testpartialsums PkPow(1.4 , 1.4) * DR * Mat(B[,1])
with varying parameters k and h.
----------------------------------------------------
One big problem with this should be explicitely mentioned, which I did
not in the previous post.
While using PPow(h) only, this means the usual Euler-summation, the results
will be either still diverging or they will be correct approximations whatever
order you choose. (The final approximation may not occur in the first (visible)
terms - but it will occur theoretically at some row index - or - it will
never converge)
Using PkPow(k,h) I cannot say such thing in generality. Here too high orders
k may lead to wrong results because the sequence of partial sums is "cut" too
early. Problems of this type are discussed under the term of "Tauberian
theorems" for summability of series and the according summation-methods.
(for which I unfortunately could not yet become an expert... ;-( )
So using PkPow gives *only heuristics* and I'd not make a "breaking news"
from a value found by this method without some crosschecking.
For series with hypergeometric growthrate I'm quite confident, that the
method is valid, however, for the fractional-tetration-series with their
higher growthrate we arrive only at some locally converging partial sums
with a certain well approximation, but I think that eventually the growth of
such series dominates any order of PkPow ("... when it begins to diverge...")
and we finally still need a conceptual strengthening of the PkPow-idea.
Gottfried
Gottfried Helms
the programs... sorry.
Am 18.09.2009 11:48 schrieb Gottfried Helms:
> F * DRH * F^-1 gives for the only relevant first column
>
> 1
> h
> h^2+h
> h^3+3*h^2+2*h
> h^4+6*h^3+11*h^2+6*h
>
> and the coefficents of h are the coefficients of the stirling matrix 1'st
> kind. So they can easily be computed when expressed as factors:
> 1
> h
> (h+1)h
> (h^2+1)(h+1)h
> (h^3+1)(h^2+1)(h+1)h
> and also the introduced factorials have to be removed.
This should be:
1
h
(h+1)h
(h+2)(h+1)h
(h+3)(h+2)(h+1)h
...sorry...
Gottfried
Gottfried Helms
Am 18.09.2009 11:48 schrieb Gottfried Helms:
> which can also be optimized when computed recursively instead of explicite
> computation of binomials, factorials and h'th, k'th powers. Then Pari/Gp
> gives this in a blink and one can experiment with some fractional power of
> the Bell matrix of tetration, or the vector of Bernoulli-numbers or ...
>
> print( PkPow(1, 1) *DR* B[,1] )
> print( PkPow(1, 2.2) *DR* B[,1] )
> print( PkPow(1.4 , 1.3)*DR* B[,1] )
print( rownorm(PkPow(1, 1)) *DR* B[,1] )
print( rownorm(PkPow(1, 2.2)) *DR* B[,1] )
print( rownorm(PkPow(1.4 , 1.3))*DR* B[,1] )
because I'm used to such functions PkPowSum(k,e), where I have the rownorm
and DR already "all-included":
PkPowSum(kord,eord,dim=n) = mrownorm(PkPow(kord,eord,dim))*matrix(dim,dim,r,c,if(r>=c,1,0))
Gottfried
mike3
<snip>
Thanks for all the answers. It does seem to work, but it's a little
tricky trying to find the parameters (k, h) for the method to get the
"right" sum.
Gottfried Helms
Yes, true.
However, there is some rule of thumb.
Eulersum (when k=1) sums geometric series, best, if h~ q, where q is
the absolute value of the quotient of consecutive terms.
So if you check the average quotients, when the consecutive quotients q_k
begin to give a smooth approximately constant sequence, then try
h ~ abs(q)
If the series has hypergeometric characteristic, say it has progression
like factorials with alternating signs, then it seems, k~1.4 to 1.7
suffices generally. If also a component of geometric growth is in the
sequence adapt the Eulerorder h as hinted above.
For summation of series with complex coefficients it seems as if in
tetration the occuring series have coefficients, whose angle in the
complex plane eventually approximate a linear progression. Then I got
best results, if the Eulerorder h is complex and reflects this angle
by something like h=exp(I*phi)*q, where q is the quotient of the absolute
values of consecutive terms, or a slight scaling of this. The idea is,
that the angle of the Euler-order rotates the entries of the complex
series asymptotically to the real axis and then euler-sums these.
But this is all handwoven, only a hint. I've no better technique myself
so far for such series. What we really need is a summation-technique
which can sum a series, whose progression is roughly q^binomial(k,2)
or at least q^binomial(k,2)/k! with the index k
I don't know, whether N�rlund-means can be configured for this. All
my approaches were insufficient so far and we have no reference values
for some divergent cases. For instance I have investigated
s(h) = e^1 + e^h + e^h^2 + e^h^3 + ...
t(h) = e^1 - e^h + e^h^2 - e^h^3 + ...
u(h) = e^(-1) + e^(-h) + e^(-h^2) + e^(-h^3) + ...
and the like for convergent and for divergent cases (h=0,5, h=2, ...)
using reordering of double-sums, decompositions to sums of geometric
series or similar approaches to that of using the bernoulli-numbers
and got some interesting results but with no conclusion so far,
mostly due to the lack of a reference value for at least one example
of a divergent case for the s(h)-type.
So let's see; perhaps there will even join some reader of sci.math
and find a way to a general solution for this...
Gottfried
mike3
> Am 19.09.2009 01:18 schrieb mike3:
>
> > On Sep 18, 4:25 am, Gottfried Helms <he...@uni-kassel.de> wrote:
> > <snip>
>
> > Thanks for all the answers. It does seem to work, but it's a little
> > tricky trying to find the parameters (k, h) for the method to get the
> > "right"sum.
>
> Yes, true.
>
> However, there is some rule of thumb.
> Eulersum (when k=1) sums geometric series, best, if h~ q, where q is
> the absolute value of the quotient of consecutive terms.
> So if you check the average quotients, when the consecutive quotients q_k
> begin to give a smooth approximately constant sequence, then try
> h ~ abs(q)
>
> If the series has hypergeometric characteristic, say it has progression
> like factorials with alternating signs, then it seems, k~1.4 to 1.7
> suffices generally. If also a component of geometric growth is in the
> sequence adapt the Eulerorder h as hinted above.
>
> tetration the occuring series have coefficients, whose angle in the
> complex plane eventually approximate a linear progression. Then I got
> best results, if the Eulerorder h is complex and reflects this angle
> by something like h=exp(I*phi)*q, where q is the quotient of the absolute
> values of consecutive terms, or a slight scaling of this. The idea is,
> that the angle of the Euler-order rotates the entries of the complex
> series asymptotically to the real axis and then euler-sums these.
> But this is all handwoven, only a hint. I've no better technique myself
> so far for such series. What we really need is asummation-technique
> or at least q^binomial(k,2)/k! with the index k
>
> my approaches were insufficient so far and we have no reference values
> for some divergent cases. For instance I have investigated
> s(h) = e^1 + e^h + e^h^2 + e^h^3 + ...
> t(h) = e^1 - e^h + e^h^2 - e^h^3 + ...
> u(h) = e^(-1) + e^(-h) + e^(-h^2) + e^(-h^3) + ...
>
> and the like for convergent and for divergent cases (h=0,5, h=2, ...)
> using reordering of double-sums, decompositions to sums of geometric
> series or similar approaches to that of using the bernoulli-numbers
> and got some interesting results but with no conclusion so far,
> mostly due to the lack of a reference value for at least one example
> of a divergent case for the s(h)-type.
>
> So let's see; perhaps there will even join some reader of sci.math
> and find a way to a general solution for this...
>
> Gottfried
What about Borel summation? You said it doesn't work for the
coefficients of
fractional iteration of the tetration (by which I assume you mean
fractional
"U tetration", i.e. iteration of exp(z) - 1 instead of exp(z)). But
what about
for the Faulhaber formula coefficients?
Gottfried Helms
>
> What about Borel summation? You said it doesn't work for the
> coefficients of
> fractional iteration of the tetration (by which I assume you mean
> fractional
> "U tetration", i.e. iteration of exp(z) - 1 instead of exp(z)). But
> what about
> for the Faulhaber formula coefficients?
until today I have not really implemented Borel-summation, due
to lack of understanding. I've read about it in K.Knopp und
G.H. Hardy but I've always difficulties with that integration-
step and the Laplace-transformation.
What I've understood is that remark about the ability of
summation of hypergeometric series. Also there are orders
of Borel-summation (in a similar sense as that for Euler-summation,
providing the ability to even sum powers of factorials), I think
I've seen this in the springer-online-encyclopedia.
I didn't try Faulhaber/Bernoulli besides of your input yet, so I
don't know more than we worked out so far. (Yes, I meant U-tetration)
do you know that I have a description for the coefficients of the
fractional iteration series for U-tetration as bivariate
polynomials? With a better analysis we may even describe an exact
growth-rate.
in http://go.helms-net.de/math/tetdocs/ContinuousfunctionalIteration.pdf
on page 21 and derived them on the pages before so one can continue that
pattern. Below I append the Pari/GP-routine for the eigendecomposition
of the Ut-matrix (note, that the U-matrix (U-tetration base e) cannot be
used because there are only units in the diagonal, take some other base)
Unfortunately I can't answer more, I'm just leaving for some
days of holiday. I'm back at 5'th of october; I'll see then.
Have a nice day -
Gottfried
Ut_EW(Ut, dim=9999, numfmt=1) = local(tu=Ut[2,2], UEW, UEWi, tt=exp(tu), tuv);
dim=max(1,min(rows(Ut),dim));
tuv=vectorv(dim,r,Ut[r,r]);
UEW=numfmt*matid(dim); \\ if used for symbolic decomposition use numfmt=1 (integer-type)
for(c=1,dim-1,
for(r=c+1,dim,
UEW[r,c]=sum(k=c,r-1,Ut[r,k]*UEW[k,c])/(tuv[c]-tuv[r])
)
);
UEWi=numfmt*matid(dim); \\ the inverse W^-1
for(r=2,dim,
forstep(c=r-1,1,-1,
UEWi[r,c]=sum(k=0,r-1-c,Ut[r-k,c]*UEWi[r,r-k])/(tuv[r]-tuv[c])
)
);
return([[tu,tt,exp(tu/tt)],UEW,tuv,UEWi]);
use:
Ut = dV(log(t))* U; \\ stirling-matrix prefixed with log of base t
UtKenn = Ut_EW(Ut) ;
then
UtKenn[2]*matdiagonal(Ut[3])*Ut[4] = Ut
and Ut[3] is vector of consecutive powers of log(t) which can be
raised to fractional powers for fractional powers of Ut:
UtPow = UtKenn[2]*dV(log(t)^h) * UtKenn[4] ;
-------------------------