[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Quickselect and the number of cycles. (VIII)
Greetings.
My previous post re. quickselect and the number of cycles contained
the right functions, but that was sheer luck. I set all indeterminates
in the solution of the DE to zero and let this be the result. Of
course the indeterminates should be computed from the initial
conditions. I am sending a Perl/MAXIMA script (the URL is
http://www.ma.utexas.edu/maxima.html) that does this.
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, $upper, $c, $len, $i, $j, @marks) = @_;
@marks = (0) x (scalar @$perm);
$c=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++;
}
$c++ if $len<=$upper;
}
}
return $c;
}
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 ? 'u' : "u^$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, $k, $n) = @ARGV;
die "$0 <n>\n" if not defined $max;
my ($initmax, @gfs) = ($max+1);
for($n=1; $n<=$initmax; $n++){
my ($f, $gf, $tgf, $i, $j, $perm, $c, $at1) =
(fact($n), [1, $initmax+1, []]);
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=cycles($perm, $max);
$tgf->[0]=$c
if not defined $tgf->[0] or $tgf->[0]>$c;
$tgf->[1]=$c
if not defined $tgf->[1] or $tgf->[1]<$c;
$tgf->[2][$c]=0
if not defined $tgf->[2][$c];
$tgf->[2][$c]++;
}
$at1 = gfat1($tgf);
$gf->[2][$j+1] = $at1->[0];
}
$gfs[$n]=$gf;
}
print "1/(1-z);\n";
print "1/(1-z)*(z";
for($k=2; $k<=$max; $k++){
print "+z^$k/$k";
}
print ");\n";
print "z*u/(1-z)/(1-z*u);\n";
print "u*SUBST(z*u, z, D1)*D3;\n";
print "u*SUBST(z*u, z, D2)*D3;\n";
print "u*SUBST(z*u, z, D1)*G(z);\n";
print "u*SUBST(z*u, z, D1)*D1;\n";
print "u*SUBST(z*u, z, D2)*D1;\n";
print "u*SUBST(z*u, z, D1)*D2;\n";
print "D1*D3;\n";
print "D2*D3;\n";
print "D1*G(z);\n";
print "FACTOR(RAT(D4+D5+D7+D8+D9+D10+D11));\n";
print "FACTOR(RAT(D6+D12));\n";
print "f*u/(1-z)/(1-z*u)*log(1/(1-z));\n";
print "f/(1-z)/(1-z*u)*log(1/(1-z*u));\n";
for($k=0; $k<=$max+1; $k++){
for($l=0; $l<=$max+1; $l++){
if($k==0 && $l==0){
print "(a0b0";
}
else{
print " + a$k" ."b$l*u^$k*z^$l";
}
}
}
print ");\n";
print "D15+D16+D17/(1-z)/(1-z*u);\n";
print "EXPAND(RAT((1-z)^2*(1-z*u)^2*DIFF(D15+D16, z)+" .
"DIFF(D17, z)*(1-z)*(1-z*u)-D17*DIFF((1-z)*(1-z*u), z)));\n";
print "EXPAND(RAT((1-z)^2*(1-z*u)^2*D13));\n";
print "RAT(SUBST((D15+D16)*(1-z)^2*(1-z*u)^2+D17*(1-z)*(1-z*u), " .
"G(z), D14));\n";
print "EXPAND(D19-D20-D21);\n";
for($k=0; $k<=$max+1; $k++){
for($l=0; $l<=$max+1; $l++){
print "COEFF(COEFF(D22, u, $k), z, $l)=0;\n";
}
}
$pos=23+($max+2)*($max+2);
my ($taylpos) = ($pos);
print "TAYTORAT(TAYLOR(D18, z, 0, " .
$initmax . "));\n";
$pos++;
for($n=1; $n<=$initmax; $n++){
print "(";
print gf2string($gfs[$n]);
print ")*z^$n/" . fact($n) . ";\n";
}
$pos+=$initmax;
print "EXPAND(D$taylpos";
for($n=1; $n<=$initmax; $n++){
print "-D" . ($taylpos+$n);
}
print ");\n";
$pos++;
my ($initial) = ($pos-1);
for($k=0; $k<=2*$initmax; $k++){
for($l=0; $l<=$initmax; $l++){
print "COEFF(COEFF(D$initial, u, $k), z, $l)=0;\n";
}
}
$pos+=(2*$initmax+1)*($initmax+1);
print "LINSOLVE([D23";
for($k=24; $k<23+($max+2)*($max+2); $k++){
print ", D$k";
}
for($k=$initial+1;
$k<$initial+1+(2*$initmax+1)*($initmax+1); $k++){
print ", D$k";
}
print "],\n[f";
for($k=0; $k<=$max+1; $k++){
for($l=0; $l<=$max+1; $l++){
print ", a$k" . "b$l";
}
}
print "]);\n";
print "SUBLIS(D$pos, D18);\n";
}