---

Computer Science Department

arXiv:math/9207221v1 [math.CA] 1 Jul 1992Convolution PolynomialsDonald E. KnuthComputer Science DepartmentStanford, California 94305–2140Abstract. The polynomials that arise as coefficients when a power series is raised to thepower x include many important special cases, which have surprising properties that arenot widely known.

This paper explains how to recognize and use such properties, and itcloses with a general result about approximating such polynomials asymptotically.A family of polynomials F0(x), F1(x), F2(x), . .

. forms a convolution family if Fn(x) has degree ≤nand if the convolution conditionFn(x + y) = Fn(x)F0(y) + Fn−1(x)F1(y) + · · · + F1(x)Fn−1(y) + F0(x)Fn(y)holds for all x and y and for all n ≥0.

Many such families are known, and they appear frequentlyin applications. For example, we can let Fn(x) = xn/n!

; the condition(x + y)nn!=nXk=0xkk!yn−k(n −k)!is equivalent to the binomial theorem for integer exponents. Or we can let Fn(x) be the binomialcoefficientxn; the corresponding identityx + yn=nXk=0xkyn −kis commonly called Vandermonde’s convolution.How special is the convolution condition?Mathematica will readily find all sequences ofpolynomials that work for, say, 0 ≤n ≤4:F[n_,x_]:=Sum[f[n,j]x^j,{j,0,n}]/n!conv[n_]:=LogicalExpand[Series[F[n,x+y],{x,0,n},{y,0,n}]==Series[Sum[F[k,x]F[n-k,y],{k,0,n}],{x,0,n},{y,0,n}]]Solve[Table[conv[n],{n,0,4}],[Flatten[Table[f[i,j],{i,0,4},{j,0,4}]]]]Mathematica replies that the F’s are either identically zero or the coefficients of Fn(x) =fn0 +fn1x + fn2x2 + · · · + fnnxn/n!

satisfyf00 = 1 ,f10 = f20 = f30 = f40 = 0 ,f22 = f 211 ,f32 = 3f11f21 ,f33 = f 311 ,f42 = 4f11f31 + 3f 221 ,f43 = 6f 211f21 ,f44 = f 411 .1

This allows us to choose f11, f21, f31, and f41 freely.Suppose we weaken the requirements by asking only that the convolution condition hold whenx = y. The definition of conv then becomes simplyconv[n_]:=LogicalExpand[Series[F[n,2x],{x,0,n}]==Series[Sum[F[k,x]F[n-k,x],{k,0,n}],{x,0,n}]]and we discover that the same solutions occur.

In other words, the weaker requirements imply thatthe strong requirements are fulfilled as well.In fact, it is not difficult to discover a simple rule that characterizes all “convolution families.”LetF(z) = 1 + F1z + F2z2 + F3z3 + · · ·be any power series with F(0) = 1. Then the polynomialsFn(x) = [zn] F(z)xform a convolution family.

Conversely, every convolution family arises in this way or is identicallyzero. (Here the notation ‘[zn] expr’ stands for what Mathematica calls Coefficient[expr,z,n].)Proof.

Let f(z) = ln F(z) = f1z + f2z2/2! + f3z3/3!

+ · · · . It is easy to verify that the coefficientof zn in F(z)x is indeed a polynomial in x of degree ≤n, because F(z)x = exf(z) = exp(xf1z +xf2z2/2!

+ xf3z3/3! + · · · ) expands to the power seriesXk1,k2,k3,...≥0xk1+k2+k3+···f k11 f k22 f k33 .

. .1!k1 k1!

2!k2 k2! 3!k3 k3!

. .

. zk1+2k2+3k3+··· ;when k1 + 2k2 + 3k3 + · · · = n the coefficient of zn is a polynomial in x with terms of degreek1 + k2 + k3 + · · · ≤n.

This construction produces a convolution family because of the rule forforming coefficients of the product F(z)x+y = F(z)xF(z)y.Conversely, suppose the polynomials Fn(x) form a convolution family. The condition F0(0) =F0(0)2 can hold only if F0(x) = 0 or F0(x) = 1.

In the former case it is easy to prove by inductionthat Fn(x) = 0 for all n. Otherwise, the condition Fn(0) = 2Fn(0) for n > 0 implies that Fn(0) = 0for n > 0. If we equate coefficients of xk on both sides ofFn(2x) = Fn(x)F0(x) + Fn−1(x)F1(x) + · · · + F1(x)Fn−1(x) + F0(x)Fn(x) ,we now find that the coefficient fnk of xk in n!

Fn(x) is forced to have certain values based onthe coefficients of F1(x), . .

. , Fn−1(x), when k > 1, because 2kfnk occurs on the left and 2fnkon the right.The coefficient fn1 can, however, be chosen freely.Any such choice must makeFn(x) = [zn] exp(xf11z + xf21z2/2!

+ xf31z3/3! + · · · ), by induction on n.Examples.

The first example mentioned above, Fn(x) = xn/n!, comes from the power seriesF(z) = ez; the second example, Fn(x) =xn, comes from F(z) = 1 + z. Several other power series2

are also known to have simple coefficients when we raise them to the power x. If F(z) = 1/(1 −z),for instance, we find[zn] (1 −z)−x =−xn(−1)n =x + n −1n.It is convenient to use the notationsxn = x(x −1) .

. .

(x −n + 1) = x!/(x −n)!xn = x(x + 1) . .

. (x + n −1) = Γ(x + n)/Γ(x)for falling factorial powers and rising factorial powers.

Sincexn= xn/n! andx+n−1n= xn/n!,our last two examples have shown that the polynomials xn/n!

and xn/n! form convolution families,corresponding to F(z) = 1 + z and F(z) = 1/(1 −z).

Similarly, the polynomialsFn(x) = x(x −s)(x −2s) . .

.x −(n −1) s)n!form a convolution family corresponding to (1 + sz)1/s when s ̸= 0.The cases F(z) = 1 + z and F(z) = 1/(1 −z) are in fact simply the cases t = 0 and t = 1 of ageneral formula for the binomial power series Bt(z), which satisfiesBt(z) = 1 + z Bt(z)t .When t is any real or complex number, exponentiation of this series is known to yield[zn] Bt(z)x =x + tnnxx + tn = x(x + tn −1) . .

. (x + tn −n + 1)n!

