Quelle  cmat.gi   Sprache: unbekannt

#############################################################################
##
#W  cmat.gi               GAP 4 package `cvec'                
##                                                            Max Neunhoeffer
##
##  Copyright (C) 2007  Max Neunhoeffer, Lehrstuhl D f. Math., RWTH Aachen
##  This file is free software, see license information at the end.
##
##  This file contains the higher levels for compact matrices over finite 
##  fields. 
##

#############################################################################
# Creation:
#############################################################################

InstallGlobalFunction( CVEC_CMatMaker_GAP, function(l,cl)
    # Makes a new CMat, given a list l with a 0 in the first place
    local greasehint,m,q,qp,ty,filts;
    if Length(l) > 0 then
        m := rec(rows := l, len := Length(l)-1, vecclass := cl,
                 scaclass := cl![CVEC_IDX_GF]);
    else
        m := rec(rows := l, len := 0, vecclass := cl,
                 scaclass := cl![CVEC_IDX_GF]);
    fi;
    m.greasehint := cl![CVEC_IDX_fieldinfo]![CVEC_IDX_bestgrease];   
         # this is the current bestgrease
    q := cl![CVEC_IDX_fieldinfo]![CVEC_IDX_q];
    qp := q^m.greasehint;
    while m.greasehint > 0 and Length(l) < qp do
        m.greasehint := m.greasehint-1;
        qp := qp/q;
    od;
    return Objectify(cl![CVEC_IDX_typecmat],m);
end );
if IsBound(CVEC_CMatMaker_C) then
    CVEC_CMatMaker := CVEC_CMatMaker_C;
else
    CVEC_CMatMaker := CVEC_CMatMaker_GAP;
fi;

InstallMethod( ConstructingFilter, "for a cmat",
  [ IsCMatRep ],
  function( m ) return IsCMatRep; end );

# `NewMatrix` was a constructor
# in GAP <= 4.12 but is a tag based operation in GAP 4.13.
Perform( [ function( filt, f, rl, l )
    local p,d,c,li,i,v;
    p := Characteristic(f);
    d := DegreeOverPrimeField(f);
    c := CVEC_NewCVecClass(p,d,rl);
    li := 0*[1..Length(l)+1];
    for i in [1..Length(l)] do
        if IsCVecRep(l[i]) and IsIdenticalObj(c,DataObj(l[i])) then
            li[i+1] := ShallowCopy(l[i]);
        elif IsVectorObj(l[i]) then
            li[i+1] := CVec(Unpack(l[i]),c);
        elif IsPlistRep(l[i]) then
            li[i+1] := CVec(l[i],c);
        else
            Error("I do not know how to handle initialization data");
            return fail;
        fi;
    od;
    return CVEC_CMatMaker(li,c);
  end ],
  function( method )
    if IS_CONSTRUCTOR( NewMatrix ) then
      # This holds in GAP <= 4.12.
      InstallMethod( NewMatrix,
        [ "IsCMatRep", "IsField and IsFinite", "IsInt", "IsList" ], method );
    elif IsBoundGlobal( "InstallTagBasedMethod" ) then
      # This should hold in GAP 4.13.
      ValueGlobal("InstallTagBasedMethod")( NewMatrix, IsCMatRep, method );
    else
      # This should not happen.
      Error( "NewMatrix is neither a constructor not a tag based operation" );
    fi;
  end );

# `NewZeroMatrix` was a constructor
# in GAP <= 4.12 but is a tag based operation in GAP 4.13.
Perform( [ function( filt, f, rows, cols )
    local p, d, c, li, i;
    p := Characteristic(f);
    d := DegreeOverPrimeField(f);
    c := CVEC_NewCVecClass(p,d,cols);
    li := 0*[1..rows+1];
    for i in [2..rows+1] do
        li[i] := CVEC_NEW(c,c![CVEC_IDX_type]);
    od;
    return CVEC_CMatMaker(li,c);
  end ],
  function( method )
    if IS_CONSTRUCTOR( NewZeroMatrix ) then
      # This holds in GAP <= 4.12.
      InstallMethod( NewZeroMatrix,
        [ "IsCMatRep", "IsField and IsFinite", "IsInt", "IsInt" ], method );
    elif IsBoundGlobal( "InstallTagBasedMethod" ) then
      # This should hold in GAP 4.13.
      ValueGlobal("InstallTagBasedMethod")( NewZeroMatrix, IsCMatRep, method );
    else
      # This should not happen.
      Error( "NewZeroMatrix is neither a constructor not a tag based operation" );
    fi;
  end );

# `NewIdentityMatrix` was a constructor
# in GAP <= 4.12 but is a tag based operation in GAP 4.13.
Perform( [ function( filt, f, rows )
    local p, d, c, li, o, i;
    p := Characteristic(f);
    d := DegreeOverPrimeField(f);
    c := CVEC_NewCVecClass(p,d,rows);
    li := 0*[1..rows+1];
    o := One(f);
    for i in [1..rows] do
        li[i+1] := CVEC_NEW(c,c![CVEC_IDX_type]);
        li[i+1,i] := o;
    od;
    return CVEC_CMatMaker(li,c);
  end ],
  function( method )
    if IS_CONSTRUCTOR( NewIdentityMatrix ) then
      # This holds in GAP <= 4.12.
      InstallMethod( NewIdentityMatrix,
        [ "IsCMatRep", "IsField and IsFinite", "IsInt" ], method );
    elif IsBoundGlobal( "InstallTagBasedMethod" ) then
      # This should hold in GAP 4.13.
      ValueGlobal("InstallTagBasedMethod")( NewIdentityMatrix, IsCMatRep, method );
    else
      # This should not happen.
      Error( "NewIdentityMatrix is neither a constructor not a tag based operation" );
    fi;
  end );

InstallMethod( CMat, "for a list of cvecs and a cvec", [IsList, IsCVecRep],
  function(l,v)
    return CMat(l,DataObj(v),true);
  end);

InstallMethod( CMat, "for a list of cvecs, a cvec, and a boolean value",
  [IsList, IsCVecRep, IsBool],
  function(l,v,checks)
    return CMat(l,DataObj(v),checks);
  end);

InstallMethod( CMat, "for a list of cvecs", [IsList],
  function(l)
    local c;
    if Length(l) = 0 or not(IsBound(l[1])) then
        Error("CMat: Cannot use one-argument version with empty list");
        return fail;
    fi;
    c := DataObj(l[1]);
    return CMat(l,c,true);
  end);

InstallMethod( CMat, "for a list of cvecs, and a boolean value", 
  [IsList, IsBool],
  function(l,checks)
    local c;
    if Length(l) = 0 or not(IsBound(l[1])) then
        Error("CMat: Cannot use two-argument version with empty list and bool");
        return fail;
    fi;
    c := DataObj(l[1]);
    return CMat(l,c,checks);
  end);

InstallMethod( CMat, "for a list of cvecs and a cvecclass", 
  [IsList, IsCVecClass],
  function(l,c)
    return CMat(l,c,true);
  end);

InstallMethod( CMat, "for a compressed GF2 matrix",
  [IsList and IsGF2MatrixRep],
  function(m)
  local c,i,l,v;
  l := 0*[1..NumberRows(m)+1];
  c := CVEC_NewCVecClass(2,1,NumberColumns(m));
  for i in [1..NumberRows(m)] do
      v := ShallowCopy(m[i]);
      PLAIN_GF2VEC(v);
      l[i+1] := CVec(v,c);
  od;
  return CVEC_CMatMaker(l,c);
end);

InstallMethod( CMat, "for a compressed 8bit matrix",
  [IsList and Is8BitMatrixRep],
  function(m)
  local c,i,l,pd,v;
  l := 0*[1..NumberRows(m)+1];
  pd := Factors(Size(DefaultFieldOfMatrix(m)));
  c := CVEC_NewCVecClass(pd[1],Length(pd),NumberColumns(m));
  for i in [1..NumberRows(m)] do
      v := ShallowCopy(m[i]);
      PLAIN_VEC8BIT(v);
      l[i+1] := CVec(v,c);
  od;
  return CVEC_CMatMaker(l,c);
end);

InstallMethod( CMat, "for a list of cvecs, a cvec class, and a boolean value", 
  [IsList,IsCVecClass,IsBool],
  function(l,cl,dochecks)
    local v,ll;
    if dochecks then
        for v in [1..Length(l)] do
            if not(IsBound(l[v])) or not(IsCVecRep(l[v])) or 
               not(IsIdenticalObj(DataObj(l[v]),cl)) then
                Error("CVEC_CMat: Not all list entries are correct vectors");
                return fail;
            fi;
        od;
    fi;
    ll := EmptyPlist(Length(l)+1);
    Add(ll,0);
    Append(ll,l);
    return CVEC_CMatMaker(ll,cl);
  end);

InstallMethod( Matrix, "for a list of cvecs, an integer, and a cmat",
  [IsList, IsInt, IsCMatRep],
  function(l,rl,m)
    local cl,i,li;
    cl := m!.vecclass;
    if rl <> cl![CVEC_IDX_len] then
        cl := CVEC_NewCVecClassSameField(cl,rl);
    fi;
    if Length(l) = 0 then
        li := [0];
        return CVEC_CMatMaker(li,cl);
    fi;
    if IsCVecRep(l[1]) then
        if not(IsIdenticalObj(cl,DataObj(l[1]))) then
            Error("Matrix: cvec not in correct class");
            return fail;
        fi;
        li := [0];
        Append(li,l);
        return CVEC_CMatMaker(li,cl);
    elif IsList(l[1]) then
        li := [0];
        for i in [1..Length(l)] do
            Add(li,CVec(l[i],cl));
        od;
        return CVEC_CMatMaker(li,cl);
    else
        ErrorNoReturn("Matrix for cmats: flat initializer not yet implemented");
    fi;
  end );

InstallMethod( Matrix, "for an empty list, an int and a GF2 matrix",
  [ IsList and IsEmpty, IsPosInt, IsGF2MatrixRep ],
  function(l,n,m) return l; end );

InstallMethod( Matrix, "for an empty list, an int and an 8bit matrix rep",
  [ IsList and IsEmpty, IsPosInt, Is8BitMatrixRep ],
  function(l,n,m) return l; end );

# Some methods to make special matrices:

InstallGlobalFunction( CVEC_ZeroMat, function(arg)
  local c,d,i,l,p,x,y;
  if Length(arg) = 2 then
      y := arg[1];
      c := arg[2];   # this must be a cvec class
      if not(IsInt(y)) and not(IsCVecClass(c)) then
          ErrorNoReturn("Usage: CVEC_ZeroMat( rows, cvecclass)");
      fi;
  elif Length(arg) = 4 then
      y := arg[1];
      x := arg[2];
      p := arg[3];
      d := arg[4];
      if not(IsInt(y) and y >= 0) or
         not(IsInt(x) and x >= 0) or
         not(IsPosInt(p) and IsPrime(p)) or
         not(IsPosInt(d) and d < CVEC_MAXDEGREE) then
          ErrorNoReturn("Usage: CVEC_ZeroMat( rows, cols, p, d )");
      fi;
      c := CVEC_NewCVecClass(p,d,x);
  else
      ErrorNoReturn("Usage: CVEC_ZeroMat( rows, [ cvecclass | cols, p, d ] )");
  fi;
  l := 0*[1..y+1];
  for i in [1..y] do
      l[i+1] := CVEC_NEW(c,c![CVEC_IDX_type]);
  od;
  return CVEC_CMatMaker(l,c);
end );

InstallMethod( ZeroMatrix, "for two integers and a cmat",
  [ IsInt, IsInt, IsCMatRep ],
  function( y, x, m )
    local c;
    if x = m!.vecclass![CVEC_IDX_len] then
        return CVEC_ZeroMat(y,m!.vecclass);
    else
        c := CVEC_NewCVecClassSameField(m!.vecclass,x);
        return CVEC_ZeroMat(y,c);
    fi;
  end );

InstallGlobalFunction( CVEC_IdentityMat, function(arg)
  local c,d,i,l,p,y;
  if Length(arg) = 1 then
      c := arg[1];   # this must be a cvec class
      if not(IsCVecClass(c)) then
          ErrorNoReturn("Usage: CVEC_IdentityMat(cvecclass)");
      fi;
      y := c![CVEC_IDX_len];
  elif Length(arg) = 3 then
      y := arg[1];
      p := arg[2];
      d := arg[3];
      if not(IsInt(y) and y >= 0) or
         not(IsPosInt(p) and IsPrime(p)) or
         not(IsPosInt(d) and d < CVEC_MAXDEGREE) then
          ErrorNoReturn("Usage: CVEC_IdentityMat( rows, p, d )");
      fi;
      c := CVEC_NewCVecClass(p,d,y);
  else
      ErrorNoReturn("Usage: CVEC_IdentityMat( [ cvecclass | rows, p, d ] )");
  fi;
  l := 0*[1..y+1];
  for i in [1..y] do
      l[i+1] := CVEC_NEW(c,c![CVEC_IDX_type]);
      l[i+1,i] := 1;   # note that this works for all fields!
  od;
  return CVEC_CMatMaker(l,c);
end );

InstallMethod( IdentityMatrix, "for an integer and a cmat",
  [ IsInt, IsCMatRep ],
  function( rows, m)
    local c;
    if rows = m!.vecclass![CVEC_IDX_len] then
        return CVEC_IdentityMat(m!.vecclass);
    else
        c := CVEC_NewCVecClassSameField(m!.vecclass,rows);
        return CVEC_IdentityMat(c);
    fi;
  end );

