A little challenge.

[Peter Luschny]

What is the general law of this rectangular array?

[1, 2, 4, ...]
[1, 3, 10, 41, 196, 1057, 6322, 41393, 293608]
[1, 4, 20, 127, 967, 8549, 85829, 962308, 11895252]
[1, 5, 33, 280, 2883, 34817, 481477, 7489454, 129259662]
[1, 6, 49, 518, 6689, 101841, 1783170, 35250562, 775700824]
[1, 7, 68, 859, 13310,243946, 5155512,123294103,3288775809]
[1, 8, 90,1321, 23851,510502,12622252,353704058,11070101343]
[...]

Tell us how you approache such a problem.

[Alois Heinz]

A(n,k)=exponential transform applied n times to identity function
evaluated at k:

# Maple code
exptr:=
proc (p) local g;
g:= proc(n) option remember;
p(n) +add (binomial (n-1, k-1) *p(k) *g(n-k), k=1..n-1)
end:
end:
A:= n-> (exptr@@n)(i->i):
seq(lprint(seq(A(n)(k), k=1..11)), n=0..11);

On Mar 3, 10:14 pm, Peter Luschny <peter.lusc...@googlemail.com>
wrote:

[Peter Luschny]

/> Alois Heinz wrote:
/> A(n,k)=exponential transform applied n times to identity function
/> evaluated at k:

Thanks Alois. Yes, this is precisely how I defined the problem.

Now my question is: How did you approach the problem? What hints did
you see? Did you look first at the rows or at the columns? Did you
apply transformations to some individual rows or columns? Or do you
have a general heuristic when approaching a problem of this kind?

In other words, please tell us more about your professional
secrets ;-)

Cheers, Peter

[Peter Luschny]

/>Alois Heinz wrote:
/> # Maple code
/> exptr:=
/> proc (p) local g;
/> g:= proc(n) option remember;
/> p(n) +add (binomial (n-1, k-1) *p(k) *g(n-k), k=1..n-1)
/> end:
/> end:
/> A:= n-> (exptr@@n)(i->i):

Here is the Sage code which I used to formulate the problem
(apart from the case for row 0, which I added later).

def h(x) :
@CachedFunction
def b(n) :
if x == 0 : return n
if n == 0 : return 1
return add(binomial(n-1,k-1)*x(k)*b(n-k) for k in (1..n))
return b

B = 9
print [h(0)(i) for i in (1..B)]
print [h(x)(i) for i in (1..B)]
print [h(h(x))(i) for i in (1..B)]
print [h(h(h(x)))(i) for i in (1..B)]
print [h(h(h(h(x))))(i) for i in (1..B)]
print [h(h(h(h(h(x)))))(i) for i in (1..B)]

It seems that Python/Sage does not offer iterated functional
composition as an operation like for example Maple does,
see "exptr@@n" in Alois' code. Any idea how to accomplish
this?

[Alois Heinz]

These are many questions, Peter. I will try to give correct answers.

I looked at columns and rows simultaneously. I remembered that
I had seen them before. One row is "Number of forests with n nodes
and height at most 1". Next row is "Natural numbers exponentiated
twice".
I have written code for this in 2008. The rest was easy and it took
me
only a few minutes.

When trying to solve a problem like this, I think that I apply
methods
similar to those described by George Pólya in his book "How to Solve
It".
So the secret is a mixture of experience (a good data base like OEIS
helps),
a set of different approaches, trial and error, and patience (if the
problem is important enough or some reward is waiting), ...

Best regards, Alois

[Peter Luschny]

Thank you for your insights, Alois.

> One row is "Number of forests with n nodes and height at most 1".

http://oeis.org/A000248

> Next row is "Natural numbers exponentiated twice".

http://oeis.org/A007550

Both sequences are from the 'Urgestein', (formerly M2857 and M3568),
which is a promising start for a generalization. So we might look for
a combinatorial interpretation in terms of set partitions.

[Peter Luschny]

Just for the records: The columns of the array are nontrivial
polynomials.

def s0(n) : return 1/factorial(0)
def s1(n) : return (n + 2)/factorial(1)
def s2(n) : return (3*n^2 + 11*n + 6)/factorial(2)
def s3(n) : return (18*n^3 + 93*n^2 + 111*n + 24)/factorial(3)
def s4(n) : return (180*n^4 + 1180*n^3 + 2160*n^2 + 1064*n + 120)/
factorial(4)
def s5(n) : return (2700*n^5 + 21225*n^4 + 51850*n^3 + 41835*n^2 +
8510*n + 720)/factorial(5)
def s6(n) : return (56700*n^6 + 515970*n^5 + 1563450*n^4 + 1770930*n^3
+
594930*n^2 + 44820*n + 5040)/factorial(6)

