Quelle  morphisms.gi   Sprache: unbekannt

#############################################################################
##
##  morphisms.gi              FinInG package
##                                                              John Bamberg
##                                                              Anton Betten
##                                                              Jan De Beule
##                                                             Philippe Cara
##                                                            Michel Lavrauw
##                                                           Max Neunhoeffer
##
##  Copyright 2018 Colorado State University
##                  Sabancı Üniversitesi
##     Università degli Studi di Padova
##     Universiteit Gent
##     University of St. Andrews
##     University of Western Australia
##                  Vrije Universiteit Brussel
##
##
##  Implementation stuff for incidence geometry morphisms
##
#############################################################################


########################################
# 20/3/14
# cmat notice: in many operations we will compute
# a matrix to initialize a form. Forms can not handle cmats, so
# Unpacking the cmats will sometimes be necessary.
########################################


############################################################
## Generic constructor operations, not intended for the user.
############################################################

## The three methods for GeometryMorphismByFunction are analogous
## with the MappingByFunction function. It behaves in exactly the same
## way except that we return an IsGeometryMorphism.

# CHECKED 27/09/11 jdb
#############################################################################
#O  GeometryMorphismByFunction( <els1>, <els2>, <fun>, <bool>, <prefun> )
##
InstallMethod( GeometryMorphismByFunction, 
 "for two times a collection of elements of any type, a function, a boolean, and a function",
 [ IsAnyElementsOfIncidenceStructure, IsAnyElementsOfIncidenceStructure,
 IsFunction, IsBool, IsFunction ],
 function( els1, els2, fun, bool, prefun )
  local morphism;
  morphism := MappingByFunction(els1, els2, fun, bool, prefun );
  SetFilterObj( morphism, IsGeometryMorphism ); 
  return morphism;
 end );  
  
# CHECKED 27/09/11 jdb
#############################################################################
#O  GeometryMorphismByFunction( <els1>, <els2>, <fun>, <prefun> )
##
InstallMethod( GeometryMorphismByFunction,
 "for two times a collection of elements of any type, and two times a function",
 [ IsAnyElementsOfIncidenceStructure, IsAnyElementsOfIncidenceStructure,
    IsFunction, IsFunction ],
 function( els1, els2, fun, inv )
  local morphism;
  morphism := MappingByFunction(els1, els2, fun, inv );
  SetFilterObj( morphism, IsGeometryMorphism ); 
  return morphism;
 end );  

# CHECKED 27/09/11 jdb
#############################################################################
#O  GeometryMorphismByFunction( <els1>, <els2>, <fun> )
##
InstallMethod( GeometryMorphismByFunction,
 "for two times a collection of elements of any type, and a function",
 [ IsAnyElementsOfIncidenceStructure, IsAnyElementsOfIncidenceStructure,
    IsFunction ],
 function( els1, els2, fun )
  local morphism;
  morphism := MappingByFunction(els1, els2, fun);
  SetFilterObj( morphism, IsGeometryMorphism ); 
  return morphism;
 end );  

#############################################################################
# Display methods
#############################################################################

InstallMethod( ViewObj, 
 "for a geometry morphism",
 [ IsGeometryMorphism ],
 function( f )
  Print("<geometry morphism from "); 
  ViewObj(Source(f));
  Print( " to " );
  ViewObj(Range(f));
  Print(">");
 end );

InstallMethod( PrintObj,
 "for a geometry morphism",
 [ IsGeometryMorphism ],
 function( f )
  Print("Geometry morphism:\n ", f, "\n");
 end );

InstallMethod( Display, 
 "for a geometry morphism",
 [ IsGeometryMorphism ],
 function( f )
  Print("Geometry morphism: ", Source(f), " -> ", Range(f), "\n");
 end );

InstallMethod( ViewObj,
 "for a geometry morphism",
 [ IsGeometryMorphism and IsMappingByFunctionWithInverseRep ],
 function( f )
  Print("<geometry morphism from "); 
  ViewObj(Source(f));
  Print( " to " );
  ViewObj(Range(f));
  Print(">");
 end );

InstallMethod( ViewObj, 
 "for a geometry morphism",
 [ IsGeometryMorphism and IsMappingByFunctionRep ],
 function( f )
  Print("<geometry morphism from "); 
  ViewObj(Source(f));
  Print( " to " );
  ViewObj(Range(f));
  Print(">");
 end );

InstallMethod( PrintObj, 
 "for a geometry morphism",
 [ IsGeometryMorphism and IsMappingByFunctionRep ],
 function( f )
  Print("Geometry morphism:\n ", f, "\n");
 end );

InstallMethod( Display,
 "for a geometry morphism",
 [ IsGeometryMorphism and IsMappingByFunctionRep ],
 function( f )
  Print("Geometry morphism: ", Source(f), " -> ", Range(f), "\n");
 end );

############################################################
## Generic methods for the Image* operations for IsGeometryMorphism. 
############################################################

# CHECKED 27/09/11 jdb
#############################################################################
#O  ImageElm( <em>, <x> )
##
InstallOtherMethod( ImageElm, 
 "for a geometry morphism and an element of an incidence structure",
 [IsGeometryMorphism, IsElementOfIncidenceStructure],
 function(em, x)
  return em!.fun(x); 
 end );

# CHECKED 27/09/11 jdb
#############################################################################
#O  \^( <x>, <em> )
##
InstallOtherMethod( \^, 
 "for an element of an incidence structure and a geometry morphism",
 [IsElementOfIncidenceStructure, IsGeometryMorphism],
 function(x, em)
  return ImageElm(em,x);
 end );


# CHECKED 27/09/11 jdb
#############################################################################
#O  ImagesSet( <em>, <x> )
##
InstallOtherMethod( ImagesSet,
 "for a geometry morphism and a collection of elements of an incidence structure",
 [IsGeometryMorphism, IsElementOfIncidenceStructureCollection],
 function(em, x)
  return List(x, t -> em!.fun(t));
 end );


# CHECKED 27/09/11 jdb
#############################################################################
#O  PreImageElm( <em>, <x> )
##
InstallOtherMethod( PreImageElm,
 "for a geometry morphism and an element of an incidence structure",
 [IsGeometryMorphism, IsElementOfIncidenceStructure],
 function(em, x)
  if IsInjective(em) then
   return em!.prefun(x); 
  else
   Error("Map is not injective");
  fi;
 end );
  
# CHECKED 27/09/11 jdb
#############################################################################
#O  PreImagesSet( <em>, <x> )
##
InstallOtherMethod( PreImagesSet,
 "for a geometry morphism and an element of an incidence structure",
 [IsGeometryMorphism, IsElementOfIncidenceStructureCollection],
 function(em, x)
  return List(x, t -> em!.prefun(t)); 
 end );

##########################################################
## User methods for the "natural geometry morphisms"
##########################################################

## The specialised operations...

# CHECKED 14/12/11 jdb + ml (and added a check that basefields of <ps1> and <ps2> are equal.)
# cmat changed 21/3/14.
#############################################################################
#O  NaturalEmbeddingBySubspace( <ps1>, <ps2>, <v> ) returns a geometry morphism
# from the projective space <ps1> into <v>, an element of the projective space
# <ps2>. <v> must be of the right type, i.e. same projective dimension as <ps1>
# and <ps1> and <ps2> must be projective spaces over the same field.
##
InstallMethod( NaturalEmbeddingBySubspace, 
 "for a projective space into another, via a specified subspace",  
 [ IsProjectiveSpace, IsProjectiveSpace, IsSubspaceOfProjectiveSpace ],
 function( geom1, geom2, v ) 
  local d1, d2, rk, f, invbasis, basis, func, pre, map, morphism, bs;
  rk := v!.type;
  basis := Unpack(v!.obj); #cmat change
  d1 := geom1!.dimension + 1;
  d2 := geom2!.dimension + 1;
  f := geom2!.basefield;
  if not v in geom2 then
   Error("Subspace is not an element of ", geom2);
  fi;
  if d2 < d1 or d1 <> rk then
   Error("Dimensions are incompatible");
  fi;
  if not f = geom1!.basefield then
   Error("Basefields must be the same");
  fi;
    ##    To find the preimage, we first find an invertible matrix C such that
    ##    [a1,..,ad]B = [a1,...,ad,0,...,0]C (where B is our d x e matrix "basis") 
    ##    is our embedding. We simply extend B to a basis for geom2.
  bs := BaseSteinitzVectors(Basis(geom2!.vectorspace), basis);
  invbasis := Inverse(Concatenation(bs!.subspace, bs!.factorspace));
  #ConvertToMatrixRep(invbasis, f); #useless now.
  func := x -> VectorSpaceToElement(geom2 , Unpack(x!.obj) * basis); #VectorSpaceToElement makes sure the 2nd arg becomes cmat.
  pre := function(y)
   local newy;
   if not y in v then
    Error("Applying preimage to an element which is not in the range");
   fi;
   newy:= Unpack(y!.obj) * invbasis; #cmat change.
   if not IsMatrix(newy) then 
    newy := newy{[1..d1]};
                #ConvertToVectorRepNC(newy, f); #useless now.
   else
    newy := newy{[1..Size(newy)]}{[1..d1]};
    #ConvertToMatrixRepNC(newy, f); #useless now.
   fi;
   return VectorSpaceToElement(geom1, newy); #VectorSpaceToElement etc.
  end;   
  morphism := GeometryMorphismByFunction(ElementsOfIncidenceStructure(geom1), 
                                           ElementsOfIncidenceStructure(geom2), 
                                           func, false, pre );
  SetIsInjective( morphism, true );  
  return morphism;
 end );

