**CONSTRAINED PERMUTATIONS AND THE PRINCIPLE OF INCLUSION-EXCLUSION**

*Philippe Flajolet, January 10, 1998*

This worksheet explores variations of some classical problems in combinatorial analysis, like the
*derangement *
problem (also called
*rencontre*
problem), or the
*menage*
problem. The descriptions of these problems as borrowed from [Comtet, 1974] are:

--
* Derangement/Rencontre*
: If guests to a party leave their hats on hooks in the cloakroom, and grab a hat at good luck when leaving, what is the probability that nobody gets back his own hat? The problem is equivalent to estimating the number of permutations without fixed point, that is to say, without singleton cycles.

--
* Menage*
: What is the number of possible ways one can arrange
married couples (=menages) around a table such that men and women alternate but no woman seats next to her husband?

These problems are in fact permutation enumeration problems from classical combinatorial analysis. The constraints considered concern, for a permutation
, its "
*succession gaps*
", that is, differences between consecutive elements,
, The derangement problem corresponds to permutations such that
, the menage problem to
. In the second case, the constraints on indices and values may also be taken to be modulo
, that is taken "cyclically".

The symbolic method in enumerative combinatorics as well as the
__Combstruct__
package that implements it are based on the concept of decomposability. Combinatorial objects defined by "constraints" are thus not normally accessible to this framework. However, as shown in this worksheet, the enumeration of various types of constrained
*permutations*
can be treated by a combination of
__Combstruct __
,
__Gfun__
, and the Maple system. One may either impose that all such gaps be
*forced*
to belong a finite set
, or dually impose that all gaps have values that

