User:Peter Luschny/SeidelTransform

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.

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})

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})

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, ~,    ~,    ~,    ~,    ~,   ~,   ~,   ~,  ~]

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

[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]

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

# 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]):

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)