# CHECKED 27/09/11 jdb + ml
# cmat changed 21/3/14.
#############################################################################
#O  NaturalEmbeddingBySubspaceNC( <ps1>, <ps2>, <v> ) 
## This operation is just like its namesake except that it 
## has no checks
##
InstallMethod( NaturalEmbeddingBySubspaceNC, 
 "for a projective space into another, via a specified subspace",  
 [ IsProjectiveSpace, IsProjectiveSpace, IsSubspaceOfProjectiveSpace ],
 function( geom1, geom2, v ) 
  local d1, f, invbasis, basis, func, pre, map, morphism, bs;
  basis := Unpack(v!.obj); #cmat change
  d1 := geom1!.dimension + 1;
  f := geom2!.basefield;
  bs := BaseSteinitzVectors(Basis(geom2!.vectorspace), basis);
  invbasis := Inverse(Concatenation(bs!.subspace, bs!.factorspace));
  # ConvertToMatrixRep(invbasis, f); #useless now.
  func := x -> VectorSpaceToElement(geom2 , x!.obj * basis); #VectorSpaceToElement etc.
  pre := function(y)
   local newy;
   newy:= Unpack(y!.obj) * invbasis; #cmat change
   if not IsMatrix(newy) then 
    newy := newy{[1..d1]}; 
    #ConvertToVectorRepNC(newy, f); useless now.
   else
    newy := newy{[1..Size(newy)]}{[1..d1]};
    #ConvertToMatrixRepNC(newy, f); useless now.
   fi;
   return VectorSpaceToElement(geom1, newy); #cfr. supra.
  end;   
  morphism := GeometryMorphismByFunction(ElementsOfIncidenceStructure(geom1), 
                                           ElementsOfIncidenceStructure(geom2), 
                                           func, false, pre );
  SetIsInjective( morphism, true );  
  return morphism;
 end );

# CHECKED 27/09/11 jdb + ml
# cmat version 20/3/14.
#############################################################################
#O  NaturalEmbeddingBySubspace( <ps1>, <ps2>, <v> ) returns a geometry morphism
# from the polar space <ps1> into <ps2>, a polar space induced as a section of
# <ps1> and <v>, the latter a subspace of the ambient projective space of <ps1>
##
InstallMethod( NaturalEmbeddingBySubspace, 
 "for a polar space into another, via a specified subspace",  
 [ IsClassicalPolarSpace, IsClassicalPolarSpace, IsSubspaceOfProjectiveSpace ], 
 function( geom1, geom2, v ) 
 local map, d1, d2, rk, f, i, basis, invbasis, bs, c1, c2, orth, tyv, quad1, quad2,
  change, perp2, formonv, ses1, ses2, ty1, ty2, newmat, func, pre, invchange;
  rk := v!.type;
  ty1 := PolarSpaceType(geom1); 
  ty2 := PolarSpaceType(geom2); 
  f := geom1!.basefield; 
  d1 := geom1!.dimension;
  d2 := geom2!.dimension;
  tyv := TypeOfSubspace( geom2, v );
    ## Check that fields are the same and v is non-degenerate 
  if geom2!.basefield <> f then
   Error("fields of both spaces must be the same");
  fi;
    #27/9/2011. jdb was wondering on the next 7 lines. Is this not just tyv = "degenerate"?
 #if not (ty2 = "parabolic" and IsEvenInt(Size(f))) then
    #  perp2 := Polarity(geom2);
    #  if perp2(v) in geom2 then
    #     Error("subspace is degenerate"); 
    #     return;
    #  fi;
    #fi;
  if rk > d2 or d1 > d2 then
   Error("dimensions are incompatible"); 
  fi; 
  if tyv = "degenerate" then
   Error("subspace is degenerate");
  elif
   tyv <> ty1 then 
   Error("non-degenerate section is not the same type as ", geom1);
  fi;  
  orth := ["parabolic", "elliptic", "hyperbolic"];
  if ty1 = ty2 or (ty1 in orth and ty2 in orth) then 
          
          ## Let B be basis of v. Then the form BM(B^T)^sigma (n.b. sigma 
          ## is a field automorphism), where M is the form for geom2,
          ## will be equivalent to the form for geom1.

   basis := Unpack(v!.obj); #cmat unpack

          ## As usual we must consider two cases, the quadratic form case
          ## and the sesquilinear form case. We then obtain a base change matrix c1.

   if HasQuadraticForm( geom2 ) and HasQuadraticForm( geom1 ) then
    quad1 := QuadraticForm(geom1); #cmat notice: these are never cmats. So keep untouched.
    quad2 := QuadraticForm(geom2);
    newmat := basis * quad2!.matrix * TransposedMat(basis);
    formonv := QuadraticFormByMatrix(newmat, f);
    c1 := BaseChangeToCanonical( quad1 );
   else
    ses1 := SesquilinearForm(geom1);
    ses2 := SesquilinearForm(geom2);
    if ty1 = "hermitian" then 
     newmat := basis * ses2!.matrix * (TransposedMat(basis))^CompanionAutomorphism(geom1);
     formonv := HermitianFormByMatrix(newmat, f);
    else
     newmat := basis * ses2!.matrix * TransposedMat(basis);
     formonv := BilinearFormByMatrix(newmat, f);
    fi;
    c1 := BaseChangeToCanonical( ses1 );
   fi;

       ## Finding isometry from geom1 to polar space defined by formofv:

   c2 := BaseChangeToCanonical( formonv );
   change := c1^-1 * c2;        
   # ConvertToMatrixRep(change, f); #became useless.
   invchange := change^-1;
   bs := BaseSteinitzVectors(Basis(geom2!.vectorspace), basis);
   invbasis := Inverse(Concatenation(bs!.subspace, bs!.factorspace));
   # ConvertToMatrixRep(invbasis, f); #same here.
   func := x -> VectorSpaceToElement( geom2, Unpack(x!.obj) * change * basis);
   pre := function(y)
    local newy;
    if not y in v then
     Error("Applying preimage to an element which is not in the range");
    fi;
    newy:= Unpack(y!.obj) * invbasis;
    if IsMatrix(newy) then 
     newy := newy{[1..Size(newy)]}{[1..rk]};
     ConvertToMatrixRepNC(newy, f);
    else
     newy := newy{[1..rk]}; 
     ConvertToVectorRepNC(newy, f);
    fi;
    return VectorSpaceToElement(geom1, newy * invchange);
   end; 
   map := GeometryMorphismByFunction(ElementsOfIncidenceStructure(geom1), 
                                         ElementsOfIncidenceStructure(geom2), 
                                         func, false, pre );
   SetIsInjective( map, true );
  else 
   Error("Polar spaces are not compatible"); 
  fi;
  return map;
 end );
  
# CHECKED 28/09/11 jdb
# cmat comments: see NaturalEmbeddingBySubspace (above).
#############################################################################
#O  NaturalEmbeddingBySubspaceNC( <ps1>, <ps2>, <v> ) 
## This operation is just like its namesake except that it 
## has no checks
##  
InstallMethod( NaturalEmbeddingBySubspaceNC, 
 "for a geometry into another, via a specified subspace, no-check version",  
 [ IsClassicalPolarSpace, IsClassicalPolarSpace, IsSubspaceOfProjectiveSpace ],
 function( geom1, geom2, v ) 
  local map, rk, f, i, basis, invbasis, bs, c1, c2, orth, tyv, quad1, quad2,
          change, perp2, formonv, ses1, ses2, ty1, ty2, newmat, func, pre, invchange;
  rk := v!.type;
  ty1 := PolarSpaceType(geom1); 
  ty2 := PolarSpaceType(geom2); 
  f := geom1!.basefield;
  tyv := TypeOfSubspace( geom2, v );
  orth := ["parabolic", "elliptic", "hyperbolic"];
  if ty1 = ty2 or (ty1 in orth and ty2 in orth) then 
          
          ## Let B be basis of v. Then the form BM(B^T)^sigma (n.b. sigma 
          ## is a field automorphism), where M is the form for geom2,
          ## will be equivalent to the form for geom1.

  basis := Unpack(v!.obj);

          ## As usual we must consider two cases, the quadratic form case
          ## and the sesquilinear form case. We then obtain a base change matrix c1.

  if HasQuadraticForm( geom2 ) and HasQuadraticForm( geom1 ) then
   quad1 := QuadraticForm(geom1);
   quad2 := QuadraticForm(geom2);
   newmat := basis * quad2!.matrix * TransposedMat(basis);
   formonv := QuadraticFormByMatrix(newmat, f);
   c1 := BaseChangeToCanonical( quad1 );
  else
   ses1 := SesquilinearForm(geom1);
   ses2 := SesquilinearForm(geom2);
   if ty1 = "hermitian" then 
    newmat := basis * ses2!.matrix * (TransposedMat(basis))^CompanionAutomorphism(geom1);
    formonv := HermitianFormByMatrix(newmat, f);
   else
    newmat := basis * ses2!.matrix * TransposedMat(basis);
    formonv := BilinearFormByMatrix(newmat, f);
   fi;
   c1 := BaseChangeToCanonical( ses1 );
  fi;

       ## Finding isometry from geom1 to polar space defined by formofv:

  c2 := BaseChangeToCanonical( formonv );
  change := c1^-1 * c2;        
  #ConvertToMatrixRepNC(change, f);
  invchange := change^-1;
  bs := BaseSteinitzVectors(Basis(geom2!.vectorspace), basis);
  invbasis := Inverse(Concatenation(bs!.subspace, bs!.factorspace));
  #ConvertToMatrixRepNC(invbasis, f);

  func := x -> VectorSpaceToElement( geom2, Unpack(x!.obj) * change * basis);
  pre := function(y)
   local newy;
   newy:= Unpack(y!.obj) * invbasis;
    if IsMatrix(newy) then 
     newy := newy{[1..Size(newy)]}{[1..rk]};
     ConvertToMatrixRepNC(newy, f);
    else
     newy := newy{[1..rk]}; 
     ConvertToVectorRepNC(newy, f);
    fi;
    return VectorSpaceToElement(geom1, newy * invchange);
   end; 
   map := GeometryMorphismByFunction(ElementsOfIncidenceStructure(geom1), 
                                         ElementsOfIncidenceStructure(geom2), 
                                         func, false, pre );
  SetIsInjective( map, true );
  else 
   Error("Polar spaces are not compatible"); 
  fi;
  return map;
 end );
  
