ANALYSIS of ALGORITHMS, Bulletin Board
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Quickselect and the number of cycles. (VI)
Greetings.
I solved the DE for the number of fixed points after a single
quickselect with the following MAXIMA script:
(The URL is http://www.ma.utexas.edu/maxima.html)
----------------------------------------------------------------------
f0*u/(1-z)/(1-z*u)*log(1/(1-z));
f1/(1-z)/(1-z*u)*log(1/(1-z*u));
(a0 + a1*z + a2*u + a3*z*u + a4*z^2 + a5*z^2*u)/(1-z)/(1-z*u);
D1+D2+D3;
EXPAND(RAT((1-z)^2*(1-z*u)^2*DIFF(D4, z)));
u + u*z + u^2*z - 3*u^2*z^2;
RAT((1 + u - 2*z*u)*(1-z)*(1-z*u)*D4);
EXPAND(D5-D6-D7);
COEFF(COEFF(D8, z, 3), u, 2)=0;
COEFF(COEFF(D8, z, 3), u, 1)=0;
COEFF(COEFF(D8, z, 2), u, 2)=0;
COEFF(COEFF(D8, z, 2), u, 1)=0;
COEFF(COEFF(D8, z, 2), u, 0)=0;
COEFF(COEFF(D8, z, 1), u, 2)=0;
COEFF(COEFF(D8, z, 1), u, 1)=0;
COEFF(COEFF(D8, z, 1), u, 0)=0;
COEFF(COEFF(D8, z, 0), u, 1)=0;
COEFF(COEFF(D8, z, 0), u, 0)=0;
LINSOLVE([D9, D10, D11, D12, D13, D14, D15, D16, D17, D18],
[f0, f1, a0, a1, a2, a3, a4, a5]);
----------------------------------------------------------------------
We have q_{n, j}'(1) = 2 H_{n+1-j} + 2 H_j - 3.
I verified this result with a Perl program that computes the
expectation from the definition and outputs the value from the formula
for the purpose of comparing the two.
----------------------------------------------------------------------
0: v 1/1 1/1
1 v
1/1
0: 2v^2 4/2 4/2
1: 2v^2 4/2 4/2
2 4v^2
8/4
0: v + 5v^3 16/6 16/6
1: 6v^3 18/6 18/6
2: v + 5v^3 16/6 16/6
3 2v + 16v^3
50/18
0: 2v + 7v^2 + 15v^4 76/24 76/24
1: 4v^2 + 20v^4 88/24 88/24
2: 4v^2 + 20v^4 88/24 88/24
3: 2v + 7v^2 + 15v^4 76/24 76/24
4 4v + 22v^2 + 70v^4
328/96
0: 9v + 18v^2 + 41v^3 + 52v^5 428/120 428/120
1: 10v^2 + 35v^3 + 75v^5 500/120 500/120
2: 6v + 28v^3 + 86v^5 520/120 520/120
3: 10v^2 + 35v^3 + 75v^5 500/120 500/120
4: 9v + 18v^2 + 41v^3 + 52v^5 428/120 428/120
5 24v + 56v^2 + 180v^3 + 340v^5
2376/600
0: 44v + 109v^2 + 128v^3 + 236v^4 + 203v^6 2808/720 2808/720
1: 54v^2 + 108v^3 + 246v^4 + 312v^6 3288/720 3288/720
2: 20v + 46v^2 + 40v^3 + 218v^4 + 396v^6 3480/720 3480/720
3: 20v + 46v^2 + 40v^3 + 218v^4 + 396v^6 3480/720 3480/720
4: 54v^2 + 108v^3 + 246v^4 + 312v^6 3288/720 3288/720
5: 44v + 109v^2 + 128v^3 + 236v^4 + 203v^6 2808/720 2808/720
6 128v + 418v^2 + 552v^3 + 1400v^4 + 1822v^6
19152/4320
----------------------------------------------------------------------
I don't know what I'll try next, the variance or the DE for the
expected number of cycles. Probably the latter.
Best regards,
Marko Riedel
----------------------------------------------------------------------
#! /usr/bin/perl -w
#
sub fact {
my ($n) = @_;
return 1 if $n==0 or $n==1;
my ($p) = (1);
$p*=$n-- while $n>1;
return $p;
}
sub qsel {
my ($perm, $lower, $upper, $index) = @_;
my ($temp);
return if $lower==$upper;
if($lower+1==$upper){
$temp=$perm->[$lower];
if($temp>$perm->[$upper]){
$perm->[$lower]=$perm->[$upper];
$perm->[$upper]=$temp;
}
return;
}
my ($pivot, $i, $j) =
($perm->[$lower], $lower+1, $upper);
while($i<=$j){
while($i<=$j && $perm->[$i]<$pivot){
$perm->[$i-1]=$perm->[$i];
$i++;
}
while($j>=$i && $perm->[$j]>$pivot){
$j--;
}
if($i<$j){
$temp=$perm->[$i];
$perm->[$i]=$perm->[$j];
$perm->[$j]=$temp;
}
}
$perm->[$i-1]=$pivot;
qsel($perm, $lower, $i-2, $index)
if $i-2>$lower && $index<=$i-2;
qsel($perm, $i, $upper, $index)
if $i<$upper && $index>=$i;
}
sub nextperm {
my ($n, $i, $j) = @_;
my (@crep, $rem, @source, @perm, $item);
for($j=1; $j<=$n; $j++){
$rem=$i%$j;
unshift @crep, $rem;
$i=($i-$rem)/$j;
}
@source=(0..$n-1);
for($j=0; $j<$n; $j++){
$item=splice @source, $crep[$j], 1;
push @perm, $item;
}
return \@perm;
}
sub cycles {
my ($perm, $c, $fp, $len, $i, $j, @marks) = @_;
@marks = (0) x (scalar @$perm);
$c=0; $fp=0;
for($i=0; $i<scalar @$perm; $i++){
if(!$marks[$i]){
$j=$i; $len=0;
while(!$marks[$j]){
$marks[$j] = 1;
$j = $perm->[$j];
$len++;
}
$fp++ if $len==1;
$c++;
}
}
return ($c, $fp);
}
sub gfadd {
my ($g1, $g2) = @_;
my ($g, $sum, $i) =
([($g1->[0]<$g2->[0] ?
$g1->[0] : $g2->[0]),
($g1->[1]>$g2->[1] ?
$g1->[1] : $g2->[1]),
[]]);
for($i=$g->[0]; $i<=$g->[1]; $i++){
$sum=0;
$sum+=$g1->[2][$i] if defined $g1->[2][$i];
$sum+=$g2->[2][$i] if defined $g2->[2][$i];
$g->[2][$i]=$sum if $sum!=0;
}
return $g;
}
sub gf2string {
my ($poly, @str, $term) = @_;
foreach my $i ($poly->[0] .. $poly->[1]){
if(defined $poly->[2][$i]){
$term=($i == 0 ? '' : ($i == 1 ? 'v' : "v^$i"));
push @str,
($poly->[2][$i] == 1 ?
($i==0 ? '1' : '') : "$poly->[2][$i]") . $term;
}
}
return join ' + ', @str;
}
sub gfat1 {
my ($poly, $total, $total1) = (@_, 0, 0);
foreach my $i ($poly->[0] .. $poly->[1]){
if(defined $poly->[2][$i]){
$total += $poly->[2][$i];
$total1 += $i * $poly->[2][$i];
}
}
return [$total1, $total];
}
sub scaledharmonic {
my ($k, $s) = @_;
my ($sum) = (0);
for(my $i=1; $i<=$k; $i++){
$sum+=$s/$i;
}
return $sum;
}
MAIN: {
my ($max, $n) = @ARGV;
die "$0 <n>\n" if not defined $max;
for($n=1; $n<=$max; $n++){
my ($f, $gf, $tgf, $i, $j, $perm, $c, $fp, $at1) =
(fact($n), [undef, undef, []]);
for($j=0; $j<$n; $j++){
$tgf=[undef, undef, []];
for($i=0; $i<$f; $i++){
$perm=nextperm($n, $i);
qsel($perm, 0, $n-1, $j);
($c, $fp)=cycles($perm);
$tgf->[0]=$fp
if not defined $tgf->[0] or $tgf->[0]>$fp;
$tgf->[1]=$fp
if not defined $tgf->[1] or $tgf->[1]<$fp;
$tgf->[2][$fp]=0
if not defined $tgf->[2][$fp];
$tgf->[2][$fp]++;
}
print "$j: " . gf2string($tgf) . " ";
$at1 = gfat1($tgf);
print "$at1->[0]/$at1->[1] ";
my ($check) = [2*scaledharmonic($j+1, $f)
+ 2*scaledharmonic($n-$j, $f)
- 3*$f, $f];
print "$check->[0]/$check->[1]\n";
$gf=($j>0 ? gfadd($gf, $tgf) : $tgf);
}
print "$n " . gf2string($gf) . "\n";
$at1 = gfat1($gf);
print "$at1->[0]/$at1->[1]\n";
}
}
----------------------------------------------------------------------
\documentclass{article}
\usepackage{latexsym,amsfonts,amssymb,amsmath,amstext,amsthm}
\begin{document}
\newcommand{\half}{\frac1{2}}
\newcommand{\thrd}{\frac1{3}}
\newcommand{\qurt}{\frac1{4}}
\newcommand{\sxth}{\frac1{6}}
\newcommand{\octa}{\frac1{8}}
\newcommand{\dode}{\frac1{12}}
\newcommand{\EVAL}[2]{\left. #1 \right|_{#2}}
\newcommand{\PAR} [1]{\left(#1\right)}
\newcommand{\MAND} {\quad \text{and } \quad}
\newcommand{\MWHERE}{\quad \text{where } \quad}
\newcommand{\DZ}{\frac{\partial}{\partial z}}
\newcommand{\DV}{\frac{\partial}{\partial v}}
\title{Quickselect and the number of fixed points}
\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 fixed points of a permutation after one of its
elements has been selected with quickselect. The more fixed points,
the more ordered the permutation.
Let $r_n(v)$ be the exponential generating function of the fixed
points of a random permutation of size $n$, where $0, 1,\ldots n-1$
are the elements being permuted. We have
\[ r_n(v) = [z^n] \frac{e^{vz - z}}{1 - z}.\]
Let $q_{n, j}(v)$ be the exponential generating function of the number
of fixed points after the element $j$, where $1 \le j \le n$ has been
selected from a random permutation of size $n$. The base cases are:
$n=1$ and $q_{1, 1}(v) = v/1 = v$ and $n=2$ and $q_{2, 1}(v) = q_{2,
2}(v) = 2v^2/2 = v^2.$ We let $r_0(v) = q_{0, j}(v) = 1.$ Let $p$
denote the pivot, where $1 \le p \le n$. The recurrence for $q_{n,
j}(v)$ is as follows.
\begin{eqnarray*}
q_{n, j}(v) & = &
\frac{v}{n}
\sum_{p=1}^{j-1} r_{p-1}(v) q_{n-p, j-p}(v) \\
& + & \frac{v}{n} r_{j-1}(v) r_{n-j}(v) \\
& + & \frac{v}{n} \sum_{p=j+1}^n q_{p-1, j}(v) r_{n-p}(v).
\end{eqnarray*}
We are interested in the asymptotics of
\[ \frac{q_n'(1)}{n}
\MWHERE
q_n(v) = \sum_{j=1}^n q_{n, j}(v).\]
\newpage
\section{Trivariate generating function}
Let
\[ R(z, v) = \frac{e^{vz - z}}{1 - z} \]
and
\[ F(z, u, v) = \sum q_{n,j}(v) u^j z^n.\]
Note that
\[ R(z, 1) = \frac1{1 - z} \]
and
\[
\PAR{ \DV R(z, v) }(z, 1) = \frac{z}{1 - z}.\]
Recall that $q_{n, l}(1) = 1$ and hence
\begin{eqnarray*}
F(z, u, 1) & = & \sum_{n \ge 1} \sum_{j=1}^n u^j z^n =
\sum_{n \ge 1} z^n
\PAR{ -1 + \frac{1 - u^{n+1}}{1 - u}} \\
& = & - \frac{z}{1 - z}
+ \frac{1}{1 - u} \frac{z}{1 - z}
- \frac{u}{1 - u} \frac{z u}{1 - z u} \\
& = & \frac{z u}{(1 - z)(1 - z u)}.
\end{eqnarray*}
We have
\begin{eqnarray*}
n \: q_{n, j}(v) u^j z^{n-1}
& = & u \: v \:
\sum_{p=1}^{j-1} r_{p-1}(v) u^{p-1} z^{p-1}
q_{n-p, j-p}(v) u^{j-p} z^{n-p} \\
& + & u \: v \:
r_{j-1}(v) u^{j-1} z^{j-1} r_{n-j}(v) z^{n-j} \\
& + & v \:
\sum_{p=j+1}^n q_{p-1, j}(v) u^j z^{p-1} r_{n-p}(v) z^{n-p}
\end{eqnarray*}
and hence
\[\DZ F(z,u,v) =
uv \: R(zu,v) F(z,u,v) + uv \: R(zu,v) R(z,v)
+ v \: R(z,v) F(z,u,v).\]
Differentiate with respect to $v$ to obtain
\begin{eqnarray*}
\DV \DZ F(z,u,v) & = &
u R(zu,v) F(z,u,v) \\ & + &
uv \DV R(zu,v) F(z,u,v) + uv R(zu,v) \DV F(z,u,v) \\ & + &
u R(zu,v) R(z,v) \\ & + &
uv \DV R(zu,v) R(z,v) + uv R(zu,v) \DV R(z,v) \\ & + &
R(z,v) F(z,u,v) \\ & + &
v \DV R(z,v) F(z,u,v) + v R(z,v) \DV F(z,u,v).
\end{eqnarray*}
Let
\[ G(z, u) = \EVAL{ \PAR{ \DV F(z,u,v) }}{v=1} \]
and set $v=1$ so that
\begin{eqnarray*}
\DZ G(z, u) & = &
u \frac1{1 - zu} \frac{z u}{(1 - z)(1 - z u)} \\
& + &
u \frac{zu}{1 - zu} \frac{z u}{(1 - z)(1 - z u)}
+ u \frac1{1 - zu} G(z, u) \\ & + &
u \frac1{1 - zu} \frac1{1 - z} \\
& + &
u \frac{zu}{1 - zu} \frac1{1 - z}
+ u \frac1{1 - zu} \frac{z}{1 - z} \\ & + &
\frac1{1 - z} \frac{z u}{(1 - z)(1 - z u)} \\
& + &
\frac{z}{1 - z} \frac{z u}{(1 - z)(1 - z u)}
+ \frac1{1 - z} G(z, u)
\end{eqnarray*}
or
\begin{eqnarray*}
\DZ G(z, u) & = &
\frac{u + u z + u^2 z - 3 u^2 z^2}
{(1 - z)^2 (1 - zu)^2}
+ \frac{1 + u - 2 z u}{(1 - z) (1 - z u)} G(z, u).
\end{eqnarray*}
Solve this to obtain
\begin{eqnarray*}
G(z,u) & = &
2 \: \frac{u}{(1 - z)(1 - z u)} \log \frac1{1 - z}
+ 2 \: \frac1{(1 - z)(1 - z u)} \log \frac1{ 1 - z u} \\
& - & 3 \: \frac{zu}{(1 - z) (1 - z u)}.
\end{eqnarray*}
This proves that
\[ q_{n, j}'(1) = 2 H_{n+1-j} + 2 H_j - 3.\]
References. Basic generating functions.
\end{document}
Date Prev |
Date Next |
Date Index |
Thread Index