BCMATH version of the Fox-Landi algorithm

<?php
/* In the pseudo-code below,  all occurrences of "\n" should be replaced by "<br>\n" */

/* Last modified 10th March 2006, 11.30am Brisbane time */

function test_absorb($j,$matrix,$n){
/* This returns 1 if state j is absorbent, 0 otherwise */
global $delete_flag;

        $flag="1";
	if($matrix[$j][$j]>="1"){
	   for($d="1";$d<=$n;$d=bcadd($d,"1")){
	       if($d!=$j && $delete_flag[$d]=="0" && $matrix[$j][$d]!="0"){
                  $flag="0";
                  break;
               }
           }
	}else{
            $flag="0";
        }
        return($flag);
}

/* Now to initially label the states:  */

function step1($matrix,$n){
global $label_count;
global $state;
global $delete_flag;
global $erg;
global $erg_count;
global $set;
global $total;

for($i="1";$i<=$n;$i=bcadd($i,"1")){
    $set[$i]["1"]=$i;
    $total[$i]="1";
}
$erg_count="0";
/* Label the absorbing states: 
   If state j is absorbing, set state[j]=1.
   Label every state i!=j with b[i][j]>=1 as transient ie. state[i]=2
*/

        $label_count="0";
  	for($j="1";$j<=$n;$j=bcadd($j,"1")){
            $state[$j]="0";
	}
        for($j="1";$j<=$n;$j=bcadd($j,"1")){
           $delete_flag[$j]="0";
        }
  	for($j="1";$j<=$n;$j=bcadd($j,"1")){
	    $result=test_absorb($j,$matrix,$n);
	    if($result=="1"){
                 $state[$j]="1";
                 print "state[$j] is absorbent
\n"; $erg_count=bcadd($erg_count,"1"); $erg[$erg_count]=$j; $label_count=bcadd($label_count,"1"); for($k="1";$k<=$n;$k=bcadd($k,"1")){ if($k!=$j && $matrix[$k][$j]>="1"){ $state[$k]="2"; /* print "state[$k] is transient
\n"; */ $label_count=bcadd($label_count,"1"); } } } } if($label_count==$n){ return("1"); }else{ return("0"); } } function step2($matrix,$n,$verbose){ global $new_matrix; global $state; global $ichain; global $r; global $label_count; global $delete_flag; global $collection; global $total; global $ergodic_state_count; global $reordered; global $erg; global $set; global $erg_count; for($d="1";$d<=$n;$d=bcadd($d,"1")){ for($e="1";$e<=$n;$e=bcadd($e,"1")){ $new_matrix[$d][$e]=$matrix[$d][$e]; } } $ergodic_state_count="0"; if($verbose=="1"){ print "about to do step1
\n"; } $t=step1($matrix,$n); if($t=="1"){ perm($matrix,$n); return; }else{ $r="0"; } while($label_count<$n){ while("1"){ if($r=="0"){/* begin a chain of states by searching the row corresponding to the first unlabelled state for a positive element. Call the row number ichain[1]. */ if($verbose){ print "starting a new chain
\n"; } for($j="1";$j<=$n;$j=bcadd($j,"1")){ if($state[$j]=="0"){ if($verbose){ print "found 1st unlabelled state[$j]
\n"; } break; } } $ichain["1"]=$j; if($verbose){ print "ichain[1]=$j
\n"; } $r="1"; } /* now extend the chain ichain[1],...,ichain[r] by searching row ichain[r] for a positive element matrix[ichain[r]][ichain[r+1]] where ichain[r+1]!=ichain[r]. */ if($verbose){ print "extending the chain "; for($j="1";$j<=$r;$j=bcadd($j,"1")){ print "ichain[$j]=$ichain[$j],"; } print "
\n"; /* It seems from the algorithm that these states are unlabelled */ } for($j="1";$j<=$n;$j=bcadd($j,"1")){ $y=$ichain[$r]; $x=$new_matrix[$y][$j]; $z=$delete_flag[$j]; if($j!=$ichain[$r] && $delete_flag[$j]=="0" && $new_matrix[$ichain[$r]][$j]>="1"){ break; } } $s=bcadd($r, "1"); $ichain[$s]=$j; /* state_ichain[r+1]=$j */ if($verbose){ print "ichain[$s]=$j
\n"; print "found positive element matrix[$ichain[$r]][$j]
\n"; } if($state[$j]=="2"){ for($t="1";$t<=$r;$t=bcadd($t,"1")){ /* the algorithm seems to indicate that these states are unlabelled. To be on the safe side, I should only increase label_count if the state was previously unlabelled */ $state[$ichain[$t]]="2"; $x=$ichain[$t]; if($verbose){ print "state[$x] is transient
\n"; } $label_count=bcadd($label_count,"1"); } if($verbose){ print "finished labelling chain transient, go to start of step 2 again
\n"; } if($label_count==$n){ if($verbose){ print "all states are now labelled
\n"; } perm($matrix,$n); return; }else{ $r="0"; break; } }else{ $flag="0"; $continue_flag="0"; for($k="1";$k<=$r;$k=bcadd($k,"1")){ if($j==$ichain[$k]){ $flag="1"; if($verbose){ print "a circuit: ichain[$s]=$ichain[$k], now do step4
\n"; } $output=step4($n,$k,$verbose); if($output=="1"){ $continue_flag="1"; } break; } } if($continue_flag=="1"){ if($verbose){ print "finished labelling states ergodic, go start of step 2
\n"; } if($label_count==$n){ if($verbose){ print "all states are now labelled
\n"; } perm($matrix,$n); return; }else{ break; } } if($flag=="0"){ if($verbose){ print "increasing r=$r by 1 and extending the chain again
\n"; } $r=bcadd($r,"1"); } } } } } function boolean_sum($u,$v){ /* Here $u and $v are >= 0. */ if($u>="1" || $v>="1"){ return("1"); }else{ return("0"); } } function boolean_product($u,$v){ /* Here $u and $v are >= 0. */ if($u>="1" && $v>="1"){ return("1"); }else{ return("0"); } } function step4($n,$k,$verbose){ /* Replace row ichain[k] by the union of rows ichain[k],...,ichain[r] Replace column ichain[k] by the union of columns ichain[k],...,ichain[r] Delete rows and columns ichain[k+1],...,ichain[r]. If state[ichain[k]] is absorbing, label states ichain[k],...,ichain[r] ergodic. Label every state h with matrix[h][ichain[k]]=1 transient. Set r=0 and return 1. if state[ichain[k]] is not absorbing, set r=k and return 0. */ global $new_matrix; global $state; global $ichain; global $r; global $label_count; global $delete_flag; global $collection; global $total; global $set; global $ergodic_state_count; global $erg; global $erg_count; global $set_union; global $set_union_card; /* replace row ichain[k] by the union of rows ichain[k],...,ichain[r] */ if($k!=$r){ if($verbose){ print "replacing row ichain[$k]=$ichain[$k] by the union of rows "; for($e=$k;$e<=$r;$e=bcadd($e,"1")){ print "ichain[$e]=$ichain[$e], "; } print"
\n"; } for($d="1";$d<=$n;$d=bcadd($d,"1")){ if($delete_flag[$d]=="0"){ $sum="0"; for($e=$k;$e<=$r;$e=bcadd($e,"1")){ $sum=boolean_sum($sum, $new_matrix[$ichain[$e]][$d]); } $new_matrix[$ichain[$k]][$d]=$sum; } } /* replace column ichain[k] by the union of columns ichain[k],...,ichain[r] */ if($verbose){ print "replacing column ichain[$k]=$ichain[$k] by the union of columns "; for($e=$k;$e<=$r;$e=bcadd($e,"1")){ print "ichain[$e]=$ichain[$e], "; } print"
\n"; } for($d="1";$d<=$n;$d=bcadd($d,"1")){ if($delete_flag[$d]=="0"){ $sum="0"; for($e=$k;$e<=$r;$e=bcadd($e,"1")){ $sum=boolean_sum($sum, $new_matrix[$d][$ichain[$e]]); } $new_matrix[$d][$ichain[$k]]=$sum; } } /*delete rows and columns $ichain[k+1],...,$ichain[r]*/ if($verbose){ print "deleting rows and columns "; } $n1=bcadd($k,"1"); if($verbose){ for($d=$n1;$d<=$r;$d=bcadd($d,"1")){ print "ichain[$d]=$ichain[$d] "; } print "
\n"; print "r=$r, k=$k
\n"; } for($d=$n1;$d<=$r;$d=bcadd($d,"1")){ $delete_flag[$ichain[$d]]="1"; } } if($verbose){ echo "new matrix=
\n"; for($d="1";$d<=$n;$d=bcadd($d,"1")){ if($delete_flag[$d]=="0"){ for($e="1";$e<=$n;$e=bcadd($e,"1")){ if($delete_flag[$e]=="0"){ print $new_matrix[$d][$e]; print " "; } } print "
\n"; } } } $temp3=bcadd($k,"1"); $temp1=$ichain[$k]; for($h=$temp3;$h<=$r;$h=bcadd($h,"1")){ $temp=$ichain[$h]; $nn=$total[$temp]; union($set[$temp1],$total[$temp1],$set[$temp],$nn); $set[$temp1]=$set_union; $total[$temp1]=$set_union_card; } if($verbose){ $t=$total[$ichain[$k]]; for($s="1";$s<=$t;$s=bcadd($s,"1")){ $temp=$set[$ichain[$k]][$s]; echo "set[$ichain[$k]][$s]=$temp, "; } echo "
\n"; } $output=test_absorb($ichain[$k],$new_matrix,$n); if($output=="1"){ $total_k=$total[$ichain[$k]]; for($s="1";$s<=$total_k;$s=bcadd($s,"1")){ $state[$set[$ichain[$k]][$s]]="1"; $erg_count=bcadd($erg_count,"1"); $erg[$erg_count]=$set[$ichain[$k]][$s]; $label_count=bcadd($label_count,"1"); $x=$set[$ichain[$k]][$s]; if($s<$total_k){ print "$x, "; }else{ print "$x "; } } $ergodic_state_count=bcadd($ergodic_state_count,"1"); if($total_k>"1"){ print "constitute ergodic chain[$ergodic_state_count]
\n"; }else{/* this can't happen, as all absorbent states were found at start */ print "constitutes ergodic chain[$ergodic_state_count]
\n"; } for($h="1";$h<=$n;$h=bcadd($h,"1")){ if($h!=$ichain[$k] && $delete_flag[$h]=="0" && $new_matrix[$h][$ichain[$k]]>="1"){ if($state[$h]=="0"){ $state[$h]="2"; if($verbose){ print "state[$h] is transient
\n"; } $label_count=bcadd($label_count,"1"); } } } $r="0"; return($output); }else{ $r=$k; $x=$ichain[$k]; if($verbose){ print "state[$x] is not absorbent
\n"; } return($output); } } function perm($matrix,$n){ global $state; global $reordered; global $erg; global $erg_count; $k="0"; for($h="1";$h<=$n;$h=bcadd($h,"1")){ if($state[$h]=="2"){ $k=bcadd($k,"1"); $trans[$k]=$h; } } for($h="1";$h<=$erg_count;$h=bcadd($h,"1")){ $reordered[$h]=$erg[$h]; } for($h="1";$h<=$k;$h=bcadd($h,"1")){ $temp=bcadd($h,$erg_count); $reordered[$temp]=$trans[$h]; if($h<$k){ echo "$trans[$h], "; }else{ echo "$trans[$h] "; } } if($k>"1"){ echo "are the transient states
\n"; } if($k=="1"){ echo "is the only transient state
\n"; } echo "
\n"; echo "row/column permutation:
\n"; echo "("; for($d="1";$d<$n;$d=bcadd($d,"1")){ echo "$reordered[$d], "; } echo "$reordered[$n])
\n "; for($d="1";$d<=$n;$d=bcadd($d,"1")){ for($e="1";$e<=$n;$e=bcadd($e,"1")){ $permuted_matrix[$d][$e]=$matrix[$reordered[$d]][$reordered[$e]]; } } return; } function boolean_matrix_product($matrix1,$matrix2,$n){ /* Calculates the boolean product of two non-negative $nx$n matrices */ for($i="1";$i<=$n;$i=bcadd($i,"1")){ for($j="1";$j<=$n;$j=bcadd($j,"1")){ $sum="0"; for($k="1";$k<=$n;$k=bcadd($k,"1")){ $temp=boolean_product($matrix1[$i][$k],$matrix2[$k][$j]); $sum=boolean_sum($sum,$temp); } $matrix3[$i][$j]=$sum; } } return($matrix3); } function boolean_matrix_power($matrix,$n,$t){ if($t=="1"){ return($matrix); } $temp=$matrix; for($k="1";$k<=$t;$k=bcadd($k,"1")){ $matrix=boolean_matrix_product($matrix,$temp,$n); } return($matrix); } function union($set1,$m,$set2,$n){ /* $set1 is composed of $m integers $set2 is composed of $n integers Their union $set_union is composed of $set_union_card integers */ global $set_union; global $set_union_card; for($k="1";$k<=$m;$k=bcadd($k,"1")){ $set_union[$k]=$set1[$k]; } $set_union_card=$m; $k=$m; for($j="1";$j<=$n;$j=bcadd($j,"1")){ $flag="0"; for($i="1";$i<=$m;$i=bcadd($i,"1")){ if($set2[$j]==$set1[$i]){ $flag="1"; break; } } if($flag=="0"){ $k=bcadd($k,"1"); $set_union[$k]=$set2[$j]; $set_union_card=$k; } } return; } ?>