# CHECKED 28/09/11 jdb (and added a check that basefields of <ps1> and <ps2> are equal.)
#cmat changed 20/3/14.
#############################################################################
#O  IsomorphismPolarSpaces( <ps1>, <ps2>, <bool> ) 
# returns the coordinate transformation from <ps1> to <ps2> (which must be
# polar spaces of the same type of course). If <bool> is true, then an intertwiner 
# is computed.
##
InstallMethod( IsomorphismPolarSpaces, 
    "for two similar polar spaces and a boolean",  
    [ IsClassicalPolarSpace, IsClassicalPolarSpace, IsBool ],
 function( ps1, ps2, computeintertwiner )
  local form1, form2, f, map, c1, c2, change, invchange, hom, ty1, ty2,
          coll1, coll2, gens1, gens2, x, y, func, inv, twinerfunc, twinerprefun, r1, r2;
  ty1 := PolarSpaceType( ps1 );
  ty2 := PolarSpaceType( ps2 );
  if ty1 = "parabolic" and ty2 = "symplectic" then
            return IsomorphismPolarSpacesProjectionFromNucleus(ps1, ps2, computeintertwiner);
        fi;
        r1 := Rank(ps1);
        r2 := Rank(ps2);
        if ty1 <> ty2 or r1 <> r2 then
   Error("Polar spaces are of different type or of different rank");
  fi;
  f := ps1!.basefield;
  if f <> ps2!.basefield then
   Error( "<ps1> and <ps2> must be polar spaces over the same field");
  fi;
  if IsEvenInt(Size(f)) and ty1 in ["parabolic", "elliptic", "hyperbolic"] then
   form1 := QuadraticForm( ps1 );
   form2 := QuadraticForm( ps2 );
  else
   form1 := SesquilinearForm( ps1 );
   form2 := SesquilinearForm( ps2 );
  fi;
  c1 := BaseChangeToCanonical( form1 );
  c2 := BaseChangeToCanonical( form2 );
  change := NewMatrix(IsCMatRep, ps1!.basefield, ps1!.dimension+1,c1^-1 * c2);  #cmat change here
  #ConvertToMatrixRep(change, f);
  invchange := Inverse(change); 
  if ty1 = "hermitian" then
   change := change^CompanionAutomorphism(ps2);
   invchange := invchange^CompanionAutomorphism(ps1);
  fi; 
  func := function(x)
   return VectorSpaceToElement(ps2, ShallowCopy(x!.obj) * change);
  end;
  inv := function(x)
   return VectorSpaceToElement(ps1, ShallowCopy(x!.obj) * invchange);
  end;
  map := GeometryMorphismByFunction(ElementsOfIncidenceStructure(ps1), 
                                      ElementsOfIncidenceStructure(ps2),
                                      func, inv);               

    # map from gens2 to gens1
  twinerprefun := function( x )
   local y;
    y := MutableCopyMat( x!.mat );
                #return ProjElWithFrob( change * y * invchange^x!.frob, x!.frob, f );  
                return ProjElWithFrob( change * y * invchange^(x!.frob^-1), x!.frob, f );  

   end;  
       
    # map from gens1 to gens2
  twinerfunc := function( y )
   local x;
    x := MutableCopyMat( y!.mat );
                #return ProjElWithFrob( invchange * x * change^y!.frob, y!.frob, f );
                return ProjElWithFrob( invchange * x * change^(y!.frob^-1), y!.frob, f ); 
   end;
  if computeintertwiner then
   if (HasIsCanonicalPolarSpace( ps1 ) and IsCanonicalPolarSpace( ps1 )) or
    HasCollineationGroup( ps1 ) then
    coll1 := CollineationGroup(ps1);
    gens1 := GeneratorsOfGroup( coll1 );  
    gens2 := List(gens1, twinerfunc);  
    coll2 := GroupWithGenerators(gens2);  
    hom := GroupHomomorphismByFunction(coll1, coll2, twinerfunc, twinerprefun); 
    SetIntertwiner( map, hom );
    UseIsomorphismRelation(coll1,coll2);                
   elif (HasIsCanonicalPolarSpace( ps2 ) and IsCanonicalPolarSpace( ps2 )) or
    HasCollineationGroup( ps2 ) then                                 
    coll2 := CollineationGroup( ps2 );  
    gens2 := GeneratorsOfGroup( coll2 );          
    gens1 := List(gens2, twinerprefun ); 
    coll1 := GroupWithGenerators(gens1);      
    hom := GroupHomomorphismByFunction(coll1, coll2, twinerfunc, twinerprefun); 
    SetIntertwiner( map, hom );      
    UseIsomorphismRelation(coll1,coll2);         
   else 
    Info(InfoFinInG, 1, "No intertwiner computed. One of the polar spaces must have a collineation group computed");
   fi;
  fi;
  return map;
 end );

# CHECKED 28/09/11 jdb
#############################################################################
#O  IsomorphismPolarSpaces( <ps1>, <ps2> ) 
# returns IsomorphismPolarSpaces( <ps1>, <ps2>, true );
##
InstallMethod( IsomorphismPolarSpaces, 
 "for two similar polar spaces",  
 [ IsClassicalPolarSpace, IsClassicalPolarSpace ],
 function( ps1, ps2 )
  return IsomorphismPolarSpaces( ps1, ps2, true );
 end );

# CHECKED 28/09/11 jdb
#############################################################################
#O  IsomorphismPolarSpacesNC( <ps1>, <ps2>, <bool> ) 
# This operation is just like its namesake except that it 
## has no checks
InstallMethod( IsomorphismPolarSpacesNC, 
    "for two similar polar spaces and a boolean",  
 [ IsClassicalPolarSpace, IsClassicalPolarSpace, IsBool ],
 function( ps1, ps2, computeintertwiner )
  local form1, form2, f, map, c1, c2, change, invchange, hom, ty1, ty2, n,
          coll1, coll2, gens1, gens2, func, inv, mono, twinerprefun, twinerfunc;
  ty1 := PolarSpaceType( ps1 );
  ty2 := PolarSpaceType( ps2 );
  f := ps1!.basefield;
  if IsEvenInt(Size(f)) and ty1 in ["parabolic", "elliptic", "hyperbolic"] then
   form1 := QuadraticForm( ps1 );
   form2 := QuadraticForm( ps2 );
  else
   form1 := SesquilinearForm( ps1 );
   form2 := SesquilinearForm( ps2 );
  fi;
  c1 := BaseChangeToCanonical( form1 );
  c2 := BaseChangeToCanonical( form2 );
  change := NewMatrix(IsCMatRep, ps1!.basefield, ps1!.dimension+1,c1^-1 * c2);  #cmat change here
  #ConvertToMatrixRepNC(change, f);
  invchange := Inverse(change);     
  func := function(x)
   return VectorSpaceToElement(ps2, ShallowCopy(x!.obj) * change);
  end;
  inv := function(x)
   return VectorSpaceToElement(ps1, ShallowCopy(x!.obj) * invchange);
  end;
  map := GeometryMorphismByFunction(ElementsOfIncidenceStructure(ps1), 
                                      ElementsOfIncidenceStructure(ps2),
                                      func, inv);               
    ## Now creating intertwiner...
 # map from gens2 to gens1
  twinerprefun := function( x )
   local y;
   y := MutableCopyMat( x!.mat );
                #return ProjElWithFrob( invchange * x * change^y!.frob, y!.frob, f );
                #return ProjElWithFrob( invchange * x * change^(y!.frob^-1), y!.frob, f );
                return ProjElWithFrob( change * y * invchange^(x!.frob^-1), x!.frob, f );
  end;  
       
    # map from gens1 to gens2
  twinerfunc := function( y )
   local x;
            x := MutableCopyMat( y!.mat );
                #return ProjElWithFrob( invchange * x * change^y!.frob, y!.frob, f );
                return ProjElWithFrob( invchange * x * change^(y!.frob^-1), y!.frob, f ); 
  end;
  if computeintertwiner then
   if (HasIsCanonicalPolarSpace( ps1 ) and IsCanonicalPolarSpace( ps1 )) or
    HasCollineationGroup( ps1 ) then
    coll1 := CollineationGroup(ps1);
    gens1 := GeneratorsOfGroup( coll1 );  
    gens2 := List(gens1, twinerfunc);  
    coll2 := GroupWithGenerators(gens2);  
    hom := GroupHomomorphismByFunction(coll1, coll2, twinerfunc, twinerprefun); 
    SetIntertwiner( map, hom );
    UseIsomorphismRelation(coll1,coll2);                
   elif (HasIsCanonicalPolarSpace( ps2 ) and IsCanonicalPolarSpace( ps2 )) or
    HasCollineationGroup( ps2 ) then                                 
    coll2 := CollineationGroup( ps2 );  
    gens2 := GeneratorsOfGroup( coll2 );          
    gens1 := List(gens2, twinerprefun ); 
    coll1 := GroupWithGenerators(gens1);      
    hom := GroupHomomorphismByFunction(coll1, coll2, twinerfunc, twinerprefun); 
    SetIntertwiner( map, hom );      
    UseIsomorphismRelation(coll1,coll2);         
   else 
    Info(InfoFinInG, 1, "No intertwiner computed. One of the polar spaces must have a collineation group computed");
   fi;
  fi;
  return map;
 end );
  
