#################################################################
 # DifSets.txt
 # Version 2.0
 # Ken W. Smith
 # Ryan L. Kaliszewski
 # Created in April 2006 from large collection of older procedures
 # Last modified June 2007.
 #     To run this program, insert the following Read statement
 # Read("DifSets.txt");
 #################################################################

 #################################################################
 # RightTranslates                                        1.1    #
 #################################################################
 Print("Reading file '(1.1) RightTranslates' with procedures: \n");
 Print("   RightTranslates( group, sublist ) \n");
 Print("   LeftTranslates( group, sublist )  \n");
 Print("   ConjugateTranslates( group, sublist ) \n");
 Print("   GroupMatrixRt( rep, group, elements ) \n");
 
 #################################################################
 # RightTranslates( group, sublist )                             #
 #################################################################
 #     RightTranslates( group, sublist ) takes a group and subset 
 # 'sublist' of the elements of the group and creates a list of all 
 # right translations of 'sublist', that is, the collection of 
 # all the images of 'sublist' under the action of right multiplication.  
 # The program returns a set of sublists of the elements of the group.
 #     A second, similar procedure, LeftTranslations( group, sublist )
 # creates the left translates of 'sublist'.
 #     A third procedure return all conjugate lists of 'sublist'.
 #
 #    If 'sublist' is a difference set in the group, then this procedure
 # will return the design created by all translates of that difference
 # set.
 #
 #     NOTE: Both the input and output of this procedure use
 # a list of group elements of a fixed group.
 #################################################################   
 RightTranslates := function( group, sublist )
 local k, elements, g, design, i, block;
    k        := Size( sublist );
    elements := Elements( group );
    design   := [ ];

    for g in elements do
       block := [ ];

       for i in [ 1..k ] do
          Add( block,  sublist[i]*g );
       od;

       Add( design, block );
    od;
  
    return design;
 end;
 #################################################################

 #################################################################
 # LeftTranslates( group, sublist )                              #
 #################################################################
 #	This procedure is a parallel to RightTranslates.  It reads
 # in a group and a subset 'sublist' of the group and creates the
 # collection of all images of 'sublist' under the action  
 # LEFT multiplication.
 #################################################################
 LeftTranslates := function( group, sublist )
 local k, elements, g, design, i, block;
    k        := Size( sublist );
    elements := Elements( group );
    design   := [ ];

    for g in elements do
       block := [ ];

       for i in [ 1..k ] do
          Add( block,  g*sublist[i] );
       od;

       Add( design, block );
    od;

    return design;
 end;
 #################################################################

 #################################################################
 # ConjugateTranslates( group, sublist )                         #
 #################################################################
 #	This procedure is a parallel to RightTranslates.  It reads
 # in a group and a subset 'sublist' of the group and creates the
 # collection of all images of 'sublist' under the action of 
 # conjugation.
 #################################################################   
 ConjugateTranslates := function( group, sublist )
 local k, elements, g, design, i, block;
    k        := Size( sublist );
    elements := Elements( group );
    design   := [ ];

    for g in elements do
       block := [ ];

       for i in [ 1..k ] do
          Add( block,  sublist[i]^g );
       od;

       Add( design, block );
    od;

    return design;
 end;
 #################################################################

 #################################################################
 # GroupMatrixRt( rep, group, elements )                         #
 #################################################################
 #   This procedure creates a group invariant matrix with with
 # representatives for the group elements.  Thus the first row
 # takes on the values of 'rep' (as a list). The group must have
 # order equalto the length of 'rep'.  'elements' is an ordered
 # list of the elements of the group.
 #################################################################
 # To make multiplication be on the left, change to:             #
 # if e[i]*e[j1]= e[j2] then block0[j1] := d[j2]; fi;            #
 #################################################################

 GroupMatrixRt := function( rep, group, elements )
 local v, mat, row, col, i, block;
    v   := Size( rep );
    mat := [  ];

    if v = Size( group ) then
       for row in [1..v] do
          block := ListWithIdenticalEntries( v, 0 );

          for col in [1..v] do
             for i in [1..v] do
                if elements[col]*elements[row]= elements[i] then 
                   block[col] := rep[i];  fi;
             od;
          od;

       Add( mat, block );
       od;

    else
       Print("The size of 'rep' is not equal "); 
       Print("to the order of the group!\n");
    fi;
    return mat;
 end;
 #################################################################




 #################################################################
 # IncidenceMatrix                                         1.2   #
 #################################################################
 Print("Reading file '(1.2) IncidenceMatrix' with procedures: \n");
 Print("   PointsToInterval( design ) \n");
 Print("   IncidenceMatrix( design ) \n");
 Print("   CharacteristicFunction( group, subset ) \n");

 # This file contains several procedures used to 
 # generate the incidence matrix of a design.
 #    PointsToInterval takes a design (a set of sets) and creates
 # an isomorphic design on the point set [1..v].
 #    SubsetToFirstRow takes subset of a group and creates from it
 # a (0,1)-list, for example, the first row of an incidence matrix
 #    IncidenceMatrix takes a design (a set of sets) on [1..v] 
 # and creates an incidence matrix. This last procedure assumes 
 # that the entries of the design are positive integers. IF
 # they are not then one MUST first run PointsToInterval( design ).
 
 #################################################################
 # PointsToInterval( design )                                    #
 #################################################################
 #   This procedure takes a design (a set of sets) and creates
 # an isomorphic design on the point set [1..v].
 # i, j, and s are iterates.
 #################################################################
 PointsToInterval := function( design )
 local blocks, v, k, points, intDesign, i, j, s;
    blocks := Size( design );
    points := Set( Union( design ) );
    v      := Size( points );

    intDesign := [ ];
    for i in [1..blocks] do 
       Add(intDesign, ShallowCopy( design[i] ) ); 
    od;

    for i in [1..blocks] do 
       k := Size( design[i] );
       for j in [1..k] do
          for s in [1..v] do 
             if design[i][j] = points[s] then 
                intDesign[i][j] := s;  fi;
          od;
       od;
    od;
    return intDesign;
 end;
 #################################################################

 #################################################################
 # IncidenceMatrix( design )                                     #
 #################################################################
 #   This procedure takes a design (a set of sets) on [1..v] 
 # and creates an incidence matrix.  This function works on 
 # designs that incorporate integers or other symbols.
 #################################################################
 IncidenceMatrix := function( design )
 local blocks, v, k, incMat, intDesign, i, j;
    intDesign := PointsToInterval(design);
    blocks    := Size( intDesign );
    v         := Size( Set( Flat( intDesign ) ) );
    k         := Size( intDesign[1] );

    incMat := NullMat(blocks,v);
    for i in [1..blocks] do 
       for j in [1..k] do
          incMat[i][intDesign[i][j]]:=1;
       od; 
    od;
    return incMat;
 end;
 #################################################################

 #################################################################
 # CharacteristicFunction( group, subset )                       #
 #################################################################
 #   This procedure takes subset of a group and creates from it
 # a (0,1)-list, for example, the first row of an incidence matrix
 #################################################################
 CharacteristicFunction := function( group, subset )
 local elements, ordSubset, ordGroup, incList, i, j;
    elements  := Elements( group );
    ordSubset := Size( subset );
    ordGroup  := Size( group );
    incList   := ListWithIdenticalEntries( ordGroup, 0 );

    for i in [1..ordSubset] do 
       for j in [1..ordGroup] do
          if subset[i] = elements[j] then 
             incList[j] := 1;  fi;
       od; 
    od;
    return incList;
 end;
 #################################################################




 #################################################################
 # ListToInteger.txt                                        1.3  #
 #################################################################
 Print("Reading file '(1.3) ListToInteger' with procedures: \n");
 Print("   IntegerToList( intIn, length, base )  \n");
 Print("   ListToInteger( vector, length, base )      \n");
 Print("   ListToMultiset( list )  \n");
 Print("   MultisetToList( multiset, listlength ) \n");

 # The procedure 'ListToInteger' turns a positive integer into a
 # list of a certain length.
 #    Using base 3, the integer 14 becomes [1,1,2] with length 3
 # or [0,0,0,0,0,1,1,2] with length 8.
 # The procedure 'IntegerToList' reverses that process.
 #################################################################
 #  IntegerToList( intIn, length, base )                         #
 #################################################################
 # This procedure turns an integer base 'base' into a list
 # Convert decimal integer to list 
 #  'intIn' is the integer we read in;
 #  'base' is the base for the computation.
 #  'length' is the length of the output vector
 #################################################################
 IntegerToList := function( intIn, length, base ) 
 local vector, i, minLength, quotient, remainder;

    minLength := LogInt( intIn, base ) + 1;
    if minLength > length then
       minLength := length;
    fi;

    vector := ListWithIdenticalEntries( length, 0 );
    for i in [1 .. minLength] do
       quotient  := QuoInt( intIn, base );
       remainder := RemInt( intIn, base );
       intIn     := quotient;
       vector[length + 1 - i] := remainder;
    od;
  
    return  vector;
 end;
 #################################################################

 #################################################################
 #  ListToInteger( vector, base )                                #
 #################################################################
 #   This procedure turns a string (as a list) of positive integers
 # into an integer.  This undoes the IntegerToList function.
 ################################################################# 
 ListToInteger := function( vector, base )
 local int_val, i, length; 
    int_val := 0;
    length  := Size( vector );

    for i in [1 .. length ] do
       int_val := int_val + vector[i] * base ^ ( length - i );
    od;

    return int_val;
 end;
 #################################################################
 
 #################################################################
 #  ListToMultiset( list )                                       #
 #################################################################
 #   This procedure turns a list of positive integers [a1, a2, ...]
 # into a set where the element 1 appears a1 times, a2 appears a2 
 # times and so on.  A (0,1)-list will be changed into the set of
 # "support" for the 1s.  (Eg. [1,1,0,1,0,0,0] becomes [1,2,4].)
 ################################################################# 
 ListToMultiset := function( list )
 local multiset, listsize, i, j; 
    multiset := [ ];
    listsize := Size( list );
 
    for i in [1..listsize] do 
       for j in [1..list[i]] do
          Add( multiset, i );
       od;
    od;

    return multiset;
 end;
 #################################################################

 #################################################################
 #  MultisetToList( multiset, listlength )                       #
 #################################################################
 #   This procedure turns a multiset of positive integers 
 # [b1, b2, ...] into a list where the number of appearances of
 # the element j in the multiset becomes the j-th coordinate of
 # the list.
 #   For example, if the list length is five then [1,3,3] will 
 # be turned into [1,0,2,0,0].  The difference set [1,2,4] in
 # C_7 will be turned into [1,1,0,1,0,0,0].
 #   This procedure is the inverse of ListToMultiset. 
 #   
 #   Known Issues:  If listlength isn't long enough to hold the
 # elements of the multiset, the list will be truncated.
 #################################################################
 MultisetToList := function( multiset, listlength )
 local list, i; 
    list := ListWithIdenticalEntries( listlength, 0 );

    for i in multiset do 
       if i <= listlength then
          list[i] := list[i] + 1;  fi;
    od;    

    return list;
 end;
 



 #################################################################
 # Convolution.txt                                        2.1    #
 #################################################################
 Print("Reading file '(2.1) Convolution' with procedures: \n");
 Print("   ConvolutionTable_f( group, elements, func ) \n");
 Print("   Convolution_AsList( CL, list1, list2 )      \n");
 Print("   Convolution_AsMultiset( CL, multset1, multset2 )  \n");

 #################################################################
 # ConvolutionTable_f( group, elements, func )                   #
 #################################################################
 #
 #   This procedure takes a group, a sorted list of the elements
 # of the group and 'func', a function previously defined on
 # the elements of the group.
 #   It returns a table whose (i,j) entry is the
 # positive integer k such that e[i] * f( e[j] ) = e[k].
 #   This is especially useful for difference sets if 
 # the function 'f' was defined by 
 #                       f:= x-> x^-1;
 #################################################################
 ConvolutionTable_f := function( group, elements, func )
 local order, CL, i0, i1, i2;
    order := Size( group );
    CL    := [ ]; 

    for i0 in [1..order] do 
       CL[i0] := [ ]; 
    od;

    for i0 in [1..order] do
       for i1 in [1..order] do 
          for i2 in [1..order] do
             if elements[i1]*func(elements[i2]) = elements[i0] then 
                CL[i1][i2] := i0;  fi;
          od; 
       od;
    od;

    return CL;
 end;
 #################################################################

 #################################################################
 # Convolution_AsList( CL, list1, list2 )                        #
 #################################################################
 #
 #   This procedure takes table CL (created, say, by 
 # ConvolutionTable_f( g, e, f ) and two lists, list1 and list2,
 # (Members of the integer group ring Zg here are represented as
 # a list of coefficients, so for example, the set e[1]+2*e[3]
 # in Z(C_5) will appear as [1,0,2,0,0].) 
 #   It then uses the convolution list CL to quickly compute 
 # list1*list2 in Zg and returns that answer as a list.
 #################################################################
 Convolution_AsList := function( CL, list1, list2 )
 local sizeCL, conprod, i1, i2;
    sizeCL  := Size( CL );
    conprod := ListWithIdenticalEntries( sizeCL, 0 );

    for i2 in [1..sizeCL] do 
       for i1 in [1..sizeCL] do  
          conprod[i2] := conprod[i2] + list1[i1]*list2[ CL[i1][i2] ];
       od; 
    od;
    return conprod;
 end;
 #################################################################

 #################################################################
 # Convolution_AsMultiset( CL, multset1, multset2 )              #
 #################################################################
 #
 #   This procedure takes table CL (created, say, by 
 # ConvolutionTable_f( g, e, f ) and two multisets, multset1 and 
 # multset2.
 # (Members of the integer group ring Zg here are represented as
 # a multiset, for example the element e[1]+2*e[3] in Z(C_5) 
 # will appear as [1,3,3].)
 #   The procedure then uses the convolution list CL to quickly 
 # compute multset1*multset2 in Zg and returns that answer 
 # as a *multiset*.
 #################################################################
 Convolution_AsMultiset := function( CL, multset1, multset2 )
 local conprod, i1, i2;
    conprod := [ ];

    for i2 in multset2 do 
       for i1 in multset1 do  
          Add( conprod, CL[i1][i2] );
       od; 
    od;
    Sort( conprod );

    return conprod;
 end;
 #################################################################




 #################################################################
 # ImagesOfFixedSize.txt                                   2.2   #
 #################################################################
 # Version 3.0
 # Ken W. Smith
 # Created in December 1999, last modified in May 2006.
 #
 # This file contains a procedure used to list homomorphic
 # images of finite groups

 Print("Reading file '(2.2) ImagesOfFixedSize' with procedure: \n");
 Print("   ImagesOfFixedSize( order, domain, im_order, range, printflag ) \n");

 #################################################################
 # ImagesOfFixedSize( order, domain, im_order, range, printflag )#
 #################################################################
 #    We look at groups of order 'order' with images of order 
 # 'im_order'.
 #    We restrict our search to groups [order, j] where j is in 
 # 'domain'.
 #    We restrict the image space to groups [im_order, i] where i 
 # is in 'range'.
 #    The procedure returns a list, 'solution', containing all j
 # where [n,j] has an image [size, i] where i in 'range'.
 #    If printflag is 0 then no details are printed during the 
 # procedure run.  Printflag values of 1 and 2 provide more details.  
 # 
 #    To receive an exhaustive list, merely set the domain (or range)
 # to be empty.  This will be interperated as making the search space
 # as large as possible:
 # IF the DOMAIN is EMPTY then we look at ALL groups of order 'order'.
 # If the RANGE is EMPTY then there are NO restrictions on the 
 # range.
 #################################################################

 ImagesOfFixedSize := function( order, domain, im_order, range, printflag )
 local ndom, nran, i, j, group, lattice, normal, is_image, solution;
    solution := [ ];
    ndom     := Size( domain );  
    nran     := Size( range );  
  
   
    if ndom = 0 then 
       domain := [1..NumberSmallGroups(order)]; 
       ndom := NumberSmallGroups(order);  
    fi;
    domain := Intersection( domain, [1..NumberSmallGroups(order)] );

    if nran = 0 then 
       range := [1..NumberSmallGroups( im_order )];  fi;
    range := Intersection(range, [1..NumberSmallGroups( im_order )]);

    if printflag >= 1 then 
       Print("We examine ", ndom , " groups of order ", order,"\n"); 
       Print("Our image space is ", range, "\n\n");
    fi;


    for j in domain do
       group := SmallGroup( order, j );
       if printflag >= 1 then 
          Print( j, ".");
          if printflag >= 2 then 
             Print(" AbelianInvariants = ", AbelianInvariants( group ));  fi;
          Print("\n");
       fi;

       lattice := NormalSubgroups( group );
       is_image := 0;
       for i in [1..Size(lattice)] do
          normal := lattice[i];
          if Size(normal) = order/im_order then 
             if IdGroup( group/normal )[2] in range then
                Add( solution, j );
                is_image := 1; 
                if printflag >= 1 then 
                   Print( "  (", IdGroup( normal )," = subgroup ",i ,") by image ");
                   Print( IdGroup(group/normal), " \n"); 
                fi;
             fi;
          fi;
       od;

       if is_image = 0 and printflag >= 1 then 
          Print("No acceptable image for group ", order,".",j,"\n");  fi;
       if printflag >= 1 then 
          Print("\n"); fi;
   od;
   solution := Set( solution );

   return solution;
 end;
 #################################################################




 #################################################################
 # RightTransversalList                                   2.3    #
 #################################################################
 Print("Reading file '(2.3) RightTransversalList' with procedures: \n");
 Print("   RightTransversalList( group, elements, normal ) \n");
 Print("   TransversalListSpecific( group, elements, normal ) \n");
 Print("   Expand( TL, d )  \n");
 Print("   InverseExpand( TL, d )  \n");

 ##################################################################
 # RightTransversalList( group, elements, normal )                #
 ##################################################################
 #   This procedure takes a group, a sorted list of elements
 # of the group, a normal subgroup, and returns a list of cosets of
 # the normal subgroup in the group.  The elements of the cosets 
 # are positive integers; if the integer i0 is listed then e[i0] 
 # is in that coset. 
 #################################################################
 RightTransversalList := function( group, elements, normal )
 local normal_elements, order, index, rightT, TL, i0, i1, i2;
    normal_elements := Elements( normal );
    order           := Size( group ); 
    index           := order/Size( normal );
    rightT          := RightTransversal( group, normal );
    TL := [ ];
    
    for i2 in [1..index] do
       Add( TL, [ ] );
       for i1 in [1..Size( normal )] do 
          for i0 in [1..order] do
             if normal_elements[i1]*rightT[i2] = elements[i0] then
                 Add( TL[i2], i0 ); 
             fi;
          od; 
       od; 
    od;

    return TL;
 end;

 ##################################################################
 # TransversalListSpecific( group, elements, normal )             #
 ##################################################################
 #   This procedure takes a group, a sorted list of elements
 # of the group, a normal subgroup, and returns a list of cosets of
 # the normal subgroup in the group.  The elements of the cosets 
 # are positive integers; if the integer i0 is listed then e[i0] 
 # is in that coset.
 #    The program has been written so that the elements of the factor
 # group are specifically tied to the elements of a group in GAP's 
 # small group library.
 #################################################################
 TransversalListSpecific := function( group, elements, normal )
 local index, factor_id, hom_grp, factor_grp, hom, factor_ele, iso, i, j, TL;
    index      := Size( group/normal );
    factor_id  := IdGroup( group/normal )[2];
    hom        := NaturalHomomorphismByNormalSubgroup(group,normal);
    hom_grp    := Images( hom );
    factor_grp := SmallGroup( index, factor_id );
    factor_ele := Elements( factor_grp );
    iso        := IsomorphismGroups( hom_grp, factor_grp );

    TL := [ ];
    for i in [1..Size(factor_grp)] do 
       Add( TL, [ ] );             
       for j in [1..Size(group)] do 
          if Image( iso, Image( hom, elements[j]) ) = factor_ele[i] then 
             Add( TL[i], j );
          fi;
       od;
    od;

    return TL;
 end;
 ##################################################################

 ##################################################################
 # Expand( TL, d )                                                #
 ##################################################################
 #   This procedure takes, TL, a list of lists, and a list d, viewed
 # as a function on TL, and expands it to a list of lists which is 
 # constant across the inside lists and when restricted to TL, is
 # consistent with d.
 #   For example, if TL = [ [1,2],[3,4], [5,6] ] and d = [ x,y,z ] then
 # the output is [ [x,x], [y,y], [z,z] ].
 #   This is useful when given an element d in Z(g/N) and expanding
 # it in Zg.
 ##################################################################
 Expand := function( TL, d )
 local i, j, output;
    output := [ ];

    for i in [1..Size(TL)] do 
       for j in [1..Size(TL[i])] do
          output[ TL[i][j] ] := d[i];
       od; 
    od;

    return output;
 end;
 ##################################################################

 ##################################################################
 # InverseExpand( TL, d )                                        #
 ##################################################################
 #   This procedure reverses Expand( TL, d ).
 # It does not check for errors; if two elements of the same coset
 # are assigned a different value, the value assigned to the coset
 # is simply the value of the first element found in its list.
 ##################################################################
 InverseExpand := function( TL, d )
 local i, j, output;
    output := [ ];

    for i in [1..Size(TL)] do 
       output[i] := d[ TL[i][1] ]; 
    od;

    return output;
 end;
 ##################################################################




 #################################################################
 # Subgroup chains                                        2.4    #
 #################################################################

 Print("Reading file '(2.4) Subgroup Chains' with procedures: \n");
 Print("   SubgroupIncidenceMatrix( order, group_id ) \n");
 Print("   MaxChains( mat ) \n");
 Print("   MaxSubgroupChains( mat, subgroup ) \n");

 #################################################################
 # SubgroupIncidenceMatrix( order, group_id )                    #            
 #################################################################
 # The following procedure creates the poset incidence matrix
 # of the lattice of normal subgroups of the group 
 # SmallGroup( order, group_id ).
 #################################################################
 SubgroupIncidenceMatrix := function( order, group_id )
 local group, elements, normals, ns_size, mat, i, j;
    group    := SmallGroup( order, group_id ); 
    elements := Elements( group );  
    normals  := NormalSubgroups( group );
    ns_size  := Size( normals );

    mat := NullMat( ns_size, ns_size );
    for i in [1..ns_size] do 
       for j in [1..ns_size] do
          if (IsSubset( normals[j], normals[i] ) and not i=j) then 
             mat[i][j] := 1;  fi;
       od;
    od;
  
   return mat;
 end;
 #################################################################

 #################################################################
 # MaxChains( mat )                                              #
 #################################################################
 #    The following procedure inputs the incidence matrix of a
 # poset and returns a list of all the maximal chains.
 #################################################################
 # Must be a (0,1) matrix of a poset.  (In particular, there cannot be
 # ones on the diagonal.)  A matrix created by SubgroupIncidenceMatrix
 # is ideal, here.
 #################################################################
 MaxChains := function( mat )
 local size, flag, chain_length, edges, i1, i2, chains, i, newchains, 
       c, e, c1;
    size         := Size(mat);
    flag         := 0;
    edges        := [ ];
    chains       := [ ];
    chain_length := 1;
    while not( mat^(chain_length+1) = NullMat( size, size ) ) do
       chain_length := chain_length + 1;
    od;

    for i1 in [1..size] do 
       for i2 in [1..size] do 
          if mat[i1][i2] = 1 then 
             Add( edges, [i1,i2] );  fi; 
       od;
    od;
    Add( chains, edges );

    for i in [2..chain_length] do
       newchains := [ ];
       for c in chains[i-1] do 
          for e in edges do
             if c[i]=e[1] then 
                c1 := ShallowCopy( c ); 
                Add( c1, e[2] ); 
                Add( newchains, c1 );
             fi;
          od; 
       od;
       Add( chains, newchains );
    od;

    return chains[ chain_length ];
 end;
 #################################################################

 #################################################################
 # MaxSubgroupChains( mat, subgroup )                            #            
 #################################################################
 #    The following procedure inputs the incidence matrix of a
 # poset and returns a list of all the maximal chains that contain
 # a specified subgroup.
 #    The procedure requires the parameter 'subgroup', which is the
 # subgroup number in NormalSubgroups( SmallGroup( order, group_id ) ).
 #    NOTE:  The procedure only finds maximal chains that contain
 # the specified subgroup.  If there is no maximal chain that
 # contains the specified subgroup then the procedure returns an
 # empty list.
 #################################################################
 # Must be a (0,1) matrix of a poset.  (In particular, there cannot be
 # ones on the diagonal.)  A matrix created by SubgroupIncidenceMatrix
 # is ideal, here.
 #################################################################
 MaxSubgroupChains := function( mat, subgroup )
 local maxchains, x, newchains;
    maxchains := MaxChains( mat );
    newchains := [ ];

    for x in maxchains do
       if subgroup in x then
          Add( newchains, x );  fi;
    od;

    return newchains;
 end;

 #################################################################
 # Partition Functions                                       2.5 #
 #################################################################
 Print("Reading file '(2.5) Partition Functions' with procedures: \n");
 Print("   PartitionRep( partition ) \n");
 Print("   IndexToRep( partition ) \n");
 Print("   QuotientGroupConvolution( CL, partition ) \n");
 Print("   MapRefinement( partition, refinement ) \n");

 #################################################################
 # PartitionRep( partition )                                     #
 #################################################################
 # The following procedure takes a partition on [1..v] and
 # creates a many-to-one map from [1..v] into [1..v] where
 # an element x is mapped to its "representative", the first 
 # element in the partition containing x.
 # 
 # The range of this function is returned, as a list.
 #################################################################
 # partition is a set of disjoint sets whose union is [1..v]
 PartitionRep := function( partition )
 local setSize, partSize, rep, i, part;
    setSize  := Size(Flat(partition));
    partSize := Size(partition);

    rep := ListWithIdenticalEntries( setSize, 0 );
    for part in partition do 
       for i in part do 
          rep[i] := part[1]; 
       od; 
    od;
    return rep;
 end;
 #################################################################

 #################################################################
 # IndexToRep( partition )                                       #
 #################################################################
 # The following procedure takes a partition on [1..v] and
 # creates a many-to-one map from [1..v] into [1..v] where
 # an element of the i-th partition is mapped to i.  (i is 
 # the "index" of the partition.
 # 
 # The range of this function is returned, as a list.
 #################################################################
 # partition is a set of disjoint sets whose union is [1..v]
 IndexToRep := function( partition )
 local parts, index_to_rep, i;
    parts := Size(partition);

    index_to_rep := ListWithIdenticalEntries( parts, 0 );
    for i in [1..parts] do 
       index_to_rep[i]:=partition[i][1]; 
    od;

    return index_to_rep;
 end;
 #################################################################

 #################################################################
 # QuotientGroupConvolution( CL, partition )                     #
 #################################################################
 # The following procedure takes a partition on [1..v] and
 # a convolution table CL on [1..v] and creates a smaller
 # convolution table on just the representatives (first
 # element of each class) of the partition.  
 #
 # If table CL is the Cayley table of a group G and the
 # partition is the cosets of a normal subgroup N, then
 # the output is a Cayley table for the group G/N.
 # 
 # The output is a Cayley table with elements which are 
 # representatives; the elements are integers but will 
 # not necessarily be 1, 2, 3, ...
 #################################################################
 QuotientGroupConvolution := function( CL, partition )
 local parts, index_to_rep, CL_small, i, rep_i, j, rep_j, partrep;
    parts        := Size( partition );
    partrep      := PartitionRep( partition );
    index_to_rep := IndexToRep( partition );
    CL_small     := NullMat( parts, parts );

    for i in [1..parts] do 
       rep_i := index_to_rep[i];
       for j in [1..parts] do
          rep_j := index_to_rep[j];
          CL_small[i][j] := partrep[CL[ rep_i ][ rep_j ]];
       od;
    od;
    return CL_small;
 end;
 #################################################################

 #################################################################
 # MapRefinement( partition, refinement )                        #
 #################################################################
 # The following procedure takes a partition on [1..v] and
 # a refinement of it (another partition on [1..v] with the
 # property that every class of this partition is a subclass
 # of the original.)  The output is a list of lists.  In 
 # the output, each class in the original partition has been
 # replaced by a list of indices corresponding to the classes
 # of the refinement which belong to that partition.
 #
 # For example, if the output is [[1,3,5], [2,4]] then the
 # original partition had two classes.  The first class
 # is the union of classes 1, 3 and 5 of the refinement.
 # The second class (of the original partition) is the union
 # of the second and fourth classes of the refinement.
 #################################################################
 MapRefinement := function( partition, refinement )
 local map, i, j, vector;
    map:=[ ];

    for j in [1..Size(partition)] do  
       vector := [ ];
       for i in [1..Size(refinement)] do 
          if IsSubset( partition[j], refinement[i]) then 
             Add( vector, i );  fi; 
       od; 
       Add( map, vector );
    od;

    return map;
 end;
 #################################################################




 #################################################################
 # Using Automorphisms                                       2.6 #
 #################################################################
 Print("Reading file '(2.6) Using Automorphisms' with procedures: \n");
 Print("   AutImages( group, elements ) \n");
 Print("   MaxValueOnSets( right_reg, aut_img, subset ) \n");
 Print("   MaxValueOnTuples( right_reg, aut_img, subset ) \n");
 Print("   ReduceSets( right_reg, aut_img, soln , printflag) \n");
 Print("   ReduceTuples( right_reg, aut_img, partition, soln, printflag ) \n");

 # For the following functions, the definitions of the parameters are
 # given as follows (where group is a group and elements is a list of the
 # elements of the group).
 #   right_reg := Image(ActionHomomorphism( group, elements, OnRight ));
 #   aut_grp   := AutomorphismGroup( group );
 #   aut_img   := Action( aut_grp, elements );


 #################################################################
 # AutImages( group, elements )                                  #
 #################################################################
 #    This procedure takes a group and its elements as parameters 
 # and returns a list with the following values:
 #   list := AutImages( group, elements );
 #
 #   list[1] = right_reg
 #   list[2] = aut_grp
 #   list[3] = aut_img
 #
 #################################################################
 AutImages := function( group, elements )
 local list;
    list := [];

    list[1] := Image(ActionHomomorphism( group, elements, OnRight ));
    list[2] := AutomorphismGroup( group );
    list[3] := Action( list[2], elements );

    return list;
 end;

 #################################################################
 # MaxValueOnSets( right_reg, aut_img, subset )                  #
 #################################################################
 #    This procedure takes right_reg and aut_img (defined above), 
 # and an integral set representing a subset of the group.
 # It then gives the element which is lexicographically least 
 # representative in the equivalence class of the subset.  
 # (Equivalence here involves translation and automorphisms.)  
 #
 # 'subset' is assumed to be a set, such as a difference set, with
 # no repeated elements.
 #################################################################
 MaxValueOnSets := function( right_reg, aut_img, subset )
 local orbit, row, mat, x, FinalList, orbit_list;
    orbit     := Orbit( right_reg, subset, OnSets);  
    mat       := [ ]; 
    FinalList := [ ];

    for x in orbit do 
       if 1 in x then 
          Add( mat, x ); fi; 
    od;

    for row in mat do
       orbit_list := Orbit( aut_img, row, OnSets );
       orbit_list := Set(orbit_list);
       Add(FinalList, orbit_list[1] );
    od;
    Sort(FinalList);

    return FinalList[1];
 end;

 #################################################################
 # MaxValueOnTuples( right_reg, aut_img, subset )                #
 #################################################################
 #    This procedure takes right_reg and aut_img (defined above), 
 # and an integral set representing a subset of the group.
 # It then gives the element which is lexicographically least 
 # representative in the equivalence class of the subset.  
 # (Equivalence here involves translation and automorphisms.)  
 #
 # 'subset' is assumed to be a multiset.
 #################################################################
 MaxValueOnTuples := function( right_reg, aut_img, subset )
 local orbit, mat, x, FinalList, row, orbit_img, sorted_img, x1;
    orbit     := Orbit( right_reg, subset, OnTuples );
    mat       := [ ]; 
    FinalList := [ ];

    for x in orbit do 
       if 1 in x then 
          Add( mat, x );  fi; 
    od;

    for row in mat do  
       orbit_img  := Orbit( aut_img, row, OnTuples );
       sorted_img := [ ]; 

       for x in orbit_img do 
          x1 := ShallowCopy(x); 
          Sort(x1); 
          Add(sorted_img, x1); 
       od;

       Sort(sorted_img);  
       Add(FinalList, sorted_img[1] );
    od;
    Sort(FinalList);
 
    return FinalList[1];
 end;
 #################################################################

 #################################################################
 # EquivalenceClass( right_reg, aut_img, subset )                #
 #################################################################
 #    This procedure takes right_reg and aut_img (defined above)
 # and then returns the entire equivalence class, along with
 # the element which is lexicographically least representative
 # in the equivalence class of 'subset'.  
 # (Equivalence here involves translation and automorphisms.)  
 #
 # subset is assumed to be a multiset
 #################################################################
 EquivalenceClass := function( right_reg, aut_img, subset )
 local orbit, mat, x, FinalList, EquivalenceClass, row, orbit_img, sorted_img, x1;
    orbit            := Orbit( right_reg, subset, OnTuples );
    mat              := [ ]; 
    EquivalenceClass := [ ];

    for x in orbit do 
       if 1 in x then 
          Add( mat, x );  fi; 
    od;

    for row in mat do  
       orbit_img  := Orbit( aut_img, row, OnTuples );
       sorted_img  := [ ]; 

       for x in orbit_img do 
          x1 := ShallowCopy(x);
          Sort(x1);
          Add( sorted_img, x1 );
       od;

       Sort(sorted_img); 
        EquivalenceClass := Union( EquivalenceClass, sorted_img);
    od;

    return EquivalenceClass;
 end;
 #################################################################

 #################################################################
 # ReduceSets( right_reg, aut_img, soln , printflag )            #
 #################################################################
 #    This procedure reduces a collection of lists (such as proposed
 # difference sets) to a single representative from each equivalence
 # class.  Thus it removes (automorphic or shift) equivalent lists
 # from the collection, leaving a single representative from the
 # equivalence class.
 #    'right_reg' and 'aut_img' are defined above.  The procedure
 # also requires  a collection, 'soln', that will be reduced.
 #    If printflag is set to 1 then this procedure will print the
 # time it takes to maximize the first set.  If it is set to 0 it
 # will not.
 #    NOTE:  This procedure assumes that the lists in soln are
 # sets with no duplicate entries (such as a difference set).
 #################################################################
 ReduceSets := function( right_reg, aut_img, soln, printflag )
 local FinalSolution, Starttime_RT, rep, x;
    FinalSolution := [ ];
    Starttime_RT  := Runtime( ); 

    for x in soln do
       rep     := MaxValueOnSets( right_reg, aut_img, x );

       Add( FinalSolution, rep );

       if printflag = 1 then 
          printflag := 0;
          Print("First set maximized in ", Runtime()-Starttime_RT," milliseconds.\n");
       fi;
    od;

    FinalSolution := Set( FinalSolution );
    return FinalSolution;
  end;

 #################################################################

 #################################################################
 # ReduceTuples( right_reg, aut_img, partition, soln, printflag )#
 #################################################################
 #    This procedure reduces a collection of lists (such as proposed
 # difference sets) to a single representative from each equivalence
 # class.  Thus it removes (automorphic or shift) equivalent lists
 # from the collection, leaving a single representative from the
 # equivalence class.
 #    'right_reg' and 'aut_img' are defined above.  The procedure
 # also requires a partition of the group and a collection, 'soln', 
 # that will be reduced.  The partition is necessary because of the
 # interpretation of the collection, 'soln'.
 #    If printflag is set to 1 then this procedure will print the
 # time it takes to maximize the first set.  If it is set to 0 it
 # will not.
 # 
 # We assume that 'soln' is a list of indices of cosets
 # (eg. [1,1,1,2,2,3,4,5,5])
 # and not a list of representatives
 # (eg. [1,1,1,8,8,13,11, 17,17])
 #################################################################
 ReduceTuples := function( right_reg, aut_img, partition, soln, printflag )
 local i, j, mult, rep, Starttime_RT, firstflag, x, FinalSolution;
    mult          := [ ];
    Starttime_RT  := Runtime( ); 
    FinalSolution := [ ];

    for i in [1..Size(soln)] do
       mult[i] := [ ];

       for j in [1..Size(soln[i])] do
          mult[i][j] := partition[soln[i][j]][1];
       od;
    od;

    for x in mult do
       rep     := MaxValueOnTuples( right_reg, aut_img, x );

       Add( FinalSolution, rep );

       if printflag = 1 then 
          printflag := 0;
          Print("First set maximized in ", Runtime()-Starttime_RT," milliseconds.\n");
       fi;
    od;

    FinalSolution := Set( FinalSolution );
    return FinalSolution;
 end;
 #################################################################




 #################################################################
 # Difference set tables                                     3.1 #
 #################################################################
 Print("Reading file '(2.7) Procedures for difference set tables': \n");
 Print("   Verify( parameters, group_id, difset ) \n");
 Print("   IncidenceMatrix_Parameters( parameters, group_id, difset ) \n");
 Print("   MultiplierGroup( parameters, group_id, difset) \n");
 
 #################################################################
 # Verify( [v,k,lambda], group_id, difset )                      #
 #################################################################
 #   This procedure takes a group ring element 'difset' 
 # with parameters [v,k,lambda], in the group SmallGroup(v,group_id) 
 # and verifies that D*D^(-1) only has two values on group 
 # elements.  
 #   The group ring element 'difset' is given as a list of 
 # positive integers {i : e[i] is in the difference set}.
 #   This returns 1 on a posititive result, and a -1 on a negative.
 #################################################################
 Verify := function( parameters, group_id, difset )
 local v, k, lambda, group, elements, inversefunc, CL, intersection_nums,
       result;
    v           := parameters[1];
    k           := parameters[2]; 
    lambda      := parameters[3];
    group       := SmallGroup(v, group_id); 
    elements    := Elements(group);
    inversefunc := x-> x^-1;
    CL          := ConvolutionTable_f( group, elements, inversefunc );

    intersection_nums := Size( Set( MultisetToList( 
                            Convolution_AsMultiset( CL, difset, difset ), v) ) );

    if intersection_nums = 2 then 
       result := 1;
    else 
       result := -1;
    fi;
    return result;
 end;

 #################################################################
 # IncidenceMatrix_Parameters( [v,k,lambda], group_id, difset )   #
 #################################################################
 #   This procedure takes a group ring element 'difset' with 
 # parameters [v,k,lambda], in the group SmallGroup(v,group_id) and 
 # constructs the incidence matrix of the design of right 
 # translates.
 #  The group ring element 'difset' is given as a list of 
 # positive integers {i : e[i] is in the difference set}.
 #################################################################
 IncidenceMatrix_Parameters := function( parameters, group_id, difset )
 local v, k, lambda, group, elements, i, ds_list, design, mat, elemdiv;
    v        := parameters[1];
    k        := parameters[2]; 
    lambda   := parameters[3];
    group    := SmallGroup(v, group_id); 
    elements := Elements(group);
    ds_list  := [ ];  

    for i in difset do 
       Add(ds_list, elements[i]); 
    od;

    design := PointsToInterval( RightTranslates( group, ds_list ) );
    mat := IncidenceMatrix( design );
    # elemdiv := Collected( ElementaryDivisorsMat( mat ) );

    return mat;
 end;
 #################################################################

 #################################################################
 # MultiplierGroup( [v,k,lambda], group_id, difset )             #
 #################################################################
 #   This procedure takes a group ring element 'difset' with
 # parameters [v,k,lambda], in the group SmallGroup(v,group_id) and
 # computes the multiplier group of the difference set.  
 # (The multiplier group is the set of automorphisms of the 
 # group which map 'difset' to a right translate.)
 #
 # The group ring element 'difset' is given as a list of 
 # positive integers {i : e[i] is in the difference set}.
 #################################################################
 MultiplierGroup := function( parameters, group_id, difset )
 local v, k, lambda, group, elements, aut_grp, aut_img, a, i, ds_list, 
       design, sorted_design, stab;
    v        := parameters[1];
    k        := parameters[2]; 
    lambda   := parameters[3];
    group    := SmallGroup(v, group_id); 
    elements := Elements(group);
    aut_grp  := AutomorphismGroup(group);
    aut_img  := Action( aut_grp, elements ); 
    ds_list  := [ ]; 

    for i in difset do 
       Add(ds_list, elements[i] ); 
    od;
    design := PointsToInterval( RightTranslates( group, ds_list ) );

    sorted_design := [ ]; 
    for i in design do  
       Add( sorted_design, Set(i) ); 
    od;

    sorted_design := Set(sorted_design);
    stab := Stabilizer( aut_img, Set( sorted_design ), OnSetsSets );
    
    return stab;
 end;
 #############################################################

 Print("\n   'DifSets.txt' procedures have been loaded into GAP.\n\n");