InstallMethod( CompanionMatrix, "for a polynomial and a cmat",
  [ IsUnivariatePolynomial, IsCMatRep ],
  function( po, m )
    local cl,i,l,ll,n;
    l := CoefficientsOfUnivariatePolynomial(po);
    n := Length(l)-1;
    if not(IsOne(l[n+1])) then
        Error("CompanionMatrix: polynomial is not monic");
        return fail;
    fi;
    cl := CVEC_NewCVecClassSameField(m!.vecclass,n);
    ll := 0*[0..n];
    ll[2] := CVEC_NEW(cl,cl![CVEC_IDX_type]);
    ll[2,n] := -l[1];
    for i in [3..n+1] do
        ll[i] := CVEC_NEW(cl,cl![CVEC_IDX_type]);
        ll[i,i-2] := 1;   # this works for all fields!
        ll[i,n] := -l[i-1];
    od;
    return CVEC_CMatMaker(ll,cl);
  end );

# `NewCompanionMatrix` was a constructor
# in GAP <= 4.12 but is a tag based operation in GAP 4.13.
Perform( [ function( ty, po, bd )
    local i,l,ll,n,one,cl;
    one := One(bd);
    l := CoefficientsOfUnivariatePolynomial(po);
    n := Length(l)-1;
    if not(IsOne(l[n+1])) then
        Error("CompanionMatrix: polynomial is not monic");
        return fail;
    fi;
    ll := NewMatrix(ty,bd,n,[]);
    cl := CVEC_NewCVecClassSameField(ll!.vecclass,n);
    Add(ll, CVEC_NEW(cl,cl![CVEC_IDX_type]));
    ll[1,n]:= -l[1];
    for i in [1..n-1] do
        Add(ll, CVEC_NEW(cl,cl![CVEC_IDX_type]));
        ll[i+1,i] := one;
        ll[i+1,n] := -l[i];
    od;
    return ll;
  end ],
  function( method )
    if IS_CONSTRUCTOR( NewCompanionMatrix ) then
      # This holds in GAP <= 4.12.
      InstallMethod( NewCompanionMatrix,
        [ "IsCMatRep", "IsUnivariatePolynomial", "IsRing" ], method );
    elif IsBoundGlobal( "InstallTagBasedMethod" ) then
      # This should hold in GAP 4.13.
      ValueGlobal("InstallTagBasedMethod")( NewCompanionMatrix, IsCMatRep, method );
    else
      # This should not happen.
      Error( "NewCompanionMatrix is neither a constructor not a tag based operation" );
    fi;
  end );

InstallGlobalFunction( CVEC_RandomMat, function(arg)
  local c,d,i,j,l,li,p,q,x,y;
  if Length(arg) = 2 then
      y := arg[1];
      c := arg[2];   # this must be a cvec class
      if not(IsInt(y)) and not(IsCVecClass(c)) then
          ErrorNoReturn("Usage: CVEC_RandomMat( rows, cvecclass)");
      fi;
      x := c![CVEC_IDX_len];
      d := c![CVEC_IDX_fieldinfo]![CVEC_IDX_d];   # used later on
      q := c![CVEC_IDX_fieldinfo]![CVEC_IDX_q];  
  elif Length(arg) = 4 then
      y := arg[1];
      x := arg[2];
      p := arg[3];
      d := arg[4];
      q := p^d;
      if not(IsInt(y) and y >= 0) or
         not(IsInt(x) and x >= 0) or
         not(IsPosInt(p) and IsPrime(p)) or
         not(IsPosInt(d) and d < CVEC_MAXDEGREE) then
          ErrorNoReturn("Usage: CVEC_RandomMat( rows, cols, p, d )");
      fi;
      c := CVEC_NewCVecClass(p,d,x);
  else
      ErrorNoReturn("Usage: CVEC_RandomMat( rows, [ cvecclass | cols, p, d ] )");
  fi;
  l := 0*[1..y+1];
  if c![CVEC_IDX_fieldinfo]![CVEC_IDX_size] <= 1 then    
      # q is an immediate integer
      li := 0*[1..x];
      for i in [1..y] do
          l[i+1] := CVEC_NEW(c,c![CVEC_IDX_type]);
          for j in [1..x] do
              li[j] := Random([0..q-1]);
          od;
          CVEC_INTREP_TO_CVEC(li,l[i+1]);
      od;
  else    # big scalars!
      li := 0*[1..x*d];
      for i in [1..y] do
          l[i+1] := CVEC_NEW(c,c![CVEC_IDX_type]);
          for j in [1..x*d] do
              li[j] := Random([0..p-1]);
          od;
          CVEC_INTREP_TO_CVEC(li,l[i+1]);
      od;
  fi;
  return CVEC_CMatMaker(l,c);
end );

InstallMethod( ChangedBaseDomain, "for a cmat and a finite field",
  [IsCMatRep,IsField and IsFinite],
  function( m, f )
    local cl,i,l;
    cl := CVEC_NewCVecClass(Characteristic(f),DegreeOverPrimeField(f),
                            m!.vecclass![CVEC_IDX_len]);
    l := [0];
    for i in [2..m!.len+1] do
        l[i] := CVec(Unpack(m!.rows[i]),cl);
    od;
    return CVEC_CMatMaker(l,cl);
  end );

InstallMethod( CompatibleVector, "for a cmat",
  [IsCMatRep],
  function( m )
    return CVEC_New(m!.vecclass);
  end );


#############################################################################
# Viewing, Printing, Displaying of cmats:
#############################################################################

InstallMethod( ViewObj, "for a cmat", [IsCMatRep],
function(m)
  local c;
  c := m!.vecclass;
  Print("<");
  if not(IsMutable(m)) then Print("immutable "); fi;
  if HasGreaseTab(m) then Print("greased "); fi;
  Print("cmat ",m!.len,"x",c![CVEC_IDX_len]," over GF(",
        c![CVEC_IDX_fieldinfo]![CVEC_IDX_p],",",
        c![CVEC_IDX_fieldinfo]![CVEC_IDX_d],")>");
end);

InstallMethod( PrintObj, "for a cmat", [IsCMatRep],
function(m)
  local c,i;
  c := m!.vecclass;
  Print("NewMatrix(IsCMatRep,GF(",
        c![CVEC_IDX_fieldinfo]![CVEC_IDX_p],",",
        c![CVEC_IDX_fieldinfo]![CVEC_IDX_d],"),",NumberColumns(m),",[");
  for i in [1..m!.len] do
      Print(Unpack(m!.rows[i+1]),",");
  od;
  Print("])");
end);
  
InstallMethod( Display, "for a cmat", 
  [IsCMatRep and IsFFECollColl],
function(m)
  local i;
  Print("[");
  for i in [1..m!.len] do
      if i <> 1 then Print(" "); fi;
      Display(m!.rows[i+1]);
  od;
  Print("]\n");
end);

InstallGlobalFunction( OverviewMat, function(M)
  local i,j,s,ts,tz,z;
  z := NumberRows(M);
  s := NumberColumns(M);
  tz := QuoInt(z+39,40);
  ts := QuoInt(s+39,40);
  for i in [1..QuoInt(z+tz-1,tz)] do
      for j in [1..QuoInt(s+ts-1,ts)] do
          if IsZero(ExtractSubMatrix(M,[1+(i-1)*tz..Minimum(i*tz,z)],
                                       [1+(j-1)*ts..Minimum(j*ts,s)])) then
              Print(".");
          else
              Print("*");
          fi;
      od;
      Print("\n");
  od;
end );

#############################################################################
# Unpacking:
#############################################################################

InstallMethod( Unpack, "for a cmat", [IsCMatRep],
  function(m)
    local mm,q,i;
    mm := [];
    for i in [2..m!.len+1] do
        Add(mm,Unpack(m!.rows[i]));
    od;
    return mm;
  end );

# (Pseudo) random matrices:

InstallMethod( Randomize, "for a cmat", [ IsCMatRep and IsMutable ],
    m -> Randomize(GlobalMersenneTwister, m) );

InstallOtherMethod( Randomize, "for a cmat and a random source",
  [ IsRandomSource, IsCMatRep and IsMutable ],
  function( m, rs )
    local i;
    for i in [2..m!.len+1] do
      Randomize(rs, m!.rows[i]);
    od;
  end );

# for compatibility with GAP < 4.11
InstallOtherMethod( Randomize, "for a cmat and a random source",
  [ IsCMatRep and IsMutable, IsRandomSource ],
    { m, rs } -> Randomize( rs, m ) );

#############################################################################
# PostMakeImmutable to make subobjects immutable:
#############################################################################

InstallMethod( PostMakeImmutable, "for a cmat", [IsCMatRep],
  function(m)
    MakeImmutable(m!.rows);
  end);

#############################################################################
# Elementary list operations for our matrices:
#############################################################################

InstallOtherMethod( Add, "for a cmat, and a cvec",
  [IsCMatRep and IsMutable, IsCVecRep],
  function(m,v)
    if not(IsIdenticalObj(DataObj(v),m!.vecclass)) then
        Error("Add: only correct cvecs allowed in this matrix");
        return fail;
    fi;
    m!.len := m!.len+1;
    m!.rows[m!.len+1] := v;
  end);
InstallOtherMethod( Add, "for a cmat, a cvec, and a position",
  [IsCMatRep and IsMutable, IsCVecRep, IsPosInt],
  function(m,v,pos)
    if not(IsIdenticalObj(DataObj(v),m!.vecclass)) then
        Error("Add: only correct cvecs allowed in this matrix");
        return fail;
    fi;
    if pos > m!.len+1 then
        Error("Add: position not possible because denseness");
    fi;
    m!.len := m!.len+1;
    Add(m!.rows,v,pos+1);
  end);

InstallOtherMethod( Remove, "for a cmat, and a position",
  [IsCMatRep and IsMutable, IsPosInt],
  function(m,pos)
    if pos < 1 or pos > m!.len then
        Error("Remove: position not possible");
        return fail;
    fi;
    m!.len := m!.len-1;
    return Remove(m!.rows,pos+1);
  end);

#InstallOtherMethod( \[\], "for a cmat, and a position", 
#  [IsCMatRep, IsPosInt],
#  function(m,pos)
#    #if pos < 1 or pos > m!.len then
#    #    Error("\\[\\]: illegal position");
#    #    return fail;
#    #fi;
#    return m!.rows[pos+1];
#  end);
InstallOtherMethod( \[\], "for a cmat, and a position",
  [IsCMatRep, IsPosInt],
  CMAT_ELM_LIST );

InstallOtherMethod( \[\]\:\=, "for a cmat, a position, and a cvec",
  [IsCMatRep and IsMutable, IsPosInt, IsCVecRep],
  function(m,pos,v)
    if pos < 1 or pos > m!.len+1 then
        Error("\\[\\]\\:\\=: illegal position");
    fi;
    if not(IsIdenticalObj(DataObj(v),m!.vecclass)) then
        Error("\\[\\]\\:\\=: can only assign cvecs to cmat");
    fi;
    if pos = m!.len+1 then
        m!.len := m!.len + 1;
    fi;
    m!.rows[pos+1] := v;
  end);

InstallMethod( MatElm, "for a cmat and two positions",
  [IsCMatRep, IsPosInt, IsPosInt],
  function( m, row, col )
    return m!.rows[row+1,col];
  end );

InstallMethod( SetMatElm, "for a cmat, two positions, and an ffe",
  [IsCMatRep and IsMutable, IsPosInt, IsPosInt, IsObject],
  function( m, row, col, el )
    m!.rows[row+1,col] := el;
  end );


# The following two method installations are useful in
# GAP 4.9 and 4.10, but won't be in GAP 4.11, where it
# suffices to install methods for MatElm and SetMatElm.
if not IsBound(ELM_MAT) then

InstallMethod( \[\], "for a cmat and two positions",
  [IsCMatRep, IsPosInt, IsPosInt],
  function( m, row, col )
    return m!.rows[row+1,col];
  end );

InstallMethod( \[\]\:\=, "for a cmat, two positions, and an ffe",
  [IsCMatRep and IsMutable, IsPosInt, IsPosInt, IsObject],
  function( m, row, col, el )
    m!.rows[row+1,col] := el;
  end );

fi;

InstallOtherMethod( \{\}, "for a cmat, and a list",
  [IsCMatRep, IsList],
  function(m,li)
    local l,what;
    what := EmptyPlist(Length(li)+1);
    Add(what,1);
    Append(what,li+1);
    l := m!.rows{what};
    return CVEC_CMatMaker(l,m!.vecclass);
  end);

InstallOtherMethod( \{\}\:\=, "for a cmat, a homogeneous list, and a cmat",
  [IsCMatRep and IsMutable, IsList, 
   IsCMatRep],
  function(m,l,n)
    local i;
    if not(IsIdenticalObj(m!.vecclass,n!.vecclass)) then
        ErrorNoReturn("{}:= : cmats not compatible");
    fi;
    for i in [1..Length(l)] do
        m!.rows[l[i]+1] := n!.rows[i+1];
    od;
  end);

# for backwards compatibility, add Length method for cmat
InstallOtherMethod( Length, "for a cmat",
  [IsCMatRep],
  function(m) return m!.len; end);

InstallOtherMethod( DimensionsMat, "for a cmat",
  [IsCMatRep and IsMatrixObj],
  function(m) return [m!.len,m!.vecclass![2]]; end );

InstallMethod( NumberRows, "for a cmat",
  [IsCMatRep],
  function(m) return m!.len; end);

InstallMethod( NumberColumns, "for a cmat",
  [IsCMatRep and IsMatrixObj],
  function(m) return m!.vecclass![2]; end );

