Quickselect and inversions. (II)

[Marko Riedel <mriedel@neuearbeit.de>]

Greetings.

I am sending a correction to the code contained in
  http://pauillac.inria.fr/algo/AofA/mailing_list/msg00432.html.

The quickselect implementation in the above message shifted the index
of the element being selected when there was no need to shift it. This 
resulted in some intermediate generating functions not being symmetric 
with respect to the location of the pivot.

Here are the correct generating functions.

0: 1
1 1
0/1
0: 2
1: 2
2 4
0/4
0: 5 + z
1: 6
2: 5 + z
3 16 + 2z
2/18
0: 15 + 6z + 2z^2 + z^3
1: 20 + 4z
2: 20 + 4z
3: 15 + 6z + 2z^2 + z^3
4 70 + 20z + 4z^2 + 2z^3
34/96
0: 52 + 33z + 15z^2 + 11z^3 + 5z^4 + 3z^5 + z^6
1: 75 + 30z + 10z^2 + 5z^3
2: 86 + 28z + 6z^2
3: 75 + 30z + 10z^2 + 5z^3
4: 52 + 33z + 15z^2 + 11z^3 + 5z^4 + 3z^5 + z^6
5 340 + 154z + 56z^2 + 32z^3 + 10z^4 + 6z^5 + 2z^6
444/600
0: 203 + 182z + 109z^2 + 81z^3 + 50z^4 + 40z^5 + 26z^6 + 15z^7 + 9z^8 + 4z^9 + z^10
1: 312 + 198z + 90z^2 + 66z^3 + 30z^4 + 18z^5 + 6z^6
2: 396 + 198z + 76z^2 + 40z^3 + 10z^4
3: 396 + 198z + 76z^2 + 40z^3 + 10z^4
4: 312 + 198z + 90z^2 + 66z^3 + 30z^4 + 18z^5 + 6z^6
5: 203 + 182z + 109z^2 + 81z^3 + 50z^4 + 40z^5 + 26z^6 + 15z^7 + 9z^8 + 4z^9 + z^10
6 1822 + 1156z + 550z^2 + 374z^3 + 180z^4 + 116z^5 + 64z^6 + 30z^7 + 18z^8 + 8z^9 + 2z^10
5508/4320
0: 877 + 1034z + 777z^2 + 631z^3 + 434z^4 + 351z^5 + 272z^6 + 206z^7 + 164z^8 + 118z^9 + 78z^10 + 49z^11 + 29z^12 + 14z^13 + 5z^14 + z^15
1: 1421 + 1274z + 763z^2 + 567z^3 + 350z^4 + 280z^5 + 182z^6 + 105z^7 + 63z^8 + 28z^9 + 7z^10
2: 1951 + 1402z + 712z^2 + 477z^3 + 255z^4 + 156z^5 + 72z^6 + 15z^7
3: 2162 + 1466z + 672z^2 + 430z^3 + 210z^4 + 80z^5 + 20z^6
4: 1951 + 1402z + 712z^2 + 477z^3 + 255z^4 + 156z^5 + 72z^6 + 15z^7
5: 1421 + 1274z + 763z^2 + 567z^3 + 350z^4 + 280z^5 + 182z^6 + 105z^7 + 63z^8 + 28z^9 + 7z^10
6: 877 + 1034z + 777z^2 + 631z^3 + 434z^4 + 351z^5 + 272z^6 + 206z^7 + 164z^8 + 118z^9 + 78z^10 + 49z^11 + 29z^12 + 14z^13 + 5z^14 + z^15
7 10660 + 8886z + 5176z^2 + 3780z^3 + 2288z^4 + 1654z^5 + 1072z^6 + 652z^7 + 454z^8 + 292z^9 + 170z^10 + 98z^11 + 58z^12 + 28z^13 + 10z^14 + 2z^15
69264/35280

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 inversions {
  my ($perm) = @_;
  my ($n, $count, $i, $j) = (scalar @$perm, 0);

  for($i=0; $i<$n; $i++){
    for($j=$i+1; $j<$n; $j++){
      $count++ if $perm->[$i]>$perm->[$j];
    }
  }

  return $count;
}

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 ? 'z' : "z^$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];
}



MAIN: {
  my ($max, $n) = @ARGV;

  die "$0 <n>\n" if not defined $max;

  for($n=1; $n<=$max; $n++){
    my ($gf, $tgf, $i, $j, $perm, $invs) =
      ([undef, undef, []]);
    for($j=0; $j<$n; $j++){
      $tgf=[undef, undef, []];
      for($i=0; $i<fact($n); $i++){
	$perm=nextperm($n, $i);
	qsel($perm, 0, $n-1, $j);
	$invs=inversions($perm);
	
	$tgf->[0]=$invs 
	  if not defined $tgf->[0] or $tgf->[0]>$invs;
	$tgf->[1]=$invs 
	  if not defined $tgf->[1] or $tgf->[1]<$invs;
	$tgf->[2][$invs]=0 
	  if not defined $tgf->[2][$invs];
	$tgf->[2][$invs]++;
      }

      print "$j: " . gf2string($tgf) . "\n";
      $gf=($j>0 ? gfadd($gf, $tgf) : $tgf);
    }
    
    print "$n " . gf2string($gf) . "\n";
    my ($at1) = gfat1($gf);
    print "$at1->[0]/$at1->[1]\n";
  }
}