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;
}
?>