InstallOtherMethod( ShallowCopy, "for a cmat",
  [IsCMatRep],
  function(m) return CVEC_CMatMaker(ShallowCopy(m!.rows),m!.vecclass); end);

InstallOtherMethod( Collected, "for a cmat",
  [IsCMatRep],
  function(m)
    return Collected(m!.rows{[2..m!.len+1]});
  end);

InstallOtherMethod( DuplicateFreeList, "for a cmat",
  [IsCMatRep],
  function(m)
    local l;
    l := DuplicateFreeList(m!.rows);
    return CVEC_CMatMaker(l,m!.vecclass);
  end);

InstallOtherMethod( Append, "for two cmats",
  [IsCMatRep and IsMutable, IsCMatRep],
  function(m1,m2)
      local i;
      if not(IsIdenticalObj(m1!.vecclass,m2!.vecclass)) then
          Error("Append: Incompatible matrices");
          return fail;
      fi;
      for i in [2..m2!.len+1] do
          Add(m1!.rows,m2!.rows[i]);
      od;
      m1!.len := m1!.len + m2!.len;
  end);

InstallOtherMethod( FilteredOp, "for a cmat and a function",
  [IsCMatRep, IsFunction],
  function(m,f)
    local l;
    l := Filtered(m!.rows{[2..m!.len+1]},f);
    Add(l,0,1);
    return CVEC_CMatMaker(l,m!.vecclass);
  end);

InstallOtherMethod( UNB_LIST, "for a cmat and a position",
  [IsCMatRep and IsMutable, IsPosInt],
  function(m,pos)
    if pos = m!.len then
        Unbind(m!.rows[m!.len+1]);
        m!.len := m!.len-1;
    else
        Error("Unbind: not possible for cmats except last entry");
    fi;
  end);


#############################################################################
# CopySubMatrix and ExtractSubMatrix:
#############################################################################

InstallGlobalFunction( CVEC_CopySubMatrix,
function(src,dst,srcli,dstli,srcpos,len,dstpos)
  local i;
  if not(IsIdenticalObj(src!.scaclass,dst!.scaclass)) then
      ErrorNoReturn("CVEC_CopySubMatrix: cmats not over common field");
  fi;
  if Length(srcli) <> Length(dstli) then
      ErrorNoReturn("CVEC_CopySubMatrix: row lists do not have equal lengths");
  fi;
  if srcpos < 1 or srcpos+len-1 > src!.vecclass![CVEC_IDX_len] or len <= 0 then
      ErrorNoReturn("CVEC_CopySubMatrix: source area not valid");
  fi;
  if dstpos < 1 or dstpos+len-1 > dst!.vecclass![CVEC_IDX_len] then
      ErrorNoReturn("CVEC_CopySubMatrix: destination area not valid");
  fi;
  if not(IsMutable(dst)) then
      ErrorNoReturn("CVEC_CopySubMatrix: destination is immutable");
  fi;
  for i in [1..Length(srcli)] do
      CVEC_SLICE(src!.rows[srcli[i]+1],dst!.rows[dstli[i]+1],
                 srcpos,len,dstpos);
  od;
end );

InstallGlobalFunction( CVEC_CopySubMatrixUgly,
function(src,dst,srows,drows,scols,dcols)
  # This handles the ugly case that scols and dcols are no ranges
  # with increment 1, we try to optimize using SLICE. We already
  # know, that they have equal nonzero length:

  local FindRuns,IntersectRuns,c,i,j,r,s;

  FindRuns := function(l)
    # l must be nonempty
    local c,i,r,s;
    r := [];
    i := 2;
    s := l[1];
    c := l[1];
    while i <= Length(l) do
        if l[i] = c+1 then
            c := c + 1;
        else
            Add(r,s);      # The start of the run
            Add(r,c-s+1);  # The length of the run
            c := l[i];
            s := l[i];
        fi;
        i := i + 1;
    od;
    Add(r,s);
    Add(r,c-s+1);
    return r;
  end;

  IntersectRuns := function(l1,l2)
    # Both are nonempty, the result are two refined runs with equal part 
    # lengths. They are given as one list of triples.
    local i1,i2,newrun,r;
    r := [];
    i1 := 1;
    i2 := 1;
    while i1 <= Length(l1) do
        # Note that i1 <= Length(l1) iff i2 <= Length(l2)
        if l1[i1+1] < l2[i2+1] then
            newrun := l1[i1+1];
            Add(r,l1[i1]);
            Add(r,newrun);
            Add(r,l2[i2]);
            l2[i2] := l2[i2] + newrun;
            l2[i2+1] := l2[i2+1] - newrun;
            i1 := i1 + 2;
        elif l1[i1+1] > l2[i2+1] then
            newrun := l2[i2+1];
            Add(r,l1[i1]);
            Add(r,newrun);
            Add(r,l2[i2]);
            l1[i1] := l1[i1] + newrun;
            l1[i1+1] := l1[i1+1] - newrun;
            i2 := i2 + 2;
        else
            newrun := l1[i1+1];
            Add(r,l1[i1]);
            Add(r,newrun);
            Add(r,l2[i2]);
            i1 := i1 + 2;
            i2 := i2 + 2;
        fi;
    od;
    return r;
  end;

  r := [FindRuns(scols),FindRuns(dcols)];
  r := IntersectRuns(r[1],r[2]);
  
  for i in [1..Length(srows)] do
      j := 1;
      while j <= Length(r) do
          CVEC_SLICE(src[srows[i]],dst[drows[i]],r[j],r[j+1],r[j+2]);
          j := j + 3;
      od;
  od;
end );

InstallOtherMethod( CopySubMatrix, "for two cmats and stuff",
  [IsCMatRep, IsCMatRep and IsMutable,
   IsList,IsList,IsList,IsList],
  function( src,dst,srows,drows,scols,dcols )
    if Length(srows) <> Length(drows) then
        Error("CVEC_CopySubMatrix: row lists must have equal length");
    fi;
    if Length(scols) = 0 or Length(srows) = 0 then return; fi;
    if not(ForAll(srows,x->x >= 1 and x <= src!.len) and
           ForAll(drows,x->x >= 1 and x <= dst!.len)) then
        Error("CVEC_CopySubMatrix: row indices out of range");
    fi;
    # These tests ensure that both matrices have at least one row!
    # Not make sure that srows and drows are plain lists:
    if IsRangeRep(srows) then
        srows := 1*srows;  # unpack it!
    fi;
    if IsRangeRep(drows) then
        drows := 1*drows;  # unpack it!
    fi;
    CVEC_COPY_SUBMATRIX(src!.rows,dst!.rows,srows,drows,scols,dcols);
  end );

InstallMethod( ExtractSubMatrix, "for a cmats and stuff",
  [IsCMatRep, IsList, IsList],
  function( mat, rows, cols )
    local i,l,res,vcl;
    vcl := mat!.vecclass;
    if cols = [1..vcl![CVEC_IDX_len]] then
        l := 0*[1..Length(rows)+1];
        for i in [1..Length(rows)] do
            l[i+1] := ShallowCopy(mat!.rows[rows[i]+1]);
        od;
        return CVEC_CMatMaker(l,vcl);
    elif Length(cols) <> vcl![CVEC_IDX_len] then
        vcl := CVEC_NewCVecClassSameField(vcl,Length(cols));
    fi;
    if not(ForAll(rows,x->x >= 1 and x <= mat!.len)) then
        ErrorNoReturn("CVEC_ExtractSubMatrix: row indices out of range");
    fi;
    # Make rows a plain list:
    res := CVEC_ZeroMat(Length(rows),vcl);
    if Length(rows) = 0 then return res; fi;
    if IsRangeRep(rows) then
        rows := 1*rows;
    fi;
    CVEC_COPY_SUBMATRIX( mat!.rows, res!.rows,rows,1*[1..Length(rows)],
                         cols, [1..Length(cols)] );
    return res;
  end );

InstallMethod( ExtractSubMatrix, "for a compressed gf2 matrix",
  [IsMatrix and IsGF2MatrixRep, IsList, IsList],
  function(m, rows, cols)
    local mm;
    mm := m{rows}{cols};
    ConvertToMatrixRep(mm,2);
    return mm;
  end );

InstallMethod( ExtractSubMatrix, "for a compressed 8bit matrix",
  [IsMatrix and Is8BitMatrixRep, IsList, IsList],
  function(m, rows, cols)
    local mm,s;
    mm := m{rows}{cols};
    s := Size(BaseDomain(m));
    ConvertToMatrixRep(mm,s);
    return mm;
  end );


#############################################################################
# Arithmetic for matrices:
#############################################################################

InstallOtherMethod( \+, "for cmats", 
  [IsCMatRep, IsCMatRep],
  function(m,n)
    local l,res,i;
    if not(IsIdenticalObj(m!.vecclass,n!.vecclass)) then
        Error("\\+: cmats not compatible");
    fi;
    if m!.len <> n!.len then
        Error("\\+: cmats do not have equal length");
    fi;
    l := 0*[1..m!.len+1];
    for i in [2..m!.len+1] do
 l[i] := ShallowCopy(m!.rows[i]);
 CVEC_ADD2(l[i],n!.rows[i],0,0);
    od;
    res := CVEC_CMatMaker(l,m!.vecclass);
    if not(IsMutable(m)) and not(IsMutable(n)) then
        MakeImmutable(res);
    fi;
    return res;
  end);

InstallOtherMethod( \-, "for cmats", 
  [IsCMatRep, IsCMatRep],
  function(m,n)
    local l,res,p,i;
    if not(IsIdenticalObj(m!.vecclass,n!.vecclass)) then
        Error("\\-: cmats not compatible");
    fi;
    if m!.len <> n!.len then
        Error("\\-: cmats do not have equal length");
    fi;
    p := m!.vecclass![CVEC_IDX_fieldinfo]![CVEC_IDX_p];
    l := 0*[1..m!.len+1];
    for i in [2..m!.len+1] do
 l[i] := ShallowCopy(m!.rows[i]);
 CVEC_ADDMUL(l[i],n!.rows[i],p-1,0,0);
    od;
    res := CVEC_CMatMaker(l,m!.vecclass);
    if not(IsMutable(m)) and not(IsMutable(n)) then
        MakeImmutable(res);
    fi;
    return res;
  end);

InstallOtherMethod( AdditiveInverseSameMutability, "for a cmat",
  [IsCMatRep],
  function(m)
    local l,res,i;
    l := 0*[1..m!.len+1];
    for i in [2..m!.len+1] do
 l[i] := -m!.rows[i];
    od;
    res := CVEC_CMatMaker(l,m!.vecclass);
    if not(IsMutable(m)) then
        MakeImmutable(res);
    fi;
    return res;
  end);
InstallOtherMethod( AdditiveInverseMutable, "for a cmat",
  [IsCMatRep],
  function(m)
    local l,i;
    l := 0*[1..m!.len+1];
    for i in [2..m!.len+1] do
 l[i] := AdditiveInverseMutable(m!.rows[i]);
    od;
    return CVEC_CMatMaker(l,m!.vecclass);
  end);

InstallOtherMethod( ZeroImmutable, "for a cmat",
  [IsCMatRep],
  function(m)
    local i,l,res,v;
    l := [0];
    v := CVEC_NEW(m!.vecclass,m!.vecclass![CVEC_IDX_type]);
    MakeImmutable(v);
    for i in [2..m!.len+1] do
        l[i] := v;
    od;
    res := CVEC_CMatMaker(l,m!.vecclass);
    MakeImmutable(res);
    return res;
  end);
InstallOtherMethod( ZeroMutable, "for a cmat",
  [IsCMatRep],
  function(m)
    local i,l;
    l := [0];
    for i in [2..m!.len+1] do
        l[i] := CVEC_NEW(m!.vecclass,m!.vecclass![CVEC_IDX_type]);
    od;
    return CVEC_CMatMaker(l,m!.vecclass);
  end);
InstallOtherMethod( ZeroSameMutability, "for a cmat",
  [IsCMatRep],
  function(m)
    if IsMutable(m) then
        return ZeroMutable(m);
    else
        return ZeroImmutable(m);
    fi;
  end);
    
InstallOtherMethod( OneMutable, "for a cmat",
  [IsCMatRep],
  function(m)
    local i,l,one,v,w;
    if m!.vecclass![CVEC_IDX_len] <> m!.len then
        #Error("OneMutable: cmat is not square");
        return fail;
    fi;
    v := CVEC_NEW(m!.vecclass,m!.vecclass![CVEC_IDX_type]);
    l := 0*[1..m!.len+1];
    one := One(m!.scaclass);
    for i in [1..m!.len] do
        w := ShallowCopy(v);
        w[i] := one;
        l[i+1] := w;
    od;
    return CVEC_CMatMaker(l,m!.vecclass);
  end);
InstallOtherMethod( OneSameMutability, "for a cmat",
  [IsCMatRep],
  function(m)
    local n;
    n := OneMutable(m);
    if not(IsMutable(m)) then
        MakeImmutable(n);
    fi;
    return n;
  end);

InstallMethod( OneImmutable, "for a matrix group generated by cmats",
  [ IsMatrixGroup ], 101,   # this has to beat looking at the family
  function ( g )
    local  gens;
    gens := GeneratorsOfGroup( g );
    if Length( gens ) > 0 then
        if IsObjWithMemory(gens[1]) or IsCMatRep(gens[1]) then
            return OneImmutable(gens[1]);
        fi;
    fi;
    TryNextMethod();
  end );

#############################################################################
# Applying Frobenius automorphisms element-wise:
#############################################################################