# CHECKED 28/09/11 jdb
#############################################################################
#O  IsomorphismPolarSpacesNC( <ps1>, <ps2> ) 
# returns IsomorphismPolarSpacesNC( <ps1>, <ps2>, true );
##
InstallMethod( IsomorphismPolarSpacesNC, 
 "for two similar polar spaces",  
 [ IsClassicalPolarSpace, IsClassicalPolarSpace ],
 function( ps1, ps2 )  
  return IsomorphismPolarSpacesNC( ps1, ps2, true );
 end );

###########################################################
## Field reduction / subfield morphisms. Helper operations.
############################################################

# CHECKED 28/09/11 jdb
#############################################################################
#O  ShrinkMat( <B>, <mat> ) 
# returns the preimage of BlownUpMat
##
InstallMethod( ShrinkMat, 
 "for a basis and a matrix",
 [ IsBasis, IsMatrix ],
 function( B, mat )
  
  ## THERE WAS A BUG IN THIS FUNCTION: (lavrauw 15/10/2010)
  ## NOT ALL MATRICES ARE BLOWNUP MATRICES! Here is a new version.
  ## Here is the explanation that could go into the documentation:
  ## The result of BlownUpMat(B,A) is the matrix, where each entry a=A_ij is replaced by
  ## the dxd matrix M_a, representing the linear map x |-> ax with respect to the basis B of
  ## the field K seen as d-dimensional vectorspace over F. This means that if 
  ## B={b_1,b_2,...,b_d}, then the first row are the coefficients of ab_1 w.r.t. B, and so on.
  ## Concersely, if we want to implement ShrinkMat, we need to check if the input is 
  ## a blown up matrix. 
  ## This can be done as follows. Let A be a dm x dn matrix. For each dxd block, say M, 
  ## we need to check that the set {b_i^(-1)\sum_{j=1}^{d}m_ij b_j: i\in \{1,..,d\}} has size one
  ## (Since this would be the a \in K represented as a linear map by M w.r.t. the basis B)
  
  ## This function is basically the inverse map of
  ## BlownUpMat.
  
 local d,vecs,m,n,bmat,blocks,bl,submat,i,j,ij,checklist,row,newmat;
 d:=Size(B);
 vecs:=BasisVectors(B);
 m:=NrRows(mat);
 n:=NrCols(mat);
  # First we check if m and n are multiples of d
 if not (IsInt(m/d) and IsInt(n/d)) then 
  Error("The matrix does not have the right dimensions");
 fi;
 blocks:=List(Cartesian([0..m/d-1],[0..n/d-1]),ij->mat{[ij[1]*d+1 .. ij[1]*d+d]}{[ij[2]*d+1 ..ij[2]*d+d]});
 newmat:=[];
 for i in [1..m/d] do
  row:=[]; 
  for j in [1..n/d] do
   #submat:=blocks[(i-1)*d+j]; #wrong line?
      submat:=blocks[(i-1)*(m/d)+j];
   checklist:=List([1..d],i->(vecs[i]^(-1)*(Sum([1..d],j->submat[i][j]*vecs[j]))));
   if Size(AsSet(checklist))<>1 then 
    Error("The matrix is not a blown up matrix");
   fi;
   Add(row,checklist[1]);
  od;
  Add(newmat,row);
 od;
 return newmat;
  # DO WE NEED TO USE ConvertToMatrixRepNC ?
 end); 
  
  
InstallMethod( ShrinkMat, 
 "for a field, a subfield and a matrix",
 [ IsField,IsField, IsMatrix ],
 function( f1,f2,mat )
 # This gives the same result as ShrinkMat(B,mat) with B the natural basis returned by
 # Basis(AsVectorSpace(f2,f1));
 local B;
 B:=Basis(AsVectorSpace(f2,f1));
 return ShrinkMat(B,mat);
end );
 

#############################################################################
#O  ShrinkVec( <f1>, <f2>, <v> ) 
# f2 is a subfield of f1 and v is vector in a vectorspace V2 over f2
# return the vector of length d2/t, where d2=dim(V2), and t=[f1:f2]
# using the natural basis Basis(AsVectorSpace(f2,f1))
##
InstallMethod( ShrinkVec, 
 "for a field, a subfield and a vector",
 [ IsField, IsField, IsVector ],
 function( f1,f2,v )
 local t,d1,d2,basvecs;
 t:=Dimension(AsVectorSpace(f2,f1));
 d1:=Length(v)/t;
 if IsInt(d1) then
  basvecs:=BasisVectors(Basis(AsVectorSpace(f2,f1)));
  return List([1..d1],i->v{[(i-1)*t+1..i*t]}*basvecs);
 else
  Error("The length of v is not divisible by the degree of the field extension");
 fi;
end );

#############################################################################
#O ShrinkVec( <f1>, <f2>, <v>, basis )
# The same operation as above, but now with a basis as extra argument
InstallMethod( ShrinkVec, 
 "for a field, a subfield and a vector",
 [ IsField, IsField, IsVector, IsBasis ],
 function( f1,f2,v,basis )
 local t,d1,d2,basvecs;
 t:=Dimension(AsVectorSpace(f2,f1));
 d1:=Length(v)/t;
 if IsInt(d1) then
  basvecs:=BasisVectors(basis);
  return List([1..d1],i->v{[(i-1)*t+1..i*t]}*basvecs);
 else
  Error("The length of v is not divisible by the degree of the field extension");
 fi;
end );



#############################################################################
#O  BlownUpSubspaceOfProjectiveSpace( <B>, <pg1> ) 
#   blows up a projective space by field reduction.
##
InstallMethod( BlownUpProjectiveSpace,
 "for a basis and a projective space",
 [ IsBasis, IsProjectiveSpace ],
 function(basis,pg1)
 local q,t,r,pg2,mat1,mat2;
  q:=basis!.q;
  t:=basis!.d;
  if not pg1!.basefield=GF(q^t) then 
   Error("The basis and the projective space are not compatible!");
  fi;
  r:=Dimension(pg1)+1;
  pg2:=PG(r*t-1,q);
  return pg2;
 end );

#############################################################################
#O  BlownUpProjectiveSpaceBySubfield( <subfield>, <pg> ) 
#   blows up a projective space by field reduction w.r.t. the canonical basis
##
InstallMethod( BlownUpProjectiveSpaceBySubfield,
 "for a field and a projective space",
 [ IsField, IsProjectiveSpace ],
 function(subfield,pg)
 local t,r,field;
        field:=LeftActingDomain(UnderlyingVectorSpace(pg));
        if not subfield in Subfields(field) then
         Error("The first argument is not a subfield of the basefield of the projective space!");
        fi;
        t:=Dimension(AsVectorSpace(subfield,field));
        r:=Dimension(pg)+1;
        return PG(r*t-1,Size(subfield));
  end );

# CHECKED 1/10/11 jdb
#############################################################################
#O  BlownUpSubspaceOfProjectiveSpace( <B>, <subspace> ) 
#  blows up a subspace of a projective space by field reduction
##
InstallMethod( BlownUpSubspaceOfProjectiveSpace,
 "for a basis and a subspace of a projective space",
 [ IsBasis, IsSubspaceOfProjectiveSpace ],
 function(basis,subspace)
 local pg1,q,t,r,pg2,mat1,mat2;
  pg1:=AmbientGeometry(subspace);
  q:=basis!.q;
  t:=basis!.d;
  if not pg1!.basefield=GF(q^t) then 
   Error("The basis and the subspace are not compatible!");
  fi;
  r:=Dimension(pg1)+1;
  pg2:=PG(r*t-1,q);
  mat1:=subspace!.obj;
  if subspace!.type = 1 then 
   mat1 := [subspace!.obj]; 
  fi; 
  mat2:=BlownUpMat(basis,mat1);
  return VectorSpaceToElement(pg2,mat2);
 end );