;see, for example, [Graham et al 1989, section 7.5, example 5], where a combinatorial proof is given.The special cases t = 2 and t = −1,B2(z) = 1 −p1 −4z2z= 1 + z + 2z2 + 5z3 + 14z4 + 42z5 + · · · ,B−1(z) = 1 +p1 + 4z2= 1 + z −z2 + 2z3 −5z4 + 14z5 −· · · ,in which the coefficients are the Catalan numbers, arise in numerous applications. For example,B2(z) is the generating function for binary trees, and B1(−z) is the reciprocal of B2(z).

We can getidentities in trigonometry by noting that B2( 12 sin θ)2= sec2(θ/2). Furthermore, if p and q areprobabilities with p + q = 1, it turns out that B2(pq) = 1/ max(p, q).

The case t = 1/2,B1/2(z) = z +√4 + z22!2= 1 + z + z22 + z323 −z527 + 2z7211 −5z9215 + 14z11219−· · · ,is another interesting series in which the Catalan numbers can be seen. The convolution polynomialsin this case are the “central factorials” x(x + n2 −1)n−1/n!

[Riordan 1968, section 6.5], also calledSteffensen polynomials [Roman and Rota 1978, example 6].3

The convolution formula corresponding to Bt(z),x + y + tnnx + yx + y + tn =nXk=0x + tkkxx + tky + t(n −k)n −kyy + t(n −k)is a rather startling generalization of Vandemonde’s convolution; it is an identity for all x, y, t,and n.The limit of Bt(z/t)t as t →∞is another important function T(z)/z; hereT(z) =Xn≥1nn−1n!zn = z + z2 + 3z32+ 8z43+ 125z524+ · · ·is called the tree function because nn−1 is the number of labeled, rooted trees. The tree functionsatisfiesT(z) = zeT (z) ,and we have the corresponding convolution family[zn]T(z)zx= [zn] exT (z) = x(x + n)n−1n!.The related power series1 + zT ′(z) =11 −T(z) =Xn≥0nnznn!= 1 + z + 2z2 + 9z32+ 32z43+ 625z524+ · · ·defines yet another convolution family of importance: We have[zn]11 −T(z)x = tn(x)n!,where tn(x) is called the tree polynomial of order n [Knuth and Pittel 1989].

The coefficients oftn(x) = tn1x + tn2x2 + · · · + tnnxn are integers with combinatorial significance; namely, tnk is thenumber of mappings of an n-element set into itself having exactly k cycles.A similar but simpler sequence arises from the coefficients of powers of ezez:n! [zn] exzez =nXk=0nkkn−kxk .The coefficient of xk is the number of idempotent mappings of an n-element set into itself, havingexactly k cycles [Harris and Schoenfeld 1967].If the reader still isn’t convinced that convolution families are worthy of detailed study, well,there’s not much hope, although another example or two might clinch the argument.

Consider thepower serieseez−1 =X bnznn!= 1 + z1! + 2z22!

+ 5z33! + 15z44!+ 52z55!+ · · · ;4

these coefficients bn are the so-called Bell numbers, the number of ways to partition sets of size ninto subsets. For example, the five partitions that make b3 = 5 are{1, 2, 3} ,{1}{2, 3} ,{1, 2}{3} ,{1, 3}{2} ,{1}{2}{3} .The corresponding convolution family is[zn] e(ez−1)x =n0+n1x +n2x2 + · · · +nnxnn!,where the Stirling numbernkis the number of partitions into exactly k subsets.Need more examples?

If the coefficients of F(z) are arbitrary nonnegative numbers with a finitesum S, then F(z)/S defines a discrete probability distribution, and the convolution polynomialFn(x) is Sx times the probability of obtaining the value n as the sum of x independent randomvariables having that distribution.A derived convolution. Every convolution family {Fn(x)} satisfies another general convolutionformula in addition to the one we began with:(x + y)nXk=0k Fk(x) Fn−k(y) = x n Fn(x + y) .For example, if Fn(x) is the convolution family corresponding to powers of Bt(z), this formula saysthat(x + y)nXk=0kx + tkkxx + tky + t(n −k)n −kyy + t(n −k) = xnx + y + tnnx + yx + y + tn ;it looks messy, but it simplifies to another amazing identity in four parameters,nXk=0x + t(n −k)n −ky + tkkyy + tk =x + y + tnnif we replace n by n + 1, k by n + 1 −k, and x by x −t + 1.

This identity has an interesting historygoing back to Rothe in 1793 see [Gould and Kauck´y 1966].The alternative convolution formula is proved by differentiating the basic identity F(z)x =Pn≥0 Fn(x)zn with respect to z and multiplying by z:xzF ′(z) F(z)x−1 =Xn≥0n Fn(x) zn .Now Pnk=0 k Fk(x) Fn−k(y) is the coefficient of zn in xzF ′(z) F(z)x+y−1, while nFn(x + y) is thecoefficient of zn in (x + y)zF ′(z) F(z)x+y−1. Q.E.D.Convolution and composition.

Once upon a time I was trying to remember the form of ageneral convolution family, so I gave Mathematica the following command:Simplify[Series[(1+Sum[A[k]z^k,{k,4}])^x,{z,0,4}]]5

The result was a surprise. Instead of presenting the coefficient of zn as a polynomial in x, Math-ematica chose another form: The coefficient of z2, for example, was 12A21x(x −1) + A2x.

In thenotation of falling factorial powers, Mathematica’s response took the form1 + A1xz + 12A21x2 + A2xz2 + 16A31x3 + A1A2x2 + A3xz3+ 124A41x4 + 12A21A2x3 +A1A3 + 12A22x2 + A4xz4 + O(z)5 .I wasn’t prepared to work with factorial powers, so I tried another tack:Simplify[Series[Exp[Sum[a[k]z^k,{k,4}]x],{z,0,4}]]This time I got ordinary polynomials in x, but—lo and behold—they were1 + a1xz + 12a21x2 + a2xz2 + 16a31x3 + a1a2x2 + a3x) z3+ 124a41x4 + 12a21a2x3 +a1a3 + 12a22x2 + a4xz4 + O(z)5 .The result was exactly the same as before, but with a’s in place of A’s, and with normal powers inplace of the factorials!So I learned a curious phenomenon: If we take any convolution family and replace each powerxk by xk, we get another convolution family. (By the way, the replacement can be done in Math-ematica by sayingExpand[F[n,x]]/.Power[x,k_]->k!Binomial[x,k];expansion is necessary in case Fn(x) has been factored.

)The proof was not difficult to find, once I psyched out how Mathematica might have come upwith its factorial-based formula: We haveexf(z) = 1 + f(z) x + f(z)22!x2 + f(z)33!x3 + · · · ,and furthermore1 + f(z)x = 1 + f(z) x + f(z)22!x2 + f(z)33!x3 + · · · .Therefore if we start with the convolution family Fn(x) corresponding to F(z) = ef(z), and replaceeach xk by xk, we get the convolution family corresponding to 1 + f(z) = 1 + ln F(z).A similar derivation shows that if we replace xk by the rising factorial power xk instead, weget the convolution family corresponding to 1/1 −f(z)= 1/1 −ln F(z). In particular, if webegin with the family Fn(x) = x(x+n)n−1/n!