InstallMethod( \^, "for a cmat and a trivial frobenius automorphism",
  [IsCMatRep and IsFFECollColl, IsMapping and IsOne],
  function( v, f )
    return v;
  end );

InstallMethod( \^, "for a mutable cmat and a trivial frobenius automorphism",
  [IsCMatRep and IsFFECollColl and IsMutable, IsMapping and IsOne],
  function( v, f )
    return MutableCopyMat(v);
  end );

InstallMethod( \^, "for a mutable cmat and a frobenius automorphism",
  [IsCMatRep and IsFFECollColl and IsMutable, IsFrobeniusAutomorphism],
  function( v, f )
    local w,i,j,rl,x,y;
    w := MutableCopyMat(v);
    rl := NumberColumns(w);
    for i in [1..NumberRows(w)] do
        x := v[i];
        y := w[i];
        for j in [1..rl] do
            y[j] := x[j]^f;
        od;
    od;
    return w;
  end );

InstallMethod( \^, "for a cmat and a frobenius automorphism",
  [IsCMatRep and IsFFECollColl, IsFrobeniusAutomorphism],
  function( v, f )
    local w,i,j,rl,x,y;
    w := MutableCopyMat(v);
    rl := NumberColumns(w);
    for i in [1..NumberRows(w)] do
        x := v[i];
        y := w[i];
        for j in [1..rl] do
            y[j] := x[j]^f;
        od;
    od;
    return MakeImmutable(w);
  end );

#############################################################################
# Multiplication with scalars:
#############################################################################

BindGlobal( "CVEC_MATRIX_TIMES_SCALAR", function(m,s)
    local i,l,res;
    l := 0*[1..m!.len+1];
    s := CVEC_HandleScalar(m!.vecclass,s);
    for i in [2..m!.len+1] do 
        l[i] := CVEC_VECTOR_TIMES_SCALAR(m!.rows[i],s); 
    od;
    res := CVEC_CMatMaker(l,m!.vecclass);
    if not(IsMutable(m)) then
        MakeImmutable(res);
    fi;
    return res;
end );
InstallOtherMethod( \*, "for a cmat", [IsCMatRep, IsInt], 
  CVEC_MATRIX_TIMES_SCALAR);
InstallOtherMethod( \*, "for a cmat", [IsCMatRep, IsFFE], 
  CVEC_MATRIX_TIMES_SCALAR);
InstallOtherMethod( \*, "for a cmat", [IsInt,IsCMatRep], 
  function(s,m) return CVEC_MATRIX_TIMES_SCALAR(m,s); end);
InstallOtherMethod( \*, "for a cmat", [IsFFE,IsCMatRep], 
  function(s,m) return CVEC_MATRIX_TIMES_SCALAR(m,s); end);


#############################################################################
# Comparison:
#############################################################################

InstallOtherMethod( \=, "for two cmats",
  [IsCMatRep, IsCMatRep],
  function(m,n)
    local i;
    if not(IsIdenticalObj(m!.vecclass,n!.vecclass)) or m!.len <> n!.len then
        return false;
    fi;
    for i in [2..m!.len+1] do
        if m!.rows[i] <> n!.rows[i] then
            return false;
        fi;
    od;
    return true;
  end);

InstallOtherMethod( \<, "for two cmats",
  [IsCMatRep, IsCMatRep],
  function(m,n)
    local i;
    if not(IsIdenticalObj(m!.vecclass,n!.vecclass)) or m!.len <> n!.len then
        return fail;
    fi;
    for i in [2..m!.len+1] do
        if m!.rows[i] < n!.rows[i] then 
            return true;
        elif n!.rows[i] < m!.rows[i] then
            return false;
        fi;
    od;
    return false;
  end);

InstallOtherMethod( IsZero, "for a cmat", [IsCMatRep],
  function(m)
    return ForAll(m!.rows,IsZero);
  end);

InstallOtherMethod( IsOne, "for a cmat", [IsCMatRep],
  function(m)
    local i,v;
    if m!.vecclass![CVEC_IDX_len] <> m!.len then
        return false;
    fi;
    for i in [1..m!.len] do
        if not(IsOne(m!.rows[i+1,i])) then
            return false;
        fi;
        v := ShallowCopy(m!.rows[i+1]);
        v[i] := 0;
        if not(IsZero(v)) then
            return false;
        fi;
    od;
    return true;
  end );


#############################################################################
# Access to the base field:
#############################################################################

InstallOtherMethod( Characteristic, "for a cmat", [IsCMatRep],
  function(m)
    return m!.vecclass![CVEC_IDX_fieldinfo]![CVEC_IDX_p];
  end);
    
InstallOtherMethod( DegreeFFE, "for a cmat", [IsCMatRep],
  function(m)
    return m!.vecclass![CVEC_IDX_fieldinfo]![CVEC_IDX_d];
  end);
    
InstallMethod( BaseDomain, "for a cmat", [IsCMatRep],
  function(m)
    local c;
    c := m!.vecclass;
    return c![CVEC_IDX_GF];
  end);
    
# compatibility with GAP <= 4.11
if IsBound(BaseField) and not IsIdenticalObj(BaseDomain, BaseField) then
InstallMethod( BaseField, "for a cmat", [IsCMatRep],
  function(m)
    local c;
    c := m!.vecclass;
    return c![CVEC_IDX_GF];
  end);
fi;

InstallMethod(FieldOfMatrixList,
  [IsListOrCollection and IsFFECollCollColl],1,
  function(l)
    local char,deg,m;
    if Length(l) = 0 then
        TryNextMethod();
    fi;
    if not(IsCMatRep(l[1])) then
        TryNextMethod();
    fi;
    deg := 1;
    char := Characteristic(l[1]);
    for m in l do
        deg := Lcm(deg,DegreeFFE(m));
        if char <> Characteristic(m) then
            Error("not all matrices over field with same characteristic");
        fi;
    od;
    return GF(char,deg);
  end);

InstallMethod(DefaultFieldOfMatrix,
  [IsCMatRep and IsFFECollColl],
  function(m)
    local f;
    return m!.vecclass![CVEC_IDX_GF];
  end);

#############################################################################
# The making of good hash functions:
#############################################################################

InstallGlobalFunction( CVEC_HashFunctionForCMats, function(x,data)
  local i,res;
  res := 0;
  for i in [2..x!.len+1] do
      res := (res * 1001 + CVEC_HashFunctionForCVecs(x!.rows[i],data)) 
             mod data[1]+1;
  od;
  return res;
end );

InstallMethod( ChooseHashFunction, "for cmats",
  [IsCMatRep,IsInt],
  function(p,hashlen)
    local bytelen,cl;
    cl := p!.vecclass;
    bytelen := cl![CVEC_IDX_wordlen] * CVEC_BYTESPERWORD;
    return rec( func := CVEC_HashFunctionForCMats,
                data := [hashlen,bytelen] );
  end );


#############################################################################
# Greasing:
#############################################################################

BindGlobal( "CVEC_SpreadTabCache", [] );

InstallGlobalFunction( CVEC_MakeSpreadTab, function(p,d,l,bitsperel)
    # Make up the spreadtab (EXTRACT values are 2^bitsperel-adic
    # expansions with digits only between 0 and p-1):
    local dim,e,i,j,k,mm,pot,spreadtab;
    if IsBound(CVEC_SpreadTabCache[p]) and
       IsBound(CVEC_SpreadTabCache[p][d]) and
       IsBound(CVEC_SpreadTabCache[p][d][l]) then
        return CVEC_SpreadTabCache[p][d][l];
    fi;
    spreadtab := [];
    dim := d*l;
    e := 0*[1..dim+1];
    j := 0;
    mm := 2^bitsperel;
    for i in [0..p^dim-1] do
        spreadtab[j+1] := i+1;
        # Now increment expansion as a p-adic expansion and modify
        # j accordingly as the value of the corresponding m-adic
        # expansion:
        k := 1;
        pot := 1;
        while true do 
            e[k] := e[k] + 1;
            j := j + pot;
            if e[k] < p then break; fi;
            e[k] := 0;
            j := j - p*pot;
            k := k + 1;
            pot := pot * mm;
        od;
    od;
    if not(IsBound(CVEC_SpreadTabCache[p])) then
        CVEC_SpreadTabCache[p] := [];
    fi;
    if not(IsBound(CVEC_SpreadTabCache[p][d])) then
        CVEC_SpreadTabCache[p][d] := [];
    fi;
    CVEC_SpreadTabCache[p][d][l] := spreadtab;
    return spreadtab;
end );

InstallOtherMethod( GreaseMat, "for a cmat",
  [IsCMatRep],
  function(m)
    if m!.vecclass![CVEC_IDX_fieldinfo]![CVEC_IDX_bestgrease] = 0 then
        Info(InfoWarning,1,"GreaseMat: bestgrease is 0, we do not grease");
        return;
    fi;
    GreaseMat(m,m!.vecclass![CVEC_IDX_fieldinfo]![CVEC_IDX_bestgrease]);
  end);

InstallMethod( GreaseMat, "for a cmat, and a level", 
  [IsCMatRep, IsInt],
  function(m,l)
    local bitsperel,d,dim,e,f,i,j,k,mm,nrblocks,p,pot,q,tablen;
    f := m!.vecclass![CVEC_IDX_fieldinfo];   # the field info
    bitsperel := f![CVEC_IDX_bitsperel];
    p := f![CVEC_IDX_p];
    d := f![CVEC_IDX_d];
    q := f![CVEC_IDX_q];
    nrblocks := QuoInt(m!.len+l-1,l);  # we do grease the last <l rows!
    tablen := q^l;  # = p^(d*l)
    m!.greaselev := l;
    m!.greaseblo := nrblocks;
    m!.greasetab := 0*[1..nrblocks];
    for i in [1..nrblocks] do
        m!.greasetab[i] := 0*[1..tablen+1+l];
        for j in [1..tablen+1+l] do
            m!.greasetab[i][j] := 
                CVEC_NEW(m!.vecclass,m!.vecclass![CVEC_IDX_type]);
        od;
        CVEC_FILL_GREASE_TAB(m!.rows,2+(i-1)*l,l,m!.greasetab[i],tablen,1);
    od;

    m!.spreadtab := CVEC_MakeSpreadTab(p,d,l,bitsperel);

    # Finally change the type:
    SetFilterObj(m,HasGreaseTab);
  end); 

InstallMethod( UnGreaseMat, "for a cmat",
  [IsCMatRep],
  function(m)
    ResetFilterObj(m,HasGreaseTab);
    Unbind(m!.greasetab);
    Unbind(m!.greaselev);
    Unbind(m!.greaseblo);
    Unbind(m!.spreadtab);
  end);

InstallGlobalFunction( CVEC_OptimizeGreaseHint, function(m,nr)
  local l,li,q;
  q := m!.vecclass![CVEC_IDX_fieldinfo]![CVEC_IDX_q];
  li := [QuoInt(nr*(q-1)*m!.len + (q-1),q)];
  l := 1;
  while l < 12 do
      li[l+1] := QuoInt(m!.len + (l-1),l)*(nr+q^l);
      if l > 1 and li[l+1] > li[l] then break; fi;
      l := l + 1;
  od;
  if li[l] < li[1] then
      m!.greasehint := l-1;
  else
      m!.greasehint := 0;
  fi;
  #Print("OptimizeGreaseHint: ",li," ==> ",m!.greasehint,"\n");
end );


#############################################################################
# Arithmetic between vectors and matrices, especially multiplication:
#############################################################################
    
InstallOtherMethod(\*, "for a cvec and a cmat, without greasing",
  [IsCVecRep, IsCMatRep],
  function(v,m)
    local i,res,vcl,s,z;
    vcl := DataObj(v);
    if not(IsIdenticalObj(vcl![CVEC_IDX_fieldinfo],
                          m!.vecclass![CVEC_IDX_fieldinfo])) then
        Error("\\*: incompatible base fields");
    fi;
    if Length(v) <> m!.len then
        Error("\\*: lengths not equal");
    fi;
    res := CVEC_NEW(m!.vecclass,m!.vecclass![CVEC_IDX_type]);  # the result
    CVEC_PROD_CVEC_CMAT_NOGREASE(res,v,m!.rows);
    if not(IsMutable(v) or IsMutable(m)) then
        MakeImmutable(res);
    fi;
    return res;
  end);
 
InstallOtherMethod(\^, "for a cvec and a cmat, without greasing",
  [IsCVecRep, IsCMatRep],
  function(v,m)
    local i,res,vcl,s,z;
    vcl := DataObj(v);
    if not(IsIdenticalObj(vcl![CVEC_IDX_fieldinfo],
                          m!.vecclass![CVEC_IDX_fieldinfo])) then
        Error("\\^: incompatible base fields");
    fi;
    if Length(v) <> m!.len then
        Error("\\^: lengths not equal");
    fi;
    res := CVEC_NEW(m!.vecclass,m!.vecclass![CVEC_IDX_type]);  # the result
    CVEC_PROD_CVEC_CMAT_NOGREASE(res,v,m!.rows);
    if not(IsMutable(v) or IsMutable(m)) then
        MakeImmutable(res);
    fi;
    return res;
  end);
 
