Quickselect and inversions. (IV)

[Marko Riedel <mriedel@neuearbeit.de>]

Greetings.

I have found the pattern that describes the "fudge factor" mentioned
in
  http://pauillac.inria.fr/algo/AofA/mailing_list/msg00434.html.
It's a binomial coefficient.  

I am sending the recurrence and some code to compute the generating
functions. I am very curious whether anyone is interested in
investigating this problem or could point me to a reference where it
is discussed.

Regards,

Marko Riedel

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

\documentclass{article}

\usepackage{latexsym,amsfonts,amssymb,amsmath,amstext,amsthm}

\begin{document}

\newcommand{\MWHERE}{\quad \text{where }  \quad}

\title{Quickselect and inversions}

\author{Marko Riedel}
\maketitle

\section{Recurrence relation}

Consider the following question. What amount of disorder is left in a
permutation after one of its elements has been selected with
quickselect?

This can be rephrased as follows. Compute the "grand average" for the
expected number of inversions left in a permutation after one of its
elements has been selected with quickselect.

Let $r_n(z)$ be the generating function of the inversions in a random
permutation of size $n$, where $0, 1,\ldots n-1$ are the elements being
permuted. We have

\[ r_n(z) = 
  \prod_{k=0}^{n-1} \sum_{l=0}^k z^l =
  \prod_{k=0}^{n-1} \frac{z^{k+1}-1}{z-1}.\]

Let $q_{n, l}(z)$ be the generating function of the number of
inversions after the element $l$, where $0 \le l < n$ has been
selected from a random permutation of size $n$. The base cases are:
$n=1$ and $q_{1, 0}(z) = 1$ and $n=2$ and $q_{2, 0}(z) = q_{2, 1}(z) =
2.$ We let $r_0(z) = q_{0, l}(z) = 1.$ Let $p$ denote the pivot. The
recurrence for $q_{n, l}(z)$ is as follows.
\begin{eqnarray*}
 q_{n, l}(z) & = &
  \sum_{p=0}^{l-1} \binom{n-1}{p} r_p(z) q_{n-1-p, l-p-1}(z) \\
  & + & \binom{n-1}{l} r_l(z) r_{n-1-l}(z) \\
  & + & \sum_{p=l+1}^{n-1} \binom{n-1}{p} q_{p, l}(z) r_{n-1-p}(z).
\end{eqnarray*}

Problem. Explain the presence of the factor $\binom{n-1}{p}.$

We are interested in the asymptotics of
\[ \frac{q_n'(1)}{q_n(1)} = \frac{q_n'(1)}{n \cdot n!}
   \MWHERE
   q_n(z) = \sum_{l=0}^{n-1} q_{n, l}(z).\]

References. Basic generating functions.

\end{document}

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


#include <stdio.h>
#include <ginac/ginac.h>

using namespace std;
using namespace GiNaC;

static symbol z("z");

ex nchoosek(ex n, ex k){
  return 
    factorial(n)/factorial(k)/factorial(n-k);
}

#define MAXN 256

static ex invgftable[MAXN];

void invgf(int k, ex *result){
  ex factor=0, term;
  int j;

  if(invgftable[k]==-1){
    for(j=0; j<k; j++){
      factor+=pow(z, j);
    }

    if(k>1){
      invgf(k-1, &term);
    }
    else{
      term=1;
    }
    term*=factor; 
    invgftable[k]=term.expand();
  }

  *result=invgftable[k];
  return;
}

static ex qselgftable[MAXN][MAXN];

void qselgfl(int n, int l, ex *result){
  ex term, prod, gf, f=factorial(n-1);
  int p;

  if(qselgftable[n][l]==-1){
    if(n==1 || n==2){
      qselgftable[n][l]=n;
    }
    else{
      term=0;
      for(p=0; p<n; p++){
	prod=1;
	if(p<l){
	  if(p>0){
	    invgf(p, &gf);
	    prod*=gf;
	  }
	  if(n-1-p>0){
	    qselgfl(n-1-p, l-p-1, &gf);
	    prod*=gf;
	  }
	}
	else if(p>l){
	  if(p>0){
	    qselgfl(p, l, &gf);
	    prod*=gf;
	  }
	  if(n-1-p>0){
	    invgf(n-1-p, &gf);
	    prod*=gf;
	  }
	}
	else{
	  if(p>0){
	    invgf(p, &gf);
	    prod*=gf;
	  }
	  if(n-1-p>0){
	    invgf(n-1-p, &gf);
	    prod*=gf;
	  }
	}
	term+=prod*nchoosek(n-1, p);
      }

      qselgftable[n][l]=term.expand();
    }
  }

  *result=qselgftable[n][l];
  return;
}

void qselgf(int n, ex *result){
  ex sum, term;
  int l;

  sum=0;
  for(l=0; l<n; l++){
    qselgfl(n, l, &term);
    sum+=term;
  }

  *result=sum;
  return;
}

void printgf(ex gf){
  int i, maxdeg=gf.degree(z);
  ex c;

  for(i=gf.ldegree(z); i<=maxdeg; i++){
    c=gf.coeff(z, i);
    if(c==0){
      continue;
    }

    if(i==0){
      cout << c;
    }
    else if(i==1){
      if(c!=1){
	cout << c << "z";
      }
      else{
	cout << "z";
      }
    }
    else{
      if(c!=1){
	cout << c << "z^" << i;
      }
      else{
	cout << "z^" << i;
      }
    }

    if(i<maxdeg){
      cout << " + ";
    }
  }
  cout << endl;

  return;
}

int main(int argc, char **argv){
  int i, j, l, n;
  ex gf, u, v, w;

  if(argc!=2){
    printf("usage: %s <n>\n", argv[0]);
    exit(1);
  }

  sscanf(argv[1], "%d", &n);
  if(n<1 || n>=MAXN){
    printf("range: 1 to %d\n", MAXN-1);
    exit(2);
  }

  for(i=0; i<MAXN; i++){
    invgftable[i]=-1;
  }

  cout << "inversion gfs" << endl;
  for(i=1; i<=n; i++){
    invgf(i, &gf);
    cout << i << " ";
    printgf(gf);
  }
  
  for(i=0; i<MAXN; i++){
    for(j=0; j<MAXN; j++){
      qselgftable[i][j]=-1;
    }
  }

  cout << "quickselect average gfs" << endl;
  for(i=1; i<=n; i++){
    for(l=0; l<i; l++){
      qselgfl(i, l, &gf);
      cout << l << ": ";
      printgf(gf);
    }
   
    qselgf(i, &gf);
    cout << i << " ";
    printgf(gf);

    u=gf.diff(z).subs(z==1);
    v=gf.subs(z==1);
    w=u/v;
    cout << u << "/" << v << " " << w << " " 
	 << evalf(u/v) << endl;
  }
}