corresponding to T(z)/z = eT (z), and if we replace xkby xk to get1n!n−1Xk=0n −1kxk+1 nn−1−k ,this must be [zn]1 −T(z)−x = tn(x)/n!, the tree polynomial.6

Indeed, we can replace each xk by k! Gk(x), where {Gk(x)} is any convolution family whatever!The previous examples, xk and xk, are merely the special cases k!xkand k!x+k−1kcorrespondingto two of the simplest and most basic families we have considered.

In general we get1 + f(z) G1(x) + f(z)22!2! G2(x) + f(z)33!3!

G3(x) + · · · ,which is none other than Gf(z)x = Gln F(z)x.For example, Gk(x) =x+2kkxx+2k = x(x+2k−1)k−1/k! is the family corresponding to B2(z).If we know the family Fn(x) corresponding to ef(z) we can replace xk by x(x + 2k −1)k−1, therebyobtaining the family that corresponds to B2f(z)=1 +p1 −4f(z)/2f(z).Convolution matrices.

I knew that such remarkable facts must have been discovered before,although they were new to me at the time. And indeed, it was not difficult to find them in books,once I knew what to look for.

(Special cases of general theorems are not always easy to recognize,because any particular formula is a special case of infinitely many generalizations, almost all ofwhich are false. )In the special case that each polynomial Fn(x) has degree exactly n, i.e., when f1 ̸= 0, thepolynomials n!

Fn(x) are said to be of binomial type [Mullin and Rota 1970]. An extensive theoryof such polynomial sequences has been developed [Rota et al 1973] [Garsia 1973] [Roman and Rota1978], based on the theory of linear operators, and the reader will find it quite interesting to comparethe instructive treatment in those papers to the related but rather different directions explored inthe present work.

A comprehensive exposition of the operator approach appears in [Roman 1984].Actually, Steffensen had defined a concept called poweroids, many years earlier [Steffensen 1941],and poweroids are almost exactly the same as sequences of binomial type; but Steffensen apparentlydid not realize that his poweroids satisfy the convolution property, which we can readily deduce(with hindsight) from equations (6) and (7) of his paper.Eri Jabotinsky introduced a nice way to understand the phenomena of convolution polynomials,by considering the infinite matrix of coefficients fnk [Jabotinsky 1947]. Let us recapitulate thenotation that was introduced informally above:exf(z) = F(z)x = 1 + F1(x) z + F2(x) z2 + · · · ;Fn(x) = (fn1x + fn2x2 + · · · + fnnxn)/n!

;f(z) = f1z + f2z2/2! + f3z3/3!

+ · · · .Then Jabotinsky’s matrix F = (fnk) is a lower triangular matrix containing the coefficients ofn! Fn(x) in the nth row.

The first few rows aref1f2f 21f33f1f2f 31f44f1f3 + 3f 226f 21 f2f 41 ,7

as we saw earlier. In general,fnk =Xn!1!k1 k1!

2!k2 k2! 3!k3 k3!

. .

. f k11 f k22 f k33 .

. .

,summed over all k1, k2, k3, . .

. ≥0 withk1 + k2 + k3 + · · · = k ,k1 + 2k2 + 3k3 + · · · = n .

(The summation is over all partitions of the integer n into k parts, where kj of the parts are equalto j.) We will call such an array a convolution matrix.If each original coefficient fj is an integer, all entries of the corresponding convolution matrixwill be integers, because the complicated quotient of factorials in the sum is an integer—it is thenumber of ways to partition a set of n elements into k subsets with exactly kj of the subsets havingsize j.

Given the first column we can compute the other columns from left to right and from topto bottom by using the recurrencefnk =n−k+1Xj=1n −1j −1fj f(n−j)(k−1) .This recurrence is based on set partitions on which the element n occurs in a subset of size j:There aren−1j−1ways to choose the other j −1 elements of the subset, and the factor f(n−j)(k−1)corresponds to partitioning the remaining n −j elements into k −1 parts.For example, if each fj = 1, the convolution matrix begins1111311761.11525101These are the numbersnkthat Mathematica calls StirlingS2[n,k]; they arose in our exampleof Bell numbers when f(z) = ez −1. Similarly, if each fj = (j −1)!, the first five rows are11123161161;245035101Mathematica calls these numbers (-1)^(n-k)StirlingS1[n,k].

In this case f(z) = ln1/(1−z),and Fn(z) =x+n−1n. The signed numbers StirlingS1[nk],1−112−31−611−6124−5035−1018

correspond to f(z) = ln(1 + z) and Fn(z) =xn. In general if we replace z by αz and x by βx,the effect is to multiply row n of the matrix by αn and to multiply column k by βk.

Thus whenβ = α−1, the net effect is to multiply fnk by αn−k. Transforming the signs by a factor (−1)n−kcorresponds to changing F(z) to 1/F(−z) and f(z) to −f(−z).

Therefore the matrix that begins1−111−31−17−611−1525−101corresponds to f(z) = 1 −e−z.Let’s look briefly at some of our other examples in matrix form. When F(z) = Bt(z), we havefj = (tj −1)j−1, which is an integer when t is an integer.

In particular, the Catalan case t = 2produces a matrix that begins131209121010718130241650335301.When t = 1/2, we can remain in an all-integer realm by replacing z by 2z and x by x/2. Thenfj = 0 when j is even, while f2j+1 = (−1)j(2j −1)!

!2:101−1010−401;90−1001If we now replace z by iz and x by x/i to eliminate the minus signs, we find that f(z) = arcsin z,because lniz+√1 −z2 = iθ when z = sin θ. Thus we can deduce a closed form for the coefficientsof ex arcsin z = B1/2(2iz)x/(2i):n!

[zn] ex arcsin z = (2i)n−1x x2i + n2 −1. .