InstallOtherMethod(\*, "for a cvec and a greased cmat",
  [IsCVecRep, IsCMatRep and HasGreaseTab],
  function(v,m)
    local i,res,vcl,l,pos,val;
    vcl := DataObj(v);
    if not(IsIdenticalObj(vcl![CVEC_IDX_fieldinfo],
                          m!.vecclass![CVEC_IDX_fieldinfo])) then
        Error("\\*: incompatible base fields");
    fi;
    if Length(v) <> m!.len then
        Error("\\*: lengths not equal");
    fi;
    res := CVEC_NEW(m!.vecclass,m!.vecclass![CVEC_IDX_type]);  # the result
    CVEC_PROD_CVEC_CMAT_GREASED(res,v,m!.greasetab,m!.spreadtab,m!.greaselev);
    if not(IsMutable(v) or IsMutable(m)) then
        MakeImmutable(res);
    fi;
    return res;
  end);
 
InstallOtherMethod(\^, "for a cvec and a greased cmat",
  [IsCVecRep, IsCMatRep and HasGreaseTab],
  function(v,m)
    local i,res,vcl,l,pos,val;
    vcl := DataObj(v);
    if not(IsIdenticalObj(vcl![CVEC_IDX_fieldinfo],
                          m!.vecclass![CVEC_IDX_fieldinfo])) then
        Error("\\^: incompatible base fields");
    fi;
    if Length(v) <> m!.len then
        Error("\\^: lengths not equal");
    fi;
    res := CVEC_NEW(m!.vecclass,m!.vecclass![CVEC_IDX_type]);  # the result
    CVEC_PROD_CVEC_CMAT_GREASED(res,v,m!.greasetab,m!.spreadtab,m!.greaselev);
    if not(IsMutable(v) or IsMutable(m)) then
        MakeImmutable(res);
    fi;
    return res;
  end);
 
InstallOtherMethod(\*, "for two cmats, second one not greased",
  [IsCMatRep, IsCMatRep],
  CVEC_PROD_CMAT_CMAT_DISPATCH );

InstallGlobalFunction( CVEC_PROD_CMAT_CMAT_BIG,
  function(m,n)
    local d,greasetab,i,j,l,lev,q,res,spreadtab,tab,tablen,vcl;
    #if not(IsIdenticalObj(m!.scaclass,n!.scaclass)) then
    #    Error("\\*: incompatible base fields");
    #fi;
    #if m!.vecclass![CVEC_IDX_len] <> n!.len then
    #    Error("\\*: lengths not matching");
    #fi;
    vcl := n!.vecclass;
    #max := Maximum(m!.len,m!.vecclass![CVEC_IDX_len],vcl![CVEC_IDX_len]);
    q := vcl![CVEC_IDX_fieldinfo]![CVEC_IDX_q];
    #if max <= 512 and q = 2 then
    #    # Make a new matrix and then go directly into the kernel:
    #    l := 0*[1..m!.len+1];
    #    for i in [2..m!.len+1] do
    #        l[i] := CVEC_NEW(vcl,vcl![CVEC_IDX_type]);
    #    od;
    #    res := CVEC_CMatMaker(l,n!.vecclass);
    #    CVEC_PROD_CMAT_CMAT_GF2_SMALL(l,m!.rows,n!.rows,max);
    #    if not(IsMutable(m) or IsMutable(n)) then
    #        MakeImmutable(res);
    #    fi;
    #    return res;
    #fi;
    if IsBound(CVEC_WinogradBounds[q]) and
        m!.len * m!.vecclass![CVEC_IDX_len] >= CVEC_WinogradBounds[q] and
        n!.len * n!.vecclass![CVEC_IDX_len] >= CVEC_WinogradBounds[q] then
        # Do the Winograd trick:
        res := CVEC_MultiplyWinograd(m,n,CVEC_WinogradBounds[q]);
        if not(IsMutable(m) or IsMutable(n)) then
            MakeImmutable(res);
        fi;
        return res;
    fi;
         
    # First make a new matrix:
    #l := 0*[1..m!.len+1];
    #for i in [2..m!.len+1] do
    #    l[i] := CVEC_NEW(vcl,vcl![CVEC_IDX_type]);
    #od;
    #res := CVEC_CMatMaker(l,n!.vecclass);
    res := CVEC_MAKE_ZERO_CMAT(m!.len,n!.vecclass);
    l := res!.rows;
    if m!.len > 0 then
        d := vcl![CVEC_IDX_fieldinfo]![CVEC_IDX_d];
        if q > CVEC.MaximumGreaseCalibration or
           not(IsBound(CVEC_CalibrationTableNoCache[q])) then
            lev := 0;
        else
            i := 1;
            if vcl![CVEC_IDX_wordlen] <= 32768 then
                tab := CVEC_CalibrationTableCache[q];
            else
                tab := CVEC_CalibrationTableNoCache[q];
            fi;
            while i <= Length(tab) and m!.len >= tab[i] do i := i + 1; od;
            lev := i-1;
        fi;
        if lev = 0 then
            # no greasing at all in this case!
            CVEC_PROD_CMAT_CMAT_NOGREASE2(l,m!.rows,n!.rows);
            # we use version 2, which unpacks full rows of m instead of
            # extracting single field entries.
        else
            spreadtab := CVEC_MakeSpreadTab(
                 vcl![CVEC_IDX_fieldinfo]![CVEC_IDX_p],
                 vcl![CVEC_IDX_fieldinfo]![CVEC_IDX_d],
                 lev, vcl![CVEC_IDX_fieldinfo]![CVEC_IDX_bitsperel]);
            tablen := vcl![CVEC_IDX_fieldinfo]![CVEC_IDX_q]^lev;
            greasetab := 0*[1..tablen+1+lev];
            for j in [1..tablen+1+lev] do
              greasetab[j] := CVEC_NEW(n!.vecclass,n!.vecclass![CVEC_IDX_type]);
            od;
            CVEC_PROD_CMAT_CMAT_WITHGREASE(l,m!.rows,n!.rows,greasetab,
                                           spreadtab,lev);
        fi;
    fi;
    if not(IsMutable(m) or IsMutable(n)) then
        MakeImmutable(res);
    fi;
    return res;
  end);

InstallOtherMethod(\*, "for two cmats, second one greased",
  [IsCMatRep, IsCMatRep and HasGreaseTab],
  function(m,n)
    local i,l,res,vcl,q;
    if not(IsIdenticalObj(m!.scaclass,n!.scaclass)) then
        Error("\\*: incompatible base fields");
    fi;
    if m!.vecclass![CVEC_IDX_len] <> n!.len then
        Error("\\*: lengths not matching");
    fi;
    vcl := n!.vecclass;
    q := vcl![CVEC_IDX_fieldinfo]![CVEC_IDX_q];
    if IsBound(CVEC_WinogradBounds[q]) and
        m!.len * m!.vecclass![CVEC_IDX_len] >= CVEC_WinogradBounds[q] and
        n!.len * n!.vecclass![CVEC_IDX_len] >= CVEC_WinogradBounds[q] then
        # Do the Winograd trick:
        res := CVEC_MultiplyWinograd(m,n,CVEC_WinogradBounds[q]);
        if not(IsMutable(m) or IsMutable(n)) then
            MakeImmutable(res);
        fi;
        return res;
    fi;
         
    # First make a new matrix:
    l := 0*[1..m!.len+1];
    for i in [2..m!.len+1] do
        l[i] := CVEC_NEW(vcl,vcl![CVEC_IDX_type]);
    od;
    res := CVEC_CMatMaker(l,n!.vecclass);
    if m!.len > 0 then
        CVEC_PROD_CMAT_CMAT_GREASED(l,m!.rows,n!.greasetab,n!.spreadtab,
                                    n!.len,n!.greaselev);
    fi;
    if not(IsMutable(m)) and not(IsMutable(n)) then
        MakeImmutable(res);
    fi;
    return res;
  end);

#############################################################################
# Inversion of matrices:
#############################################################################

BindGlobal( "CVEC_INVERT_FFE",function(helper)
  helper[1] := helper[1]^-1;
end );

InstallGlobalFunction( CVEC_InverseWithoutGrease, function(m)
    # Now make a new identity matrix:
    local helper,i,l,mc,mi,vcl;
    vcl := m!.vecclass;
    l := [0];
    for i in [m!.len+1,m!.len..2] do
        l[i] := CVEC_NEW(vcl,vcl![CVEC_IDX_type]);
        l[i,i-1] := 1;   # note that this works for all fields!
    od;
    mi := CVEC_CMatMaker(l,vcl);
    # Now make a copy of the matrix:
    mc := MutableCopyMat(m);

    # Now do the real work:
    helper := CVEC_New(vcl);
    i := CVEC_CMAT_INVERSE(mi!.rows,mc!.rows,CVEC_INVERT_FFE,helper);

    if i <> fail then
        return mi;
    else
        return fail;
    fi;
  end );

InstallGlobalFunction (CVEC_InverseWithGrease,
  function(m,lev)
    local greasetab1,greasetab2,helper,i,l,mc,mi,spreadtab,tablen,vcl;
    vcl := m!.vecclass;
    if m!.len <> vcl![CVEC_IDX_len] then return fail; fi;
    if m!.len = 0 then return fail; fi;
    if m!.len = 1 then
        l := [0,CVEC_New(vcl)];
        i := m!.rows[2,1]^-1;
        if i = fail then
            return fail;
        fi;
        l[2,1] := i;
        return CVEC_CMatMaker(l,m!.vecclass);
    fi;
    # Now make a new identity matrix:
    l := [0];
    for i in [m!.len+1,m!.len..2] do
        l[i] := CVEC_NEW(vcl,vcl![CVEC_IDX_type]);
        l[i,i-1] := 1;   # note that this works for all fields!
    od;
    mi := CVEC_CMatMaker(l,vcl);
    # Now make a copy of the matrix:
    mc := MutableCopyMat(m);

    # Prepare to grease:
    spreadtab := CVEC_MakeSpreadTab(
         vcl![CVEC_IDX_fieldinfo]![CVEC_IDX_p],
         vcl![CVEC_IDX_fieldinfo]![CVEC_IDX_d],
         lev, vcl![CVEC_IDX_fieldinfo]![CVEC_IDX_bitsperel]);
    tablen := vcl![CVEC_IDX_fieldinfo]![CVEC_IDX_q]^lev;
    greasetab1 := 0*[1..tablen+1+lev];
    greasetab2 := 0*[1..tablen+1+lev];
    for i in [1..tablen+1+lev] do
      greasetab1[i] := CVEC_NEW(vcl,vcl![CVEC_IDX_type]);
      greasetab2[i] := CVEC_NEW(vcl,vcl![CVEC_IDX_type]);
    od;

    # Now do the real work:
    helper := CVEC_New(vcl);
    i := CVEC_CMAT_INVERSE_GREASE(mi!.rows,mc!.rows,CVEC_INVERT_FFE,helper,
                                  [greasetab1,greasetab2,spreadtab,lev,tablen]);

    if i <> fail then
        return mi;
    else
        return fail;
    fi;
  end );

InstallOtherMethod( InverseMutable, "for a square cmat",
  [IsCMatRep],
  function(m)
    local i,l,vcl;
    vcl := m!.vecclass;
    if m!.len <> vcl![CVEC_IDX_len] then return fail; fi;
    if m!.len = 0 then return fail; fi;
    if m!.len = 1 then
        l := [0,CVEC_New(vcl)];
        i := m!.rows[2,1]^-1;
        if i = fail then
            return fail;
        fi;
        l[2,1] := i;
        return CVEC_CMatMaker(l,m!.vecclass);
    fi;
    if m!.greasehint = 0 or m!.len < 100 then
        return CVEC_InverseWithoutGrease(m);
    else
        return CVEC_InverseWithGrease(m,m!.greasehint);
    fi;
  end );

InstallOtherMethod( InverseSameMutability, "for a square cmat",
  [IsCMatRep],
  function(m)
    local mi;
    mi := InverseMutable(m);
    if mi <> fail and not(IsMutable(m)) then
        MakeImmutable(mi);
    fi;
    return mi;
  end );


#############################################################################
# Transposition:
#############################################################################

InstallOtherMethod( TransposedMatOp, "for a cmat",
  [IsCMatRep],
  function(m)
    # First make a new matrix:
    local c,ct,i,l,mt,newlen;
    c := m!.vecclass;
    ct := CVEC_NewCVecClassSameField(c,m!.len);
    newlen := c![CVEC_IDX_len];
    l := 0*[1..newlen+1];
    for i in [2..newlen+1] do
        l[i] := CVEC_NEW(ct,ct![CVEC_IDX_type]);
    od;
    mt := CVEC_CMatMaker(l,ct);
    if m!.len > 0 and mt!.len > 0 then
        CVEC_TRANSPOSED_MAT(m!.rows,mt!.rows);
    fi;
    return mt;
  end);

InstallOtherMethod( TransposedMat, "for a cmat",
  [IsCMatRep],
  function(m)
    local mt;
    mt := TransposedMatOp(m);
    MakeImmutable(mt);
    return mt;
  end);

#############################################################################
# I/O for Matrices:
#############################################################################

BindGlobal( "CVEC_64BIT_NUMBER_TO_STRING_LITTLE_ENDIAN", function(n)
  local i,s;
  s := "        ";
  for i in [1..8] do
      s[i] := CHAR_INT(RemInt(n,256));
      n := QuoInt(n,256);
  od;
  return s;
end );