#############################################################################
#O  BlownUpSubspaceOfProjectiveSpaceBySubfield( <subfield>, <subspace> )
# blows up a subspace of projective space by field reduction
# This is w.r.t. to the canonical basis of the field over the subfield.
##
InstallMethod( BlownUpSubspaceOfProjectiveSpaceBySubfield,
 "for a field and a subspace of a projective space",
 [ IsField, IsSubspaceOfProjectiveSpace],
 function(subfield,subspace)
 local pg,field,basis;
  pg:=AmbientGeometry(subspace);
  field:=LeftActingDomain(UnderlyingVectorSpace(pg));
  if not subfield in Subfields(field) then
   Error("The first argument is not a subfield of the basefield of the projective space!");
  fi;
  basis:=Basis(AsVectorSpace(subfield,field));
  return BlownUpSubspaceOfProjectiveSpace(basis,subspace);
 end );

#############################################################################
#O  IsDesarguesianSpreadElement( <B>, <subspace> ) 
# checks if a subspace is a blown up point using field reduction, w.r.t. a basis
##
InstallMethod( IsDesarguesianSpreadElement, 
 "for a basis and a subspace of a projective space",
 [ IsBasis, IsSubspaceOfProjectiveSpace ],
 function(basis,subspace)
  local flag,q,t,pg1,pg2,em,basvecs,v,v1,i,mat,rt,r;
  flag:=true;
  q:=basis!.q;
  t:=basis!.d;
  pg2:=AmbientGeometry(subspace);
  rt:=Dimension(pg2)+1;
  if not (Dimension(subspace)+1 = t and IsInt(rt/t)) then 
   flag:=false;
  else
   r:=rt/t;
   pg1:=PG(r-1,q^t);
   basvecs:=BasisVectors(basis);
   v:=Coordinates(Random(Points(subspace)));
   v1:=List([1..r],i->v{[(i-1)*t+1..i*t]}*basvecs);
   mat:=BlownUpMat(basis,[v1]);
   flag:=subspace = VectorSpaceToElement(pg2,mat);
  fi;
  return flag;
 end );
   
# To be done: check this method! (but it is currently never used).   
#############################################################################
#O  IsBlownUpSubspaceOfProjectiveSpace( <B>, <mat> ) 
# checks if a subspace is blown up using field reduction, w.r.t. a basis
##
InstallMethod( IsBlownUpSubspaceOfProjectiveSpace, 
 "for a basis and a subspace of a projective space",
 [ IsBasis, IsSubspaceOfProjectiveSpace ],
 # It's important that we include a basis in the arguments, 
 # since for every subspace, provided it has the right dimension, there exists some
 # basis with respect to which the subspace is blown up.
 function(basis,subspace)
  local flag,q,F,t,K,pg2,rt,r,pg1,basvecs,mat1,span,x,v,v1,i;
  flag:=true;
  q:=basis!.q;
  F:=GF(q);
  t:=basis!.d;
  K:=GF(q^t);
  pg2:=AmbientGeometry(subspace);
  rt:=Dimension(pg2)+1;
  if not (IsInt(rt/t) and IsInt((Dimension(subspace)+1)/t)) then 
   flag:=false;
  else
   r:=rt/t;
   pg1:=PG(r-1,q^t);
   basvecs:=BasisVectors(basis);
   mat1:=[];
   span:=[];
   repeat
    repeat x:=Random(Points(subspace)); until not x in span;
     v:=Coordinates(x);
     v1:=List([1..r],i->v{[(i-1)*t+1..i*t]}*basvecs);
     Add(mat1,v1);
     span:=VectorSpaceToElement(pg2,BlownUpMat(basis,mat1));
    until Dimension(span)=Dimension(subspace);
    if not span = subspace then 
     flag:= false;
    fi;
  fi;
  return flag;
 end );   
  
#############################################################################
#  PROJECTIVE SPACES
#############################################################################

#
#  A recent reference on field reduction of projective spaces is
#  "Field reduction and linear sets in finite geomety" by Lavrauw
#   and Van de Voorde (preprint 2013)

#############################################################################
#O  NaturalEmbeddingByFieldReduction( <geom1>, <field>, <basis> ) 
# <geom2> is a projective space over a field K, <geom1> is a projective space
# over a field extension L, and considering L as a vector space over K, yields 
# that <geom1> and <geom2> have the same ambient vectorspace over K, then this
# operation returns the natural embedding, i.e. with relation to the given basis
# of L over K.
# Important note: the intertwiner has as source the *Projectivity group* only, not
# the full collineation group!
##
InstallMethod( NaturalEmbeddingByFieldReduction, 
 "for two projective spaces and a basis",
     [ IsProjectiveSpace, IsField, IsBasis ],
 function( pg1, f2, basis )
  
  ## This morphism contains a func and prefunc with built-in check.
 
  local map, f1, d1, d2, t, fun, prefun, g1, gens, newgens, g2, twiner, hom, hominv, q1, q2, pg2;
  f1 := pg1!.basefield; 
  q1 := Size(f1);
  q2 := Size(f2);
  t := LogInt(q1,q2);
  if not q2^t = q1 then
   Error( "<f2> is not a subfield of the base field of <pg1> ");
  fi;
  d1 := pg1!.dimension + 1;
  d2 := d1*t;
  pg2 := ProjectiveSpace(d2-1,q2);
  if not (IsBasis(basis) and f1=GF((basis!.q)^basis!.d) and f2=GF(basis!.q) and d1*(basis!.d)=d2) then
   Error("The basis is not a basis or is not compatible with the basefields of the geometries");
  fi;

  fun := function( subspace ) # This map blows up a subspace of geom1 to a subspace of geom2
            local mat1,mat2;
            #return BlownUpSubspaceOfProjectiveSpace(basis,x);
            mat1 := subspace!.obj;
            if subspace!.type = 1 then
                mat1 := [subspace!.obj];
            fi;
            mat2 := BlownUpMat(basis,mat1);
            return VectorSpaceToElement(pg2,mat2);
  end;

        prefun := function( subspace ) # This map is the inverse of func and returns an error, or a subspace of geom1
            local flag,basvecs,mat1,span,x,v,v1,i,vecs;
            flag:=true;
            if not subspace in pg2 then
                Error("The input is not in the range of the field reduction map!");
            fi;
            if not IsInt((Dimension(subspace)+1)/t) then
                flag := false;
            else
                basvecs:=BasisVectors(basis);
                vecs := Unpack(UnderlyingObject(subspace));
                mat1 := List(vecs,x->List([1..d1],i->x{[(i-1)*t+1..i*t]}*basvecs));
                span := VectorSpaceToElement(pg2,BlownUpMat(basis,mat1));
                if not span = subspace then
                    flag := false;
                fi;
            fi;
            if flag= false then
                Error("The input is not in the range of the field reduction map!");
            fi;
            return VectorSpaceToElement(pg1,mat1);
        end;

     #prefun := function( x ) #I want to find out why this is not correct...
  # return VectorSpaceToElement(pg1,ShrinkMat(basis,x!.obj));
  #end;
       
  map := GeometryMorphismByFunction(ElementsOfIncidenceStructure(pg1),
                                         ElementsOfIncidenceStructure(pg2),
                                            fun, false, prefun);
  
  SetIsInjective( map, true );
        
        ## Now creating intertwiner
        ## 12/6/14: added check to see whether m is really a projectivity
  
        hom := function( m )
   local image;
            if not IsOne(m!.frob) then
                return fail;
                Info(InfoFinInG, 1, "<el> is not a projectivity");
            fi;
   image := BlownUpMat(basis, m!.mat); 
   #ConvertToMatrixRepNC( image, f2 );
   return CollineationOfProjectiveSpace(image, f2);
  end;

  hominv := function( m )
   local preimage;
            if not IsOne(m!.frob) then
                return fail;
                Info(InfoFinInG, 1, "<el> is not a projectivity");
            fi;
   preimage := ShrinkMat(basis, Unpack(m!.mat));
   #ConvertToMatrixRepNC( preimage, f1 );       
   return CollineationOfProjectiveSpace(preimage, f1);
  end;

  g1 := ProjectivityGroup( pg1 );
  gens := GeneratorsOfGroup( g1 );
  newgens := List(gens, hom);
  g2 := Group( newgens );
  SetSize(g2, Size(g1));
  twiner := GroupHomomorphismByFunction(g1, g2, hom, hominv);
  SetIntertwiner( map, twiner);
  return map;
 end );

#############################################################################
#O  NaturalEmbeddingByFieldReduction( <geom1>, <field> ) 
# <field> is a field K, <geom1> is a projective space
# over a field extension L, and considering L as a vector space over K, yields 
# that a projective space geom2 that has the same ambient vectorspace over K. This
# operation returns the natural embedding, i.e. with relation to the standard basis of
# L over K.
##
InstallMethod( NaturalEmbeddingByFieldReduction, 
 "for a projective space and a field",
 [ IsProjectiveSpace, IsField ],
 function( pg1, f2 )
  local basis;
  basis:=Basis(AsVectorSpace(f2,pg1!.basefield));
  return NaturalEmbeddingByFieldReduction(pg1,f2,basis);
 end );