[Arie Groeneveld]

After understanding how things work (exponential transform and the recursive Maple and Sage programs) here's an iterative program in Haskell:

import Data.List
import Control.Arrow
import Control.Monad

zippWith = zipWith. zipWith
pascalTriangle = iterate (ap (zipWith (+). (0:)) (++[0])) [1]
 
seqGen n =
  iterate (tail. foldl (\cs xs -> cs ++[sum $ zipWith (*) cs (reverse xs)] ) [1]
           . zippWith (*) pascalTriangle. tail. inits) [1..n]

Some output executed in the GHC interpreter GHCi:

*Main> mapM_ print $ take 11 $ seqGen 9
[1,2,3,4,5,6,7,8,9]
[1,3,10,41,196,1057,6322,41393,293608]
[1,4,20,127,967,8549,85829,962308,11895252]
[1,5,33,280,2883,34817,481477,7489454,129259662]
[1,6,49,518,6689,101841,1783170,35250562,775700824]
[1,7,68,859,13310,243946,5155512,123294103,3288775809]
[1,8,90,1321,23851,510502,12622252,353704058,11070101343]
[1,9,115,1922,39597,968624,27407429,879897348,31581220687]
[1,10,143,2680,62013,1705872,54333217,1965191524,79524291987]
[1,11,174,3613,92744,2832951,100274470,4033230317,181603522854]
[1,12,208,4739,133615,4486411,174669967,7731854648,383320429534]

[Peter Luschny]

Welcome Arie! Very neat!

Some day I will have to dip deeper into Haskell.
Peter

[Arie Groeneveld]

Perhaps the following finds are of some interest too.

The sum of the coefficients of the column polynomials gives you the second row.

Looking at the numerators of the coefficients of the polynomials I noticed that

the coefficients of the highest  exponents can be find in A006472

1 1 3 18 180 2700 56700 1587600 57153600

the numerators of the coefficients of the lowest exponents are factorial (n+1)

1 2 6 24 120 720 5040 40320 362880

[Paul D. Hanna]

As you knowm, when trying to determine the general law behind an
array,
one of the most common transforms to attempt first is the inverse
binomial transform.
When we apply that transform on the first row, and find something
promising.

OBSERVATION.

The inverse binomial transform of row 1 yields
A138911 = [1, 2, 5, 19, 81, 401, 2233, 13721, 91969, ...].

We notice that the e.g.f. of A138911, A(x), satisfies:
[x^n/n!] A(x)*exp(-n*x) = 1 for n>=0.

APPROACH.
A138911 gives us a hunch to investigate: generate a table
of coefficients in R(n,x)*exp(-n*x) for each row e.g.f. R(n,x)
and analyze the terms found along the main diagonal.