InstallGlobalFunction( CVEC_WriteMat, function(f,m)
  local buf,c,chead,dhead,header,i,magic,phead,rhead;
  if not(IsFile(f)) then
      Error("CVEC_WriteMat: first argument must be a file");
      return fail;
  fi;
  if not(IsCMatRep(m)) then
      Error("CVEC_WriteMat: second argument must be a cmat");
      return fail;
  fi;
  c := m!.vecclass;
  Info(InfoCVec,2,"CVEC_WriteMat: Writing ",m!.len,"x",
       c![CVEC_IDX_len]," matrix over GF(",
       c![CVEC_IDX_fieldinfo]![CVEC_IDX_p],",",
       c![CVEC_IDX_fieldinfo]![CVEC_IDX_d],")");
  # Make the header:
  magic := "GAPCMat1";
  phead := CVEC_64BIT_NUMBER_TO_STRING_LITTLE_ENDIAN(
           c![CVEC_IDX_fieldinfo]![CVEC_IDX_p]);
  dhead := CVEC_64BIT_NUMBER_TO_STRING_LITTLE_ENDIAN(
           c![CVEC_IDX_fieldinfo]![CVEC_IDX_d]);
  rhead := CVEC_64BIT_NUMBER_TO_STRING_LITTLE_ENDIAN(m!.len);
  chead := CVEC_64BIT_NUMBER_TO_STRING_LITTLE_ENDIAN(
           c![CVEC_IDX_len]);
  header := Concatenation(magic,phead,dhead,rhead,chead);
  if IO_Write(f,header) <> 40 then
      Info(InfoCVec,1,"CVEC_WriteMat: Write error during writing of header");
      return fail;
  fi;
  buf := "";  # will be made longer automatically
  for i in [1..m!.len] do
      CVEC_CVEC_TO_EXTREP(m!.rows[i+1],buf);
      if IO_Write(f,buf) <> Length(buf) then
          Info(InfoCVec,1,"CVEC_WriteMat: Write error");
          return fail;
      fi;
  od;
  return true;
end );

InstallGlobalFunction( CVEC_WriteMatToFile, function(fn,m)
  local f;
  f := IO_File(fn,"w");
  if f = fail then
      Info(InfoCVec,1,"CVEC_WriteMatToFile: Cannot create file");
      return fail;
  fi;
  if CVEC_WriteMat(f,m) = fail then return fail; fi;
  if IO_Close(f) = fail then
      Info(InfoCVec,1,"CVEC_WriteMatToFile: Cannot close file");
      return fail;
  fi;
  return true;
end );

InstallGlobalFunction( CVEC_WriteMatsToFile, function(fnpref,l)
  local i;
  if not(IsString(fnpref)) then
      Error("CVEC_WriteMatsToFile: fnpref must be a string");
      return fail;
  fi;
  if not(IsList(l)) then
      Error("CVEC_WriteMatsToFile: l must be list");
      return fail;
  fi;
  for i in [1..Length(l)] do
      if CVEC_WriteMatToFile(Concatenation(fnpref,String(i)),l[i]) = fail then
          return fail;
      fi;
  od;
  return true;
end );

BindGlobal( "CVEC_STRING_LITTLE_ENDIAN_TO_64BIT_NUMBER", function(s)
  local i,n;
  n := 0;
  for i in [8,7..1] do
      n := n * 256 + INT_CHAR(s[i]);
  od;
  return n;
end );

InstallGlobalFunction( CVEC_ReadMat, function(f)
  local buf,c,cols,d,header,i,len,m,p,rows;
  if not(IsFile(f)) then
      Error("CVEC_ReadMat: first argument must be a file");
      return fail;
  fi;
  header := IO_ReadBlock(f,40);
  if Length(header) < 40 then
      Info(InfoCVec,1,"CVEC_ReadMat: Cannot read header");
      return fail;
  fi;

  # Check and process the header:
  if header{[1..8]} <> "GAPCMat1" then
      Info(InfoCVec,1,"CVEC_ReadMat: Magic of header incorrect");
      return fail;
  fi;
  p := CVEC_STRING_LITTLE_ENDIAN_TO_64BIT_NUMBER(header{[9..16]});
  d := CVEC_STRING_LITTLE_ENDIAN_TO_64BIT_NUMBER(header{[17..24]});
  rows := CVEC_STRING_LITTLE_ENDIAN_TO_64BIT_NUMBER(header{[25..32]});
  cols := CVEC_STRING_LITTLE_ENDIAN_TO_64BIT_NUMBER(header{[33..40]});
  Info(InfoCVec,2,"CVEC_ReadMat: Reading ",rows,"x",cols," matrix over GF(",
       p,",",d,")");
   
  c := CVEC_NewCVecClass(p,d,cols);
  m := CVEC_ZeroMat(rows,c);
  buf := "";  # will be made longer automatically
  if rows > 0 then
      CVEC_CVEC_TO_EXTREP(m!.rows[2],buf);   # to get the length right
      len := Length(buf);
  else
      len := 0;
  fi;

  for i in [1..rows] do
      buf := IO_ReadBlock(f,len);
      if len <> Length(buf) then
          Info(InfoCVec,1,"CVEC_ReadMat: Read error");
          Error();
          return fail;
      fi;
      CVEC_EXTREP_TO_CVEC(buf,m!.rows[i+1]);
  od;
  return m;
end );

InstallGlobalFunction( CVEC_ReadMatFromFile, function(fn)
  local f,m;
  f := IO_File(fn,"r");
  if f = fail then
      Info(InfoCVec,1,"CVEC_ReadMatFromFile: Cannot open file");
      return fail;
  fi;
  m := CVEC_ReadMat(f);
  if m = fail then return fail; fi;
  IO_Close(f);
  return m;
end );

InstallGlobalFunction( CVEC_ReadMatsFromFile, function(fnpref)
  local f,i,l,m;
  if not(IsString(fnpref)) then
      Error("CVEC_ReadMatsFromFile: fnpref must be a string");
      return fail;
  fi;
  f := IO_File(Concatenation(fnpref,"1"),"r");
  if f = fail then
      Error("CVEC_ReadMatsFromFile: no matrices there");
      return fail;
  else
      IO_Close(f);
  fi;
  l := [];
  i := 1;
  while true do
      f := IO_File(Concatenation(fnpref,String(i)),"r");
      if f = fail then break; fi;
      IO_Close(f);
      m := CVEC_ReadMatFromFile(Concatenation(fnpref,String(i)));
      if m = fail then
          return fail;
      else
          Add(l,m);
          i := i + 1;
      fi;
  od;
  return l;
end );

#############################################################################
# Further handling of matrices: 
#############################################################################

InstallOtherMethod( PositionNonZero, "for a cmat",
  [ IsCMatRep and IsMatrixObj ],
  function(m)
    local i;
    i := 1;
    while i <= m!.len and IsZero(m!.rows[i+1]) do i := i + 1; od;
    return i;
  end );

InstallOtherMethod( PositionNonZero, "for a cmat, and an integer",
  [ IsCMatRep, IsInt ],
  function(m,j)
    local i;
    i := Maximum(j+1,1);
    while i <= m!.len and IsZero(m!.rows[i+1]) do i := i + 1; od;
    if i > m!.len then
        return m!.len + 1;
    else
        return i;
    fi;
  end );

InstallOtherMethod( PositionLastNonZero, "for a cmat",
  [ IsCMatRep and IsMatrixObj ],
  function(m)
    local i;
    i := m!.len;;
    while i >= 1 and IsZero(m!.rows[i+1]) do i := i - 1; od;
    return i;
  end );

InstallOtherMethod( PositionLastNonZero, "for a cmat, and an integer",
  [ IsCMatRep and IsMatrixObj, IsInt ],
  function(m,j)
    local i;
    i := Minimum(j-1,m!.len);
    while i >= 1 and IsZero(m!.rows[i+1]) do i := i - 1; od;
    return i;
  end );

InstallOtherMethod( IsDiagonalMat, "for a cmat", [IsCMatRep],
  function(m)
    local i,mi;
    mi := Minimum(m!.len,m!.vecclass![CVEC_IDX_len]);
    i := 1;
    while i <= mi do
        if PositionNonZero(m!.rows[i+1]) < i or
           PositionLastNonZero(m!.rows[i+1]) > i then
            return false;
        fi;
        i := i + 1;
    od;
    while i <= m!.len do
        if not(IsZero(m!.rows[i+1])) then
            return false;
        fi;
        i := i + 1;
    od;
    return true;
  end);

InstallOtherMethod( IsUpperTriangularMat, "for a cmat", 
  [IsCMatRep],
  function(m)
    local i,mi;
    mi := Minimum(m!.len,m!.vecclass![CVEC_IDX_len]);
    i := 1;
    while i <= mi do
        if PositionNonZero(m!.rows[i+1]) < i then
            return false;
        fi;
        i := i + 1;
    od;
    while i <= m!.len do
        if not(IsZero(m!.rows[i+1])) then
            return false;
        fi;
        i := i + 1;
    od;
    return true;
  end);

InstallOtherMethod( IsLowerTriangularMat, "for a cmat", 
  [IsCMatRep],
  function(m)
    local i,mi;
    mi := Minimum(m!.len,m!.vecclass![CVEC_IDX_len]);
    i := 1;
    while i <= mi do
        if PositionLastNonZero(m!.rows[i+1]) > i then
            return false;
        fi;
        i := i + 1;
    od;
    while i <= m!.len do
        if not(IsZero(m!.rows[i+1])) then
            return false;
        fi;
        i := i + 1;
    od;
    return true;
  end);

#############################################################################
# Copying of matrices:
#############################################################################

InstallOtherMethod( MutableCopyMat, "for a cmat", 
  [IsCMatRep and IsMatrixObj],
  function(m)
    local l,i;
    l := 0*[1..m!.len+1];
    for i in [2..m!.len+1] do
        l[i] := ShallowCopy(m!.rows[i]);
    od;
    Unbind(l[1]);
    return CVEC_CMatMaker(l,m!.vecclass);
  end);

#############################################################################
# KroneckerProduct:
#############################################################################

# The generic method suffices, however, for compatibility with
# GAP < 4.5 we install it here anew:

InstallOtherMethod( KroneckerProduct, "for cmats", 
               [ IsCMatRep and IsMatrixObj, IsCMatRep and IsMatrixObj ],
  function( A, B )
    local rowsA, rowsB, colsA, colsB, newclass, AxB, i, j;
      rowsA := NumberRows(A);
      colsA := NumberColumns(A);
      rowsB := NumberRows(B);
      colsB := NumberColumns(B);

      AxB := ZeroMatrix( rowsA * rowsB, colsA * colsB, A );

      # Cache matrices
      # not implemented yet

      for i in [1..rowsA] do
 for j in [1..colsA] do
   CopySubMatrix( A[i,j] * B, AxB, 
    [ 1 .. rowsB ], [ rowsB * (i-1) + 1 .. rowsB * i ],
    [ 1 .. colsB ], [ (j-1) * colsB + 1 .. j * colsB ] );
 od;
      od;
      return AxB;
    end );

#############################################################################
# (Un-)Pickling:
#############################################################################

InstallMethod( IO_Pickle, "for cmats",
  [IsFile, IsCMatRep and IsList],
  function( f, m )
    local tag;
    if IsMutable(m) then tag := "MCMA";
    else tag := "ICMA"; fi;
    if IO_Write(f,tag) = fail then return IO_Error; fi;
    if CVEC_WriteMat( f, m ) = fail then return IO_Error; fi;
    return IO_OK;
  end );

IO_Unpicklers.MCMA :=
  function( f )
    local m;
    m := CVEC_ReadMat( f );
    if m = fail then return IO_Error; fi;
    return m;
  end;

IO_Unpicklers.ICMA :=
  function( f )
    local m;
    m := CVEC_ReadMat( f );
    if m = fail then return IO_Error; fi;
    MakeImmutable(m);
    return m;
  end;


#############################################################################
# Memory usage information:
#############################################################################

InstallOtherMethod( Memory, "for a cmat", [ IsCMatRep ],
  function( m )
    local bpw,bpb;
    # Bytes per word:
    bpw := GAPInfo.BytesPerVariable;
    # Bytes per bag (in addition to content):
    bpb := 8 + 2*bpw;   # this counts the header and the master pointer!
    if NumberRows(m) = 0 then
        return 2*bpb + SIZE_OBJ(m) + SIZE_OBJ(m!.rows);
    else
        return 2*bpb + SIZE_OBJ(m) + SIZE_OBJ(m!.rows)
               + NumberRows(m) * Memory(m!.rows[2]);
    fi;
    # FIXME: this does not include greased data!
  end );


#############################################################################
# Grease calibration:
#############################################################################

CVEC.MaximumGreaseCalibration := 1024;

InstallGlobalFunction( CVEC_ComputeVectorLengthsForCalibration,
  function()
    local bpw,cl,epw,f,i,l,le,lencache,lennocache;
    l := Filtered([2..CVEC.MaximumGreaseCalibration],IsPrimePowerInt);
    le := Length(l);
    for i in [1..le] do
        f := Factors(l[i]);
        l[i] := [f[1],Length(f),l[i]];
    od;
    cl := List([1..le],i->CVecClass(l[i][1],l[i][2],1));
    lencache := 0*[1..le];
    lennocache := 0*[1..le];
    bpw := GAPInfo.BytesPerVariable;
    for i in [1..le] do
        epw := cl[i]![CVEC_IDX_fieldinfo]![CVEC_IDX_elsperword];
        # 64kB is surely in the cache:
        lencache[i] := QuoInt(65536 / bpw * epw,
                              cl[i]![CVEC_IDX_fieldinfo]![CVEC_IDX_d]);
        # 16MB is supposedly out of the cache:
        lennocache[i] := QuoInt(16*1024*1024 / bpw * epw,
                                cl[i]![CVEC_IDX_fieldinfo]![CVEC_IDX_d]);
    od;
    return rec( pdq := l, le := le, lencache := lencache, 
                lennocache := lennocache );
  end );