#############################################################################
#O  NaturalEmbeddingByFieldReduction( <geom1>, <geom2>, <basis> ) 
# <geom2> is a projective space over a field K, <geom1> is a projective space
# over a field extension L, and considering L as a vector space over K, yields 
# that <geom1> and <geom2> have the same ambient vectorspace over K, then this
# operation returns the natural embedding, i.e. with relation to the given basis
# of L over K.
##
InstallMethod( NaturalEmbeddingByFieldReduction, 
 "for two projective spaces and a basis",
 [ IsProjectiveSpace, IsProjectiveSpace, IsBasis ],
 function( pg1, pg2, basis )
  return NaturalEmbeddingByFieldReduction(pg1,pg2!.basefield,basis);
 end);

# CHECKED 28/09/11 jdb
#############################################################################
#O  NaturalEmbeddingByFieldReduction( <geom1>, <geom2> ) 
# <geom2> is a projective space over a field K, <geom1> is a projective space
# over a field extension L, and considering L as a vector space over K, yields 
# that <geom1> and <geom2> have the same ambient vectorspace over K, then this
# operation returns the natural embedding, i.e. with relation to the standard basis of
# L over K.
##
InstallMethod( NaturalEmbeddingByFieldReduction, 
 "for two projective spaces",
 [ IsProjectiveSpace, IsProjectiveSpace ],
 function( pg1, pg2 )
  local basis;
  basis:=Basis(AsVectorSpace(pg2!.basefield,pg1!.basefield));
  return NaturalEmbeddingByFieldReduction(pg1,pg2!.basefield,basis);
 end );










############################################################################
#  POLAR SPACES
##############################################################################


#
#  Two references on field reduction of polar spaces:
#  [1] from group theoretic point of view: "Polar spaces and embeddings of classical groups" by Nick Gill
#  (New Zealand J. Math)
#  [2] from finite geometry point of view: "Field reduction and linear sets in finite geomety" by Lavrauw
#   and Van de Voorde (preprint 2013)
#
#  We will use [2] as it is more convenient to us.


####
#
# First the field reduction of bilinear forms, quadratic forms and hermitian forms
#
#####


#############################################################################
#O  BilinearFormFieldReduction( <bil1>, <f2>, <alpha>, <basis> ) 
# <bil1> is a bilinear form over <f1>, <f2> is a subfield of <f1>. 
# <alpha> determines the linear map from f1 to <f2>, f1 the field extension of <f2>
# with basis <basis>. This operation then
# returns the bilinear form L(bil1(.,.)), L(x) = T(alpha*x), T the trace from <f1> to <f2>.
##
# REMARK (ml 31/10/13) The fourth argument is a basis of L over K. This just gives
# extra functionality to the user, and it is used as third argument in the 
# operation ShrinkVec. 
# There is also a version of this field reduction operation, which does not use
# a basis as fourth argument. In this case the standard GAP basis is used.

InstallMethod( BilinearFormFieldReduction,
 "for a bilinear form and a field",
 [ IsBilinearForm, IsField, IsFFE, IsBasis ],
 function(bil1,f2,alpha,basis)
 # f2 is a subfield of the basefield of the bilinear form bil1
 local mat1,f1,d1,t,d2,V2,V1,phi,b2,b2vecs,mat2,i,row,j,bil2;
 mat1:=bil1!.matrix;
 f1:=bil1!.basefield;
 d1:=NrRows(mat1);
 t:=Dimension(AsVectorSpace(f2,f1));
 d2:=d1*t;
 V2:=f2^d2;
 V1:=f1^d1;
 b2:=Basis(V2);
 b2vecs:=BasisVectors(b2);
 mat2:=[];
 for i in [1..d2] do 
  row:=[];
  for j in [1..d2] do
   Add(row,Trace(f1,f2,alpha*(ShrinkVec(f1,f2,b2vecs[i],basis)*mat1*ShrinkVec(f1,f2,b2vecs[j],basis))));
  od;
  Add(mat2,row);
 od;
 bil2:=BilinearFormByMatrix(mat2,f2);
 return bil2;
end );

#############################################################################
#O  BilinearFormFieldReduction( <bil1>, <f2>, <alpha> ) 
# The same as above, but without specified basis. In this case the canonical basis is used.
##
InstallMethod( BilinearFormFieldReduction,
 "for a bilinear form and a field",
 [ IsBilinearForm, IsField, IsFFE ],
 function(bil1,f2,alpha)
 # f2 is a subfield of the basefield of the bilinear form bil1
 return BilinearFormFieldReduction(bil1,f2,alpha,Basis(AsVectorSpace(f2,bil1!.basefield)));
end );

#############################################################################
#O  QuadraticFormFieldReduction( <qf1>, <f2>, <alpha>, <basis> ) 
# <bil1> is a bilinear form over <f1>, <f2> is a subfield of <f1>. This operation
# returns the bilinear form T(alpha*qf1(.,.)), T the trace from <f1> to <f2>.
##
# REMARK (ml 31/10/13) The fourth argument is a basis of L over K. This just gives
# extra functionality to the user, and it is used as third argument in the 
# operation ShrinkVec. 
# There is also a version of this field reduction operation, which does not use
# a basis as fourth argument. In this case the standard GAP basis is used.

InstallMethod( QuadraticFormFieldReduction,
 "for a quadratic form and a field",
 [ IsQuadraticForm, IsField, IsFFE, IsBasis ],
 function(qf1,f2,alpha,basis)
 # f2 is a subfield of the basefield for the quadratic form q1
 local f1,basvecs,d1,t,d2,V2,V1,phi,b2,b2vecs,mat2,i,bil1,j,qf2;
 f1:=qf1!.basefield;
 basvecs:=BasisVectors(basis);
 d1:=NrRows(qf1!.matrix);
 t:=Dimension(AsVectorSpace(f2,f1));
 d2:=d1*t;
 V2:=f2^d2;
 V1:=f1^d1;
 b2:=Basis(V2);
 b2vecs:=BasisVectors(b2);
 mat2:=IdentityMat(d2,f2);
 for i in [1..d2] do
  mat2[i,i]:=Trace(f1,f2,alpha*((ShrinkVec(f1,f2,b2vecs[i],basis))^qf1));
 od;
 for i in [1..d2-1] do
  for j in [i+1..d2] do
   mat2[i,j]:=Trace(f1,f2,alpha*((ShrinkVec(f1,f2,b2vecs[i]+b2vecs[j],basis))^qf1
       -(ShrinkVec(f1,f2,b2vecs[i],basis))^qf1-(ShrinkVec(f1,f2,b2vecs[j],basis))^qf1));
   #mat2[j,i]:=mat2[i,j]; THESE entries need to be zero
  od;
 od;
 qf2:=QuadraticFormByMatrix(mat2,f2);
 return qf2;
end );

#############################################################################
# #O  QuadraticFormFieldReduction( <qf1>, <f2>, <alpha> )
# The same as above, but without specified basis. In this case the canonical basis is used.
##
InstallMethod( QuadraticFormFieldReduction,
 "for a quadratic form and a field",
 [ IsQuadraticForm, IsField, IsFFE ],
 function(qf1,f2,alpha)
 local basis,qf2;
 # f2 is a subfield of the basefield for the quadratic form q1
 basis:=Basis(AsVectorSpace(f2,qf1!.basefield));
 qf2:=QuadraticFormFieldReduction(qf1,f2,alpha,basis);
 return qf2;
end );


#############################################################################
#O  HermitianFormFieldReduction( <hf1>, <f2>, <alpha>, <basis> ) 
# <hf1> is a hermitian form over <f1>, <f2> is a subfield of <f1>. This operation
# returns the form T(alpha*bil1(.,.)). Depending on the degree of the field extension, this
# yields either a bilinear form or a hermitian form.
##
# REMARK (ml 31/10/13) The fourth argument is a basis of L over K. This just gives
# extra functionality to the user, and it is used as third argument in the 
# operation ShrinkVec. 
# There is also a version of this field reduction operation, which does not use
# a basis as fourth argument. In this case the standard GAP basis is used.

InstallMethod( HermitianFormFieldReduction,
 "for a hermitian form and a field",
 [ IsHermitianForm, IsField, IsFFE, IsBasis ],
 function(hf1,f2,alpha,basis)
 # f2 is a subfield of the basefield for the hermitian form hf1
 # here the basefield is always a square
 local f1,d1,t,d2,V2,V1,phi,b2,b2vecs,mat2,i,row,j,hf2;
 f1:=hf1!.basefield;
 d1:=NrRows(hf1!.matrix);
 t:=Dimension(AsVectorSpace(f2,f1));
 d2:=d1*t;
 V2:=f2^d2;
 V1:=f1^d1;
 b2:=Basis(V2);
 b2vecs:=BasisVectors(b2);
 mat2:=[];
 for i in [1..d2] do 
  row:=[];
  for j in [1..d2] do
   Add(row,Trace(f1,f2,alpha*([ShrinkVec(f1,f2,b2vecs[i],basis),ShrinkVec(f1,f2,b2vecs[j],basis)]^hf1)));
  od;
  Add(mat2,row);
 od;
 # checking parity of t is indeed sufficient to decide, see table in manual.
 if IsOddInt(t) then 
  hf2:=HermitianFormByMatrix(mat2,f2);
 else
  hf2:=BilinearFormByMatrix(mat2,f2);
 fi;
 return hf2;
end );


