#############################################################################
##
#M  DeterminantMatDivFree( <M> )                            Timothy Debaillie
##                                                               Robert Morse
##                                                             Marcus Wassmer
##
##  Division free method. This is an alternative to the fraction free
##  version found in the library when division of matrix entries is 
##  expensive or not possible.
## 
##  This algorithm implements a division free method found in 
## 
##  "Determinant: Combinatorics, Algorithms, and Complexity", M. Mahajan, 
##  V. Vinay, Chicago Journal of Theoretical Computer Science, 1997 
##
##  The run time is $O(n^4)$
##  Auxillary storage size $n^2+n + C$
##
##  Our implementation has two runtime optimizations (both noted 
##  by Mahajan and Vinay)
##    1. Partial monomial sums, subtractions, and products are done at
##       each level. 
##    2. Prefix property is maintained allowing for a pruning of many
##       vertices at each level
##
##  and two auxillary storage size optimizations
##    1. only the upper triangular and diagonal portion of the 
##       auxillary storage is used.
##    2. Level information storage is reused (2 levels).
##
##
## $Id: divfreedet.g,v 1.9 2001/11/29 13:26:39 unialg Exp $
##
##
DeclareOperation("DeterminantMatDivFree", [IsMatrix]);

InstallMethod( DeterminantMatDivFree,
    "Division-free method",
    true,
    [ IsMatrix and IsRectangularTable], 0,
    function ( M )

        local u,v,w,i,   ## indices
              a,b,c,x,y, ## temp indices
              temp,      ## temp variable
              nlevel,    ## next level 
              clevel,    ## current level
              pmone,     ## plus or minus one
              zero,      ## zero of the ring
              n,         ## size of the matrix
              Vs,        ## final sum
              V;         ## graph

        # check that the argument is a square matrix and set the size
        n := Length(M);
        if n <> Length(M[1])  then
            Error("DeterminantMatDivFree: <mat> must be a square matrix");
        fi;

        ## initialze the final sum, the vertex set, initial parity
        ## and level indexes
        ##
        zero := Zero(M[1][1]);
        Vs := zero;
        V := [];
        pmone := (-One(M[1][1]))^((n mod 2)+1);
        clevel := 1; nlevel := 2;
     
        ##  Vertices are indexed [u,v,i] holding the (partial) monomials 
        ##  whose sums will form the determinant
        ##    where i = depth in the tree (current and next reusing
        ##              level storage)
        ##          u,v indices in the matrix
        ##
        ##  Only the upper triangular portion of the storage space is
        ##  needed. It is easier to create lower triangular data type
        ##  which we do here and index via index arithmetic.
        ##
        for u in [1..n] do
            Add(V,[]);
            for v in [1..u] do
                Add(V[u],[zero,zero]);
            od;
            ## Initialize the level 0 nodes with +/- one, depending on 
            ## the initial parity determined by the size of the matrix
            ##
            V[u][u][clevel] := pmone;
        od;

        ##  Here are the $O(n^4)$ edges labeled by the elements of
        ##  the matrix $M$. We build up products of the labels which form
        ##  the monomials which make up the determinant.
        ##    
        ##  1. Parity of monomials are maintained implicitly.
        ##  2. Partial sums for some vertices are not part of the final 
        ##     answer and can be pruned.
        ##
        for i in [0..n-2] do
            for u in [1..i+2] do  ## <---- pruning of vertices
                for v in [u..n] do         ## (maintains the prefix property)
                    for w in [u+1..n] do

                        ## translate indices to lower triangluar coordinates
                        ## 
                        a := n-u+1; b := n-w+1; c := n-v+1;
                        V[a][b][nlevel]:= V[a][b][nlevel]+
                            V[a][c][clevel]*M[v][w];
                        V[b][b][nlevel]:= V[b][b][nlevel]-
                            V[a][c][clevel]*M[v][u];
                    od;
                od;
            od;
 
            ## set the new current and next level. The new next level 
            ## is intialized to zero 
            ##
            temp   := nlevel; nlevel := clevel; clevel := temp;
            for x in [1..n] do
                for y in [1..x] do
                    V[x][y][nlevel] := zero;
                od;
            od;
        od; 

        ##  with the final level, we form the last monomial product and then
        ##  sum these monomials (parity has been accounted for)
        ##  to find the determinant.
        ##
        for u in [1..n] do
            for v in [u..n] do
                Vs := Vs + V[n-u+1][n-v+1][clevel]*M[v][u]; 
            od;
        od;
     
        ##  Return the final sum
        ##
        return Vs;

    end);