The general idea of counting by inclusion-exclusion is to enumerate by generating functions (GF's) objects with a number of
*distinguished exceptions*
to a set of constraints. The principle is as follows: If
is the bivariate generating function of such objects, where
records size and
records the number of distinguished exceptions of some type
, then inclusion-exclusion [e.g., Comtet, 1974] provides
as the univariate generating function of
-free configurations.

Here, we show cases where so-called "
*permutation templates*
" can be described by finite-state mechanisms, so that a multivariate generating function of templates is directly accessible to combstruct. (A template describes a class of permutations by specifying what happens at some places while others are free, which is rendered by a "dont-care" symbol.) A specialization of the multivariate GF that involves a sign-change (for inclusion-exclusion) and an integral transformation effected on an auxiliary variable (for "filling" the free positions in templates and transforming them into permutations) yields a counting GF for the original problem. Generating functions, recurrences, and numerical values can be obtained automatically in this framework.

This demonstration worksheet was originally inspired by works of Kostas Hatzis (Patras) relative to edge-disjoint paths in random graphs and of Bruno Codenotti (Pisa) relative to computing permanents of circulants matrices.

**References**
.

[Comtet, 1974] : L. Comtet,
*Advanced Combinatorics*
, Reidel, 1974.

[EIS]: N. Sloane and S. Plouffe,
*The Encyclopedia of Integer Sequences*
, Academic Press, 1995.

This Maple worksheet is based on the current versions of the Maple packages
__Combstruct__
and
__Gfun__
* *
(for version V.4) that can be found under

`> `
**restart;**

`> `
**with(combstruct);**

`> `
**with(gfun);**

**1. Permutations with forbidden "position gaps" and the principle of Inclusion-Exclusion**

This section is devoted to the combinatorial basis on which this whole worksheet relies. We define the notion of
__templates__** **
that are simple combinatorial objects out of which constrained permutations can be built. An elementary integral transformation implements the inclusion-exclusion principle and leads to generating functions of constrained permutations.

The set
is a subset of
with
being its maximal element. We propose to count the number
of permutations
such that
* none*
of the position gaps
belongs to
. By inclusion-exclusion, we first enumerate permutations with a certain number of distinguished exceptions.

For this purpose, one needs to determine

= the number of permutations such that of the position gaps are distinguished and forced to belong to , while the remaining positions may be occupied by arbitrary elements, as long as the permutation property is satisfied.

**Templates**
. The diagram above is the beginning of a permutation template (read from top to bottom) such that the sequence of position gaps is

2, 0, *, 0, 2, *, 1.

There, the don't care symbol "*" means that no constraint is imposed at the corresponding position. Thus, any permutation of 1,2,3,etc, that starts like

3, 2, *, 4, 7, *, 8, etc

will match the template. As they specify classes of permutations, the templates are such that there is at most one position occupied (represented by a cross) in each line and each column. Legal templates can thus be viewed as a language over a finite alphabet of dominos where each domino corresponds to an element of . Clearly, the constraints that each column is occupied by at most one cross can be tested with finite memory, i.e., they can be described by a finite state model.

The generating functions that are defined in this way, , are such that

= the coefficient in represents the number of templates (domino placements) of size that have don't-care symbols (hence distinguished occurrences of a position gap lying in ).

Then, by filling the don't-care symbols in templates, one has

.

Finally, by inclusion-exclusion, the number of permutations without any position gap in is

.

This is also expressible as an integral transform

.

since .

**Plotting routines archive**

The diagram above was produced by the following simple routines:

`> `
**domino:=proc(x,y,l) local j; seq([[x+j,y],[x+j+1,y],[x+j+1,y+1],[x+j,y+1],[x+j,y]],j=0..l-1); end;**

`> `
**domino2:=proc(x,y,l,s) domino(x,y,l),[[x+s,y],[x+s+1,y+1]],[[x+s+1,y],[x+s,y+1]]; end;**

`> `
**plot([domino2(1,-1,3,2),domino2(2,-2,3,0),domino(3,-3,3),domino2(4,-4,3,0),domino2(5,-5,3,2),domino(6,-6,3),domino2(7,-7,3,1)],scaling=constrained,axes=none,color=blue,labels=[`A permutation template`,``],thickness=3):**

**2. Templates and permutations with **
**-exceptions (programme)**

Templates are decomposable objects and thus they can be specified in the Combstruct language. As it turns out, they are described by a finite state model that translates all the allowed transitions between adjacent positions. This section provides the main routines that compute the specification of templates associated with a given set of position gaps.

We consider templates with some positions marked that are elements of (exceptions) and "don't-care" positions.

Templates can be described by a finite state device. For instance, if
, we have a linear version of the classical
*menage problem*
: If the choice
of
in
has been made at stage
, that is
has been chosen, then the permutation property implies that one must have
, that is to say, the choice
is forbidden (at stage
) in a template immediately after (at stage
) a choice
. In other words, the language of templates is such that a domino of type "
" cannot be followed by a domino of type "
".

More generally, a
*state*
is any subset
of
whose meaning is that the values in
are unavailable as current position gaps, due to previous "commitments". Only don't cares or values of
in
that are compatible with
can be then be taken, and there is a transition to a new state that records the necessary information about occupied positions. The finite state system thus comprises
states, and the cardinality of the alphabet is equal to the cardinality of
plus one (for the don't care symbol).

This can be expressed by combstruct specifications. Let be a letter that represents a position gap equal to . An template of length is thus described by a word of length over the alphabet formed with the for all in and a don't care symbol.

First build the alphabet:

`> `
**build_alphabet:=proc(Omega) local omega; seq(a[omega]=Atom,omega=Omega),dontcare=Prod(Atom,marked), marked=Epsilon end;**

`> `
**build_alphabet({0,1,3});**

Next, build the transitions. The auxiliary procedure decreases all the elements of a set by .

`> `
**MinusOne:=proc(S) map(proc(x) x-1 end,S) minus {-1} end; **

The transitions (in a raw form) are thus given by the compatibility relations between each state and the symbols.

`> `
**build_transition:=proc(state,Omega) local x,t; if state={} then t:=Epsilon else t:=NULL fi;s[state]=Union(t,Prod(dontcare,s[MinusOne(state)]),seq(Prod(a[x],s[MinusOne(state union {x})]),x=Omega minus state)) end;
**

`> `
**build_transition({},{0,1},1);build_transition({0},{0,1},1);**

The grammar is then obtained by collecting transitions and taking into account the don't care symbol. The initial state is the empty set and it is also the final state (since no domino can "protrude").

`> `
**build_grammar:=proc(Omega) [s[{}],{build_alphabet(Omega)} union map(build_transition,combstruct[allstructs](Subset({$0..max(op(Omega))-1})),Omega),unlabelled] end;**

Here is an example:

`> `
**build_grammar({0,1});**

In general, the grammar comprises nonterminal, which corresponds to a finite-stae device with as many states. The template GF's are thus rational functions of degree , where is the maximal element of .

**3. Special cases of permutations with forbidden "position gaps"**

In this section, we enumerate various classes of permutations associated with a given finite set
of forbidden position gaps. The process is in three steps: (i) Generate the grammar by means of the procedures of Section 2; (ii) compute automatically a bivariate GF of templates by means of
__Combstruct[gfsolve]__
; (iii) transform the corresponding GF by means of the integral transform of Section 1 that implements the inclusion-exclusion principle.All cases can be treated automatically, and we detail here the derangement and menage problems. Several generalizations of the menage problem are also tabulated. The ordinary GF's turn out to have hypergeometric forms and satisfy simple (but combinatorially nonobvious) linear recurrences of the holonomic type that can once more be derived automatically.

**Derangements (or Rencontres)**

A
* derangement*
[Comtet, 1974] is a permutation without fixed point, that is to say, the set
. of position gaps is forbidden. A derangement is thus a permutation without singleton cycles.

`> `
**G0:=build_grammar({0},0);**

The language of templates is thus just the set of all words over the alphabet where is shorthand for a don't-care symbol.

`> `
**draw(G0,size=10);**

`> `
**seq(count(G0,size=n),n=0..10);**

The bivariate GF is obtained by combstruct[gfsolve]:

`> `
**gfsolve(op(2,G0),op(3,G0),z,[[u,marked]]);**

In particular,

`> `
**g0:=subs(",s[{}](z,u));**

The counting sequence is obtained by the transformation corresponding to the inclusion-exclusion principle.

`> `
**series(subs([z=-z,u=-t],g0),z=0,11):map(proc(x) int(x*exp(-t),t=0..infinity) end,");**

The sequence
, is of course classical. It is
**M1937**
of [EIS] where it is described as "subfactorial or rencontre numbers".

The ordinary generating function (OGF) is also accessible in closed-form as

`> `
**P0:=int(exp(-t)*subs([z=-z,u=-t],g0),t=0..infinity);**

This involves the exponential integral , where the exponential integral symbol means

The OGF is purely divergent and P0 is to be interpreted as asymptotic to the right-hand side as tends to .

`> `
**map(simplify,asympt(subs(z=-1/y,P0),y,10));**

**Verification with Combstruct**
. Derangements are characterized by their cycle decomposition as labelled sets of cycles with length at least 1. They are thus specifiable in combstruct.

`> `
**der:=[S,{S=Set(Cycle(Z,card>1))},labelled];**

`> `
**seq(count(der,size=n),n=0..10);**

`> `
**gfsolve(op(2,der),op(3,der),z);**

The exponential generating function (EGF) is well-known

`> `
**der_egf:=simplify(subs(",S(z)));**

`> `
**series(der_egf,z=0,13);**

We check that the coefficients in the EGF and in the OGF are the same. (The two GF's are related by the Laplace transform.)

**Asymptotics**
. The number of derangements of size
satisfies the well-known asymptotic formula:

,

a fact which is obvious given the singularity of the EGF at 1:

`> `
**series(der_egf,z=1,3);**

**Menage numbers**

The
*menage numbers*
are defined by the excluded set
. In nineteenth century terminology, this is phrased as follows: In how many ways can one place couples with men and women alternating in such a way that no husband can seat next to his wife? The problem is usually posed in terms of a round table (=cyclic permutations); see [Comtet, 1974]. In this section, the linear version of the problem is treated.

`> `
**G01:=build_grammar({0,1},1);**

`> `
**g01:=subs(gfsolve(op(2,G01),op(3,G01),z,[[u,marked]]),s[{}](z,u));**

The linear "menage numbers" immediately result by the general transformation.

`> `
**series(subs([z=-z,u=-t],g01),z=0,17):ser01:=map(proc(x) int(x*exp(-t),t=0..infinity) end,");**

This is sequence
**M3020**
of [
*EIS*
], defined there as "sums of menage numbers".

**Checking**
. Combstruct is also useful for checking and debugging purposes. Consistency of templates can be checked for small values of
by expansions and by
* combstruct[allstructs]*
.

`> `
**series(g01,z=0,10):map(expand,");**

`> `
**allstructs(G01,size=2);allstructs(G01,size=3);**

**Holonomic forms**
. A recurrence can be heuristically guessed by
, proved by the hypergeometric method detailed here (this is however specific to the problem at hand), and also obtained automatically by
. Problems of this type a priori belong to the so-called
* holonomic class*
, introduced by Zeilberger, that is to say, sequences that satisfy linear recurrences with polynomial coefficients. Alternatively, the GF's satisfy linear differential equations with polynomial coefficients.

As a first approach, the procedure
__gfun[listtodiffeq]__
provides a plausible differential equation that accounts for a given number sequence.

`> `
**l01:=[seq(coeff(ser01,z,i),i=0..16)];**

`> `
**ode01:=listtodiffeq(l01,Y(z));**

This suggests a formal representation of the OGF of the basic menage problem.

`> `
**dsolve(op(1,op(1,ode01)),Y(z));**

This guess can be validated by direct integration:

`> `
**P01:=int(subs([z=-z,u=-t],g01)*exp(-t),t=0..infinity);**

**Hypergeometrics**
. Alternatively any expression with the exponential integral
can be expressed in terms of a divergent hypergeometrics (that is also the OGF of permutations). One has

`> `
**y*exp(y)*Ei(1,y)=asympt(y*exp(y)*Ei(1,y),y,10);**

that is also expressible as a -hypergeometric:

`> `
**hypergeom([1,1],[],z)=series(hypergeom([1,1],[],z),z=0,12);**

Thus the menage OGF has an hypergeometric form:

`> `
**hg_menage:=1/(1+z)*hypergeom([1,1],[],z/(1+z)^2);**

`> `
**series(hg_menage,z=0,10);**

This gives a single alternating sum for the coefficients that is the counterpart of the one of [Comtet, 1974] for the circular menage problem.

This process gives rise to a general procedure for conversion from exponential integral to hypergeometric form:

`> `
**convert_hypergeom:=proc(e) subs(Ei=proc(a,b) if a<>1 then ERROR(`unable to convert`) else exp(-b)/b*hypergeom([1,1],[],-1/b) fi end,e); simplify("); end;**

`> `
**convert_hypergeom(P01);**

**Asymptotics**
. Experimentally, we have

`> `
**seq(op(i,l01)/(i-1)!*1.0,i=1..16);**

`> `
**exp(-2)=exp(-2.);**

Now, two values are forbidden at each position. Empirically, the number of menage numbers is seen to satisfy

.

This is to be compared to the asymptotics of derangements. The asymptotic estimate can be established directly from the alternating sum expression of coefficients.

**A simplified menage problem**

The case is in fact a variant of the derangement problem.

`> `
**G1:=build_grammar({1},1):**

`> `
**g1:=subs(gfsolve(op(2,G1),op(3,G1),z,[[u,marked]]),s[{}](z,u));**

`> `
**series(subs([z=-z,u=-t],g1),z=0,11):map(proc(x) int(x*exp(-t),t=0..infinity) end,");**

The sequence appears as
**M2905**
in [EIS] where the EGF is given:

`> `
**exp(-z)/(1-z)^2=series(exp(-z)/(1-z)^2,z=0,8);**

**2-Menage numbers**

This is a generalization of the menage problem where : In order to avoid family fights, couples are now at distance at least 3 (!).

`> `
**G012:=build_grammar({0,1,2},2);**

`> `
**gfsolve(op(2,G012),op(3,G012),z,[[u,marked]]); g012:=subs(",s[{}](z,u));**

`> `
**series(subs([z=-z,u=-t],g012),z=0,13):P012_ser:=map(proc(x) int(x*exp(-t),t=0..infinity) end,");**

Only the beginning of this sequence (till
included) appears in [EIS] as sequence
**M3970. **
The computation above makes it possible to extend the values of [EIS] to an arbitrary order.

**Holonomic forms**
. Since the denominator of
has degree 1 in
, the Laplace-like integral applied to it must be expressible in terms of the exponential integral
. This, in turn is equivalent to a divergent hypergeometric form.

`> `
**convert(subs([z=-z,u=-t],g012),parfrac,u);h012:=int("*exp(-t),t=0..infinity);**

`> `
**hh012:=convert_hypergeom(h012);**

`> `
**series(",z=0,12);**

One can then get the holonomic forms by .

`> `
**holexprtodiffeq(hh012,y(z)): subs(_C[0]=1,");**

`> `
**diffeqtorec(",y(z),u(n));**

This gives a fast procedure for computing the numbers.

`> `
**s012:=rectoproc(",u(n),remember);**

`> `
**seq(s012(j),j=0..30);**

**Asymptotics**
. This sequence obeys a general asymptotic pattern that was alluded to before. Verification to high orders is easy thanks to the holonomic recurrences.

`> `
**seq(1.0*s012(m)/m!,m=0..50);**

Such a numerical computation supports the assumption that asymptotically, the number of such permutations should grows like

.

In effect, one has:

`> `
**exp(-3)=exp(-3.0);**

(See the conclusion section for a proof of this observation.)

**Variants of 2-menage numbers**

The case is close to the classical menage problem.

`> `
**G12:=build_grammar({1,2},2):gfsolve(op(2,G12),op(3,G12),z,[[u,marked]]): g12:=subs(",s[{}](z,u));**

`> `
**series(subs([z=-z,u=-t],g12),z=0,11):map(proc(x) int(x*exp(-t),t=0..infinity) end,");**

**3-Menage numbers**

The basic version is defined by .

`> `
**G0123:=build_grammar({0,1,2,3},3):**

`> `
**gfsolve(op(2,G0123),op(3,G0123),z,[[u,marked]]): g0123:=subs(",s[{}](z,u));**

`> `
**series(subs([z=-z,u=-t],g0123),z=0,11):map(proc(x) int(x*exp(-t),t=0..infinity) end,");**

This sequence is not to be found in [EIS].

**Variants of the 3-Menage problem**

We just offer here a few random samples of sequences.

First, the case .

`> `
**G03:=build_grammar({0,3},3):**

`> `
**gfsolve(op(2,G03),op(3,G03),z,[[u,marked]]): g03:=subs(",s[{}](z,u));**

`> `
**series(subs([z=-z,u=-t],g03),z=0,11):map(proc(x) int(x*exp(-t),t=0..infinity) end,");**

Next, the case .

`> `
**G013:=build_grammar({0,1,3},3):gfsolve(op(2,"),op(3,"),z,[[u,marked]]): g013:=subs(",s[{}](z,u));series(subs([z=-z,u=-t],g013),z=0,11):map(proc(x) int(x*exp(-t),t=0..infinity) end,");**

And the case .

`> `
**G023:=build_grammar({0,2,3},3):gfsolve(op(2,"),op(3,"),z,[[u,marked]]): g023:=subs(",s[{}](z,u));series(subs([z=-z,u=-t],g023),z=0,11):map(proc(x) int(x*exp(-t),t=0..infinity) end,");**

Finally, the case .

`> `
**G13:=build_grammar({1,3},3):gfsolve(op(2,"),op(3,"),z,[[u,marked]]): g13:=subs(",s[{}](z,u));series(subs([z=-z,u=-t],g13),z=0,11):map(proc(x) int(x*exp(-t),t=0..infinity) end,");**

**4-Menage numbers**

This and the next example serve to test the system on models of a fairly large size. The rational GF for the 4-menage proble is of degree 16, and the computations only take a few seconds.

`> `
**G01234:=build_grammar({0,1,2,3,4},4):gfsolve(op(2,G01234),op(3,G01234),z,[[u,marked]]):**

`> `
**g01234:=subs(",s[{}](z,u));
**

`> `
**series(subs([z=-z,u=-t],g01234),z=0,11):map(proc(x) int(x*exp(-t),t=0..infinity) end,");**

`> `
**convert(g01234,parfrac,u);**

**5-Menage numbers**

The computation time is now of the order of a minute (on a slow machine with about cycles per second). The denominator is now of degree 32 in and of degree 6 in .

`> `
**G012345:=build_grammar({0,1,2,3,4,5},5): gfsolve(op(2,G012345),op(3,G012345),z,[[u,marked]]):**

`> `
**g012345:=subs(",s[{}](z,u));
**

`> `
**series(subs([z=-z,u=-t],g012345),z=0,13):map(proc(x) int(x*exp(-t),t=0..infinity) end,");**

`> `
**denom(g012345):series(",u,8):**

`> `
**map(factor,");**

Only partial patterns are apparent.

**4. Cyclically constrained position gaps**

The problems we consider in this section are of so-called "cyclic" or "toroidal type". In the menage formulation, one may for instace have a round table (the classical formulation of the menage problem is of this type). Cyclic templates are in fact easily derivable from the linear ones that were considered earlier. This gives directly the enumeration of permutations with forced position gaps (in the cyclic sense) and the generating functions are, like for templates, rational functions. The inclusion-exclusion argument still aplies in the form already encountered and this provides the enumeration of permutations with forbidden position gaps (in the cyclic sense). In this way, we solve automatically the classical menage problem and its generalizations while completing in passing some of the entries of the
*Encyclopedia of Integer Sequences*
[EIS].

**Forced position gaps**

One first seeks to enumerate the permutations such that all position gaps are
*forced*
to lie in some fixed set
. For
a subset of the nonnegative integers, the problem only makes sense if one considers the constraints to be "

The system of transitions is the same as before. However, in the graph of transitions, one must restrict attention to paths such that the initial and the final states are the same: this corresponds to the fact that constraints from the end of the permutations are propagated back to the beginning. In this basic enumeration problem, no inclusion-exclusion argument is needed, and the problem is directly modelled by a finite-state device related to the earlier ones.

`> `
**build_transition:=proc(state,Omega,initial_state) local x,t; if state=initial_state then t:=Epsilon else t:=NULL fi;s[state]=Union(t,Prod(dontcare,s[MinusOne(state)]),seq(Prod(a[x],s[MinusOne(state union {x})]),x=Omega minus state)) end;**

`> `
**build_transition({1,2},{0,1,2,3},{});**

The following procedure builds the specification that corresponds to a loop that starts and ends with in the transition graph. It is a simple modification of the previous procedure .

`> `
**build_loop:=proc(Omega,initial_state) [s[{}],{build_alphabet(Omega)} union map(build_transition,combstruct[allstructs](Subset({$0..max(op(Omega))-1})),Omega,initial_state),unlabelled] end;**

`> `
**build_loop({0},{});build_loop({0,1},{0});**

We could construct a composite grammar for all loops, it is however simpler to solve for the GF of each loop.

`> `
**gf_loops:=proc(Omega) local is;
normal(add(subs(gfsolve(op(2,build_loop(Omega,is)),unlabelled,z,[[u,marked]]),s[is](z,u)),is=combstruct[allstructs](Subset({$0..max(op(Omega))-1})))) end;**

The bivariate GF counts configurations with possible exceptions (it is useful for the analysis of configurations with forbidden gaps).

In the case where only singletons are allowed, then there is only 1 possibility for each size .

`> `
**l0:=gf_loops({0},0); subs(u=0,");series(",z=0,10);**

In the case , the first element (either 0 or 1) conditions all the rest.

`> `
**l01:=gf_loops({0,1},1);subs(u=0,");series(",z=0,10);**

The case is the first one that is nontrivial.

`> `
**l012:=gf_loops({0,1,2},2);subs(u=0,");series(",z=0,10);**

The sequence
,etc, is sequence
**M2396**
of [
*EIS*
].

`> `
**l0123:=gf_loops({0,1,2,3},3);subs(u=0,");series(",z=0,12);**

**Generating functions**
. These are clearly rational. Let
be the GF that arises in each case. Then, a general result concerning generating functions of loops in graphs [Biggs,
*Algebraic Graph Theory*
, 1993] implies that
, where
is the dimension of the state space, is a logarithmic derivative,

For instance:

`> `
**L0123:=subs(u=0,l0123);**

`> `
**int((L0123-8)/z,z);**

This is also true of the bivariate GF:

`> `
**int((l0123-8)/z,z);**

**Recurrences**
. The GF's for forced position gaps are always
*rational*
. Linear recurrences (with constant coefficients) are then available by gfun.

`> `
**series(L0123,z=0,12);**

`> `
**holexprtodiffeq(L0123,y(z));**

`> `
**recl0123:=diffeqtorec(",y(z),u(n));**

High index values can be obtained quickly by the binary powering method that is implemented in
__gfun[rectoproc]__
. The recurrence needs to be made homogeneous, first.

`> `
**op(1,recl0123)-subs(n=n+1,op(1,recl0123));**

`> `
**sl0123:=rectoproc({",u(0)=8,u(1)=4,u(2)=8,u(3)=16},u(n));**

`> `
**sl0123(1000);**

The complexity of this computation is thus of O(log(n)) arithmetic operations for the th term:

`> `
**for m in [1000,2000,4000,8000,16000,32000] do st:=time():sl0123(m):print(time()-st); od:**

The computation of the value of index 32000 takes typically 1 second of computer time though the number has 8470 digits:

`> `
**sl0123(32000)*1.0;**

(Such fast algorithms are useful for instance when analyzing patterns in strings, like in DNA sequences.)

**Forbidden position gaps**

We are dealing here with menage problems in the classical sense of a circular table. The same integral transformation as in the noncyclic/nontoroidal case works for obtaining the GF of permutations with excluded patterns,
*i.e*
., no position gap in
.

The case
gives back the derangement numbers (again, this is
**M1937**
of [EIS]).

`> `
**series(subs([z=-z,u=-t],l0),z=0,13):map(proc(x) int(x*exp(-t),t=0..infinity) end,");**

Next comes the case .

`> `
**series(subs([z=-z,u=-t],l01),z=0,13):map(proc(x) int(x*exp(-t),t=0..infinity) end,");**

This, apart from the first two values (see the next section for possible adjustments) is exactly the sequence of
*Menage Numbers*
that is listed as
**M2062 **
in [EIS].

`> `
**int(x*exp(-t)*subs([z=-z,u=-t],l01-2),t=0..infinity);**

`> `
**convert_hypergeom(");**

This gives the ordinary GF for the menage problem. This GF is closely related to that of the linear menage problem.

The generalized menage numbers that correspond to
appear as
**M2121**
in [EIS] under the name "discordant permutations of length
".

`> `
**series(subs([z=-z,u=-t],l012),z=0,13):map(proc(x) int(x*exp(-t),t=0..infinity) end,");**

The next sequence called "hit polynomials" is given as
**M2154**
in [EIS], where it stops at the value 369321.

`> `
**series(subs([z=-z,u=-t],l0123),z=0,13):map(proc(x) int(x*exp(-t),t=0..infinity) end,");**

It is possible to complete the table to arbitrary orders, obtaining exact forms of the GF, recurrences, etc. For instance:

`> `
**int(subs([z=-z,u=-t],l0123)*exp(-t),t=0..infinity):**

`> `
**hl0123:=convert_hypergeom(");**

`> `
**series(hl0123,z=0,33);**

**Conclusions**

The combinatorial packages may be used to facilitate the enumeration of permutations with contraints on "positions". In essence, a finite-state model is solved that leads to a bivariate generating function and the GF's of interest appear to be integral transforms of . In particular these GF are all holonomic. Thus, the corresponding counting sequences can be determined in arithmetic operations.

If partial fraction decompositions are effected on
, then the GF's appear to be systematically expressible as compositions of the divergent hypergeometric series
with algebraic functions. The simplest of these forms have been made explicit above in some cases, like derangements, menage, and successions. The same method also allows for counting
-exceptions in
*all *
permutations by means of bivariate generating functions. A typical problem in this range is the statistics of the number of "conflicts" in a random assignment of seats in the menage problem.

The fast algorithms make it possible to perform numerical observations and conjecture various asymptotic patterns. One of these, which is confirmed in the particular cases treated above is [see, I. Vardi,
*Computational Recreations with Mathematica*
, Addison-Wesley, 1991; p. 123]:

**Property**
.
*The fraction of permutations with gaps that do not belong to a finite set *
* is asymptotic to*

* *
*,*

*where *
*. More generally, the probability that a random permutation of size *
* has *
* exceptions of type *
* (with *
* fixed) is given by a Poisson law of parameter *
*,*

* *
*.*

A weaker form of this property (for position gaps where
is an initial segment of the integers) is established by probabilistic arguments in [Barbour, Holst, Svanson,
*Poisson Approximation*
, Oxford, 1992].

The approach via hypergeometrics points to a unified derivation via generating functions and divergent series that encompasses more general constraints (like succession gaps and the case where is not an initial segment of the integers). The principle of the generating function method is that

,

provided the argument of the hypergeometric is a function that is analytic at the origin. For instance, in the case of the menage problem, one has

`> `
**P01=convert_hypergeom(P01);**

`> `
**series(1/(1+2*z+z^2),z=0,6);**

so that the asymptotic proportion of menage permutations is
*provedly*
equal to
. The method works for any of the examples discussed in this worksheet.

A companion worksheet shows that a similar approach applies to "succession gaps", by which we means constraints on the values of successive elements in a permutation. The corresponding solutions serve in the analysis of multiconnectivity in random interconnection graphs.