#############################################################################
# #O  HermitianFormFieldReduction( <hf11>, <f2>, <alpha> ) 
# The same as above, but without specified basis. In this case the canonical basis is used.
##

InstallMethod( HermitianFormFieldReduction,
 "for a hermitian form and a field",
 [ IsHermitianForm, IsField, IsFFE ],
 function(hf1,f2,alpha)
 # f2 is a subfield of the basefield for the hermitian form hf1
 # here the basefield is always a square
 local basis;
 
 basis:=Basis(AsVectorSpace(f2,hf1!.basefield));
 return HermitianFormFieldReduction(hf1,f2,alpha,basis);
end );



##########################################################
#
# Next the embeddings for polar spaces based on the field reduction of the forms above
# 
#############################################################################

# master version for the user: handle all parameters.
#############################################################################
#O  NaturalEmbeddingByFieldReduction( <ps1>, <f2>, <alpha>, <basis>, <bool> ) 
# returns a geometry morphism, described below. If <bool> is true, then an intertwiner
# is computed.
##
InstallMethod (NaturalEmbeddingByFieldReduction,
 "for a polar space, a field, a basis, a finite field element, and a boolean",
 [IsClassicalPolarSpace, IsField, IsFFE, IsBasis, IsBool],
 function(ps1,f2,alpha,basis,computeintertwiner)
 # f2 must be a subfield of the basefield of the classical polar space
 local type1,qf1,qf2,ps2,bil1,bil2,hf1,hf2,fun,prefun,f1,
     map,hom,hominv,g1,gens,newgens,g2,twiner, t,q1,q2,d1,d2;
 type1 := PolarSpaceType(ps1);
 f1:=ps1!.basefield;
 q1 := Size(f1);
 q2 := Size(f2);
 t := LogInt(q1,q2);
 d1 := ps1!.dimension + 1;
 d2 := d1*t;

 # 1. the polar space is of orthogonal type
 if type1 in ["hyperbolic","elliptic","parabolic"] then
  qf1:=QuadraticForm(ps1);
  qf2:=QuadraticFormFieldReduction(qf1,f2,alpha,basis);
  
  if IsSingularForm(qf2) then 
   Error("The field reduction does not yield a natural embedding");
  else ps2:=PolarSpace(qf2);
  fi;
 # 2. the polar space is of symplectic type
 elif type1 in ["symplectic"] then
  bil1:=SesquilinearForm(ps1);
  bil2:=BilinearFormFieldReduction(bil1,f2,alpha,basis);
  ps2:=PolarSpace(bil2);
 # 3. the polar space is of hermitian type 
 elif type1 in ["hermitian"] then
  hf1:=SesquilinearForm(ps1);
  hf2:=HermitianFormFieldReduction(hf1,f2,alpha,basis);
  ps2:=PolarSpace(hf2);
 fi;
 
 #em:=NaturalEmbeddingByFieldReduction(AmbientSpace(ps1),AmbientSpace(ps2));
 # this is the field reduction for projective spaces PG(r-1,q^t) -> PG(rt-1,q)
 
 #fun:=function(x)
 # local projfun;
 # projfun:=em!.fun;
 # return VectorSpaceToElement(ps2,projfun(x)!.obj);
 #end;
 
    fun := function( subspace )
        local mat1,mat2;
        mat1 := subspace!.obj;
        if subspace!.type = 1 then
            mat1 := [subspace!.obj];
        fi;
        mat2 := BlownUpMat(basis,mat1);
        return VectorSpaceToElement(ps2,mat2);
    end;


 #prefun:=function(x)
 # local projprefun;
 # projprefun:=em!.prefun;
 # return VectorSpaceToElement(ps1,projprefun(x)!.obj);
 #end;

    # new version (2/7/14, much faster since we avoid the repeat loops and Random calls.
    prefun := function( subspace ) # This map is the inverse of func and returns an error, or a subspace of geom1
        local flag,basvecs,mat1,span,x,v,v1,i,vecs;
        flag:=true;
        if not subspace in ps2 then
            Error("The input is not in the range of the field reduction map!");
        fi;
        if not IsInt((Dimension(subspace)+1)/t) then
            flag:=false;
        else
            basvecs:=BasisVectors(basis);
            mat1:=[];
            span:=[];
            vecs := Unpack(UnderlyingObject(subspace));
            mat1 := List(vecs,x->List([1..d1],i->x{[(i-1)*t+1..i*t]}*basvecs));
            span := VectorSpaceToElement(AmbientSpace(ps2),BlownUpMat(basis,mat1)); 
                #the ambientspace makes sure that there is no error message here yet.
            if not span=subspace then
                flag := false;
            fi;
        fi;
        if flag= false then
            Error("The input is not in the range of the field reduction map!");
        fi;
        return VectorSpaceToElement(ps1,mat1);
    end;

 map := GeometryMorphismByFunction(ElementsOfIncidenceStructure(ps1), ElementsOfIncidenceStructure(ps2),
                                           fun, false, prefun);
  
 SetIsInjective( map, true );
 
 ## Now creating intertwiner
 
 if computeintertwiner then
        if (HasIsCanonicalPolarSpace( ps1 ) and IsCanonicalPolarSpace( ps1 )) or
            HasCollineationGroup( ps1 ) then
            
            hom := function( m )
                local image;
                image := BlownUpMat(basis, m!.mat);
                ConvertToMatrixRepNC( image, f2 );
                return CollineationOfProjectiveSpace(image, f2);
            end;
            hominv := function( m )
                local preimage;
                preimage := ShrinkMat(basis, Unpack(m!.mat));
                ConvertToMatrixRepNC( preimage, f1 );
                return CollineationOfProjectiveSpace(preimage, f1);
            end;
            g1 := IsometryGroup( ps1 );
            gens := GeneratorsOfGroup( g1 );
            newgens := List(gens, hom);
            g2 := Group( newgens );
            SetSize(g2, Size(g1));
            twiner := GroupHomomorphismByFunction(g1, g2, hom, hominv);
            SetIntertwiner( map, twiner);
        fi;
 fi;

    return map;
 end );


#first version: user wants a particular alpha, and basis, agrees with bool=true.
#############################################################################
#O  NaturalEmbeddingByFieldReduction( <geom1>, <f2>, <alpha>, <basis> )
# returns NaturalEmbeddingByFieldReduction( <geom1>, <f2>, <alpha>, <basis>, true )
#
InstallMethod (NaturalEmbeddingByFieldReduction,
 "for a polar space and a field, a finite field element, and a basis",
 [IsClassicalPolarSpace, IsField, IsFFE, IsBasis],
 function(ps1,f2,alpha,basis)
  return NaturalEmbeddingByFieldReduction(ps1,f2,alpha,basis,true);
 end );

#second particular version: user wants a particular alpha, and bool, agrees with basis.
#############################################################################
#O  NaturalEmbeddingByFieldReduction( <geom1>, <f2>, <alpha>, <bool> )
# returns NaturalEmbeddingByFieldReduction( <geom1>, <f2>, <alpha>, <basis>, true )
# where <basis> is the canonical basis of BaseField(geom1) over <f2>.
#
InstallMethod (NaturalEmbeddingByFieldReduction,
 "for a polar space and a field, a finite field element, and a boolean",
 [IsClassicalPolarSpace, IsField, IsFFE, IsBool],
 function(ps1,f2,alpha,bool)
  return NaturalEmbeddingByFieldReduction(ps1,f2,alpha,Basis(AsVectorSpace(f2,BaseField(ps1))),bool);
 end );

#third particular version: user wants a particular alpha, agrees with basis and bool.
#############################################################################
#O  NaturalEmbeddingByFieldReduction( <geom1>, <f2>, <alpha> )
# returns NaturalEmbeddingByFieldReduction( <geom1>, <f2>, <alpha>, <basis>, true )
# where <basis> is the canonical basis of BaseField(geom1) over <f2>.
#
InstallMethod (NaturalEmbeddingByFieldReduction,
 "for a polar space, a field, and a finite field element",
 [IsClassicalPolarSpace, IsField, IsFFE],
 function(ps1,f2,alpha)
  return NaturalEmbeddingByFieldReduction(ps1,f2,alpha,Basis(AsVectorSpace(f2,BaseField(ps1))),true);
 end );

# fourth particular version: user agrees but wants to control intertwiner
#############################################################################
#O  NaturalEmbeddingByFieldReduction( <geom1>, <f2>, <bool> ) 
# returns NaturalEmbeddingByFieldReduction( <geom1>, <f2>, <one>, <basis>, <bool> )
# where <basis> is the canonical basis of BaseField(geom1) over <f2>, and
# <one> is One(BaseField(geom1)). This might be usefull if you agree with the defaults, but doesn't 
# want the intertwiner.
#
InstallMethod (NaturalEmbeddingByFieldReduction,
 "for a polar space, a field, and a boolean",
 [IsClassicalPolarSpace, IsField, IsBool],
 function(ps1,f2,bool)
  return NaturalEmbeddingByFieldReduction(ps1,f2,One(f2),Basis(AsVectorSpace(f2,BaseField(ps1))),bool);
 end );