InstallGlobalFunction( CVEC_FastFill,
  function( v )
    local cl,d,e,i,l,p,q,x;
    cl := CVecClass(v);
    p := cl![CVEC_IDX_fieldinfo]![CVEC_IDX_p];
    d := cl![CVEC_IDX_fieldinfo]![CVEC_IDX_d];
    q := cl![CVEC_IDX_fieldinfo]![CVEC_IDX_q];
    x := 1;
    l := Length(v);
    for i in [1..Minimum(Length(v),1000)] do
        v[i] := x;
        x := x + 1;
        if x = q then x := 1; fi;
    od;
    i := 1;
    while 1000*i+1 <= Length(v) do
        e := Minimum(Length(v),1000*(i+1));
        CopySubVector(v,v,[1..e-1000*i],[1000*i+1..e]);
        i := i + 1;
    od;
  end );

InstallGlobalFunction( GreaseCalibration,
  function()
    local cl,d,gd,i,info,inters,j,l,limits,minpos,mult,p,q,sc,t,t1,t2,t3,
          v1,v2,v3;

    info := CVEC_ComputeVectorLengthsForCalibration();
    gd := rec();
    gd.info := info;
    gd.tfCache := 0*[1..info.le];
    gd.tfNoCache := 0*[1..info.le];
    gd.tfPrimRootCache := 0*[1..info.le];
    gd.tfPrimRootNoCache := 0*[1..info.le];
    CVEC_CalibrationTableCache := [];
    CVEC_CalibrationTableNoCache := [];
    for i in [1..info.le] do
        p := info.pdq[i][1];
        d := info.pdq[i][2];
        q := info.pdq[i][3];

        # Make a short vector:
        cl := CVecClass(p,d,info.lencache[i]);

        v1 := CVEC_New(cl);
        CVEC_FastFill(v1);
        v2 := ShallowCopy(v1);
        v3 := CVEC_New(cl);

        # Calibrate mult:
        repeat    # a garbage collection can happen during calibration:
            mult := 0;
            t := Runtime();
            while Runtime()-t < 10 do
                mult := mult + 1;
                CVEC_ADD3(v3,v1,v2);
            od;

            t := Runtime();
            for j in [1..mult] do CVEC_ADD3(v3,v1,v2); od;
            t1 := Runtime() - t;
        until t1 > 0;

        if d = 1 then
            t := Runtime();
            sc := 2^Log2Int(p)-1;
            for j in [1..mult] do CVEC_MUL2(v3,v1,sc); od;
            t2 := Runtime() - t;
            gd.tfCache[i] := Maximum(1,QuoInt(t2,t1));
        else
            t := Runtime();
            for j in [1..mult] do CVEC_MUL2(v3,v1,p); od;
            t2 := Runtime() - t;
            t := Runtime();
            for j in [1..mult] do CVEC_ADDMUL(v3,v1,q-1,0,0); od;
            t3 := Runtime()-t;
            gd.tfPrimRootCache[i] := Maximum(1,QuoInt(t2,t1));
            gd.tfCache[i] := Maximum(1,QuoInt(t3,t1));
        fi;

        # Make long vectors:
        cl := CVecClass(p,d,info.lennocache[i]);

        v1 := CVEC_New(cl);
        CVEC_FastFill(v1);
        v2 := ShallowCopy(v1);
        v3 := CVEC_New(cl);

        # Calibrate mult:
        repeat    # a garbage collection can happen during calibration:
            mult := 0;
            t := Runtime();
            while Runtime()-t < 10 do
                mult := mult + 1;
                CVEC_ADD3(v3,v1,v2);
            od;

            t := Runtime();
            for j in [1..mult] do CVEC_ADD3(v3,v1,v2); od;
            t1 := Runtime() - t;
        until t1 > 0;

        if d = 1 then
            t := Runtime();
            sc := 2^Log2Int(p)-1;
            for j in [1..mult] do CVEC_MUL2(v3,v1,sc); od;
            t2 := Runtime() - t;
            gd.tfNoCache[i] := Maximum(1,QuoInt(t2,t1));
        else
            t := Runtime();
            for j in [1..mult] do CVEC_MUL2(v3,v1,p); od;
            t2 := Runtime() - t;
            t := Runtime();
            for j in [1..mult] do CVEC_ADDMUL(v3,v1,q-1,0,0); od;
            t3 := Runtime()-t;
            gd.tfPrimRootNoCache[i] := Maximum(1,QuoInt(t2,t1));
            gd.tfNoCache[i] := Maximum(1,QuoInt(t3,t1));
        fi;

        # Now we can determine the best grease levels:

        # First for in-cache:
        # Compute intersections of cost lines for levels 1..8:
        limits := List([1..7],l->q^l*(l*q-l-1));
        limits := Filtered(limits,x->x <= 1000000);
        # Now compute intersection of each cost line with level 0:
        inters := [];
        for l in [1..Length(limits)+1] do
            if not IsZero( (q-1)/q*(gd.tfCache[i]+1)-1/l ) then
                inters[l] := (q^l/l-d+(d-1)*gd.tfPrimRootCache[i])/
                                 ((q-1)/q*(gd.tfCache[i]+1)-1/l);
            else
                inters[l] := infinity;
            fi;
        od;
        minpos := 1;
        for j in [2..Length(inters)] do
            if inters[j] < inters[minpos] then minpos := j; fi;
        od;
        # Do not use grease levels < minpos:
        Add(limits,Int(inters[minpos]),1);
        for j in [2..minpos] do limits[j] := limits[1]; od;
        CVEC_CalibrationTableCache[q] := limits;

        # Now out of cache:
        # Compute intersections of cost lines for levels 1..8:
        limits := List([1..7],l->q^l*(l*q-l-1));
        limits := Filtered(limits,x->x <= 1000000);
        # Now compute intersection of each cost line with level 0:
        inters := [];
        for l in [1..Length(limits)+1] do
            if not IsZero( (q-1)/q*(gd.tfCache[i]+1)-1/l ) then
                inters[l] := (q^l/l-d+(d-1)*gd.tfPrimRootNoCache[i])/
                                 ((q-1)/q*(gd.tfNoCache[i]+1)-1/l);
            else
                inters[l] := infinity;
            fi;
        od;
        minpos := 1;
        for j in [2..Length(inters)] do
            if inters[j] < inters[minpos] then minpos := j; fi;
        od;
        # Do not use grease levels < minpos:
        Add(limits,Int(inters[minpos]),1);
        for j in [2..minpos] do limits[j] := limits[1]; od;
        CVEC_CalibrationTableNoCache[q] := limits;

        Print(i,"/",info.le,"\r");
    od;
    Print("\n");

    CVEC.greasedata := gd;
    return gd;
  end );

InstallGlobalFunction( CVEC_StoreGreaseCalibration,
  function()
    # Prints the current grease calibration data into a file from where it can
    # read after later startups.
    local MyStringList,f,i;
    MyStringList := function(l)
        return JoinStringsWithSeparator(List(l,String),",");
    end;
    f := IO_File("greasecalibration.g","w");
    IO_Write(f,"# This is grease calibration data:\n");
    IO_Write(f,"CVEC_CalibrationTableCache := [");
    for i in [1..Length(CVEC_CalibrationTableCache)] do
        if IsBound(CVEC_CalibrationTableCache[i]) then
            IO_Write(f,"\n# q=",i,":\n");
            IO_Write(f,"[",MyStringList(CVEC_CalibrationTableCache[i]),"],");
        else
            IO_Write(f,",");
        fi;
    od;
    IO_Write(f,"];\n\n");
    IO_Write(f,"CVEC_CalibrationTableNoCache := [");
    for i in [1..Length(CVEC_CalibrationTableNoCache)] do
        if IsBound(CVEC_CalibrationTableNoCache[i]) then
          IO_Write(f,"\n# q=",i,":\n");
          IO_Write(f,"[",MyStringList(CVEC_CalibrationTableNoCache[i]),"],");
        else
          IO_Write(f,",");
        fi;
    od;
    IO_Write(f,"];\n\n");
    IO_Close(f);
  end );

InstallGlobalFunction(CVEC_AddMat, function(a,b,mult)
  local i;
  for i in [1..Length(a)] do
      AddRowVector(a[i],b[i],mult);
  od;
end);

InstallGlobalFunction(CVEC_MulMat, function(a,mult)
  local i;
  for i in [1..Length(a)] do
      MultVector(a[i],mult);
  od;
end);


InstallGlobalFunction( CVEC_MultiplyWinograd, function(M,N,limit)
  # Call with matrices M and N of fitting dimensions for multiplication.
  # R must be either a matrix of the dimensions of M*N or false. limit
  # is an integer. Does M*N. A new matrix with the result is
  # returned. For matrices whose
  # product of the dimensions is smaller or equal to limit the standard
  # matrix multiplication is used, otherwise the Winograd trick is done.
  local a11,a12,a21,a22,b11,b12,b21,b22,l,l2,lhi,lhilo,llo,m,m1,m2,m3,m4,
        m5,m6,m7,mhi,mhilo,mlo,mo,n,n2,nhi,nhilo,nlo,o,r,s1,s2,s3,s4,s5,s6,
        s7,s8,t,t1,t2,ze,R;
  if NumberRows(M) * NumberColumns(M) < limit or
     NumberRows(N) * NumberColumns(N) < limit then
      return M*N; 
  fi;
  if Memory(M) >= 5000000 then 
      return CVEC_MultiplyWinogradMemory(M,N,limit);
  fi;
  #Print("Wino ",NumberRows(M)," ",NumberColumns(M)," ",NumberRows(N)," ",NumberColumns(N),"\n");
  t := Runtime();
  # From now on we have a result matrix:
  l := NumberRows(M);
  m := NumberColumns(M);   # = NumberRows(N)
  n := NumberColumns(N);
  l2 := QuoInt(l+1,2);
  m2 := QuoInt(m+1,2);
  n2 := QuoInt(n+1,2);
  llo := [1..l2];
  lhi := [l2+1..l];
  if IsOddInt(l) then
      lhilo := [1..l2-1];
  else
      lhilo := [1..l2];
  fi;
  mlo := [1..m2];
  mhi := [m2+1..m];
  if IsOddInt(m) then
      mhilo := [1..m2-1];
  else
      mhilo := [1..m2];
  fi;
  nlo := [1..n2];
  nhi := [n2+1..n];
  if IsOddInt(n) then
      nhilo := [1..n2-1];
  else
      nhilo := [1..n2];
  fi;

  o := One(BaseDomain(M));
  mo := -o;
  ze := Zero(BaseDomain(M));

  a11 := ExtractSubMatrix(M,llo,mlo);
  a21 := ZeroMutable(a11);
  CopySubMatrix(M,a21,lhi,lhilo,mlo,mlo);
  a22 := ZeroMutable(a11);
  CopySubMatrix(M,a22,lhi,lhilo,mhi,mhilo);
  s1 := a21+a22;
  s2 := s1-a11;
  s3 := a11-a21;
  a12 := ZeroMutable(a11);
  CopySubMatrix(M,a12,llo,llo,mhi,mhilo);
  s4 := a12-s2;

  b11 := ExtractSubMatrix(N,mlo,nlo);
  b12 := ZeroMutable(b11);
  CopySubMatrix(N,b12,mlo,mlo,nhi,nhilo);
  s5 := b12-b11;
  b22 := ZeroMutable(b11);
  CopySubMatrix(N,b22,mhi,mhilo,nhi,nhilo);
  s6 := b22-s5;
  s7 := b22-b12;
  b21 := ZeroMutable(b11);
  CopySubMatrix(N,b21,mhi,mhilo,nlo,nlo);
  s8 := s6-b21;
  #Print(Runtime()-t,"\n"); t := Runtime();

  # To save allocations:
  Unbind(a21);
  Unbind(b12);

  # Now the multiplications:
  m1 := CVEC_MultiplyWinograd(s2,s6,limit);
  #Print(Runtime()-t,"\n"); t := Runtime();
  m2 := CVEC_MultiplyWinograd(a11,b11,limit);
  #Print(Runtime()-t,"\n"); t := Runtime();
  m3 := CVEC_MultiplyWinograd(a12,b21,limit);
  #Print(Runtime()-t,"\n"); t := Runtime();
  m4 := CVEC_MultiplyWinograd(s3,s7,limit);
  #Print(Runtime()-t,"\n"); t := Runtime();
  m5 := CVEC_MultiplyWinograd(s1,s5,limit);
  #Print(Runtime()-t,"\n"); t := Runtime();
  m6 := CVEC_MultiplyWinograd(s4,b22,limit);
  #Print(Runtime()-t,"\n"); t := Runtime();
  m7 := CVEC_MultiplyWinograd(a22,s8,limit);
  #Print(Runtime()-t,"\n"); t := Runtime();

  Unbind(s1); Unbind(s2); Unbind(s3); Unbind(s4);
  Unbind(s5); Unbind(s6); Unbind(s7); Unbind(s8);

  t1 := m1;
  CVEC_AddMat(t1,m2,o);
  t2 := m4;
  CVEC_AddMat(t2,t1,o);
  #Print(Runtime()-t,"\n"); t := Runtime();
  
  Unbind(m1);

  # Create a new result matrix if necessary:
  R := ZeroMatrix(l,n,M);

  # Put together the result:
  CVEC_AddMat(m2,m3,o);
  CopySubMatrix(m2,R,llo,llo,nlo,nlo);
  CVEC_AddMat(t1,m5,o);
  CVEC_AddMat(t1,m6,o);
  CopySubMatrix(t1,R,llo,llo,nhilo,nhi);
  CopySubMatrix(t2-m7,R,lhilo,lhi,nlo,nlo);
  CVEC_AddMat(t2,m5,o);
  CopySubMatrix(t2,R,lhilo,lhi,nhilo,nhi);
  #Print(Runtime()-t,"\n"); t := Runtime();

  return R;
end);