. x2i −n2 + 1=(x2(x2 + 22) .

. .x2 + (n −2)2,n even;x(x2 + 12)(x2 + 32) .

. .x2 +(n −2)2,n odd.This remarkable formula is equivalent to the theorem of [Gomes Teixeira 1896].If fj = 21−j when j is odd but fj = 0 when j is even, we get the convolution matrix corre-sponding to e2x sinh(z/2):1011401010111605201010501164091160354010102101401.9

Again we could stay in an all-integer realm if we replaced z by 2z and x by x/2; but the surprisingthing in this case is that the entries in even-numbered rows and columns are all integers before wemake any such replacement. The reason is that the entries satisfy fnk = k2f(n−2)k/4 + f(n−2)(k−2).

(See [Riordan 1968, pages 213–217], where the notation T(n, k) is used for these “central factorialnumbers” fnk. )We can complete our listing of noteworthy examples by setting fj = Pnk=1 nn−k−1nk; thenwe get the coefficients of the tree polynomials:13117911429518115691220305301.The sum of the entries in row n is nn.Composition and iteration.

Jabotinski’s main reason for defining things as he did was hisobservation that the product of convolution matrices is a convolution matrix. Indeed, if F and G arethe convolution matrices corresponding to the functions exf(z) and exg(z) we have the vector/matrixidentitiesexf(z) −1 = (z, z2/2!, z3/3!, .

. . ) F (x, x2, x3, .

. .

)Texg(z) −1 = (z, z2/2!, z3/3!, . .

. ) G (x, x2, x3, .

. .

)TIf we now replace xk in exf(z) by k! Gk(x), as in our earlier discussion, we get(z, z2/2!, z3/3!, .

. . ) FG1(x), 2!

G(x), 3! G3(x), .

. .T=(z, z2/2!, z3/3!, .

. . ) FG (x, x2, x3, .

. .

)T=f(z), f(z)2/2!, f(z)3/3!, . .

.G (x, x2, x3, . .

. )T=exg(f(z)) −1 .Multiplication of convolution matrices corresponds to composition of the functions in the exponent.Why did the function corresponding to FG turn out to be gf(z)instead of fg(z)?

Jabotin-sky, in fact, defined his matrices as the transposes of those given here.The rows of his (up-per triangular) matrices were the power series f(z)k, while the columns were the polynomialsFn(x) = [zn] exf(z); with those conventions the product of his matrices F TGT corresponded tofg(z). (In fact, he defined a considerably more general representation, in which the matrix Fcould be U −1FU for any nonsingular matrix U.) However, when our interest is focussed on thepolynomials n!

Fn(x), as when we study Stirling numbers or tree polynomials or the Stirling poly-nomials to be discussed below, it is more natural to work with lower triangular matrices and toinsert factorial coefficients, as Comtet did [Comtet 1970, section 3.7]. The two conventions are iso-morphic.

Without the factorials, convolution matrices are sometimes called renewal arrays [Rogers1978]. We would get a non-reversed order if we had been accustomed to using postfix notation10

(z)f for functions, as we do for operations such as squaring or taking transposes or factorials; thengf(z)would be(z)fg.Recall that the Stirling numbersnkcorrespond to f(z) = ez −1, and the Stirling numbersnkcorrespond to g(z) = ln1/(1 −z). Therefore if we multiply Stirling’s triangles we get theconvolution matrix121FG =661,2636121150250120201which corresponds to gf(z)= ln1/(2 −ez).

Voila! These convolution polynomials representthe coefficients of (2 −ez)−x.

[Cayley 1859] showed that (2 −ez)−1 is the exponential generatingfunction for the sequence 1, 3, 13, 75, 541, . .

. , which counts preferential arrangements of n objects,i.e., different outcomes of sorting when equality is possible as well as inequality.

The coefficient(fg)nk is the number of preferential arrangements in which the “current minimum” changes k timeswhen we examine the elements one by one in some fixed order. (See [Graham et al 1989, exercise7.44].

)Similarly, the reverse matrix product yields the so-called Lah numbers [Lah 1955],121GF =661;2436121120240120201here fj = j! and the rows represent the coefficients of expxfg(z)= exp(xz+xz2+xz3+· · · ).

In-deed, the convolution polynomials in this case are the generalized Laguerre polynomials L(−1)n(−x),which Mathematica calls LaguerreL[n,-1,-x].These polynomials can also be expressed asLn(−x) −Ln−1(−x); or as LaguerreL[n,-x]-LaguerreL[n-1,-x] if we sayUnprotect[LaguerreL]; LaguerreL[-1,x_]:=0; Protect[LaguerreL]first. The row sums 1, 3, 13, 73, 501, .

. .

of GF enumerate “sets of lists” [Motzkin 1971]; the coeffi-cients are (GF)nk = n! [zn] fg(z)k/k!

=nkn−1k−1(n −k)! [Riordan 1968, exercise 5.7].Since convolution matrices are closed under multiplication, they are also closed under exponen-tiation, i.e., under taking of powers.

The qth power F q of a convolution matrix then correspondsto q-fold iteration of the function ln F = f. Let us denote ff(z)by f [2](z); in general, the qthiterate f [q](z) is defined to be ff [q−1](z), where f [0](z) = z. This is Mathematica’s Nest[f,z,q].The qth iterate can be obtained by doing O(log q) matrix multiplications, but in the interestingcase f ′(0) = f1 = 1 we can also compute the coefficients of f [q](z) by using formulas in which qis simply a numerical parameter.

Namely, as suggested by [Jabotinsky 1947], we can express thematrix power F q asI + (F −I)q = I +q1(F −I) +q2(F −I)2 +q3(F −I)3 + · · · .11

This infinite series converges, because the entry in row n and column k of (F −I)j is zero for allj > n −k. When q is any positive integer, the result defined in this way is a convolution matrix.Furthermore, the matrix entries are all polynomials in q.

Therefore the matrix obtained by thisinfinite series is a convolution matrix for all values of q.Another formula for the entries of F q was presented in [Jabotinsky 1963]. Let f (q)nk be theelement in row n and column k; thenf (q)nk =mXl=0ql(F −I)lnk=mXj=0f (j)nkmXl=jqllj(−1)l−j=mXj=0f (j)nkqj mXl=jq −jl −j(−1)l−j=mXj=0f (j)nkqjq −j −1m −j(−1)m−j ,for any m ≥n−k.

Indeed, we have p(q) = Pmj=0 p(j)qjq−j−1m−j(−1)m−j whenever p is a polynomialof degree ≤m; this is a special case of Lagrange interpolation.It is interesting to set q = 1/2 and compute convolution square roots of the Stirling numbermatrices. We have11/211/83/2105/4311/325/85512=111131176111525101;11/215/83/215/413/431109/3275/810512=11123161161245035101.The function z+z2/4+z3/48+z5/3840−7z6/92160+· · · therefore lies “halfway” beween z and ez−1 =z +z2/2!+z3/3!+· · · , and the function z +z2/4+5z3/48+5z4/96+109z5/3840+497z6/30720+· · ·lies halfway between z and ln 1/(1 −z) = z + z2/2 + z3/3 + · · · .

These half-iterates are unfamiliarfunctions; but it is not difficult to prove that z/(1 −z/2) = z + z2/2 + z3/4 + · · · is halfwaybetween z and z/(1 −z) = z + z2 + z3 + · · · . In general when f(z) = z/(1 −czk)1/k we havef [q](z) = z/(1 −qczk)1/k.It seems natural to conjecture that the coefficients of f [q](z) are positive for q > 0 whenf(z) = ln 1/(1 −z); but this conjecture turns out to be false, because Mathematica reports that12

[z8] f [q](z) = −11q/241920 + O(q2). Is there a simple necessary and sufficient condition on f thatcharacterizes when all coefficients of f [q] are nonnegative for nonnegative q?

This will happen ifand only if the entries in the first column ofln F = (F −I) −12(F −I)2 + 13(F −I)3 −· · ·are nonnegative. (See [Kuczma 1968] for iteration theory and an extensive bibliography.)Reversion.

The case q = −1 of iteration is often called reversion of series, although Mathematicauses the more proper name InverseSeries. Given f(z) = f1z + f2z2/2!

+ · · · , we seek g(z) =f [−1](z) such that gf(z)= z. This is clearly equivalent to finding the first column of the inverseof the convolution matrix.The inverse does not exist when f1 = 0, because the diagonal of F is zero in that case.Otherwise we can assume that f1 = 1, because f1gf(z/f1)= z when g reverts the power seriesf(z/f1).When f1 = 1 we can obtain the inverse by setting q = −1 in our general formula for iteration.But Lagrange’s celebrated inversion theorem for power series tells us that there is another, moreinformative, way to compute the function g = f [−1].

Let us set bF(z) = f(z)/z = 1 + f2z/2! +f3z2/3!

+ · · · . Then Lagrange’s theorem states that the elements of the matrix G = F −1 aregnk = (n −1)!

(k −1)!bFn−k(−n) ,where bFn(x) denotes the convolution family corresponding to bF(z).There is a surprisingly simple way to prove Lagrange’s theorem, using our knowledge of con-volution families. Note first thatfnk = n!

[znxk] exf(z) = n!k! [zn] f(z)k = n!k!

[zn−k] bF(z)k ;thereforefnk = n!k!bFn−k(k) .Now we need only verify that the matrix product GF is the identity, by computing its element inrow n and column m:nXk=mgnkfkm =nXk=m(n −1)! (k −1)!bFn−k (−n) k!m!bFk−m(m) .When m = n the sum is obviously 1.

When m = n −p for p > 0 it is (n −1)!/(n −p)! timesnXk=n−pk bFn−k(−n) bFk−n+p(n −p) =pXk=0(n −k) bFk(−n) bFp−k(n −p)= npXk=0bFk(−n) bFp−k(n −p) −pXk=0k bFk(−n) bFp−k(n −p)= n bFp(−p) −n bFp(−p) = 013

by the original convolution formula and the one we derived from it. The proof is complete.Extending the matrix.

The simple formula for fnk that we used to prove Lagrange’s theoremwhen f1 = 1 can be written in another suggestive form, if we replace k by n −k:fn(n−k) = nk bFk(n −k) .For every fixed k, this is a polynomial in n, of degree ≤2k. Therefore we can define the quantityfy(y−k) for all real or complex y to be yk bFk(y −k); and in particular we can define fnk in thismanner for all integers n and k, letting fnk = 0 when k > n. For example, in the case of Stirlingnumbers this analysis establishes the well-known fact that yy−kand yy−kare polynomials in yof degree 2k, and that these polynomials are multiples of yk+1 = y(y −1) .

. .

(y −k) when k > 0.The two flavors of Stirling numbers are related in two important ways. First, their matricesare inverse to each other if we attach the signs (−1)n−k to the elements in one matrix:nXk=0nk km(−1)n−k =mXk=0nk km(−1)n−k = δmn .This follows since the numbersnkcorrespond to f(z) = ez −1 and the numbersnk(−1)n−kcorrespond to g(z) = ln(1 + z), as mentioned earlier, and we have gf(z)= z.The other important relationship beweennkandnkis the striking identitynk=−k−n,which holds for all integers n and k when we use the polynomial extension method.

We can provein fact, that the analogous relationfnk = (−1)k−ng(−k)(−n)holds in the extended matrices F and G that correspond to any pair of inverse functions gf(z)= z,when f ′(0) = 1. For we have(−1)n−kg(−k)(−n) = (−1)n−k(−k −1)(−k −2) .

. .

(−n) bFn−k(k) = n!k!bFn−k(k) = fnkin the formulas above. (The interesting history of the identitynk=−k−nis traced in [Knuth1992].

The fact that the analogous formula holds in any convolution matrix was pointed out by IraGessel after he had read a draft of that paper. See also [Jabotinski 1953]; [Carlitz 1978]; [Romanand Rota 1978, section 10].

)Suppose we denote the Lah numbersnkn−1k−1(n −k)! by nk.

The extended matrix in thatcase has a pleasantly symmetrical propertynk =−k−n ,14

because the corresponding function f(z) = z/(1 −z) satisfies f−f(−z)= z. (Compare [Mullinand Rota 1969, section 9].) Near the origin n = k = 0, the nonzero entries look like this:.

. .1. .

.121. . .3661. .

.2462111216612436121.........Still more convolutions. Our proof of Lagrange’s theorem yields yet another corollary.

Supposegf(z)= z and f ′(0) = 1, and let bF(z) = f(z)/z, bG(z) = g(z)/z. Then the equationgnk = n!k!bGn−k(k) = (n −1)!

(k −1)!bFn−k(−n)tell us, after replacing n by n + k, that the identityn + kkbGn(k) = bFn(−n −k)holds for all positive integers k. Thus the polynomials bGn(x) and bFn(x) must be related by theformula(x + n) bGn(x) = x bFn(−x −n) .Now bFn(x) is an arbitrary convolution family, and bFn(−x) is another. We can conclude that if{Fn(x)} is any convolution family, then so is the set of polynomials {xFn(x + n)/(x + n)}.

Indeed,if Fn(x) corresponds to the coefficients of F(z)x, our argument proves that the coefficients of G(z)xare x Fn(x + n)/(x + n), where zG(z) is the inverse of the power series z/F(z):G(z) = FzG(z),Gz/F(z)= F(z) .The case F(z) = 1 + z and G(z) = 1/(1 −z) provides a simple example, where we know thatFn(x) =xnand Gn(x) =x+n−1n= xFn(x + n)/(x + n).A more interesting example arises when F(z) = zez/(ez −1) = z + z/(ez −1) = 1 + z/2 +B2z2/2!+B4z4/4!+· · · ; then F(−z) is the exponential generating function for the Bernoulli numbers.The convolution family for F(z)x is Fn(x) = xσn(x), where σn(x) is called a Stirling polynomial. (Actually σ0(x) = 1/x, but σn(x) is a genuine polynomial when n ≥1.) The function G suchthat Gz/F(z)= F(z) is G(z) = z−1 ln1/(1 −z); therefore the convolution family for G(z)x isGn(x) = xFn(x + n)/(x + n) = xσn(x + n).15

In this example the convolution family for exzG(z) = (1 −z)−x isx + n −1n= 1n!n0+n1x + · · · +nnxn;thereforenn −k=n! (n −k)!

Gk(n −k) =n! (n −k)!

(n −k) σk(n) = n(n −1) . .

. (n −k) σk(n) .We also havenn −k=k −n−n= (k −n)(k −1 −n) .

. .

(−n) σk(k −n) .These formulas, which are polynomials in n of degree 2k for every fixed k, explain why the σ func-tions are called Stirling polynomials. Notice that σn(1) = (−1)nBn/n!

; it can also be shown thatσn(0) = −Bn/(n · n! ).The process of going from Fn(x) to xFn(x + n)/(x + n) can be iterated: Another replacementgives xFn(x+2n)/(x+2n), and after t iterations we discover that the polynomials xFn(x+tn)/(x+tn) also form a convolution family.This holds for all positive integers t, and the convolutioncondition is expressible as a set of polynomial relations in t; therefore xFn(x + tn)/(x + tn) is aconvolution family for all complex numbers t. If Fn(x) = [zn] F(z)x, then xFn(x + tn)/(x + tn) =[zn] Ft(z)x, where Ft(z) is defined implicitly by the equationFt(z) = FzFt(z)t.In particular, we could have deduced the convolution properties of the coefficients of Bt(z)x in thisway.Let us restate what we have just proved, combining it with the “derived convolution formula”obtained earlier:Theorem.

Let Fn(x) be any family of polynomials in x such that Fn(x) has degree ≤n. IfFn(2x) =nXk=0Fk(x) Fn−k(x)holds for all n and x, then the following identities hold for all n, x, y, and t:(x + y) Fn(x + y + tn)x + y + tn=nXk=0x Fk(x + tk)x + tky Fn−ky + t(n −k)y + t(n −k);n Fn(x + y + tn)x + y + tn=nXk=1k Fk(x + tk)x + tky Fn−ky + t(n −k)y + t(n −k).Additional constructions.

We have considered several ways to create new convolution familiesfrom given ones, by multiplication or exponentiation of the associated convolution matrices, or by16

replacing Fn(x) by x Fn(x + tn)/(x + tn). It is also clear that the polynomials αnFn(βx) form aconvolution family whenever the polynomials Fn(x) do.One further operation deserves to be mentioned: If Fn(x) and Gn(x) are convolution families,then so is the family Hn(x) defined byHn(x) =nXk=0Fk(x) Gn−k(x) .This is obvious, since Hn(x) = [zn] F(z)xG(z)x.

The corresponding operation on matrices F =(fnk), G = (gnk), H = (hnk) ishnk =Xi,jnjfji g(n−j)(k−i) .If we denote this binary operation by H = F ◦G, it is interesting to observe that the associativelaw holds: (E ◦F) ◦G = E ◦(F ◦G) is true for all matrices E, F, G, not just for convolutionmatrices. A convolution matrix is characterized by the special property F ◦F = F diag(2, 4, 8, .

. .

).The construction just mentioned is merely a special case of the one-parameter familyH(t)n (x) =nXk=0Fk(x) Gn−k(x + tk) .Again, {H(t)n (x)} turns out to be a convolution family, for arbitrary t: We havenXk=0H(t)n (x)zn =Xn≥k≥0Fk(x) Gn−k(x + tk)zn =Xn,k≥0Fk(x) Gn(x + tk)zn+k=Xk≥0Fk(x)zkG(z)x+tk = G(z)xFzG(z)tx,so Hn(x) = [zn]G(z)FzG(z)tx.Applications. What’s the use of all this?

Well, we have shown that many interesting convolutionfamilies exist, and that we can deduce nonobvious facts with comparatively little effort once weknow that we’re dealing with a convolution family.One moral to be drawn is therefore the following. Whenever you encounter a triangular patternof numbers that you haven’t seen before, check to see if the first three rows have the formaba2c3aba3for some a, b, c. (You may have to multiply or divide the nth row by n!

first, and/or reflect itsentries left to right.) If so, and if the problem you are investigating is mathematically “clean,”chances are good that the fourth row will look liked4ac + 3b26a2ba4 .17

And if so, chances are excellent that you are dealing with a convolution family. And if so, you maywell be able to solve your problem.In fact, exactly that scenario has helped the author on several occasions.Asymptotics.

Once you have identified a convolution family Fn(x), you may well want to knowthe approximate value of Fn(x) when n and x are large. The remainder of this paper discussesa remarkable general power series expansion, discovered with the help of Mathematica, whichaccounts for the behavior of Fn(x) when n/x stays bounded and reasonably small as x →∞,although n may also vary as a function of x.

We will assume that Fn(x) is the coefficient of zn inF(z)x, where F(0) = F ′(0) = 1.Our starting point is the classical “saddle point method,” which shows that in many cases thecoefficient of zn in a power series P(z) can be approximated by considering the value of P at a pointwhere the derivative of P(z)/zn is zero. (See [Good 1957].) In our case we have P(z) = exf(z),where f(z) = ln F(z) = z + f2z2/2!

+ · · · ; and the derivative is zero when x f ′(z) = n/z. Let thissaddle point occur at z = s; thus, we haves f ′(s) = n/x .Near s we have f(z) = f(s) + (z −s)f ′(s) + O(z −s)2; so we will base our approximation on theassumption that the O(z −s)2contribution is zero.

The approximation to Fn(x) will be eFn(x),whereeFn(x) = [zn] expx f(s) + x (z −s) f ′(s)= ex(f(s)−sf′(s))n!xnf ′(s)n = F(s)xn! nesn.First let’s look at some examples; later we will show that the ratio Fn(x)/ eFn(x) is well behavedas a formal power series.

Throughout this discussion we will lety = n/x ;our goal, remember, is to find approximations that are valid when y is not too large, as x andpossibly n go to ∞.The simplest example is, of course, F(z) = ez and f(z) = z; but we needn’t sneeze at itbecause it will give us some useful calibration. In this case f ′′(z) = 0, so our approximation willbe exact.

We have s = y, henceeFn(x) = exyn! neyn= enn!xen= xnn!

= Fn(x) .Next let’s consider the case F(z) = T(z)/z, f(z) = T(z), when we know that Fn(x) =x(x + n)n−1/n!. In this case z T ′(z) = T(z)/1 −T(z), so we have T(s)/1 −T(s)= y orT(s) =y1 + y ,s =y1 + y e−y/(1+y)18

because T(z) = zeT (z). ThereforeeFn(x) = exy/(1+y)n!n(1 + y)ey e−y/(1+y)n= (x + n)nn!

;the ratio Fn(x)/ eFn(x) = x/(x + n) = 1/(1 + y) is indeed near 1 when y is small.If F(z) = 1 + z we find, similarly, s = y/(1 −y) andn! eFn(x) =11 −yx n(1 −y)eyn=xxe−n(x −n)x−n ;by Stirling’s approximation we also haven!

Fn(x) =x! (x −n)!

=xxe−n(x −n)x−n (1 −y)−1/21 + O(x−1).Again the ratio Fn(x)/ eFn(x) is near 1. In general if F(z) = Bt(z) the saddle point s turns out tobe y1 + (t −1)yt−1/(1 + ty)t, andn!

eFn(x) =(x + tn)x+tne−nx + (t −1)nx+(t−1)n ;a similar analysis shows that this approximation is quite good, for any fixed t.We know thatFn(x) = xnn!1 + fn(n−1)x+ fn(n−2)x2+ · · ·and that fn(n−k) is always a polynomial in n of degree ≤2k. Therefore if n2/x →0 as x →∞, wecan simply use the approximation Fn(x) = (xn/n!

)1 + O(n2/x). But there are many applicationswhere we need a good estimate of Fn(x) when n2/x →∞while n/x →0; for example, x might ben log n. In such cases eFn(x) is close to Fn(x) but xn/n!

is not.We can express s/y as a power series in y by inverting the power series expression sf ′(s) = y:s/y = 1 −f2y + (4f 22 −f3)y2/2 + (15f2f3 −30f 32 −f4)y3/6 + · · · .From this we can get a formal series for eFn(x),eFn(x) = xnn!expn(s/y)(1 + f2s/2! + f3s2/3!

+ · · · ) −n(s/y)n= xnn!1 + nf22y + 3n2f 22 −12nf 22 + 4nf324y2 + O(n3y3).We can also use the formulafn(n−k) =Xnk+k2+k3+···2!k2 k2! 3!k3 k3!

. .

. f k22 f k33.

. .

,19

where the sum is over all nonnegative k2, k3, . .

. with k2 + 2k3 + · · · = k, to writeFn(x) = xnn!1 + nf2 −f2 + O(x−1)2y+ 3n2f 22 −18nf 22 + 4nf3 + 33f 22 −12f3 + O(x−1)24y2 + O(n3y3).These series are not useful asymptotically unless ny = n2/x is small.

But the approximation eFn(x)itself is excellent, because amazing cancellations occur when we compute the ratio:Fn(x)eFn(x)= 1 −f22 y + 11f 22 −4f38y2 + O(y3) + O(x−1) .Theorem. When F(z) = exp(z + f2z2/2!

+ f3z3/3! + · · · ) and the functions Fn(x) and eFn(x) aredefined as above, the ratio Fn(x)/ eFn(x) can be written as a formal power series Pi,j≥0 cijyix−j,where y = n/x and the coefficients cij are polynomials in f2, f3, .

. .

.The derivation just given shows that we can write Fn(x)/ eFn(x) as a formal power series of theform Pi,j≥0 aijnix−j, where aij = 0 when i > 2j; the surprising thing is that we also have aij = 0whenever i > j. Therefore we can let cij = ai(i+j).To prove the theorem, we let R(z) = 1+ R1z + R2z2 + · · · stand for the terms neglected in ourapproximation:F(z)x = exf(s)−xsf′(s)1 + nsz1!

+ n2s2z22! + n3s3z33!

+ · · ·R(z) .The coefficient of zn isFn(x) = eFn(x)1 + R1s + n −1nR2s2 + (n −1)(n −2)n2R3s3 + · · ·;so the ratio Fn(x)/ eFn(x) is equal toXk≥0nknk Rksk =Xj,k≥0(−n)−jkk −jRksk =Xj(−n)−jPj ,where Pj = Pk kk−jRksk is a certain power series in s and x. The coefficients Rk are themselvespower series in s and x, because we haveR(z) = expx(z −s)2 f ′′(s)2!+ x(z −s)3 f ′′′(s)3!+ · · ·.We know from the discussion above thatkk −j= k(k −1) .

. .

(k −j) σj(k)20

is a polynomial in k. Therefore we can writePj =ϑϑ −jR(z)z=s,where ϑ is the operator that takes zk 7→k zk for all k; i.e., ϑG(z) = z G′(z) for all power seriesG(z). The theorem will be proved if we can show that Pj/nj is a formal power series in y and x−1,and if the sum of these formal power series over all j is also such a series.Consider, for example, the simplest case P0 = R(s); obviously P0 = 1.

The next simplest caseis P1 = ϑϑ−1R(z)z=s = 12ϑ(ϑ −1)R(z)z=s. It is easy to see thatϑj = zjDj ,where D is the differentiation operator D G(z) = G′(z), because zjDj takes zk into kjzk.

ThereforeP1 = 12s2R′′(s) = 12 xs2f ′′(s) .It follows that P1/n = 12(s/y)sf ′′(s) is a power series in y; it begins 12f2y + 12(f3 −f 22 )y2 + · · · .Now let’s consider Pj in general. We will use the fact that the Stirling numbers kk−jcan berepresented in the formkk −j= pj1kj + 1+ pj2kj + 2+ · · · + pjj k2j,where the coefficients pji are the positive integers in the following triangular array:123620152413021010512092423802520945.This array is clearly not a convolution matrix; but the theory developed above implies that thenumbers j!

pji/(i + j)! , namely1/22/31/43/211/8,24/513/311/16202285/125/61/32do form the convolution matrix for the powers of exp(z/2 + z2/3 + z3/4 + · · · ).

The expression kk−j= Pji=1 pji kj+iwas independently discovered by [Appell 1880], [Jordan 1933], and [Ward1934]. The number of permutations of i + j elements having no fixed points and exactly i cyclesis pji, an “associated Stirling number of the first kind” [Riordan 1958, section 4.4] [Comtet 1970,exercise 6.7].It follows thatPj = pj1sj+1 R(j+1)(s)(j + 1)!

+ pj2sj+2 R(j+2)(s)(j + 2)! + · · · + pjjs2j R(2j)(s)(2j)!.21

Now R(z) is a sum of terms having the formailxi(z −s)l ,where l ≥2i and where ail is a power series in s. Such a term contributes ailxislpj(l−j) to Pj; soit contributes ail(s/y)jsl−jxi−jpj(l−j) to Pj/nj. This contribution is nonzero only if j < l ≤2j.Since l ≥2i, we have i ≤j; so Pj/nj is a power series in y and x−1.For a fixed value of j −i, the smallest power of y that can occur in Pj/nj is y2i−j = yj−2(j−i).Therefore only a finite number of terms of Pj Pj/(−n)j contribute to any given power of y and x−1.This completes the proof.A careful analysis of the proof, and a bit of Mathematica hacking, yields the more preciseresultFn(x)eFn(x)=1(1 + s2y−1d2)1/2 +(s/y)3Ax(1 + s2y−1d2)7/2 + O(x−2) ,where A =112s3y−1d32 −34sd22 −12s2d2d3 −524s3d23 + 13yd3 + 18s3d2d4 + 18syd4 and dk = f (k)(s).Acknowledgment.

I wish to thank Ira Gessel and Svante Janson for stimulating my interest inthis subject and for their helpful comments on the first draft. Ira Gessel and Richard Brent alsointroduced me to several relevant references.ReferencesAppell, P.1880.“D´eveloppement en s´erie enti`ere de (1 + ax)1/x.” Archiv der Mathematik undPhysik 65: 171–175.Bell, E. T.1934.“Exponential numbers.” American Mathematical Monthly 41: 411–419.Carlitz, L.1978.“Generalized Stirling and related numbers.” Rivista di Matematica della Uni-versit`a di Parma, serie 4, 4: 79–99.Cayley, A.1859.“On the analytical forms called trees.

Second part.” Philosophical Magazine18: 371–378. Reprinted in Cayley’s Collected Mathematical Papers 4: 112–115.Comtet, Louis.1970.Analyse Combinatoire.

Presses Universitaires de France. (English trans-lation, Advanced Combinatorics, D. Reidel, Dordrecht, 1974.

)Garsia, Adriano M.1973.“An expos´e of the Mullin-Rota theory of polynomials of binomialtype.” Linear and Multilinear Algebra 1: 47–65.Gomes Teixeira, F.1896.“Sur le d´eveloppement de xk en s´erie ordonn´ee suivant les puissancesdu sinus de la variable.” Nouvelles Annales de Math´ematiques, s´erie 3, 15: 270–274.Good, I. J.1957.“Saddle-point methods for the multinomial distribution.” Annals of Mathe-matical Statistics 28: 861–881.Gould, H. W., and Kauck´y, I.1966.“Evaluation of a class of binomial coefficient summations.”Journal of Combinatorial Theory 1: 233–248.Graham, Ronald L., Knuth, Donald E., and Patashnik, Oren.1989.Concrete Mathematics.Addison-Wesley, Reading, Massachusetts.22

Harris, Bernard, and Schoenfeld, Lowell.1967.“The number of idempotent elements in sym-metric semigroups.” Journal of Combinatorial Theory 3: 122–135.Jabotinsky, Eri.1947.“Sur la repr´esentation de la composition de fonctions par un produit dematrices. Applicaton `a l’it´eration de ex et de ex −1.” Comptes Rendus Hebdomadaires desSciences de L’Academie des Sciences, 224: 323–324.Jabotinsky, Eri.1953.“Representation of functions by matrices.

Application to Faber polyno-mials.” Proceedings of the American Mathematical Society 4: 546–553.Jabotinsky, Eri.1963.“Analytic iteration.” Transactions of the American Mathematical Society108: 457–477.Jordan, Charles.1933.“On Stirling’s numbers.” Tˆohoku Mathematical Journal 37: 254–278.Knuth, Donald E.1992.“Two notes on notation.” American Mathematical Monthly 99: 403–422.Knuth, Donald E. and Pittel, Boris.1989.“A recurrence related to trees.” Proceedings of theAmerican Mathematical Society 105: 335–349.Kuczma, Marek.1968.Functional Equations in a Single Variable. Polish Scientific Publishers,Warsaw.Lah, I.1955.“Eine neue Art von Zahlen, ihre Eigenschaften und Anwendung in der mathemat-ischen Statistik,” Mitteilungsblatt f¨ur Mathematische Statistik 7: 203–212.Motzkin, T. S.1971.“Sorting numbers for cylinders and other classification numbers.” Pro-ceedings of Symposia in Pure Mathematics 19: 167–176.Mullin, Ronald, and Rota, Gian-Carlo.1970.“On the foundations of combinatorial theory.

III.Theory of binomial enumeration.” In Graph Theory and Its Applications, edited by BernardHarris (Academic Press, 1970), 167–213.Riordan, John.1958.An Introduction to Combinatorial Analysis. John Wiley & Sons, NewYork.Riordan, John.1968.Combinatorial Identities.

John Wiley & Sons, New York.Rogers, D. G.1978.“Pascal triangles, Catalan numbers and renewal arrays.” Discrete Mathe-matics 22: 301–310.Roman, Steven.1978.The Umbral Calculus. Academic Press, Orlando.Roman, Steven M., and Rota, Gian-Carlo.1978.“The umbral calculus.” Advances in Mathe-matics 27: 95–188.Rota, Gian-Carlo, Kahaner, D., and Odlyzko, A.1973.“On the foundations of combinatorialtheory.

VIII. Finite operator calculus.” Journal of Mathematical Analysis and Applications42: 884–760.Reprinted in Rota, Gian-Carlo, Finite Operator Calculus (Academic Press,1975), 7–82.Steffensen, J. F.1941.“The poweroid, an extension of the mathematical notion of power.” ActaMathematica 73: 333–366.Ward, Morgan.1934.“The representation of Stirling’s numbers and Stirling’s polynomials assums of factorials.” American Journal of Mathematics 56: 87–95.23

---

출처: arXiv:9207.221 • 원문 보기