# fifth particular version: user agrees with everything
#############################################################################
#O  NaturalEmbeddingByFieldReduction( <geom1>, <f2> ) 
# returns NaturalEmbeddingByFieldReduction( <geom1>, <f2>, <alpha>, <basis>, true )
# where <basis> is the canonical basis of BaseField(geom1) over <f2>, and
# <alpha> is One(f2).
InstallMethod (NaturalEmbeddingByFieldReduction,
 "for a polar space and a field",
 [IsClassicalPolarSpace, IsField],
 function(ps1,f2)
  return NaturalEmbeddingByFieldReduction(ps1,f2,One(f2),Basis(AsVectorSpace(f2,BaseField(ps1))),true);
 end );

# CHECKED 28/09/11 jdb
#############################################################################
#O  NaturalEmbeddingByFieldReduction( <geom1>, <geom2> ) 
# returns NaturalEmbeddingByFieldReduction(<geom1>, <geom2> )
##
InstallMethod( NaturalEmbeddingByFieldReduction, 
 "for two polar spaces",
 [ IsClassicalPolarSpace, IsClassicalPolarSpace ],
 function( geom1, geom2 )
  return NaturalEmbeddingByFieldReduction( geom1, geom2, true );
 end );

# ml 31/10/2013
# changed with new prefun 2/7/14 jdb
#####################################################################################
# O  NaturalEmbeddingByFieldReduction
# This NaturalEmbeddingByFieldReduction takes two polar spaces and returns an embedding
# obtained by field reduction if such an embedding exists, otherwise it returns an error.
##
InstallMethod (NaturalEmbeddingByFieldReduction,
 "for two classical polar spaces",
 [IsClassicalPolarSpace, IsClassicalPolarSpace, IsBool],
 function(ps1, ps2, computeintertwiner)

 #####################################################################
 #  List of possible embeddings of ps1 in PG(r-1, q^t) -> ps2 in PG(rt-1, q)
 #
 #  Hermitian
 #   H -> H (t odd)
 #   H -> W (t even)
 #   H -> Q+ (t even, r even)
 #   H -> Q- (t even, r odd)
 #
 #  Symplectic
 #   W -> W (always)
 #
 #  Orthogonal 
 # r even   
 # Q+ -> Q+ (always)
 # Q- -> Q- (always)
 # r odd
 # Q -> Q (t odd, q odd)
 #   Q -> Q+ (t even, q odd)
 #   Q -> Q- (t even, q odd)
 #####################################################################

    local t, r, type1, type2, form1, f1, f2, newps, sigma, map,
    gamma, bil1, q, n, alpha, basis, is_square, iso, form2, fun, prefun,
          c1, c2, cps2, cps2inv, d1, hom, hominv, g1, g2, gens, newgens, twiner;
     
 type1 := PolarSpaceType(ps1);
 type2 := PolarSpaceType(ps2);
 if type1 in ["hyperbolic", "parabolic", "elliptic"] then 
    form1 := QuadraticForm(ps1);
 else
    form1 := SesquilinearForm(ps1);
 fi;
    d1 := ps1!.dimension+1;
 f1 := ps1!.basefield;
 f2 := ps2!.basefield;
 q := Size(f2);
 
 if not IsSubset(f1,f2) then
    Error("Fields are incompatible");
 fi;
 t := Log(Size(f1),q);
 r := ProjectiveDimension(ps1)+1;
 basis := CanonicalBasis(AsVectorSpace(f2,f1));
    is_square := function(x, f) return IsZero(x) or IsEvenInt(LogFFE(x,PrimitiveRoot(f))); end;
      
 if IsSesquilinearForm(form1) then
        if type1 = "hermitian" and type2 = "hermitian" and IsOddInt(t) then
   Info(InfoFinInG, 1, "These polar spaces are suitable for field reduction");
    # need alpha to be fixed by sigma:
   alpha := One(f1);
   #map := NaturalEmbeddingByFieldReduction( ps1, f2, alpha, basis );
   form2 := HermitianFormFieldReduction(form1,f2,alpha,basis);
     elif type1 = "hermitian" and type2 = "symplectic" and IsEvenInt(t) then
   Info(InfoFinInG, 1, "These polar spaces are suitable for field reduction");
   if IsEvenInt(q) then
    # need alpha to be fixed by sigma:
    alpha := One(f1);    
   else
    # need sigma(alpha)=-alpha
    sigma := Sqrt(Size(f1));
    alpha := First(f1, x->not IsZero(x) and x^sigma = -x);    
   fi; 
   #map := NaturalEmbeddingByFieldReduction( ps1, f2, alpha, basis );
            form2 := HermitianFormFieldReduction(form1,f2,alpha,basis);
  
     elif type1 = "hermitian" and type2 = "hyperbolic" and IsEvenInt(t) and IsEvenInt(r) 
    and IsOddInt(q) then
   Info(InfoFinInG, 1, "These polar spaces are suitable for field reduction");
   alpha := One(f1); 
   #map := NaturalEmbeddingByFieldReduction( ps1, f2, alpha, basis );
            form2 := HermitianFormFieldReduction(form1,f2,alpha,basis);
      
     elif type1 = "hermitian" and type2 = "elliptic" and IsEvenInt(t) and IsOddInt(r) 
    and IsOddInt(q) then
   Info(InfoFinInG, 1, "These polar spaces are suitable for field reduction");
   alpha := One(f1); 
   #map := NaturalEmbeddingByFieldReduction( ps1, f2, alpha, basis );
            form2 := HermitianFormFieldReduction(form1,f2,alpha,basis);
   
     elif type1 = "symplectic" and type2 = "symplectic" then
   Info(InfoFinInG, 1, "These polar spaces are suitable for field reduction");
   alpha := One(f1);    
   #map := NaturalEmbeddingByFieldReduction( ps1, f2, alpha, basis );
            form2 := BilinearFormFieldReduction(form1,f2,alpha,basis);
        else
    Error("These polar spaces are not suitable for field reduction.");
  fi;
 else  # form1 is a quadratic form
  if type1 = "hyperbolic" and type2 = "hyperbolic" then
   Info(InfoFinInG, 1, "These polar spaces are suitable for field reduction");
   alpha := One(f1);
   #map := NaturalEmbeddingByFieldReduction( ps1, f2, alpha, basis );
            form2 := QuadraticFormFieldReduction(form1,f2,alpha,basis);
  elif type1 = "elliptic" and type2 = "elliptic" then
   Info(InfoFinInG, 1, "These polar spaces are suitable for field reduction");
   alpha := One(f1);
   #map := NaturalEmbeddingByFieldReduction( ps1, f2, alpha, basis );
            form2 := QuadraticFormFieldReduction(form1,f2,alpha,basis);
  elif type1 = "parabolic" and type2 = "parabolic" and IsOddInt(t) and IsOddInt(q) then
   Info(InfoFinInG, 1, "These polar spaces are suitable for field reduction");
   alpha := One(f1);
   #map := NaturalEmbeddingByFieldReduction( ps1, f2, alpha, basis );
            form2 := QuadraticFormFieldReduction(form1,f2,alpha,basis);

  # If ps1 is parabolic and ps2 is hyperbolic or elliptic, 
  # then the choice of alpha that we will use for the embedding
  # depends on the isometry type of ps1. The isometry type is determined by
  # gamma being a square or not, where gamma is
  # (-1)^(r-1)/2 times the determinant of the bilinear form, divided by two.
  # The conditions we use here are taken from [Lavrauw - Van de Voorde: "Field
  # reduction and linear sets in finite geometry", preprint], they are
  # somewhat easier to work with than the conditions in Gill's paper.

  elif type1 = "parabolic" and type2 = "hyperbolic" and IsEvenInt(t) and IsOddInt(q) then
   Info(InfoFinInG, 1, "These polar spaces are suitable for field reduction: Parabolic -> Hyperbolic");
   bil1:=AssociatedBilinearForm(form1); # important to use this instead of 
   # BilinearFormByQuadraticForm (see documentation of the Forms package).
   gamma:=((-1)^((r-1)/2)) * Determinant(bil1!.matrix)/2;
   if q^(t/2) mod 4 = 1 then # the product of alpha and gamma must be a nonsquare
     alpha := First(f1, a -> not IsZero(a) and not is_square( a*gamma, f1 ) );
   else # the product of alpha and gamma must be a square
    alpha := First(f1, a -> not IsZero(a) and is_square( a*gamma, f1 ) );
   fi;
   #map := NaturalEmbeddingByFieldReduction( ps1, f2, alpha, basis );
            form2 := QuadraticFormFieldReduction(form1,f2,alpha,basis);

  elif type1 = "parabolic" and type2 = "elliptic" and IsEvenInt(t) and IsOddInt(q) then
   Info(InfoFinInG, 1, "These polar spaces are suitable for field reduction: Parabolic -> Elliptic");
   bil1:=AssociatedBilinearForm(form1); # important to use this instead of 
   # BilinearFormByQuadraticForm (see documentation of the Forms package).
--> --------------------

--> maximum size reached

--> --------------------