InstallGlobalFunction( CVEC_MultiplyWinogradMemory, function(M,N,limit)
  # Memory efficient variant.
  # Call with matrices M and N of fitting dimensions for multiplication.
  # limit is an integer. Does M*N. A new matrix with the result is
  # returned. For matrices whose
  # product of the dimensions is smaller or equal to limit the standard
  # matrix multiplication is used, otherwise the Winograd trick is done.
  # This needs in each step at most 6 additional quarter-matrices,
  # thus the whole recursive procedure needs at most 2 additional
  # matrices to do A=M*N.

  local ClearLastColumn,ClearLastRow,R,R11,R12,R21,R22,l,l2,lhi,lhilo,llo,
        m,m2,mhi,mhilo,mlo,mo,n,n2,nhi,nhilo,nlo,o,t,x,xx,y,yy,z,ze,zz;

  #Print("WinoMem ",NumberRows(M)," ",NumberColumns(M)," ",NumberRows(N)," ",
  #      NumberColumns(N),"\n");

  ClearLastRow := function(m)
      MultVector(m[NumberRows(m)],Zero(BaseDomain(m)));
  end;
  ClearLastColumn := function(m)
      local i,r,z;
      r := NumberColumns(m);
      z := Zero(BaseDomain(m));
      for i in [1..NumberRows(m)] do
          m[i,r] := z;
      od;
  end;

  t := Runtime();
  # From now on we have a result matrix:
  l := NumberRows(M);
  m := NumberColumns(M);   # = NumberRows(N)
  n := NumberColumns(N);
  l2 := QuoInt(l+1,2);
  m2 := QuoInt(m+1,2);
  n2 := QuoInt(n+1,2);
  llo := [1..l2];
  lhi := [l2+1..l];
  if IsOddInt(l) then
      lhilo := [1..l2-1];
  else
      lhilo := [1..l2];
  fi;
  mlo := [1..m2];
  mhi := [m2+1..m];
  if IsOddInt(m) then
      mhilo := [1..m2-1];
  else
      mhilo := [1..m2];
  fi;
  nlo := [1..n2];
  nhi := [n2+1..n];
  if IsOddInt(n) then
      nhilo := [1..n2-1];
  else
      nhilo := [1..n2];
  fi;

  o := One(BaseDomain(M));
  mo := -o;
  ze := Zero(BaseDomain(M));

  #Print(Runtime()-t,"\n"); t := Runtime();
  # Start with computing R11:
  x := ZeroMatrix(l2,m2,M); 
  CopySubMatrix(M,x,llo,llo,mhi,mhilo);     # b
  y := ZeroMatrix(m2,n2,N);
  CopySubMatrix(N,y,mhi,mhilo,nlo,nlo);     # B
  R11 := CVEC_MultiplyWinograd(x,y,limit);  # bB
  #Print(Runtime()-t,"\n"); t := Runtime();
  CopySubMatrix(M,x,llo,llo,mlo,mlo);       # a
  CopySubMatrix(N,y,mlo,mlo,nlo,nlo);       # A
  z := CVEC_MultiplyWinograd(x,y,limit);    # aA
  CVEC_AddMat(R11,z,o);                     # aA+bB
  #Print(Runtime()-t,"\n"); t := Runtime();
  # Now R11 is ready, z contains aA, x contains a and y contains A

  # We now compute u, v, and w:
  
  # First w, to this end we need (a-c-d) and (A+D-C)
  xx := ZeroMatrix(l2,m2,M);
  CopySubMatrix(M,xx,lhi,lhilo,mhi,mhilo);  # d
  CVEC_AddMat(x,xx,-o);                     # a-d
  CopySubMatrix(M,xx,lhi,lhilo,mlo,mlo);    # c
  CVEC_AddMat(x,xx,-o);                     # a-c-d
  yy := ZeroMatrix(m2,n2,N);
  CopySubMatrix(N,yy,mhi,mhilo,nhi,nhilo);  # D
  CVEC_AddMat(y,yy,o);                      # A+D
  CopySubMatrix(N,yy,mlo,mlo,nhi,nhilo);    # C
  CVEC_AddMat(y,yy,-o);                     # A+D-C
  zz := CVEC_MultiplyWinograd(x,y,limit);   # (a-c-d)*(A+D-C)
  CVEC_AddMat(z,zz,-o);                     # w = aA + (c+d-a)*(A+D-C)
  # Now z holds w, x holds a-c-d, y holds A+D-C, xx holds c, yy holds C
  # zz is free but allocated l2 x n2
  Unbind(zz);
  #Print(Runtime()-t,"\n"); t := Runtime();
  # Here we have allocated 6 smaller matrices (plus R11) plus recursion!
  
  # Now compute (a+b-c-d)D since we have a-c-d at hand:
  CopySubMatrix(M,xx,llo,llo,mhi,mhilo);    # b
  if IsOddInt(m) then ClearLastColumn(xx); fi;
  CVEC_AddMat(x,xx,o);                      # a+b-c-d
  CopySubMatrix(N,yy,mhi,mhilo,nhi,nhilo);  # D
  if IsOddInt(n) then ClearLastColumn(yy); fi;
  if IsOddInt(m) then ClearLastRow(yy);fi;
  R12 := CVEC_MultiplyWinograd(x,yy,limit);  # (a+b-c-d)*D
  #Print(Runtime()-t,"\n"); t := Runtime();

  # Now compute d(B+C-A-D) since we have A+D-C at hand:
  CopySubMatrix(N,yy,mhi,mhilo,nlo,nlo);    # B
  if IsOddInt(m) then ClearLastRow(yy); fi;
  CVEC_AddMat(yy,y,-o);                     # B+C-A-D
  CopySubMatrix(M,xx,lhi,lhilo,mhi,mhilo);  # d
  if IsOddInt(l) then ClearLastRow(xx); fi;
  if IsOddInt(m) then ClearLastColumn(xx); fi;
  R21 := CVEC_MultiplyWinograd(xx,yy,limit); # d*(B+C-A-D)
  #Print(Runtime()-t,"\n"); t := Runtime();

  # Now add w to all three results:
  R22 := z;                                 # w
  CVEC_AddMat(R12,z,o);                     # w+(a+b-c-d)*D
  CVEC_AddMat(R21,z,o);                     # w+d*(B+C-A-D)
  # Now we only need u and v, but x, y, xx, yy are all allocated
  # R11 is ready, R12, R21, R22 are all allocated
  # z and zz are free since R22 inherited z and zz was freed

  # Now compute u and add it to R21 and R22:
  CopySubMatrix(M,x,lhi,lhilo,mlo,mlo);     # c
  if IsOddInt(l) then ClearLastRow(x); fi;
  CopySubMatrix(M,xx,llo,llo,mlo,mlo);      # a
  CVEC_AddMat(xx,x,-o);                     # a-c
  CopySubMatrix(N,y,mlo,mlo,nhi,nhilo);     # C
  if IsOddInt(n) then ClearLastColumn(y); fi;
  CopySubMatrix(N,yy,mhi,mhilo,nhi,nhilo);  # D
  if IsOddInt(m) then ClearLastRow(yy); fi;
  if IsOddInt(n) then ClearLastColumn(yy); fi;
  CVEC_AddMat(yy,y,-o);                     # D-C
  z := CVEC_MultiplyWinograd(xx,yy,limit);  # (c-a)*(C-D)=(a-c)*(D-C)
  CVEC_AddMat(R21,z,o);                     # R21 ready
  CVEC_AddMat(R22,z,o);                     # w+u
  Unbind(z);
  # x is still c and y is still C
  #Print(Runtime()-t,"\n"); t := Runtime();

  # Now compute v and add it to R12 and R22:
  CopySubMatrix(M,xx,lhi,lhilo,mhi,mhilo);  # d
  if IsOddInt(l) then ClearLastRow(xx); fi;
  if IsOddInt(m) then ClearLastColumn(xx); fi;
  CVEC_AddMat(x,xx,o);                      # c+d
  CopySubMatrix(N,yy,mlo,mlo,nlo,nlo);      # A
  CVEC_AddMat(y,yy,-o);                     # C-A
  z := CVEC_MultiplyWinograd(x,y,limit);    # (c+d)*(C-A) = v
  CVEC_AddMat(R12,z,o);                     # R12 ready
  CVEC_AddMat(R22,z,o);                     # R22 ready
  #Print(Runtime()-t,"\n"); t := Runtime();

  Unbind(x); Unbind(y); Unbind(z); Unbind(xx); Unbind(yy);
  # Now only R11, R12, R21, R22 are still allocated

  # Put together the result:
  R := ZeroMatrix(l,n,M);
  CopySubMatrix(R11,R,llo,llo,nlo,nlo);
  CopySubMatrix(R12,R,llo,llo,nhilo,nhi);
  CopySubMatrix(R21,R,lhilo,lhi,nlo,nlo);
  CopySubMatrix(R22,R,lhilo,lhi,nhilo,nhi);
  #Print(Runtime()-t,"\n"); t := Runtime();
  #Print("Out WinoMem\n");
  return R;
end);

InstallGlobalFunction( CVEC_ValueLaurentPoly,
  function( f, x )
    local i,j,n,o,s,val;
    f := CoefficientsOfLaurentPolynomial( f );
    i := Length( f[1] );
    if i = 0 then
        return ZeroSameMutability(x);
    elif i = 1 then
        if f[2] <> 0 then return f[1][1] * x^f[2]; fi;
        if IsOne(f[1][1]) then   # return one
            return x^0;
        fi;
        val := ZeroMutable(x);
        s := f[1][1];
        for j in [1..Length(x)] do val[j][j] := s; od;
        if not(IsMutable(x)) then
            MakeImmutable(val);
        fi;
        return val;
    fi;
    # Now we really have to do something:
    val := MutableCopyMat(x);
    n := Length(x);
    s := f[1][i];
    if not(IsOne(s)) then
        for j in [1..n] do MultVector(val[j],s); od;
    fi;
    o := One(BaseDomain(x));
    while true do
        i := i - 1;
        # Add a multiple of the identity:
        if not(IsZero(f[1][i])) then
            s := f[1][i];
            for j in [1..n] do val[j,j] := s + val[j,j]; od;
        fi;
        if i = 1 then break; fi;
        val := x*val;   # this is mutable!
    od;
    if f[2] <> 0 then val := val * x^f[2]; fi;
    if not(IsMutable(x)) then MakeImmutable(val); fi;
    return val;
  end );

InstallMethod( Value, "for a univariate Laurent polynomial, a cmat, and one",
  [ IsLaurentPolynomial, IsRingElement and IsCMatRep, 
    IsRingElement and IsCMatRep ],
  function( p, x, one ) return CVEC_ValueLaurentPoly(p,x); end );

InstallOtherMethod( Value, "for a univariate Laurent polynomial, and a cmat", 
  [ IsLaurentPolynomial, IsRingElement and IsCMatRep ],
  CVEC_ValueLaurentPoly );

InstallMethod( ScalarProductsRows, "for two matrices, and a positive integer",
  [ IsMatrixObj, IsMatrixObj, IsPosInt ],
  function( m, n, l )
    local i,sum;
    sum := ScalarProduct( m[1], n[1] );
    for i in [2..l] do
        sum := sum + ScalarProduct( m[i], n[i] );
    od;
    return sum;
  end );

InstallMethod( ScalarProductsRows, "for two cmats, and a positive integer",
  [ IsCMatRep, IsCMatRep, IsPosInt ],
  function( m, n, l )
    local i,rm,rn,sum;
    rm := m!.rows;
    rn := n!.rows;
    sum := ScalarProduct( rm[2], rn[2] );
    for i in [3..l+1] do
        sum := sum + ScalarProduct( rm[i], rn[i] );
    od;
    return sum;
  end );

InstallMethod( ScalarProductsRows, 
  "for two cmats over a small prime field and a pos. integer",
  [ IsCMatRep and IsCVecRepOverSmallField and IsCVecRepOverPrimeField,
    IsCMatRep and IsCVecRepOverSmallField and IsCVecRepOverPrimeField,
    IsPosInt ],
  CMATS_SCALAR_PRODUCTS_ROWS );

InstallMethod( EntryOfMatrixProduct, "for cmats over prime fields",
  [ IsCMatRep and IsCVecRepOverSmallField and IsCVecRepOverPrimeField,
    IsCMatRep and IsCVecRepOverSmallField and IsCVecRepOverPrimeField,
    IsPosInt, IsPosInt ],
  CMAT_ENTRY_OF_MAT_PROD );

InstallMethod( EntryOfMatrixProduct, "generic method",
  [ IsMatrixObj, IsMatrixObj, IsPosInt, IsPosInt ],
  function( m, n, i, j )
    local f,k,res;
    f := BaseDomain(m);
    res := Zero(f);
    for k in [1..NumberColumns(m)] do
        res := res + m[i,k] * n[k,j];
    od;
    return res;
  end );

##
##  This program is free software; you can redistribute it and/or modify
##  it under the terms of the GNU General Public License as published by
##  the Free Software Foundation; either version 2 of the License,
##  or (at your option) any later version.
##
##  This program is distributed in the hope that it will be useful,
##  but WITHOUT ANY WARRANTY; without even the implied warranty of
##  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
##  GNU General Public License for more details.
##
##  You should have received a copy of the GNU General Public License
##  along with this program; if not, write to the Free Software
##  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
##