RESULTS.
--------------------------------------------------
ROW 1.
ROW=[1, 3, 10, 41, 196, 1057, 6322, 41393, 293608]
EGF=Ser(vector(#ROW,n,ROW[n]/(n-1)!))
for(n=1,#EGF,print(Vec(serlaplace(EGF*exp(-n*x)))))

[1, 2, 5, 19, 81, 401, 2233, 13721, 91969]
[1, 1, 2, 9, 28, 145, 726, 4249, 27000]
[1, 0, 1, 5, 1, 79, 121, 1511, 6721]
[1, -1, 2, 1, -12, 113, -422, 2441, -6584]
[1, -2, 5, -9, 1, 157, -1263, 8173, -45087]
[1, -3, 10, -31, 76, 1, -1922, 19841, -153896]
[1, -4, 17, -71, 273, -805, 1, 29339, -359135]
[1, -5, 26, -135, 676, -3071, 10626, 1, -525144]
[1, -6, 37, -229, 1393, -8087, 42313, -167839, 1]

Diagonal: All 1's - E.g.f.: exp(x)-1.

PARI CODE: to get row 1 from diagonal:
DIAG=[1, 1, 1, 1, 1, 1, 1, 1, 1]
{a(n)=DIAG[n+1] + sum(k=0, n-1, binomial(n, k) * (n+1)*(k+1)^(n-
k-1)*DIAG[k+1] )}
for(n=0,#DIAG-1,print1(a(n),","))

--------------------------------------------------
ROW 2.
ROW=[1, 4, 20, 127, 967, 8549, 85829, 962308, 11895252]
EGF=Ser(vector(#ROW,n,ROW[n]/(n-1)!))
for(n=1,#EGF,print(Vec(serlaplace(EGF*exp(-n*x)))))

[1, 3, 13, 78, 564, 4803, 46777, 511241, 6182363]
[1, 2, 8, 47, 319, 2647, 25037, 267402, 3168676]
[1, 1, 5, 28, 172, 1451, 13109, 138055, 1602171]
[1, 0, 4, 15, 87, 825, 6493, 71624, 794132]
[1, -1, 5, 2, 52, 499, 2609, 40893, 362107]
[1, -2, 8, -17, 79, 203, 437, 31246, 83748]
[1, -3, 13, -48, 204, -453, 877, 28907, -158629]
[1, -4, 20, -97, 487, -2099, 7829, 4140, -322604]
[1, -5, 29, -170, 1012, -5725, 29993, -115631, 21147]

Diagonal: A000110 - E.g.f.: exp(exp(x)-1)-1.
Bell or exponential numbers
[1, 2, 5, 15, 52, 203, 877, 4140, 21147, ...]

PARI CODE: to get row 2 from diagonal:
DIAG=[1, 2, 5, 15, 52, 203, 877, 4140, 21147]
{a(n)=DIAG[n+1] + sum(k=0, n-1, binomial(n, k) * (n+1)*(k+1)^(n-
k-1)*DIAG[k+1] )}
for(n=0,#DIAG-1,print1(a(n),","))

--------------------------------------------------
ROW 3.
ROW=[1, 5, 33, 280, 2883, 34817, 481477, 7489454, 129259662]
EGF=Ser(vector(#ROW,n,ROW[n]/(n-1)!))
for(n=1,#EGF,print(Vec(serlaplace(EGF*exp(-n*x)))))

[1, 4, 24, 195, 1942, 22896, 310686, 4758508, 81062649]
[1, 3, 17, 134, 1291, 14915, 198877, 3002900, 50537278]
[1, 2, 12, 91, 846, 9644, 126310, 1883144, 31330713]
[1, 1, 9, 60, 547, 6213, 79485, 1174834, 19316622]
[1, 0, 8, 35, 358, 3992, 49342, 731684, 11830777]
[1, -1, 9, 10, 267, 2471, 30181, 458648, 7158654]
[1, -2, 12, -21, 286, 1140, 19302, 290080, 4214553]
[1, -3, 17, -64, 451, -631, 17365, 167894, 2392078]
[1, -4, 24, -125, 822, -3712, 29470, 14684, 1606137]

Diagonal: A000258 - E.g.f.: exp(exp(exp(x)-1)-1)-1.
[1, 3, 12, 60, 358, 2471, 19302, 167894, 1606137, ...]

PARI CODE: to get row 3 from diagonal:
DIAG=[1, 3, 12, 60, 358, 2471, 19302, 167894, 1606137]
{a(n)=DIAG[n+1] + sum(k=0, n-1, binomial(n, k) * (n+1)*(k+1)^(n-
k-1)*DIAG[k+1] )}
for(n=0,#DIAG-1,print1(a(n),","))

--------------------------------------------------
ROW 4.
ROW=[1, 6, 49, 518, 6689, 101841, 1783170, 35250562, 775700824]
EGF=Ser(vector(#ROW,n,ROW[n]/(n-1)!))
for(n=1,#EGF,print(Vec(serlaplace(EGF*exp(-n*x)))))

[1, 5, 38, 388, 4888, 73115, 1262799, 24690060, 538362539]
[1, 4, 29, 288, 3545, 52199, 890210, 17227618, 372427448]
[1, 3, 22, 212, 2552, 37083, 624843, 11977750, 256842043]
[1, 2, 17, 154, 1825, 26237, 436698, 8300250, 176606360]
[1, 1, 14, 108, 1304, 18491, 303815, 5735512, 121082219]
[1, 0, 13, 68, 953, 12915, 210474, 3954970, 82755064]
[1, -1, 14, 28, 760, 8699, 146115, 2721618, 56348603]
[1, -2, 17, -18, 737, 5033, 104978, 1855570, 38231768]
[1, -3, 22, -76, 920, 987, 86463, 1199620, 26097835]

Diagonal: A000307 - E.g.f.: exp(exp(exp(exp(x)-1)-1)-1)-1.
[1, 4, 22, 154, 1304, 12915, 146115, 1855570, 26097835, ...]

PARI CODE: to get row 4 from diagonal:
DIAG=[1, 4, 22, 154, 1304, 12915, 146115, 1855570, 26097835]
{a(n)=DIAG[n+1] + sum(k=0, n-1, binomial(n, k) * (n+1)*(k+1)^(n-
k-1)*DIAG[k+1] )}
for(n=0,#DIAG-1,print1(a(n),","))

--------------------------------------------------

From the above results, we formulate the general law for all the
rows.

GENERAL LAW.
Let us denote T as the table to be defined.

Given R(n,x) is the e.g.f. of row n of table T, we find that:
[x^k/k!] R(n,x)*exp(-k*x) = [x^k/k!] {n-th iteration of exp(x)-1}
which forms a diagonal in the array of coefficients in the
iterated inverse binomial transforms of row n of table T.

We now have the following

FORMULA.
Let E^n(k) denote the coefficient of x^k/k! in the n-th iteration of
exp(x)-1,
then

T(n,k) = E^n(k+1) + Sum_{j=0..k-1} binomial(k,j)*(k+1)*(j+1)^(k-j-1) *
E^n(j+1)

where

E^n(k) = [x^k/k!] R(n,x)*exp(-k*x) and

R(n,x) = Sum_{k>=0} T(n,k)*x^k/k! is the e.g.f. of the n-th row of T.

PROGRAM:
(PARI)
/* Coefficient of x^k/k! in the n-th Iteration of exp(x)-1: */
{EXPIT(n,k)=local(Exp1=exp(x
+x*O(x^k))-1,G=x);for(i=1,n,G=subst(Exp1,x,G));k!*polcoeff(G,k)}

/* Coefficient of x^k/k! in the e.g.f. of the n-th row: */
{T(n,k)=EXPIT(n,k+1) + sum(j=0, k-1, binomial(k, j) * (k+1)*(j+1)^(k-
j-1)*EXPIT(n,j+1) )}

/* Print the table: */
for(n=0,8,for(k=0,8,print1(T(n,k),", "));print(""))

OUTPUT:
1, 2, 3, 4, 5, 6, 7, 8, 9,
1, 3, 10, 41, 196, 1057, 6322, 41393, 293608,
1, 4, 20, 127, 967, 8549, 85829, 962308, 11895252,
1, 5, 33, 280, 2883, 34817, 481477, 7489454, 129259662,
1, 6, 49, 518, 6689, 101841, 1783170, 35250562, 775700824,
1, 7, 68, 859, 13310, 243946, 5155512, 123294103, 3288775809,
1, 8, 90, 1321, 23851, 510502, 12622252, 353704058, 11070101343,
1, 9, 115, 1922, 39597, 968624, 27407429, 879897348, 31581220687,
1, 10, 143, 2680, 62013, 1705872, 54333217, 1965191524, 79524291987,
...

Paul D. Hanna

[Peter Luschny]

Thanks for the very nice explanation of your method Paul.

(Just for the record: Paul send me his solution by pm only a few
hours after I had posted this challenge on the seqfan list.)

Let me sum up a little bit.

Our problem is very similar to a famous problem once
investigated by E. T. Bell (The Iterated Exponential Integers,
Annals of Mathematics, 39 (1938), 539-557).

Indeed, if we start from Alois' elegant solution

exptr := proc(p) local g;
g := proc(n) option remember; local k;
if n = 0 then 1 else
add(binomial (n-1, k-1)*p(k)*g(n-k), k=1..n) fi
end:
end:

then

B := n-> (exptr@@n)(i->1):
seq(lprint(seq(B(n)(k), k=0..6)), n=0..6);

leads to Alois' A144150 (rows and cols interchanged) wheras

A := n -> (exptr@@n)(i->i):
seq(lprint(seq(A(n)(k), k=0..6)), n=0..6);

gives the solution to our problem.

Thus it is exactly the same procedure, once applied to the
sequence i->1 (Bell) and next to the sequence i->i (our case).

Therefore, does it make sense to introduce a 'Bell-transformation'
which maps an arbitrary sequence s to (exptr@@n)(s)?

In any case, the problem of a nice combinatorial interpretation
for our case seems to be still open. Any suggestions?

Peter

P.S. Note that the square array which I will finally submit to OEIS
will start at k=0 and thus will have an additional first column.

[Peter Luschny]

Thanks to all which took up this little challenge.

The sequence number is A209631. Please add your formulas and programs
to the draft. (Note that the case k=0 might require adjustments.)

Peter