User:Peter Luschny/SeidelTransform

An old operation on sequences:
the Seidel transform

The Seidel-Euler triangle (open in a new window to enlarge).

The Seidel-Euler triangle

The picture above shows a triangle which illustrates an algorithm which Seidel describes:

"Abgesehen von der 1 an der Spitze, und von den Nullen, mit welchen die Zeilen abwechselnd links und rechts beginnen, ist jede Zahl der Tafel die Summe aus der neben ihr stehenden kleineren, und der gerade über der letzteren befindlichen."

Apart from the 1 at the tip, and of the zeros with which the lines start alternately left and right, any number in the table is the sum of the smaller standing next to it, and the one just over the latter.

Seidel describes the course of the computation "... always starts on the side where the 0 is, and alternately calculates from left to right and from right to left (boustrophedon).") ("... stets auf der Seite beginnt, wo in ihnen die 0 steht, und also abwechselnd von links nach rechts und von rechts nach links rechnet (βoνστρoφηδóν).") Because of this characteristic the algorithm is also called 'boustrophedon transform'.

The triangle above computes the secant U_2n and the tangent T_2n+1 numbers (n≥0) in unparalleled simplicity and elegance just starting from the 1. An implementation of the triangle with Sage (Python):

def SeidelEulerTriangle(n) :
    A = {-1:0, 0:1}
    k = 0; e = 1
    for i in (0..n) :
        Am = 0
        A[k + e] = 0
        e = -e
        for j in (0..i) :
            Am += A[k]
            A[k] = Am
            k += e
        print [A[z] for z in (-i//2..i//2)]

SeidelEulerTriangle(11)

[1]
[0, 1]
[1, 1, 0]
[0, 1, 2, 2]
[5, 5, 4, 2, 0]
[0, 5, 10, 14, 16, 16]
[61, 61, 56, 46, 32, 16, 0]
[0, 61, 122, 178, 224, 256, 272, 272]
[1385, 1385, 1324, 1202, 1024, 800, 544, 272, 0]
[0, 1385, 2770, 4094, 5296, 6320, 7120, 7664, 7936, 7936]
[50521, 50521, 49136, 46366, 42272, 36976, 30656, 23536, 15872, 7936, 0]
[0,50521,101042,150178,196544,238816,275792,306448,329984,345856,353792,353792]

The computation of the Euler numbers

Number of alternating permutations (Secant + Tangent)

To extract the number of alternating permutations we can simply collect the nonzero diagonal elements alternately from the right and the left hand side of Seidel's triangle.

def A000111_list(n) :
    R = []; A = {-1:0, 0:1}; k = 0; e = 1
    for i in (0..n) :
        Am = 0; A[k + e] = 0; e = -e
        for j in (0..i) : Am += A[k]; A[k] = Am; k += e
        R.append(A[-i//2] if e < 0 else A[i//2])
    return R
A000111_list(22) 

Secant numbers (Euler numbers)

The elements of the left diagonal (the secant numbers) are the signless Euler numbers, a name given by Raabe in 1851 to the secant numbers. Sylvester followed Raabe's lead in 1861 in "Note on the numbers of Bernoulli and Euler ...". Ever since this name is an established convention in mathematical writings. All major computer algebra systems (Maple, Mathematica, Maxima, Sage) implement Euler numbers according to this definition. In the second half of the last century the Handbook of Mathematical Functions and nowadays the NIST Digital Library of Mathematical Functions define the Euler numbers in this way.

For some reason R. P. Stanley does not follow this traditional convention. He calls the number of alternating permutations Euler numbers. The present author does not think that Stanley's naming deserves adoption. It generates only confusion with Raabe's definition. A much better name for the number of alternating permutations seems to be 'André numbers' in honor of Désiré André. André studied 2-alternating permutations in 1879 ("Développements de sec x et tg x") and in 1881 ("Mémoire sur les permutations alternées") . Moreover, when generalizing the notion of alternating permutations to n-alternating permutations a good identifier is needed which is different from the name for generalized secant numbers (generalized Euler numbers). Table 1 and Table 2 below show the first few values of these generalizations.

Without the zero entries the Euler numbers can be computed as

# n -> [a(0), a(1), ..., a(n-1)] for n > 0.
def A000364_list(len) :
    R = []; A = {-1:0, 0:1}; k = 0; e = 1
    for i in (0..2*len-1) :
        Am = 0; A[k + e] = 0; e = -e
        for j in (0..i) : Am += A[k]; A[k] = Am; k += e
        if e < 0 : R.append(A[-i//2])
    return R

Tangent numbers

The elements of the right diagonal are the tangent numbers. Without the zero entries they can be computed as

# n -> [a(1), ..., a(n)] for n >= 1.
def A000182_list(len) :
    R = []; A = {-1:0, 0:1}; k = 0; e = 1
    for i in (0..2*len-1) :
        Am = 0; A[k + e] = 0; e = -e
        for j in (0..i) : Am += A[k]; A[k] = Am; k += e
        if e > 0 : R.append(A[i//2])
    return R

Median Euler numbers

These are the elements in the even numbered rows of the central column of Seidel's triangle. They can be computed as

def A000657(n) :
    R = []; A = {-1:0, 0:1}; k = 0; e = 1
    for i in (0..n) :
        Am = 0; A[k + e] = 0; e = -e
        for j in (0..i) :
            Am += A[k]; A[k] = Am; k += e
        if e < 0 :
            R.append(A[0])
    return R
    
A000657(9)  
1, 1, 4, 46, 1024, 36976, 1965664, 144361456, 13997185024

Generalized Euler numbers

Euler gave the expansion of sec in 1755 ("Institutiones Calculi Differentialis" ) as well as their connection to permutations. Following a proposal of Raabe these numbers are nowadays called Euler numbers by almost all major mathematicians.

Trivially the Euler numbers are a subsequence of the expansion of sec + tan and thus count alternating permutations. But it is a special kind of permutations: the permutations which act on sets whose cardinality is a multiple of 2. (For example E₄=5 count the permutations {2143, 3142, 3241, 4132, 4231}).

As soon as one introduces the general notion of n-alternating permutations and observes that the ordinary alternating permutations are the same as the 2-alternating permutations one sees that more goes on than an even-odd glasperlenspiel. The generalized Euler numbers count the n-alternating permutations in the symmetric group of permutations of [m] where n divides m. Such permutations can be expected to have more interesting properties then the bulk of permutations which act on a [m] where m is not a multiple of n.

Tabel 1.  Euler Numbers     A181985

m\n	0	1	2	3	4	5
1	1	1	1	1	1	1
2	1	1	5	61	1385	50521	A000364
3	1	1	19	1513	315523	136085041	A002115
4	1	1	69	33661	60376809	288294050521	A211212
5	1	1	251	750751	11593285251	613498040952501
6	1	1	923	17116009	2301250545971	1364944703949044401
			A030662	A211213	A181991

A very nice feature of Seidel's algorithm is that it carries over from the classical Euler-secant-2-alternating-permutation case to the general Euler-n-alternating-permutation case with essentially no modification.

EulerNumbers := proc(n, len) local E, dim, i, k;
dim := n*(len-1); E := array(0..dim, 0..dim); E[0, 0] := 1;

for i from 1 to dim do
if i mod n = 0 then E[i, 0] := 0 ;
   for k from i-1 by -1 to 0 do E[k, i-k] := E[k+1, i-k-1] + E[k, i-k-1] od;
else E[0, i] := 0;
   for k from 1 by 1 to i do E[k, i-k] := E[k-1, i-k+1] + E[k-1, i-k] od;
fi 
od; seq(E[0, n*k], k=0..len-1) end:

for n from 1 to 6 do print(EulerNumbers(n, 6)) od; # gives the table above.

The computation of the André numbers

For an integer n>0, a permutation s = s₁ ... s_k is an n-alternating permutation if it has the property that s_i < s_i+1 if and only if n divides i. The André numbers count the n-alternating permutations.

Tabel 2.  André Numbers     A181937

n\k	0	1	2	3	4	5	6	7	8	9	10	11
1	1	1	1	1	1	1	1	1	1	1	1	1
2	1	1	1	2	5	16	61	272	1385	7936	50521	353792	A000111
3	1	1	1	1	3	9	19	99	477	1513	11259	74601	A178963
4	1	1	1	1	1	4	14	34	69	496	2896	11056	A178964
5	1	1	1	1	1	1	5	20	55	125	251	2300	A181936
6	1	1	1	1	1	1	1	6	27	83	209	461

As said before, the Seidel's algorithm carries over to the general case of n-alternating permutations with virtually no modification. Convince yourself:

AndreNumbers := proc(modulus, dim) local E,DIM,n,k;
DIM := dim-1; E := array(0..DIM, 0..DIM); E[0,0] := 1;
for n from 1 to DIM do
if n mod modulus = 0 then E[n,0] := 0 ;
   for k from n-1 by -1 to 0 do E[k,n-k] := E[k+1,n-k-1] + E[k,n-k-1] od;
else E[0,n] := 0;
   for k from 1 by 1 to n do E[k,n-k] := E[k-1,n-k+1] + E[k-1,n-k] od;
fi od; [E[0,0],seq(E[k,0]+E[0,k],k=1..DIM)] end:

André numbers of 3- and 4-alternating permutations

A178963_list := dim -> AndreNumbers(3, dim):
A178963_list(12); 1, 1, 1, 1, 3, 9, 19, 99, 477, 1513, 11259, 74601

According to Mendes and Remmel the exponential generating function is

(3+2*sqrt(3)*exp(x/2)*sin(sqrt(3)*x/2))/(exp(-x)+2*exp(x/2)*cos(sqrt(3)*x/2)).

Hackers might call this a good obfuscation of the mathematical meaning.

In the case n = 4:

A178964_list := dim -> AndreNumbers(4, dim):
A178964_list(12); 1, 1, 1, 1, 1, 4, 14, 34, 69, 496, 2896, 11056

According to Mendes and Remmel the exponential generating function is in this case

(1+sqrt(2)*sin(x/sqrt(2))*cosh(x/sqrt(2))+sin(x/sqrt(2))*sinh(x/sqrt(2)))
/(cos(x/sqrt(2))*cosh(x/sqrt(2))).

Hackers might call this a perfect obfuscation of the mathematical meaning.

André numbers of n-alternating permutations

The generating function of 23-alternating permutations has not yet been published to the author's knowledge. In view of the e. g. f. in the case n = 4 it might well be the case that the generating function in terms of elementary functions (if such a formula exists at all) is long and prone to error. What is the answer? Arguable a short and general algorithm. Indeed, it is easy to compute the André numbers also in the case n = 23 with the Seidel transform.

AndreNumbers(23, 34);
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
23, 299, 2599, 17549, 98279, 475019, 2035799, 7888724, 28048799, 92561039]

André numbers and Pascal's triangle minus 1.

A more compact display of the André numbers can be achieved if the rows are shifted k places to the left (thus disregarding k leading 1's). In table 3 below additionally the rows and columns are interchanged.

Tabel 3.  André Numbers by diagonals

n\k	1	2	3	4	5	6	7	8	9	10
1	1	2	3	4	5	6	7	8	9	10	A000027 (k≥1)
2	1	5	9	14	20	27	35	44	54	65	A000096 (k≥2)
3	1	16	19	34	55	83	119	164	219	285	A062748 (k≥3)
4	1	61	99	69	125	209	329	494	714	1000	A063258 (k≥4)
5	1	272	477	496	251	461	791	1286	2001	3002	A062988 (k≥5)
6	1	1385	1513	2896	2300	923	1715	3002	5004	8007	A124089 (k≥6)

The diagonal of this array are the central binomial coefficients minus 1, A030662. More is true: all the values on the diagonal and to the right of the diagonal (color blue) obey a simple formula: T(n, k) = binomial(n + k, k) - 1. In other words, the upper triangle of table 3 is the upper triangle of Pascal's triangle minus 1 indexed as a square array (see the table of A014473).

Milan Janjic pointed out the combinatorial interpretation of T(n, k) = binomial(n + k, k) - 1. If X is a (n+k)-set and Y a fixed k-subset of X then T(n, k) is equal to the number of n-subsets of X intersecting Y.

A little bit of Sage hacking confirms this interpretation (the proof is not hard, though).

def intersectsubsets(n, k) :
    X = Set(range(n))
    Y = Subsets(X, k).first()
    L = Subsets(X, n-k).list()
    F = filter(lambda x: x.intersection(Y) <> Set([]), L)
    return len(F)

for n in(0..6): print [n], [intersectsubsets(n+k, k) for k in (0..9)]

To sum up: For k ≥ n André(k, n + k) = binomial(n + k, k) - 1.

André polynomials.

Let $u=\left(s^{2}-2\right)^{1/2}$ and

$F(s,x)=ux-2\log \left({\frac {1}{2}}\left(e^{ux}\left(1-{\frac {s}{u}}\right)+{\frac {s}{u}}+1\right)\right)$

then for n ≥ 0 the sequence of polynomials

$p_{n}(s)=(n+2)![x^{n+2}]F(s,x)$

are the André polynomials. (For a different definition see A094503.)

def AndrePoly(n) :
    s = var('s'); u = sqrt(s^2-2)
    egf = u*x-2*ln((exp(u*x)*(1-s/u)+s/u+1)/2)
    return factorial(n+2)*egf.series(x,n+4).coeff(x,n+2)

def AndreCoeffs(n) : return [AndrePoly(n).coeff(s,i) for i in (0..n)]

n	André coefficients											André numbers
0	1											1
1	0	1										1
2	1	0	1									2
3	0	4	0	1								5
4	4	0	11	0	1							16
5	0	34	0	26	0	1						61
6	34	0	180	0	57	0	1					272
7	0	496	0	768	0	120	0	1				1385
8	496	0	4288	0	2904	0	247	0	1			7936
9	0	11056	0	28768	0	10194	0	502	0	1		50521
10	11056	0	141584	0	166042	0	34096	0	1013	0	1	353792

Observe that the (nonzero) subdiagonal is the same as the subdiagonal of A008292, the Eulerian numbers A000295. The leftmost column shows the aerated reduced tangent numbers A002105, which, by a comment of Wenjin Woan, count increasing complete binary trees with exponential generating function sec^2(x/sqrt(2)).

Without zeros and in row-reversed order the triangle A094503 is generated by

def A094503_row(n) : return [AndrePoly(n).coeff(s,n-2*i) for i in (0..n//2)]

However, probably the best way to define the André polynomials is via a simple recursion of their coefficients.

def AndreCoeff(n, k) :
    if n==k : return 1 
    if k<1 or k>n : return 0 
    return ((n-k)//2+1)*AndreCoeff(n-1,k-1)+k*AndreCoeff(n-1,k+1)

for n in (0..8): print "p_",n,"(x) =",
                       add(AndreCoeff(n,k)*x^k for k in (0..n))

p_0(x) = 1
p_1(x) = x
p_2(x) = x^2
p_3(x) = x+x^3
p_4(x) = 4*x^2+x^4
p_5(x) = 4*x+11*x^3+x^5
p_6(x) = 34*x^2+26*x^4+x^6
p_7(x) = 34*x+180*x^3+57*x^5+x^7
p_8(x) = 496*x^2+768*x^4+120*x^6+x^8

In this setup the polynomials fit the enumeration of the André numbers better than A094503 as A000111(n) = 1 for n=0,1,2.

The Swiss-Knife polynomials as a generalization of the Seidel-Euler triangle.

Let us now study how the Swiss-Knife polynomials fit into the picture. I have written several times about these amazing polynomials in this blog and for instance on my homepage.

The Swiss-Knife polynomials can be defined by the exponential generating function 2*exp(x*t)/(exp(t)+exp(-t)) and computed with Sage as

def sigma(n, x) : 
    def A(k) : return 0 if k % 4 == 3 else 2^((1-k)//2)*(-1)^((k+1)//4) 
    S = add(A(k)*(-1)^v*binomial(k,v)*(v+x+1)^n for k in (0..n) for v in (0..k))
    return expand(S) 

Next let us build the difference table of the Swiss-Knife polynomials:

We immediately see (for instance from our discussion of the Swiss-Knife polynomials) that the sequence of polynomials in the first row of the table can be simply calculated by:

def tau(n, x) : return sigma(n, x-1)

The exponential generating function of τ_n(x) is 2*exp(t*(x-2))/(1+exp(-2*t)).

g := t -> 2*exp(t*(x-2))/(1+exp(-2*t)):
s := n -> series(g(t),t,n+2): tau := n -> n!*coeff(s(n),t,n):
for i from 0 to 9 do T[i] := unapply(sort(expand(tau(i))),x) od;

A nice feature of this difference table can be seen if evaluated at x = 0 and x = 1. Evaluated at x = 0 we find Seidel's 'Doppel-Treppen-Schema' from the top of this page in a signed version.

1	-1	0	2	0	-16	0
0	-1	2	2	-16	-16
-1	1	4	-14	-32
0	5	-10	-46
5	-5	-56
0	-61
-61

If we evaluate the difference table at x = 1 we get the transpose of the 'Doppel-Treppen-Schema' albeit with a different sign scheme.

1	0	-1	0	5	0	-61
1	-1	-1	5	5	-61
0	-2	4	10	-56
-2	2	14	-46
0	16	-32
16	-16
0

The investigative reader might also whish to know what happens if the difference table is evaluated at x = 1/2 or at x = -1.

The Swiss-Knife median polynomials

The Swiss-Knife median polynomials are given by the main diagonal of the difference table of the Swiss-Knife polynomials.

M(0, x)	1
M(1, x)	x^2-x-1
M(2, x)	x^4-2*x^3-5*x^2+6*x+4
M(3, x)	x^6-3*x^5-12*x^4+29*x^3+57*x^2-72*x-46
M(4, x)	x^8-4*x^7-22*x^6+80*x^5+261*x^4-660*x^3-1264*x^2+1608*x+1024

                          1
                      1, -1, -1
                   1, -2, -5, 6, 4
             1, -3, -12, 29, 57, -72, -46
     1, -4, -22, 80, 261, -660, -1264, 1608, 1024

We note that evaluating the Swiss-Knife median polynomials at x = 1 or at x = 0 gives the signed Euler median numbers A099023 (the unsigned variant is listed in A000657, see the implementation above).

1, -1, 4, -46, 1024, -36976, 1965664, -144361456, 13997185024, ...

A closer look at Seidel's algorithm.

We will first give an overview in four steps, each step slightly more general than the preceding one. The scripts are only fragments which focus on the parts that change. The entire script will be given in the appendix.

The starting point is Seidel's triangle as shown in the facsimile on the top of this page.

The classical Seidel triangle.

The boustrophedon algorithm (Sec-Tan coefficients).

Seidel := proc(    # Input:
  dim::integer,    # length of sequence
)

plow = proc(n)
   if is_even(n) then # right to left backward
      E[n,0] = 0
      ....
   else               # left to right forward
      E[0,n] = 0
      ....
   end

E = array(0..dim-1, 0..dim-1); E[0,0] = 1
for n from 1 to dim-1 do plow(n) od
return E

The boolean Seidel transform.

Recall that the characteristic of the Seidel transform was that the flow of computation is turning like oxen in ploughing. This means you 'plow' across the triangle from right to left, then turn around and 'plow' from left to right, and so on.

The next step introduces a small but crucial change. Now we have the freedom to start a new line of computation at whatever side of the triangle we want. The monotone zig-zag is now replaced by zig-zag-zag or zag-zig-zag-zag or whatever. Oxen cannot fly from one side of the acre to the other but boolean values can be easily changed.

The most simple example for this is to change in the Seidel algorithm the condition 'n mod 2 = 0' by 'n mod k = 0'. And voilá, the classical boustrophedon algorithm has turned into a boolophedon algorithm and the computation of the coefficients of sec + tan is generalized to the computation of the André numbers. This is exactly what we have done above in our implementations of the generalized Euler and André numbers.

More general, in the implementation below the boolean control sequence has become a parameter of the procedure, a function of type integer -> boolean.

Seidel = proc(
  dim::integer,     # length of sequence
  contr::procedure, # integer -> boolean
)

plow = proc(n, b)
   if b then # right to left backward
      E[n,0] = 0
      ....
   else # left to right forward
      E[0,n] = 0
      ....
   end

E = array(0..dim-1, 0..dim-1); E[0,0] = 1
for n from 1 to dim-1 do plow(n, contr(n)) od
return E

As a test case let us compute something bizarre. What is the Seidel triangle of the prime numbers? Here the control function is n -> is_prime(n).

The triangle starts:

[  1   0   1   2   0  11   0 181   0]
[  1   1   1   2  11  11 181 181   ~]
[  0   0   3   9  22 170 362   ~   ~]
[  0   3   6  31 148 532   ~   ~   ~]
[  3   3  37 117 680   ~   ~   ~   ~]
[  0  40  80 797   ~   ~   ~   ~   ~]
[ 40  40 877   ~   ~   ~   ~   ~   ~]
[  0 917   ~   ~   ~   ~   ~   ~   ~]
[917   ~   ~   ~   ~   ~   ~   ~   ~]

Adding the first row and the first column we get

1, 1, 1, 2, 3, 11, 40, 181, 917, 5263, 19144, 54924, 524073, ..

More generally, to every sequence A(n) there is associated the boolean indicator function defined as

 is_an_A(n)  = true if there exists an integer k such that n = A(k), and false otherwise.

Taking this indicator function as the control function of the Seidel triangle the sum of the first row and the first column is the boolean Seidel transform of A. For instance the classical André sequence 1, 1, 1, 2, 5, 16, 61, 272,... is the boolean Seidel transform of the even numbers 0, 2, 4, 6, 8,...

Note that the terms of the sequence A(n) themselves do not enter into the computation of the triangle. This will be different in the next steps of the generalization where the initial values of the first row and the first column are initialized by values passed to the transform.

To sum up: the boolean Seidel transform has the function type (signature)

boolean sequence -> integer sequence

The two basic examples which we have considered were

[n -> n mod 2 = 0]  ->  [sec + tan]
[n -> is_prime(n)]    ->  [ no name ].

Interpreting the control function as the indicator function of a sequence we might equivalently say:

[0, 2, 4, 6, 8,10, ...]  ->  [1, 1, 1, 2, 5, 16, 61, 272,...]
[2, 3, 5, 7, 11, 13,...]  ->  [ 1, 1, 1, 2, 3, 11, 40, 181, ...].

The boolean Seidel transform is the most basic and immediate generalization of the Seidel triangle.

The generalized Seidel transform.

The next step in generalizing Seidel's algorithm is to provide additionally to the control function a sequence for the first term of those rows which are marked 'true' by the control function. These were set to 0 in the original setup by Seidel.

Seidel = proc(
  dim::integer,     # length of sequence
  contr::procedure, # integer -> boolean
  S::list,          # integer list of length dim
)

plow = proc(n, b)
   if b then # right to left backward
      E[n,0] = S[n]
      ....
   else # left to right forward
      E[0,n] = 0
      ....
   end

E = array(0..dim-1, 0..dim-1); E[0,0] = 1
for n from 1 to dim-1 do plow(n, contr(n)) od
return the pair of sequences
   seq(E[i,0] for i in {contr(i) = false})
   seq(E[0,i] for i in {contr(i) = true})

Let us give an example. Let the control sequence be n -> n mod 2 = 0 as in the original Seidel algorithm and S the constant sequence 1 (instead of the constant sequence 0 in the original case). Then we get the Salie numbers A000795 on the right hand side of the triangle.

Because of the symmetry in the computation there is also a second sequence which qualifies as an output: those numbers on the left hand side which are not predetermined by the input sequence. In our example we get the sequence A003719, the odd part of the expansion of tan(x)*cosh(x).

To sum up: the generalized Seidel transform has the function type (signature)

boolean sequence, integer sequence -> integer sequence, integer sequence

The two basic examples which we have considered were

[n -> n mod 2 = 0,  n -> 0]  ->  [sec, tan]
[n -> n mod 2 = 0,  n -> 1]  ->  [cosh*tan, cosh/cos].

Of course the two resulting generating functions can be recombined if desirable. Adding sec and tan leads to the classical André numbers, adding cosh*tan and cosh/cos gives cosh*(sec+tan) and is listed as Frank Ellermann's A062272.

The Seidel operation.

The last step of our generalization is quite obvious. We can also predetermine those initial values on the left hand side of the triangle for which the control sequence is false. An additional input sequence takes care for that. Now the setup is more symmetric.

Seidel := proc(
  dim::integer,     # length of sequence
  contr::procedure, # integer -> boolean
  signed::boolean,  # signed version? yes/no
  S::list,          # integer list of length dim
  T::list           # integer list of length dim
)

plow = proc(n, b)
if b then # right to left backward
   E[n,0] = S[n] 
   ...
else # left to right forward
   E[0,n] = T[n];
   ...

E = array(0..dim-1, 0..dim-1); E[0,0] = 1
for n from 1 to dim-1 do plow(n, contr(n)) od
return the pair of sequences
   seq(E[i,0] for i in {contr(i) = false})
   seq(E[0,i] for i in {contr(i) = true})

As an example we might choose S = T = (-1)^n. This leads to Frank Ellermann's A062162, which is the inverse binomial transform of the classical André numbers by a comment of Paul Barry and has the exponential generating function exp(-x)(sec+tan) as observed by Vladeta Jovovic.  

contr := n -> n mod 2 = 0:
T := [seq((-1)^k,k=1..dim)]:
Seidel(10, contr, T, T):

Inverse binomial transform of the classical André numbers.
[1, 0, 0, 1, 0, 5, 10, 61, 280, 1665]

[1,   -1,    0,   -1,    0,   -1,  10,  -1, 280, -1] 
[0,    1,   -1,    1,   -1,   11,   9, 281, 279,  ~] 
[1,    0,    2,    0,   12,   20, 272, 560,   ~,  ~] 
[1,    2,    2,   12,   32,  252, 832,   ~,   ~,  ~] 
[1,    4,   10,   44,  220, 1084,   ~,   ~,   ~,  ~] 
[5,    6,   54,  176, 1304,    ~,   ~,   ~,   ~,  ~] 
[1,   60,  122, 1480,    ~,    ~,   ~,   ~,   ~,  ~] 
[61,  62, 1602,    ~,    ~,    ~,   ~,   ~,   ~,  ~] 
[1, 1664,    ~,    ~,    ~,    ~,   ~,   ~,   ~,  ~] 
[1665, ~,    ~,    ~,    ~,    ~,   ~,   ~,   ~,  ~]

Note that the Seidel product is not commutative. (This is sometimes expressed somewhat cloudy by saying that there are two kinds of Seidel transforms.) For instance with

contr := n -> n mod 2 = 0:
S := [seq(0,k=1..dim)]:
T := [seq(1,k=1..dim)]:

one finds for Seidel(dim, contr, T, S)

[1, 5, 41, 685]                         [A003719]
[1, 2, 12, 152, 3472]                   [A000795]
[1, 2, 10, 116, 2572]                   // median
[1, 1, 2, 5, 12, 41, 152, 685, 3472]    [A062272]

whereas Seidel(dim, contr, S, T) gives

[2, 6, 52, 896]
[1, 3, 17, 203, 4577]
[1, 2, 12, 150, 3376]  // median
[1, 2, 3, 6, 17, 52, 203, 896, 4577]

To sum up: the Seidel operation has the function type

integer sequence, integer sequence -> integer sequence

Here the control function n -> n mod 2 = 0 is tacitly assumed and the two resulting sequences are combined by Return[k] = E[0, k] if control(k) is true else E[k, 0].

The two basic examples which we have considered were

[n -> 0, n -> 1]  ->  [A062272]
[n -> 1, n -> 0]  ->  [1, 2, 3, 6, 17, 52, 203, 896, 4577,..].

The generalized Seidel operation has the function type

boolean sequence, integer sequence, integer sequence
    -> integer sequence, integer sequence

Of course the two resulting generating functions can be combined if desirable.

Appendix

Maple implementation

# Boolophedon = Generalized Seidel Algorithm
# Copyright(c) 2012 Peter Luschny

Seidel := proc(
   dim::integer,     # length of sequence
   contr::procedure, # f: integer -> boolean
   S::list,          # integer list of length dim
   T::list           # integer list of length dim
)

local E, M, B, U, V, DIM, plow, mat, n, k, r, c;

mat := proc(E)   
# Prints triangle as a matrix.
# For readability only: put `~` in non-used entries.
for n from 0 to DIM do for k from 0 to DIM do
if not type(E[n,k], integer) then E[n,k] := `~` fi od od;
M := convert(E, matrix); print(M);
M end;

plow := proc(n, cond) local k;

if cond then # right to left backward
   E[n,0] := S[n];
   for k from n-1 by -1 to 0 do
       E[k,n-k] := E[k+1,n-k-1] + E[k,n-k-1]
   od
else # left to right forward
   E[0,n] := T[n];
   for k from 1 by 1 to n do
       E[k,n-k] := E[k-1,n-k+1] + E[k-1,n-k]
   od
fi end:

DIM := dim-1; E := array(0..DIM, 0..DIM); E[0,0] := 1;
for n from 1 to DIM do plow(n, contr(n)) od;

B := select(k->evalb(not contr(k)),[$0..DIM]);
r := [seq(E[k,0],k=B)];        # first row   
B := select(k->evalb(contr(k)),[$0..DIM]);
c := [seq(E[0,k],k=B)];        # first column
V := [seq(E[k,k],k=0..DIM/2)]; # median (diagonal)
U := [seq(`if`(contr(k),E[0,k],E[k,0]),k=0..DIM)];

print(r); print(c); print(V); print(U); mat(E);
end:

# Sample Usage:
Seidel(9, n -> n mod 2 = 0, [0$9], [0$9]):
Seidel(9, n -> n mod 2 = 0, [1$9], [1$9]):
Seidel(9, n -> isprime(n),  [0$9], [0$9]):

Sage implementation

def Seidel(
   dim,    # integer, length of sequence
   contr,  # function, integer -> boolean
   S,      # integer list of length dim
   T       # integer list of length dim
):
    E = matrix(ZZ, dim); E[0,0] = 1
        
    def plow(n, cond) :
        if cond : # right to left backward
            E[n,0] = S[n] 
            for k in range(n-1,-1,-1) : 
                E[k,n-k] = E[k+1,n-k-1] + E[k,n-k-1]
        else : # left to right forward
            E[0,n] = T[n]
            for k in range(1,n+1,1) : 
                E[k,n-k] = E[k-1,n-k+1] + E[k-1,n-k]

    for n in (1..dim-1) : plow(n, contr(n))
    
    print [E[k,0] for k in range(dim) if not contr(k)] # first row   
    print [E[0,k] for k in range(dim) if contr(k)]     # first column
    print [E[k,k] for k in range((dim+1)//2)]          # median
    print [E[0,k] if contr(k) else E[k,0] for k in range(dim)] # all

    return E  

# Sample Usage
print Seidel(9, lambda x: x%2 == 0, [0]*9, [0]*9)
print Seidel(9, lambda x: x%2 == 0, [1]*9, [1]*9)
print Seidel(9, is_prime, [0]*9, [0]*9)    

Further references

Seidel's algorithm and its combinatorial interpetation was studied among others by Dumont and Viennot , Dumont , Arnold , Millar, Sloane and Young , and Zeng and Zhou .

Next month we will look at a second algorithm proposed by Seidel in the same paper which can be seen as an optimization of the one given here which can be applied for instance in the computation of the Bernoulli numbers (now here) and the Catalan triangle (now here).

References