%invmaj.tex: The Number of Inversions and the Major Index of Permutations are Asymptotically Joint-Independently Normal
%%a Plain TeX file by Andrew Baxter and Doron Zeilberger (xx pages)

%begin macros
\def\L{{\cal L}}
\def\Z{{\cal Z}}
\def\N{{\cal N}}
\baselineskip=14pt
\parskip=10pt
\def\halmos{\hbox{\vrule height0.15cm width0.01cm\vbox{\hrule height
 0.01cm width0.2cm \vskip0.15cm \hrule height 0.01cm width0.2cm}\vrule
 height0.15cm width 0.01cm}}
\font\eightrm=cmr8 \font\sixrm=cmr6
\font\eighttt=cmtt8
\magnification=\magstephalf
\def\P{{\cal P}}
\def\Q{{\cal Q}}
\def\1{{\overline{1}}}
\def\2{{\overline{2}}}
\parindent=0pt
\overfullrule=0in
\def\Tilde{\char126\relax}
\def\frac#1#2{{#1 \over #2}}
%\headline={\rm  \ifodd\pageno  \RightHead  \else  \LeftHead  \fi}
%\def\RightHead{\centerline{
%Title
%}}
%\def\LeftHead{ \centerline{Doron Zeilberger}}
%end macros
\bf
\centerline
{
The Number of Inversions and the Major Index of Permutations
}
\centerline
{
are Asymptotically Joint-Independently-Normal
}
\rm
\bigskip
\centerline{
{\it Andrew BAXTER}{$^1$} and {\it Doron ZEILBERGER}\footnote{$^1$}
{\eightrm  \raggedright
Department of Mathematics, Rutgers University (New Brunswick),
Hill Center-Busch Campus, 110 Frelinghuysen Rd., Piscataway,
NJ 08854-8019, USA.
%\break
{\eighttt [baxter,zeilberg]  at math dot rutgers dot edu} ,
\hfill \break
{\eighttt http://www.math.rutgers.edu/\~{}[baxter,zeilberg]} .
April 5, 2010.
Accompanied by the Maple package {\eighttt InvMaj}
downloadable from
\hfill\break
{\eighttt http://www.math.rutgers.edu/\~{}zeilberg/mamarim/mamarimhtml/invmaj.html} .
\hfill\break
The work of both authors was supported in part by the USA National Science Foundation.
}
}

{\bf Abstract:} We use {\it recurrences} (alias {\it difference equations}) to prove the
longstanding conjecture that the two most important {\it permutation statistics},
namely the number of inversions and the major index, are asymptotically joint-independently-normal.
We even derive  more-precise-than-needed asymptotic formulas for the (normalized) mixed moments.

{\bf Human Statistics}

{\it Human} statistics are numerical attributes defined on {\it humans}, for example,
{\it longevity}, {\it height}, {\it weight}, {\it IQ}, and it is well-known, at least empirically, that
these are, each separately, {\it asymptotically normal}, which means that if you draw a histogram
with the statistical data, it would look like a {\it bell-curve}. It is also true that
they are usually joint-asymptotically-normal, but usually not {\it independently} so.
But if you compute empirically the {\it correlation matrix}, you would get, asymptotically
(i.e. for ``large'' populations) that they are close to being distributed according
to a multivariate (generalized) {\it Gaussian} $exp (-Q(x_1, x_2, ...))$ with $Q(x_1,x_2, \dots)$ a
certain {\it quadratic form} that can be deduced from the correlation matrix.


{\bf Permutation Statistics}

Let our population be the set of {\it permutations} of $\{1,2, \dots, n \}$.
They too, can be assigned {\it numerical attributes}, and the great
classical combinatorialist, Dominique Foata, (who got his 3rd-cycle doctorate in statistics!)
coined the term {\it permutation statistics} for them.

The most important permutation statistic is the {\bf number of inversions}, $inv(\pi)$,
that counts the number of pairs $1 \leq i < j \leq n$ such that $\pi_i>\pi_j$
(and ranges from $0$ to $n(n-1)/2$). For example, $inv(314625)=5$,
corresponding to the set of pairs $\{(1,2),(1,5),(3,5),(4,5),(4,6)\}$.
It features in the definition of the determinant, and Netto
proved that their {\it probability generating function} 
(the polynomial in $q$ such that its coefficient of $q^i$ is the probability that
a uniformly-at-random $n$-permutation has $i$ inversions) is given by
$$
{
{(1)(1+q)(1+q+q^2) \cdots (1+q+q^2+ \dots + q^{n-1} )}
\over
{n!}
}=
{
{\prod_{i=1}^{n} (1-q^i)}
\over
{n!(1-q)^n}
} \quad .
$$


The {\it second most} important permutation statistic is the {\it major index}, $maj(\pi)$, that is the
sum of the places $i$, where $\pi_{i}>\pi_{i+1}$
For example, $maj(314625)=1+4=5$, because at $i=1$ and $i=4$ we have  descents.
Major Percy Alexander MacMahon[M] famously proved that the probability generating function
for the major index is also given by that very same formula.
In other words the permutation statistics $inv$ and $maj$ are {\it equidistributed}.
Dominique Foata[Fo] gave a seminal lovely {\it bijective} proof that proved the stronger statement
that $inv$ and $maj$ are equidistributed also when restricted to permutations ending with a given integer.

William Feller ([Fe], p.257) proved that the number of inversions (and hence also the major index) is
{\it asymptotically normal} in the following sense.  Feller easily computed the {\it expectation},
$$
E[inv]=m_n=n(n-1)/4 \quad ,
$$
and the {\it variance},
$$
(\sigma_n)^2={{2n^3+3n^2-5n} \over {72}} \quad .
$$
If we denote by $\Z_n$ the {\it centralized} and {\it normalized} random variable
$$
\Z_n ={{inv-m_n} \over {\sigma_n}} \quad ,
$$
then $\Z_n \rightarrow \N$ in distribution as $n \rightarrow \infty$,
where $\N$ is the Gaussian distribution whose {\it probability density function} is
$e^{-x^2/2}/\sqrt{2\pi}$.


A {\it computer-generated} proof, which gives much more detail regarding the precise asymptotics, can be obtained using Zeilberger's 
Maple package {\tt http://www.math.rutgers.edu/\Tilde zeilberg/tokhniot/AsymptoticMoments}
that accompanies the article [Z].

Hence both $inv$ and $maj$ are {\it individually} asymptotically normal, but what about
their {\it interaction}, in other words, what can you say about the
limit of 
$$
\frac{1}{da \, db} Pr( \sigma_n a  \leq inv(\pi)-m_n \leq  \sigma_n (a+da) \quad AND \quad
 \sigma_n b  \leq maj(\pi)-m_n \leq  \sigma_n (b+db) \quad) \quad ,
$$
as $n \rightarrow \infty$ and $da,db \rightarrow 0$?

In this article, we prove that this limit exists and equals 
$(2 \pi)^{-1} e^{-a^2/2-b^2/2}$. In  other words,
$inv$ and $maj$ are asymptotically {\it joint-independently-normal}.


{\bf A Brief History}

It all started when the great Swedish probabilist Svante Janson (member of the Swedish Academy of Science,
that awards the Nobel prizes) asked Donald Knuth (one of the greatest computer scientists of all time,
winner of the Turing and Kyoto prizes, among many other honors) about
the {\it asymptotic covariance} of $inv$ and $maj$. Neither of these luminaries knew the answer,
so Don Knuth asked one of us (DZ). DZ didn't know the answer either, so he asked his
beloved servant, Shalosh B. Ekhad, who {\it immediately} ([E]) produced, not just the asymptotics,
but the {\it exact answer}!. It turned out to be $n(n-1)/8$. In particular, the
{\it correlation}, $Cov(inv,maj)/\sigma_n^2= \frac{\frac{n(n-1)}{8}}{\frac{2n^3+3n^2-5n}{72}}=\frac{9}{2n}+O(1/n^2)$
tends to zero as $n$ goes to infinity. It followed that for large $n$, $inv$ and $maj$ are practically {\it uncorrelated}.

But there are lots of pairs of random variables that are uncorrelated yet not independent. 
A convenient way to 
prove that $X_n:=(inv-m_n)/\sigma_n$ and $Y_n:=(maj-m_n)/\sigma_n$ are asymptotically independent (we already know that they
are both normal) is to use the {method of moments}, and to prove that the {\it mixed moments}
$$
M_{r,s}(n):=E[(X_n)^{r} (Y_n)^{s}] \quad,
$$
tend to the mixed moments of $\N \times \N$ as $n \rightarrow \infty$.  In other words, for $r,s \geq 1$:
$$
\lim_{n \rightarrow \infty} M_{2r,2s}(n)={{(2r)!} \over {2^r r!}}{{(2s)!} \over {2^s s!}} \quad ,
\eqno(EE)
$$
$$
\lim_{n \rightarrow \infty} M_{2r-1,2s}(n) \, = \, 0 \quad ,
\eqno(OE)
$$
$$
\lim_{n \rightarrow \infty} M_{2r,2s-1}(n) \, = \, 0 \quad ,
\eqno(EO)
$$
$$
\lim_{n \rightarrow \infty} M_{2r-1,2s-1}(n) \, = \, 0 \quad .
\eqno(OO)
$$

Ekhad's brilliant approach to the Janson-Knuth question merely settled the case $r=1,s=1$ of
$(OO)$. Of course, because of symmetry $(OE)$ and $(EO)$ are trivially true (before taking the limits!,
i.e. $M_{2r,2s-1}(n) \equiv 0 $ and $M_{2r-1,2s}(n) \equiv 0$).

One natural approach would be to extend Ekhad's brilliant answer for $M_{1,1}(n)$ to the general case,
and try to derive closed-form expressions for $M_{r,s}(n)$ for larger $r$ and $s$.
Since Ekhad's proof is so brief, we can cite it here in full.


``Svante Janson asked Don Knuth,
who asked me, about the covariance of $inv$ and $maj$.
The answer is ${{n} \choose {2}}/4$. To prove it, I asked
Shalosh to compute the average of the quantity 
$(inv ( \pi) - E(inv))(maj( \pi) - E(maj)) $ over all permutations of
a given length $n$, and it gave me,
for $n=1,2,3,4,5$, the values $0,1/4,3/4,3/2,5/2$, respectively. Since
we know a priori\footnote{$^2$}
{\eightrm This is the old trick to compute moments
of combinatorial `statistics', described nicely [GKP], section 8.2,
by changing the order of summation. It applies equally well
to covariance. 
Rather than actually carrying out the
gory details, we observe that this is always a polynomial whose
degree is trivial to bound.}
that this is a polynomial of degree $\leq 4$, this must be it! \halmos''.

Obviously this brute-brute-force approach would be hopeless
for deriving polynomial expressions for the moments $M_{r,s}(n)$ for
larger $r$ and $s$. As we will soon see, the degree of
the polynomial $M_{2r,2s}(n)$ is $3(r+s)$, so for example,
in order to (rigorously) guess  $M_{10,10}(n)$,
we would need $31$ data points, requiring the computer to
examine more than $31! > 0.822 \cdot 10^{34}$ permutations.

It is perhaps surprising, then, that an {\it  inspired}, ``empirical'' (yet {\bf fully rigorous}),
``brute-force'' approach {\it does} indeed work.
The first step is to have a more efficient way to
compute the moments $M_{r,s}(n)$, for {\it specific} $n$ and {\it specific}
$r$ and $s$. We will do this by first designing an efficient way
to generate the {\it probability generating function}, let's call it $G(n)(p,q)$, for
the pair of statistics $(inv, maj)$.  Above we see Netto's  beautiful
closed-form expressions for $G(n)(1,q)$ and $G(n)(p,1)$, but no such closed-form expression seems to exist for
the bi-variate generating function, so the best that we
can hope for is to find a {\it recurrence scheme}.

{\bf A Combinatorial Interlude}

Setting aside probability for a few moments, and we focus on a fast
algorithm for computing
$$
H(n)(p,q):=\sum_{\pi \in S_n} p^{inv(\pi)}q^{maj(\pi)} \quad ,
$$
for $n$ up to, say, $50$.

Define the {\it weight} of a permutation $\pi$ to be $p^{inv(\pi)}q^{maj(\pi)}$.
Suppose that $\pi \in S_n$ ends with $i$, so we can write
$\pi=\pi' i$, where $\pi'$ is a permutation of 
$\{1, \dots , i-1, i+1, \dots n \}$.

When you chop-off $i$ from $\pi$ you always lose $n-i$ inversions (that is, $\pi'$ has $n-i$ fewer inversions than $\pi$).  
The major index, however, decreases by $n-1$ if the last letter of $\pi'$, let's call it $j$, is larger than $i$, and if $j<i$
 the major index does not decrease at all.
So writing $\pi=\pi'' j i$, we have
$$
inv(\pi'' ji)=inv(\pi''j)+n-i \quad,
$$
$$
maj(\pi'' ji) = \cases{
maj(\pi''j) ,& if $j<i$ ;\cr
maj(\pi''j)+n-1 ,&  if $j>i$.\cr}   
$$
Combining, we have
$$
weight(\pi'' ji) = \cases{
p^{n-i}weight(\pi''j) ,& if $j<i$ ;\cr
p^{n-i}q^{n-1} weight(\pi''j),&  if $j>i$.\cr}  
$$
So in order to compute $H(n)(p,q)$, we need to introduce the
more general weight-enumerators of those permutations in $S_n$
that end with an $i$. Let's call these $F(n,i)(p,q)$.
In symbols:
$$
F(n,i)(p,q):=\sum_
{
{{\pi \in S_n} \atop {\pi_n=i}}
} p^{inv(\pi)}q^{maj(\pi)} \quad .
$$
It follows that (let's omit the arguments $(p,q)$ from now on):
$$
F(n,i)=
p^{n-i} \sum_{j=1}^{i-1} F(n-1,j)+p^{n-i}q^{n-1} \sum_{j=i}^{n-1} F(n-1,j) 
\quad .
\eqno(Fni)
$$
Note that when we chop off the last entry, $i$, from
$\pi=\pi'' j i$, we see $\pi'' j$ is a permutation of
$\{1, \dots, i-1, i+1, n \}$.  We then ``reduce''
$\pi'' j$ to a permutation of $\{1, \dots, n-1\}$ by
reducing all the entries larger than $i$ by $1$, and so
the summation ranges from $j=i$ to $j=n-1$ rather than
from $j=i+1$ to $j=n$.

Replacing $i$ by $i+1$ in the above equation, we have:
$$
F(n,i+1)=
p^{n-i-1} \sum_{j=1}^{i} F(n-1,j)+p^{n-i-1}q^{n-1} 
\sum_{j=i+1}^{n-1} F(n-1,j) \quad .
$$
Subtracting the former equation from $p$ times the latter we get
$$
F(n,i)-pF(n,i+1)=
$$
$$
p^{n-i} \sum_{j=1}^{i-1} F(n-1,j)+p^{n-i}q^{n-1} \sum_{j=i}^{n-1} F(n-1,j) 
$$
$$
-p^{n-i} \sum_{j=1}^{i} F(n-1,j)-p^{n-i}q^{n-1} \sum_{j=i+1}^{n-1} F(n-1,j)
$$
$$
=p^{n-i} \left ( \, \sum_{j=1}^{i-1} F(n-1,j) - \sum_{j=1}^{i} F(n-1,j) \, \right )
+p^{n-i}q^{n-1} \left ( \, \sum_{j=i}^{n-1} F(n-1,j)  - \sum_{j=i+1}^{n-1} F(n-1,j) \,  \right )
$$
$$
=-p^{n-i}F(n-1,i)+p^{n-i}q^{n-1}F(n-1,i)=p^{n-i}(q^{n-1}-1)F(n-1,i) \quad.
$$
Rearranging, we get:
$$
F(n,i)=pF(n,i+1)+p^{n-i}(q^{n-1}-1)F(n-1,i) \quad .
\eqno(RecF)
$$
This enables us to compute $F(n,i)$ for $i<n$ starting
with $F(n,n)$. We still need to specify $F(n,n)$, and for this
we do need the $\sum$ symbol, namely we use Eq. $(Fni)$ with
$i=n$:
$$
F(n,n)= \sum_{j=1}^{n-1} F(n-1,j)
\quad .
\eqno(Fnn)
$$
The recurrence $(RecF)$ together with the final condition $(Fnn)$,
and the trivial {\it initial condition} $F(1,i)=\delta_{i,1}$,
enables us to efficiently compute $F(n,i)$ for numeric $(n,i)$, 
for $\{(n,i) \vert 1 \leq i \leq n \leq N\}$ for any finite $N$ (not too large, but not too small either:
e.g., $N=100$ is still plausible).
In particular, we can compile a table of $H(n)(p,q)=F(n+1,n+1)(p,q)$, for $n \leq N-1$.

{\bf A crash course in multivariable enumerative probability}

Suppose that you have a finite set of objects $S$ and
several {\it statistics} $f_1(s), \dots , f_r(s)$.
The {\it multi-variable generating function}
(weight-enumerator under the weight $x_1^{f_1(s)} \cdots x_r^{f_r(s)}$)
is defined to be:
$$
\sum_{s \in S} x_1^{f_1(s)} \cdots x_r^{f_r(s)} \quad .
$$
Suppose that you pick an element $s \in S$ {\it uniformly at random}
and you want the multivariable generating function such that
the coefficient of $x_1^{a_1} \cdots x_r^{a_r}$ would give you the
probability that $f_1(s)=a_1, \dots, f_{r-1}(s)=a_{r-1}, and f_r(s)=a_r$.
It is given by:
$$
P(x_1, \dots, x_r)=
{{1} \over {|S|}}\sum_{s \in S} x_1^{f_1(s)} \cdots x_r^{f_r(s)} \quad .
$$
The {\it expectations}, $\bar f_1, \dots, \bar f_r$ are simply
$$
\bar f_i =\left ( {{\partial} \over {\partial x_i}} P \right) (1, \dots, 1) 
\quad.
$$
The {\it centralized} probability generating function, which has mean $0$, is
$$
\tilde P(x_1, \dots, x_r)=\frac{P(x_1, \dots, x_r)} 
{x_1^{\bar f_1}  \cdots x_r^{\bar f_r}} \quad .
$$
The mixed moments (about the mean) of the statistics
$f_1(s), \dots, f_r(s)$ are defined by
$$
Mom[a_1, \dots, a_r] :=
{{1} \over {|S|}}
\sum_{s \in S} (f_1(s)-\bar f_1)^{a_1} \cdots (f_r(s)-\bar f_r)^{a_r} \quad .
$$
Often it is more convenient to consider the
mixed {\it factorial moments}, using the combinatorial ``powers''
$z^{(r)}:=z(z-1) \cdots (z-r+1)$, better known as the
{\it falling-factorials}. The mixed factorial moments
are defined analogously by
$$
FM[a_1, \dots, a_r]=
{{1} \over {|S|}}
\sum_{s \in S} (f_1(s)-\bar f_1)^{(a_1)} \cdots (f_r(s)-\bar f_r)^{(a_r)} 
\quad .
$$
Once you know the $FM$'s for all $a_1, \dots, a_r \leq M$, you can
easily figure out the $Mom$'s, using Stirling numbers of the
second kind (see [Z]). For our present purposes the sets $S$ lie in the 
family of {\it finite} sets $S_n$, and so the leading terms
of the $FM$'s and the $Mom$'s are the same.  Hence for the leading
asymptotics it suffices to consider the easier $FM$'s.

The best way to compute $FM[a_1, \dots, a_r]$
for all $0 \leq a_1, \dots, a_r \leq N$ for some fixed positive integer
$N$ is via the {\it Taylor expansion}  of $\tilde P(x_1, \dots, x_r)$
around $(x_1, \dots, x_r)=(1,\dots, 1)$, or equivalently,
the {\it Maclaurin expansion} of 
$\tilde P(1+x_1, \dots, 1+x_r)$
$$
\tilde P(1+x_1, \dots, 1+x_r)=\sum_{\alpha_1, \dots, \alpha_r \geq 0}
\frac{FM(\alpha_1, \dots, \alpha_r)}{\alpha_1 ! \cdots \alpha_r !}
 x_1^{\alpha_1} \cdots x_r^{\alpha_r} \quad .
$$

{\bf Back to inv-maj}

In the present approach, we need to tolerate the
more general discrete function $F(n,i)(p,q)$ even though
ultimately we are only interested in $H(n)(p,q)=F(n+1,n+1)(p,q)$.
We will prove the stronger statement that even if you restrict
attention to those $(n-1)!$ permutations that end with a specific
$i$, it is still asymptotically-joint-independently normal.

Since the averages of both $inv$ and $maj$ over the permutations that end in $i$ is $n-i+(n-1)(n-2)/4$,
the {\it centralized probability generating function} corresponding
to $F(n,i)(p,q)$ is:
$$
G(n,i)(p,q):=
\frac{F(n,i)(p,q)}{(n-1)!(pq)^{n-i+(n-1)(n-2)/4}} \quad .
$$
The recurrence $(RecF)$ becomes
$$
G(n,i)={{1} \over {q}}G(n,i+1)+
\frac{p^{n-i}(q^{n-1}-1)}{(pq)^{n/2} (n-1)}G(n-1,i) \quad ,
\eqno(RecG)
$$
and the {\it final condition} becomes
$$
G(n,n)= \frac{1}{n-1} \sum_{j=1}^{n-1} (pq)^{n/2-j} G(n-1,j)\quad .
\eqno(Gnn)
$$
We also need the obvious {\it initial condition}
$G(1,i)=\delta_{i,1}$.

{\bf Guessing Polynomial Expressions for the Factorial Moments}

Equipped with these efficient recurrences, our computer computes $G(n,i)(p,q)$ for many values of $n$ and $i$.
 Then for each of these it computes the
factorial moments $FM(r,s)(n,i)$, for many small numeric $(r,s)$
by computing the (intial terms of the) taylor series centered at $(p,q)=(1,1)$.
Then it fixes numeric values of $(r,s)$ and uses {\it polynomial interpolation} to {\it guess}
explicit polynomial expressions for $FM(r,s)(n,i)$ as polynomials
in $(n,i)$, for that pair $(r,s)$.  The process is repeated for many pairs $(r,s)$. We note
that it is obvious, both from the combinatorics and from the
recurrences, that $FM(r,s)(n,i)$ are always polynomials in $(n,i)$ for
any fixed numeric $r$ and $s$.

It would have been nice if we could guess {\it closed-form} expressions
for $FM(r,s)(n,i)$ for {\it symbolic} $(r,s)$, but no such closed-form
exists as far as we know, and besides, it is too much to ask for and more than we need.
{\it But} to prove {\it asymptotic normality} we only need the
{\it leading terms}. Viewing the leading terms, our beloved computer
easily conjectures the following expressions, for integers $r,s \geq 0$
$$
FM(2r,2s)(n,i)=\frac{(2r)!}{2^r r!}\frac{(2s)!}{2^s s!}
\left ( \frac{1}{36} \right )^{r+s} n^{3r+3s}
+(lower-degree-terms-in-(n,i)) \quad ,
$$
$$
FM(2r,2s-1)(n,i)=\frac{(2r)!}{2^r r!}\frac{(2s)!}{2^s s!}
\left ( \frac{1}{36} \right )^{r+s-1} n^{3r+3s-6}
[-(s-1)n^3-6rn^2i+18rni^2-12ri^3]
$$
$$	
+(lower-degree-terms-in-(n,i)) \quad ,
$$
$$
FM(2r-1,2s)(n,i)=-\frac{(2r)!}{2^r r!}\frac{(2s)!}{2^s s!}
\left ( \frac{1}{36} \right )^{r+s-1} (r-1)n^{3r+3s-3}
+(lower-degree-terms-in-(n,i)) \quad ,
$$
$$
FM(2r-1,2s-1)(n,i)=\frac{(2r)!}{2^r r!}\frac{(2s)!}{2^s s!}
\left ( \frac{1}{36} \right )^{r+s-1} \frac{9}{2} n^{3r+3s-6} (n-2i)^2
+(lower-degree-terms-in-(n,i)) \quad .
$$

{\bf Nice conjectures but what about proofs?}

While we prefer the empirical approach of guessing, an alternative
approach to finding many $FM(r,s)(n,i)$'s is to first use
$(RecG)$ and $(Gnn)$. Write $G(n,i)(1+p,1+q)$ as an {\it infinite}
{\bf generic} Taylor series around $(0,0)$, and write-down the
implied {\it infinite-order} recurrence expressing
$FM(r,s)(n,i)$ in terms of $FM(r',s')$ with $r'+s' <r+s$.
Note that in order to generate $FM(r,s)$ we only need finitely
many terms. Of course, as we have already commented, there is
no hope for finding a general expression for 
$FM(2r,2s)(n,i)$ , $FM(2r,2s-1)(n,i)$ , $FM(2r-1,2s)(n,i)$  and
$FM(2r-1,2s-1)(n,i)$, depending explicitly on $r$ and $s$ as well as on $n$ 
and $i$, but to prove by induction on $r$ and $s$,
that the above {\it leading terms}  are valid, all we need is
verify that the leading terms of the implied recurrences for the
$FM(r,s)$'s (easily derivable by hand, although we used the computer)
are consistent with the above explicit expressions.

The implication of $(RecG)$ is
$$
FM(r,s)(n,i)-FM(r,s)(n,i+1)=
sFM(r,s-1)(n,i+1)-sFM(r,s-1)(n-1,i)+ (lower-order-terms) \quad ,
$$
while the implication of $Gnn$ starts with
$$
FM(r,s)(n,n)=
\frac{1}{n-1} \sum_{i=1}^{n-1} FM(r,s)(n-1,i) -
\frac{s}{2(n-1)} \sum_{i=1}^{n-1} (2i-n)FM(r,s-1)(n-1,i) -
$$
$$
\frac{r}{2(n-1)} \sum_{i=1}^{n-1} (2i-n)FM(r-1,s)(n-1,i) 
+
\frac{rs}{4(n-1)} \sum_{i=1}^{n-1} (2i-n)^2 FM(r-1,s-1)(n-1,i) 
$$
$$
+ (lower-order-terms) \quad .
$$
The next step is to refine these recurrences into the four cases where $(r,s)$ are
$(even,even)$, $(even,odd)$, $(odd,even)$, and $(odd,odd)$.  
Once we have these four recurrences, we (or better, our computer) plugs in the above conjectured
expressions for the leading terms of $FM(r',s')(n,i)$ for $r'+s'<r+s$.  When the dust settles,
we are left with an expression with leading terms matching those in our conjecture for $FM(r,s)(n,i)$.
The proof follows from a routine yet intricate proof by induction on $r+s$.

{\bf The Maple package InvMaj}

All the nitty-gritty calculations described above, that
constitutes  {\bf a fully rigorous} proof, may be found
in the Maple package {\tt InvMaj} accompanying this article.
This package is available from the webpage of the present article:

{\tt http://www.math.rutgers.edu/\Tilde zeilberg/mamarim/mamarimhtml/invmaj.html} \quad ,

where the reader can also find some sample input and output.
The direct url of the package is: \hfill \break {\tt http://www.math.rutgers.edu/\Tilde zeilberg/tokhniot/InvMaj} \quad .

{\bf La Grande Finale}

The special case $r=1,s=0$ and $r=0,s=1$ give
$$
FM(2,0)(n,i)=\frac{1}{36}n^3+ O(n^2) \quad
$$
$$
FM(0,2)(n,i)=\frac{1}{36}n^3+ O(n^2) \quad .
$$
So 
$$
\frac{FM(2r,2s)}{FM(2,0)^rFM(0,2)^s}=
\frac{(2r)!}{2^r r!}\frac{(2s)!}{2^s s!}+O(1/n) \quad ,
$$
$$
\frac{FM(2r,2s-1)}{FM(2,0)^rFM(0,2)^{s-1/2}}=
o(1/n) \quad ,
$$
$$
\frac{FM(2r-1,2s)}{FM(2,0)^{r-1/2}FM(0,2)^s}=
o(1/n) \quad ,
$$
$$
\frac{FM(2r-1,2s-1)}{FM(2,0)^{r-1/2}FM(0,2)^{s-1/2}}=
O(1/n) \quad .
$$
And we see that as $n \rightarrow \infty$ these indeed converge
to the mixed moments of the famous mixed moments
of the bivariate independent normal distribution
$e^{-a^2/2-b^2/2}/(2 \pi)$ \halmos.

{\bf Encore: A more refined asymptotics for the (Normalized) Mixed Moments}

With more effort, we (or rather, our computer) can guess-and-prove
the following asymptotics for the case of interest
$(n+1,n+1)$, i.e. the asymptotic expressions for the
(genuine, not factorial) normalized mixed-moments
divided by the appropriate powers of the variances
(also known as mixed-alpha coefficients), let's call
them $\alpha(r,s)(n)$ . Indeed, according to S. B. Ekhad, we have:
$$
\alpha (2r,2s)(n)=\frac{(2r)!}{2^r r!}\frac{(2s)!}{2^s s!}
\left (  1- \frac{9(r^2+s^2-r-s)}{25} \cdot \frac{1}{n} +O(\frac{1}{n^2}) \right )
$$
$$
\alpha (2r-1,2s-1)(n)=\frac{(2r)!}{2^r r!}\frac{(2s)!}{2^s s!}
\left ( \frac{9}{2n}+
(\frac{-81}{50}(r^2+s^2)+\frac{243}{50}(r+s)-\frac{1773}{100})
\frac{1}{n^2}
+O(\frac{1}{n^3}) \right ) \quad .
$$
Of course, by {\it symmetry} $\alpha(2r,2s-1)$ and
$\alpha(2r-1,2s)$ are {\it identically} (not just asymptotically!)
zero. 

{\bf References}

[E] Shalosh B. Ekhad, {\it The joy of brute force: the covariance
of the major index and the number of inversions}, Personal Journal
of S.B.Ekhad and D. Zeilberger, 
{\tt http://www.math.rutgers.edu/~Tilde zeilberg/pj.html}, ca. 1995.

[Fe] William Feller, ``{\it An Introduction to Probability
Theory and Its Application}'', volume 1, three editions.
John Wiley and sons. First edition: 1950. Second edition:
1957. Third edition: 1968.

[Fo] Dominique Foata, {\it On the Netto inversion number of a sequence}, Proc. Amer. Math. Soc. {\bf 19} (1968), 236-240.

[GKP] Ronald Graham, Oren Patashnik, 
and Donald E. Knuth, {\it Concrete Mathematics: A Foundation for
Computer Science}, Addison Wesley, Reading, 1989.

[M] Percy A. MacMahon, {\it The indices of permutations}, Amer. J. Math. {\bf 35} (1913), 281-322.

[Z] Doron Zeilberger, {\it The automatic central limit theorems
generator (and much more!)}, in: ``Advances in Combinatorial
Mathematics'' (in honor of Georgy P. Egorychev), Ilias S. Kotsireas
and Eugene V. Zima, eds., Springer, 2009.

\end










