text:
InstallGlobalFunction(MitM_SCSCPServerHandshake,
function(in_stream, out_stream)
    local msg, get, pi, client_versions, versions, version;
    # Send connection initiation message (5.1.1)
    msg := Concatenation("<?scscp service_name=\"{}\" service_version=\"{}\" ",
                         "service_id=\"{}\" scscp_versions=\"{}\" ?>");
    msg := StringFormatted(msg,
                           "gap", "4.10dev", String(IO_getpid()),
                           JoinStringsWithSeparator(MitM_SCSCPVersions, " "));
    WriteLine(out_stream, msg);
    # TODO: we are required to support SCSCP v1.0

    # Version negotiation (5.1.2)
    get := MitM_ReadToPI(in_stream);
    if get.success = false then
        Info(InfoMitMSCSCP, 2,
             "Version negotiation failed: no version request from client");
        return fail;
    fi;
    pi := Concatenation("<", get.pi, ">");
    pi := GetSTag(pi, 2);
    if not IsBound(pi.attributes.version) then
        Info(InfoMitMSCSCP, 2,
             "Version negotiation failed: no version request from client");
        return fail;
    fi;
    client_versions := SplitString(pi.attributes.version, "", " \n\r\t");
    versions := Filtered(client_versions, v -> v in MitM_SCSCPVersions);
    if IsEmpty(versions) then
        Info(InfoMitMSCSCP, 2,
             Concatenation("Version negotiation failed: client requested ",
                           Concatenation(client_versions),
                           ", none of which is supported"));
        WriteLine(out_stream,
                  "<?scscp quit reason=\"not supported version\" ?>");
        return fail;
    fi;
    version := Maximum(versions); # choose the highest lexicographically
    WriteLine(out_stream, Concatenation("<?scscp version=\"",
                                        version, "\" ?>"));
    # Successful handshake, now ready for procedure calls
    return version;
end);

InstallGlobalFunction(MitM_SCSCPClientHandshake,
function(in_stream, out_stream)
    local get, pi, start, finish;
    # Get connection initiation message (5.1.1)
    get := MitM_ReadToPI(in_stream);
    if not get.success then
        Info(InfoMitMSCSCP, 2,
             "Initiation failed: no processing instruction received");
        return false;
    fi;
    # Use XML parser to get at the information
    pi := Concatenation("<", get.pi, ">");
    pi := GetSTag(pi, 2);
    if pi.name <> "scscp" then
        Info(InfoMitMSCSCP, 2,
             "Initiation failed: no SCSCP instruction received");
        return false;
    elif Set(RecNames(pi.attributes)) <> ["scscp_versions",
                                          "service_id",
                                          "service_name",
                                          "service_version"] then
        Info(InfoMitMSCSCP, 2,
             "Initiation failed: bad connection initiation message received");
        return false;
    elif not "1.3" in SplitString(pi.attributes.scscp_versions, "", " ") then
        # TODO: we are required to support SCSCP v1.0
        Info(InfoMitMSCSCP, 2,
             "Initiation failed: server offers no SCSCP version we know");
        return false;
    fi;

    # Version negotiation (5.1.2)
    WriteLine(out_stream, "<?scscp version=\"1.3\" ?>");
    get := MitM_ReadToPI(in_stream);
    if not get.success then
        Info(InfoMitMSCSCP, 2,
             "Version negotiation failed: no valid response from server");
        return false;
    elif PositionSublist(get.pi, "quit") <> fail then
        Info(InfoMitMSCSCP, 2,
             "Version negotiation failed: server quit with following message:");
        start := Position(get.pi, '\"');
        finish := Position(get.pi, '\"', start);
        Info(InfoMitMSCSCP, 2, get.pi{[start..finish]});
        return false;
    fi;
    pi := Concatenation("<", get.pi, ">");
    pi := GetSTag(pi, 2);
    if not IsBound(pi.attributes.version) then
        Info(InfoMitMSCSCP, 2,
             "Version negotiation failed: server sent no version number");
        return false;
    elif pi.attributes.version <> "1.3" then
        Info(InfoMitMSCSCP, 2,
             Concatenation("Version negotiation failed: server chose version ",
                           pi.attributes.version,
                           ", which is not supported"));
        return false;
    fi;

    # Successful handshake, now ready for procedure calls
    return true;
end);

InstallGlobalFunction(MitM_ReadSCSCP,
function(stream)
    local get, pi, r;
    # Get first processing instruction
    get := MitM_ReadToPI(stream);
    # It should be an SCSCP start instruction
    if not get.success then
        return MitM_Error("no processing instruction found");
    fi;
    pi := get.pi;
    pi := SplitString(pi, "", " \n\r\t");
    if pi[1] <> "scscp" or pi[2] <> "start" or Length(pi) > 2 then
        return MitM_Error("no SCSCP start instruction found");
    fi;
    # Read to the end instruction
    get := MitM_ReadToPI(stream);
    if not get.success then
        return MitM_Error("only one processing instruction found");
    fi;
    pi := get.pi;
    pi := SplitString(pi, "", " \n\r\t");
    if pi[1] <> "scscp" or pi[2] <> "end" or Length(pi) > 2 then
        return MitM_Error("no SCSCP end instruction found");
    fi;
    # The in-between stuff should be an OMOBJ
    r := MitM_XMLToOMRec(get.pre);
    if MitM_Tag(r) <> "OMOBJ" then
        return MitM_Error("no OMOBJ object found");
    fi;
    return MitM_Result(r);
end);

InstallGlobalFunction(MitM_ReadToPI,
function(stream)
    local pre, len, char, byte, error, lastchar, pi;
    pre := "";
    len := 0;
    char := fail;
    repeat
        byte := ReadByte(stream);
        if byte = fail then
            return rec(success := false, pre := pre,
                       error := "stream ended before an XML PI was found");
        fi;
        lastchar := char;
        char := CharInt(byte);
        len := len + 1;
        pre[len] := char;
    until lastchar = '<' and char = '?';
    Remove(pre);
    Remove(pre);
    pi := "";
    len := 0;
    repeat
        byte := ReadByte(stream);
        if byte = fail then
            return rec(success := false, pre := pre, pi := pi,
                       error := "stream ended in middle of XML PI");
        elif byte = 60 then
            return rec(success := false, pre := pre, pi := pi,
                       error := "'<' character illegal inside XML PI");
        fi;
        lastchar := char;
        char := CharInt(byte);
        len := len + 1;
        if len > 4092 then # max length 4094 including "<?"
            return rec(success := false, pre := pre, pi := pi,
                       error := "XML PI cannot be longer than 4094 characters");
        fi;
        pi[len] := char;
    until lastchar = '?' and char = '>';
    Remove(pi);
    Remove(pi);
    return rec(success := true, pre := pre, pi := pi);
end);

max_stars_repo_name: markuspf/MathsInTheMiddle

max_issues_repo_name: markuspf/MathsInTheMiddle

max_forks_repo_name: markuspf/MathsInTheMiddle


text:
# Copyright (c) 2018-2021, Carnegie Mellon University
# See LICENSE for details

##    extractP3(P, k) takes a bit matrix P and a boundary k, and returns
##    submatrix P3 where P is partitioned as
##
##    P4 | P3
##    -------
##    P2 | P1
##
##    If P is of size N by N, P3 is size N-k by k.
##
extractP3 := meth(P, k)
   local dim, row, i, j, res;
   dim := Length(P);
   res := [];
   for i in [1..dim-k] do
      row := [];
      for j in [dim-k+1..dim] do
         Append(row, [P[i][j]]);
      od;
      Append(res, [row]);
   od;
   return res;
end;

##    extractP1(P, k) takes a bit matrix P and a boundary k, and returns
##    submatrix P1 where P is partitioned as
##
##    P4 | P3
##    -------
##    P2 | P1
##
##    P1 is size k by k.
##
extractP1 := meth(P, k)
   local dim, row, i, j, res;
   dim := Length(P);
   res := [];
   for i in [dim-k+1..dim] do
      row := [];
      for j in [dim-k+1..dim] do
         Append(row, [P[i][j]]);
      od;
      Append(res, [row]);
   od;
   return res;
end;

##    NonZeroRows(P) determines which rows of P are not all-zero.  It returns a
##       vector (indexed starting from 1) that indicates these rows.
##
NonZeroRows := meth(P2)
   local i, l, t, res;
   res := [];

   for i in [1..Length(P2)] do
      t := 0;

      for l in [1..Length(P2[i])] do
         if (P2[i][l] <> 0 * GF(2).one) then t := t + 1; fi;
      od;

      if (t <> 0) then Append(res, [i]); fi;

   od;
   return res;
end;

##    ZeroRows(P) determines which rows of P are all zero.  It returns a
##       vector (indexed starting from 1) that indicates these rows.
##
ZeroRows := meth(P1)
   local i, l, t, res;
   res := [];

   for i in [1..Length(P1)] do
      t := 0;

      for l in [1..Length(P1[i])] do
         if (P1[i][l] <> 0 * GF(2).one) then t := t + 1; fi;
      od;

      if (t = 0) then Append(res, [i]); fi;

   od;
   return res;
end;

##   buildH2(i, j, n, k) builds matrix H2 given i, j, n, and k as in
##      Algorithm 5.2, Case 1, Line 3, in the Streaming Permutation
##      paper.
##
buildH2 := meth(i, j, n, k)
   local res, l, c;
   res := [];
   c := 1;
   for l in [1..n-k] do
      if l in j then
         Append(res, [BasisVec(k, i[c]-1) * GF(2).one]);
	 c := c+1;
      else
	 Append(res, [ List([1..k], x -> 0 * GF(2).one) ]);
      fi;
   od;

   return TransposedMat(res);
end;

##   buildH(n, k, H2) constructs matrix H given size n, partition
##      location k, and submatrix H2.
##
##      I(n-k) |  0
##      -------|-----
##        H2   | I(k)
##
buildH := meth(n, k, H2)
   local res, i;
   res := [];
   for i in [1..n-k] do
      Append(res, [BasisVec(n, i-1) * GF(2).one]);
   od;
   for i in [n-k+1..n] do
      Append(res, [Concatenation(H2[i-(n-k)], BasisVec(k, i-(n-k)-1) * GF(2).one )]);
   od;

   return res;
end;

##   GetTopLeftOnes(P, k) counts the number of ones in the top
##      left corner.
##
GetTopLeftOnes := meth(P, k)

#    local count, x, size;

#    size := Length(P);
#    count := 1;

#    repeat x := P[count][count]; count := count+1; until x <> GF(2).one;

#    return Min2(count-2, size-k);

    local size, count, x, y, Pt;

    Pt := TransposedMat(P);
    size := Length(P);
    count := 0;

    repeat
	count := count+1;
        x := P[count];
	y := Pt[count];
    until ((x <> BasisVec(size, count-1) * GF(2).one) or (y <> BasisVec(size, count-1) * GF(2).one) or (count = size));

    return Min2(count-1, size-k);

end;

##   AddTopLeftOnes(P, l) adds l ones to the top left corner.
##
AddTopLeftOnes := meth(P, l)
   local size, res, i, t;
   size := Length(P);
   res := [];

   for i in [1..l] do
      Append(res, [BasisVec(size+l, i-1)]*GF(2).one);
   od;

   t := List([1..l], e->0*GF(2).one);

   for i in [1..size] do
      Append(res, [Concatenation(t, P[i])]);
   od;

   return res;

end;

##   TransposeVector(x) transposes 1d 'horizontal' vector x into
##      a 2d 'vertical' vector.
##
TransposeVector := meth(x)
   local i, res;
   res := [];
   if IsList(x[1]) then
       for i in [1 .. Length(x)] do res[i] := x[i][1]; od;
   else
      for i in [1 .. Length(x)] do res[i] := [x[i]]; od;
   fi;
   return res;
end;

##   CheckPerm(P) returns 1 if P is a permutation matrix,
##      and 0 otherwise.
##
CheckPerm := meth(P)
    local i, j, t, row;
    P := BitMatrixToInts(P);

    for i in [1 .. Length(P)] do
       if ( Sum( P[i]{ [ 1..Length(P[i]) ] } ) <> 1) then
	   return 0;
       fi;
    od;

    return 1;
end;

switchCols := meth(P, a, b)
   local ptrans, row;
   ptrans := TransposedMat(P);
   row := ptrans[a];
   ptrans[a] := ptrans[b];
   ptrans[b] := row;

   return TransposedMat(ptrans);
end;

Declare(BRAMPerm);
Declare(STensor);
Declare(MyTriangulizeMat);
#Declare(BRAMPermStreamOne);
Declare(PermMatrixRowVals);

# ===============================================
# A container for a streaming permutation
# StreamPermBit(<SPL>, <q>)
#   where <SPL> is the permutation as an SPL object, and <q>
#   is the streaming width.
# ===============================================
Class(StreamPermBit, BaseMat, SumsBase, rec(
	abbrevs := [(s, q) -> [s, q]],

	new := (self, prm, ss) >> SPL(WithBases(self, rec(
	     _children := [prm],
	     streamSize := ss,

	     dimensions := prm.dims(),
	     ))),

	rChildren := self >> self._children,
	rSetChild := meth(self, n, what) self._children[n] := what; end,
	child := (self, n) >> self._children[1],
	children := self >> self._children,

	createCode := meth(self)
	   local res, p3, p1, r, k, M, N, P, x, i, j, l, K, size, count, ptrans, z, G, q1, q3, q3b, p1b, G;

       # if (self.streamSize = 1) then
#           return BRAMPermStreamOne(PermMatrixRowVals(MatSPL(self.child(1))),
#                       self.child(1));
#       fi;

	   k := Log2Int(self.streamSize);
	   P := self._children[1]._children[1];
	   size := Log2Int(self.dimensions[1]);

       # If our streaming width is 1, then we are done.  P = I*P.
       #    We no longer do this, because we use the general version above.
	   if (self.streamSize = 1) then
	      #return BRAMPerm(P, MatSPL(I(size))*GF(2).one, 1);
           return BRAMPerm(MatSPL(I(size))*GF(2).one, P, 1);
	   fi;

	   # if P is the bit form of I tensor P',
	   # we can remove the top-left corner bits
	   # (up to (size-k) bits, and replace them
	   # with a streaming tensor product of a
	   # smaller permutation
	   count := GetTopLeftOnes(P,k);
	   P := P{[count+1..Length(P)]}{[count+1..Length(P)]};

	   size := size-count;

	   p3 := extractP3(P, k);
	   p1 := extractP1(P, k);

	   r := Rank(p1);

	   if (r = k) then
	       M := P;
	       N := MatSPL(I(Length(M))) * GF(2).one;
	   else
           if (CheckPerm(P) = 1) then
	           j := NonZeroRows(p3);
	           i := ZeroRows(p1);
	       else
		       p1b := TransposedMat(Copy(p1));
		       G := MyTriangulizeMat(p1b);
		       G := TransposedMat(G);
		       q1 := p1 * G;
		       q3 := p3 * G;

 		       z := LinearIndependentColumns(TransposedMat(q1));
 		       i := [];

 		       for l in [1 .. k] do
   		       if (l in z) then
 		               ;
               else
 		           Add(i, l);
 		       fi;
 		       od;

 		       q3b := [];

 		       for l in [1 .. k] do
 		       Append(q3b, q3{[r+1..k]});
 		       od;

 		       j := LinearIndependentColumns(TransposedMat(q3b));
           fi;

	       K := buildH2(i, j, size, k);

	       N := buildH(size, k, K) * GF(2).one;
	       M := N*P;
	   fi;

	   if (Maximum(Maximum(N * M - P)) <> 0 * GF(2).one) then
	       Print("ERROR: Could not factor permutation matrix.\n");
	   fi;

	   if (count = 0) then
	       return BRAMPerm(N, M, self.streamSize);
	   else
	       return STensor(BRAMPerm(N,M,self.streamSize), 2^count, self.streamSize);
	   fi;
        end,

	print := (self,i,is) >> Print(self.name, "(", self._children[1], ", ", self.streamSize, ")"),

	toAMat := self >> AMatMat(BitsToPermMatrix(self._children[1]._children[1])),
	dims := self >> self._children[1].dims(),

	sums := self >> self,

));

# ===============================================
# A streaming permutation that has been factored into reading and writing
# into a dual-ported RAM
# BRAMPerm(M, N)
#   where M and N are bit representations of the factorized perm
#
# (see Pueschel, Milder, and Hoe.  "Permuting Streaming Data Using RAMs.")
# ===============================================
Class(BRAMPerm, BaseMat, SumsBase, rec(
	abbrevs := [(m,n,q)-> [m,n,q]],

	new := (self, m, n,q) >> SPL(WithBases(self, rec(
		    _children  := [m,n],
		    dimensions := [2^Length(m), 2^Length(m)],
		    streamSize := q
	     ))),

	rChildren := self >> self._children,
	rSetChild := meth(self, n, what) self._children[n] := what; end,
	child     := (self, n) >> self._children[n],
	children  := self >> self._children,

	createCode := self >> self,
        print      := (self,i,is) >> Print(self.name, "(", BitMatrixToInts(self._children[1]), ", ", BitMatrixToInts(self._children[2]), ", ", self.streamSize, ")"),

	toAMat := self >> AMatMat(BitsToPermMatrix(self.child(1)) * BitsToPermMatrix(self.child(2))),
	dims   := self >> [2^(Length(self._children[1])), 2^(Length(self._children[1]))],

	sums   := self >> self
));

## PermMatrixRowVals(p) takes permutation matrix p and returns a list
##   of the '1' location in each row of the matrix.
PermMatrixRowVals := meth(x)
   local length, i, res;

   res := [];
   length := Length(x);
   for i in [1 .. length] do
      res[i] := BasisVecToNum(x[i]);
   od;
   return res;
end;

#############################################################################
##
#F  MyTriangulizeMat( <mat> ) . . . . . bring a matrix in upper triangular form
##
##    This is a modified version of TriangulizeMat included with Gap.
##       I have modified it to keep track of the elementary row
##       operations performed, and return their product.
MyTriangulizeMat := function ( mat )
    local   m, n, i, j, k, row, row2, zero, G;

    if mat <> [] then

       # get the size of the matrix
       m := Length(mat);
       n := Length(mat[1]);
       zero := 0*mat[1][1];
       G := MatSPL(I(Length(mat)));

       # run through all columns of the matrix
       i := 0;
       for k  in [1..n]  do

           # find a nonzero entry in this column
           j := i + 1;
           while j <= m and mat[j][k] = zero  do j := j + 1;  od;

           # if there is a nonzero entry
           if j <= m  then

               # increment the rank
               InfoMatrix2(k," \c");
               i := i + 1;

               # make its row the current row and normalize it
               row := mat[j];  mat[j] := mat[i];  mat[i] := row[k]^-1 * row;
	       row2 := G[j]; G[j] := G[i]; G[i] := row[k]^-1 * row2;
	       # Print("aG = ", G, "\n");
	       # Print("amat = ", mat, "\n");

               # clear all entries in this column
               for j  in [1..m] do
                   if  i <> j  and mat[j][k] <> zero  then
		       G[j] := G[j] - mat[j][k] * G[i];
                       mat[j] := mat[j] - mat[j][k] * mat[i];
		       # Print("bG = ", G, "\n");
		       # Print("G[i] = ", G[i], "\n");
		       # Print("mat[j][k] = ", mat[j][k], "\n");
		       # Print("bmat = ", mat, "\n");
                   fi;
               od;

           fi;

       od;

    fi;

    return G;

end;

ext: gi


text:
#
# cls 15: A_2
#
handle15char2:=function()
    local tmp1,tmp2,o,info,rep,cent,u;
    o:=AllClasses(orbs)[15];
    info:=infos[15];

    rep:=Representative(o);
    #pr[7]=[0,0,1,1]-->[0,1,1,0]
    #pr[13]=[0,1,2,1]-->[1,1,2,0]
    tmp1:=ApplyRootsReflections(rep,[7,13]);
    u:=Unipotent(algebraicU(orbs),[[4,1]]);
    tmp1:=Conj(tmp1,u);

    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[7,13]);
    tmp2:=Conj(tmp2,u);

    return [tmp1,tmp2];
end;

#
# cls 14: \tilde A_2
#
handle14char2:=function()
    local tmp1,tmp2,o,info,rep,cent,u;
    o:=AllClasses(orbs)[14];
    info:=infos[14];

    rep:=Representative(o);
    #pr[8]=[1,0,1,1]-->[0,1,1,1]
    #pr[10]=[0,0,2,1]-->[0,1,2,0]
    # merge si [13,11],[19,20]
    tmp1:=ApplyRootsReflections(rep,[8,10]);
    u:=Unipotent(chevalleyAdj(rep),[[3,1]],Ordering(rep));
    tmp1:=Conj(tmp1,u);

    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[8,10]);
    u:=Unipotent(chevalleyAdj(rep),[[3,1]],Ordering(cent));
    tmp2:=Conj(tmp2,u);

    return [tmp1,tmp2];
end;

#
# cls 12: B_2
#
handle12char2:=function()
    local tmp1,tmp2,o,info,rep,cent;
    o:=AllClasses(orbs)[12];
    info:=infos[12];

    rep:=Representative(o);
    tmp1:=ApplyRootsReflections(rep,[9]);#pr[9]=[0,1,1,1]-->[1,1,1,0]
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[9]);#pr[9]=[0,1,1,1]-->[1,1,1,0]
    return [tmp1,tmp2];
end;

#
# cls 11: B_2^{(2)} [6,12,22]
#
# here there is no component group
#
handle11char2:=function()
    local rep,cc,zz,uu;

    rep:=Unipotent(alg,[[6,1],[12,1],[22,1]]);
    cc:=UnipotentCentralizer(rep);
    zz:=Zentrum(cc[1]);
    uu:=Unipotent(alg,[[1,1]],Ordering(rep));# A_1 the reductive part of the Levi

    return [rep,cc,zz,uu,TorusWeights(zz,[1])];
end;

#
# cls 8: C_3(a_1)
#
# here there is no component group
#
handle8char2:=function()
    local o,rep,cc,ccc,zz,zzz,uu,dist_par;

    o:=AllClasses(orbs)[8];
    rep:=BorelRep(o);
    vals:=ShallowCopy(avars);
    vals[1]:=Zero(APR);
    vals[2]:=Zero(APR);
    vals[5]:=Zero(APR);
    vals[6]:=Zero(APR);
    vals[8]:=Zero(APR);
    vals[12]:=vals[11];
    vals[15]:=vals[9];
    vals[18]:=vals[13];
    vals[20]:=vals[16];
    cc:=Value(generic,avars,vals);
    # a_9+a_9^2+a_11^2=0
    #ccc:=Value(cc,[avars[9]],[Zero(APR)]); NU: a_9=0 <=> a_11=0

    zz:=Zentrum(infos[8][1]); # centrul centralizatorului calculat inainte
    dist_par:=ToPositiveBorel(o,0Level(DistinguishedParabolic(o)));
    zzz:=Subword(zz,dist_par);
#    uu:=Unipotent(alg,[[1,1]],Ordering(rep));# A_1 the reductive part of the Levi

    return [cc,zz,zzz];
end;

##
## cls 7: C_3(a_1)
##
#handle7char3:=function()
#    local tmp1,tmp2,o,info,rep,cent;
#    o:=AllClasses(orbs)[7];
#    info:=infos[7];
#
#    rep:=Representative(o);
#    tmp1:=ConjugateByTorus(rep,1,-1);# h_1 --> h_4
#    cent:=FromPositiveBorel(o,info[5]);
#    tmp2:=ConjugateByTorus(cent,1,-1);# h_1 --> h_4
#    return [tmp1,tmp2];
#end;

#
# cls 6: F_4(a_3) [4,6,10,18] (finished?)
#
handle6char2:=function()
    local tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,o,info,rep,cent,uu;
    o:=AllClasses(orbs)[6];
    info:=infos[6];

    rep:=Representative(o);
    tmp1:=ApplyRootsReflections(rep,[3]);#pr[3]=[0,0,1,0]-->[0,0,1,0]
    uu:=Unipotent(chevalleyAdj(rep),[[2,-1],[6,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);
    #tmp1:=ConjugateByTori(tmp1,[1,2,3,4],[1,-1,aaa,aaa^2]);

    tmp3:=ApplyRootsReflections(rep,[5,1]);
    uu:=Unipotent(chevalleyAdj(rep),[[2,1]],Ordering(rep));#[1..24]);
    tmp3:=Conj(tmp3,uu);
    tmp3:=ConjNegUa(tmp3,2,1);
    uu:=Unipotent(chevalleyAdj(rep),[[6,1],[4,1]],Ordering(rep));
    tmp3:=Conj(tmp3,uu);

    #tmp5:=ConjugateByTori(rep,[3],[-1]);

    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[3]);#pr[3]=[0,0,1,0]-->[0,0,1,0]
    uu:=Unipotent(chevalleyAdj(cent),[[2,-1],[6,1]]);
    tmp2:=Conj(tmp2,uu);
    #tmp2:=ConjugateByTori(tmp2,[1,2,3,4],[1,-1,aaa,aaa^2]);

    tmp4:=ApplyRootsReflections(cent,[5,1]);
    uu:=Unipotent(chevalleyAdj(rep),[[2,1]],Ordering(cent));#[1..24]);
    tmp4:=Conj(tmp4,uu);
    tmp4:=ConjNegUa(tmp4,2,1);
    uu:=Unipotent(chevalleyAdj(rep),[[6,1],[4,1]],Ordering(cent));
    tmp4:=Conj(tmp4,uu);
    #tmp4:=ConjugateByTori(tmp4,[3],[aaa]);

    #tmp6:=ConjugateByTori(cent,[3],[-1]);

    return [tmp1,tmp3,tmp2,tmp4];
end;

#
# cls 3: F_4(a_2) (finished?)
#
handle3char2:=function()
    local tmp1,tmp2,o,info,rep,cent,uu,
    avars,vals,comp1,comp2;
    o:=AllClasses(orbs)[3];
    info:=infos[3];

    rep:=Representative(o);
    tmp1:=ApplyRootsReflections(rep,[2,3]);
    uu:=Unipotent(chevalleyAdj(rep),[[1,1],[10,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);

    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[2,3]);
    uu:=Unipotent(chevalleyAdj(cent),[[1,1],[10,1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);

    avars:=IndeterminatesOfPolynomialRing(ring(chevalleyAdj(rep)));
    vals:=ShallowCopy(avars);
    vals[1]:=One(APR);
    vals[2]:=Zero(APR);
    vals[3]:=Zero(APR);
    vals[4]:=vals[7];
    vals[6]:=vals[7]+One(APR);
    vals[9]:=vals[7];
    vals[5]:=vals[6]+vals[9];
    vals[10]:=One(APR)+vals[7];
    vals[11]:=vals[16]+vals[7]+vals[7]^2;
    vals[12]:=vals[7]*vals[9]+vals[7]+vals[16];
    vals[13]:=vals[9];
    vals[15]:=vals[16]+vals[7]^2;
    vals[17]:=vals[19]+vals[16]+vals[7]+vals[7]^2;
    vals[18]:=vals[16]+vals[7]+vals[7]^2;
    vals[20]:=vals[16]^2+vals[7]+vals[7]^4;
    vals[22]:=vals[7]^4+vals[7]^2+vals[16]^2+vals[7]+vals[16];
    comp1:=Value(generic,avars,vals);

    vals:=ShallowCopy(avars);
    vals[1]:=Zero(APR);
    vals[2]:=Zero(APR);
    vals[3]:=Zero(APR);
    vals[7]:=vals[4];
    vals[6]:=vals[7];
    vals[9]:=vals[7];
    vals[5]:=vals[6]+vals[9];
    vals[10]:=vals[7];
    vals[13]:=vals[9];
    vals[12]:=vals[7]*vals[9]+vals[16];
    vals[11]:=vals[16]+vals[7]^2;
    vals[15]:=vals[16]+vals[7]^2;
    vals[17]:=vals[19]+vals[16]+vals[7]^2;
    vals[18]:=vals[16]+vals[7]^2;
    vals[20]:=vals[16]^2+vals[7]^4;
    vals[22]:=vals[7]^4+vals[7]^2+vals[16]^2+vals[16];
    comp2:=Value(generic,avars,vals);

    return [comp1,comp2,tmp1,tmp2];
end;
#vals:=ShallowCopy(avars);
#vals[1]:=One(APR);
#vals[2]:=Zero(APR);
#vals[3]:=Zero(APR);
#vals[4]:=vals[7];
#vals[6]:=vals[7]+One(APR);
#vals[9]:=vals[7];
#vals[5]:=vals[6]+vals[9];
#vals[10]:=One(APR)+vals[7];
#vals[11]:=vals[16]+vals[7]+vals[7]^2;
#vals[12]:=vals[7]*vals[9]+vals[7]+vals[16];
#vals[13]:=vals[9];
#vals[15]:=vals[16]+vals[7]^2;
#vals[17]:=vals[19]+vals[16]+vals[7]+vals[7]^2;
#vals[18]:=vals[16]+vals[7]+vals[7]^2;
#vals[20]:=vals[16]^2+vals[7]+vals[7]^4;
#vals[22]:=vals[7]^4+vals[7]^2+vals[16]^2+vals[7]+vals[16];
#
#xx:=Value(generic,avars,vals);
#vals:=List([1..Length(avars)],i->Zero(APR));
#vals[7]:=One(APR);
#xx0:=Value(xx,avars,vals);
#xxx:=Value(cc[1],[avars[1],avars[2],avars[3],avars[4],avars[5]],[Zero(APR),Zero(APR),Zero(APR),avars[7]^2+avars[8]*avars[7]+a#vars[7]*avars[14]+One(APR),One(APR)]);
#
#vals:=ShallowCopy(avars);
#vals[1]:=Zero(APR);
#vals[2]:=Zero(APR);
#vals[3]:=Zero(APR);
#vals[7]:=vals[4];
#vals[6]:=vals[7];
#vals[9]:=vals[7];
#vals[5]:=vals[6]+vals[9];
#vals[10]:=vals[7];
#vals[13]:=vals[9];
#vals[12]:=vals[7]*vals[9]+vals[16];
#vals[11]:=vals[16]+vals[7]^2;
#vals[15]:=vals[16]+vals[7]^2;
#vals[17]:=vals[19]+vals[16]+vals[7]^2;
#vals[18]:=vals[16]+vals[7]^2;
#vals[20]:=vals[16]^2+vals[7]^4;
#vals[22]:=vals[7]^4+vals[7]^2+vals[16]^2+vals[16];
##Print("ok");
#yy:=Value(generic,avars,vals);
#vals[8]:=avars[1];#
#vals[16]:=avars[2];#
#vals[14]:=avars[3];#
#vals[19]:=avars[5];#
#vals[21]:=avars[6];#
#vals[23]:=avars[9];#
#vals[24]:=avars[10];#
##vals:=ShallowCopy(avars);
#yy0:=Value(yy,avars,vals);
#yyy:=Value(cc[1],[avars[1],avars[2],avars[3],avars[4],avars[5]],[Zero(APR),Zero(APR),Zero(APR),avars[8]*avars[7]+avars[7]*avars[14],Zero(APR)]);
#
#vals:=List([1..Length(avars)],i->Zero(APR));
#vals[4]:=One(APR);
#vals[16]:=One(APR);
#yy1:=Value(yy,avars,vals);

#
# cls 2: F_4(a_1)
#
handle2char2:=function()
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=AllClasses(orbs)[2];
    info:=infos[2];

    rep:=Representative(o);
    tmp1:=ConjugateByTori(rep,[3,1],[-1,-1]);# h_1 --> h_4, h_3 --> h_3
    uu:=Unipotent(chevalleyAdj(rep),[[3,-1]],Ordering(rep));#pr[3]=[0,0,1,0]-->[0,0,1,0]
    tmp1:=Conj(tmp1,uu);

    cent:=FromPositiveBorel(o,info[5]);
    tmp2:=ConjugateByTori(cent,[3,1],[-1,-1]);# h_1 --> h_4, h_3 --> h_3
    uu:=Unipotent(chevalleyAdj(cent),[[3,-1]],Ordering(cent));#pr[3]=[0,0,1,0]-->[0,0,1,0]
    tmp2:=Conj(tmp2,uu);

    return [tmp1,tmp2];
end;


text:
#############################################################################
##
#W  automsg.gi               automgrp package                  Yevgen Muntyan
#W                                                             Dmytro Savchuk
##
#Y  Copyright (C) 2003 - 2018 Yevgen Muntyan, Dmytro Savchuk
##

###############################################################################
##
#M  AutomatonSemigroup(<list>)
##
InstallMethod(AutomatonSemigroup, "for [IsList]", [IsList],
function (list)
  return AutomatonSemigroup(list, false);
end);

###############################################################################
##
#M  AutomatonSemigroup(<list>, <bind_vars>)
##
InstallMethod(AutomatonSemigroup, "for [IsList, IsBool]", [IsList, IsBool],
function (list, bind_vars)
  if not AG_IsCorrectAutomatonList(list, false) then
    Error("in AutomatonSemigroup(IsList):\n",
          "  given list is not a correct list representing automaton\n");
  fi;

  # XXX
  return SemigroupOfAutomFamily(AutomFamily(list, bind_vars));
end);

###############################################################################
##
#M  AutomatonSemigroup(<list>, <names>)
##
InstallMethod(AutomatonSemigroup, "for [IsList, IsList]", [IsList, IsList],
function (list, names)
  if not AG_IsCorrectAutomatonList(list, false) then
    Error("error in AutomatonSemigroup(IsList, IsList):\n",
          "  given list is not a correct list representing automaton\n");
  fi;

  # XXX
  return SemigroupOfAutomFamily(AutomFamily(list, names));
end);

###############################################################################
##
#M  AutomatonSemigroup(<list>, <names>, <bind_vars>)
##
InstallMethod(AutomatonSemigroup, "for [IsList, IsList, IsBool]",
              [IsList, IsList, IsBool],
function (list, names, bind_vars)
  if not AG_IsCorrectAutomatonList(list, false) then
    Error("error in AutomatonSemigroup(IsList):\n",
          "  given list is not a correct list representing automaton\n");
  fi;

  #XXX
  return SemigroupOfAutomFamily(AutomFamily(list, names, bind_vars));
end);

###############################################################################
##
#M  AutomatonSemigroup(<string>)
#M  AutomatonSemigroup(<string>, <bind_vars>)
##
InstallMethod(AutomatonSemigroup, "for [IsString]", [IsString],
function(string)
    return AutomatonSemigroup(string, AG_Globals.bind_vars_autom_family);
end);
InstallMethod(AutomatonSemigroup, "AutomatonSemigroup(IsString, IsBool]", [IsString, IsBool],
function(string, bind_vars)
    local s;
    s := AG_ParseAutomatonString(string);
    return AutomatonSemigroup(s[2], s[1], bind_vars);
end);

###############################################################################
##
#M  AutomatonSemigroup(<A>)
#M  AutomatonSemigroup(<A>, <bind_vars>)
##
InstallMethod(AutomatonSemigroup, "for [IsMealyAutomaton]", [IsMealyAutomaton],
function(A)
  return AutomatonSemigroup(AutomatonList(A), A!.states);
end);

InstallMethod(AutomatonSemigroup, "for [IsMealyAutomaton, IsBool]", [IsMealyAutomaton, IsBool],
function(A, bind_vars)
  return AutomatonSemigroup(AutomatonList(A), A!.states, bind_vars);
end);

###############################################################################
##
#M  IsSelfSimilar(<G>)
##
InstallMethod(IsSelfSimilar, "for [IsAutomSemigroup]",
              [IsAutomSemigroup],
function(G)
  local g, i, res;
  res := true;
  for g in GeneratorsOfSemigroup(G) do
    for i in [1..UnderlyingAutomFamily(G)!.deg] do
      res := Section(g, i) in G;
      if res = fail then
        TryNextMethod();
      elif not res then
        return false;
      fi;
    od;
  od;
  return true;
end);

###############################################################################
##
#M  UnderlyingAutomFamily(<G>)
##
InstallMethod(UnderlyingAutomFamily, "for [IsAutomSemigroup]",
              [IsAutomSemigroup],
function(G)
  return FamilyObj(GeneratorsOfSemigroup(G)[1]);
end);

# ###############################################################################
# ##
# #M  UseSubsetRelation(<G>)
# ##
# InstallMethod(UseSubsetRelation,
#               "for [IsAutomSemigroup, IsAutomSemigroup]",
#               [IsAutomSemigroup, IsAutomSemigroup],
# function(super, sub)
#   ## the full group is self similar, so if <super> is smaller than the full
#   ##  group then sub is smaller either
#   if HasIsGroupOfAutomFamily(super) then
#     if not IsGroupOfAutomFamily(super) then
#       SetIsGroupOfAutomFamily(sub, false); fi; fi;
#   TryNextMethod();
# end);

# ###############################################################################
# ##
# #M  __AG_SubgroupOnLevel(<G>, <gens>, <level>)
# ##
# InstallMethod(__AG_SubgroupOnLevel, [IsAutomGroup,
#                                  IsList and IsTreeAutomorphismCollection,
#                                  IsPosInt],
# function(G, gens, level)
#   local overgroup;
#
#   if IsEmpty(gens) or (Length(gens) = 1 and IsOne(gens[1])) then
#     return TrivialSubgroup(G);
#   fi;
#
#   if HasIsGroupOfAutomFamily(G) and IsGroupOfAutomFamily(G) then
#     overgroup := G;
#   else
#     overgroup := GroupOfAutomFamily(UnderlyingAutomFamily(G));
#   fi;
#
#   return SubgroupNC(overgroup, gens);
# end);
#
# InstallOtherMethod(__AG_SubgroupOnLevel, [IsAutomGroup, IsList and IsEmpty, IsPosInt],
# function(G, gens, level)
#   return TrivialSubgroup(G);
# end);
#
# InstallMethod(__AG_SubgroupOnLevel, [IsTreeAutomorphismGroup,
#                                  IsList and IsAutomCollection,
#                                  IsPosInt],
# function(G, gens, level)
#   local overgroup;
#
#   overgroup := GroupOfAutomFamily(FamilyObj(gens[1]));
#
#   if Length(gens) = 1 and IsOne(gens[1]) then
#     return TrivialSubgroup(overgroup);
#   fi;
#
#   return SubgroupNC(overgroup, gens);
# end);
#
# InstallMethod(__AG_SimplifyGroupGenerators, [IsList and IsAutomCollection],
# function(gens)
#   local words, fam;
#
#   if IsEmpty(gens) then
#     return [];
#   fi;
#
#   fam := FamilyObj(gens[1]);
#   words := FreeGeneratorsOfGroup(Group(List(gens, a -> a!.word)));
#
#   if fam!.use_rws and not IsEmpty(words) then
#     words := AG_ReducedForm(fam!.rws, words);
#     words := FreeGeneratorsOfGroup(Group(words));
#   fi;
#
#   return List(words, w -> Autom(w, fam));
# end);

###############################################################################
##
#M  DegreeOfTree(<G>)
##
InstallMethod(DegreeOfTree, "for [IsAutomSemigroup]",
              [IsAutomSemigroup],
function(G)
  return DegreeOfTree(UnderlyingAutomFamily(G));
end);

InstallMethod(TopDegreeOfTree, "for [IsAutomSemigroup]",
              [IsAutomSemigroup],
function(G)
  return DegreeOfTree(UnderlyingAutomFamily(G));
end);

###############################################################################
##
#M  PrintObj(<G>)
##
InstallMethod(PrintObj, "for [IsAutomatonSemigroup]",
              [IsAutomatonSemigroup],
function(G)
  Print("AutomatonSemigroup(\"", String(G), "\")");
end);

#############################################################################
##
#M  String(<G>)
##
InstallMethod(String, "for [IsAutomSemigroup]", [IsAutomSemigroup],
function(G)
  local i, gens, formatone, s;

  formatone := function(a)
    return Concatenation(String(a), " = ", String(Decompose(a)));
  end;

  if IsMonoid(G) then
    gens := GeneratorsOfMonoid(G);
  else
    gens := GeneratorsOfSemigroup(G);
  fi;

  s := "";
  for i in [1..Length(gens)] do
    Append(s, formatone(gens[i]));
    if i <> Length(gens) then
      Append(s, ", ");
    fi;
  od;

  return s;
end);

###############################################################################
##
#M  Display(<G>)
##
InstallMethod(Display, "for [IsAutomSemigroup]",
              [IsAutomSemigroup],
function(G)
  local i, gens, printone;

  printone := function(a)
    Print(a, " = ", Decompose(a));
  end;

  gens := GeneratorsOfSemigroup(G);
  if gens = [] then Print("< >"); fi;
  if Length(gens) = 1 then
    Print("< "); printone(gens[1]); Print(" >");
  else
    Print("< "); printone(gens[1]); Print(", \n");
    for i in [2..Length(gens)-1] do
      Print("  "); printone(gens[i]); Print(", \n");
    od;
    Print("  "); printone(gens[Length(gens)]); Print(" >");
  fi;
end);

###############################################################################
##
#M  ViewObj(<G>)
##
InstallMethod(ViewObj, "for [IsAutomSemigroup]",
              [IsAutomSemigroup],
function(G)
  local i, gens;
  gens := List(GeneratorsOfSemigroup(G), g -> Word(g));
  if gens = [] then Print("< >"); fi;
  Print("< ");
  for i in [1..Length(gens)-1] do
    if IsOne(gens[i]) then
      Print(AG_Globals.identity_symbol, ", ");
    else
      Print(gens[i], ", ");
    fi;
  od;
  if IsOne(gens[Length(gens)]) then
    Print(AG_Globals.identity_symbol, " >");
  else
    Print(gens[Length(gens)], " >");
  fi;
end);

###############################################################################
##
#M  IsTrivial(G)
##
InstallMethod(IsTrivial, "for [IsAutomSemigroup]", [IsAutomSemigroup],
function (G)
  local g;
  for g in GeneratorsOfSemigroup(G) do
    if not IsOne(g) then return false; fi;
  od;
  return true;
end);

###############################################################################
##
#M  Size(G)
##
InstallMethod(Size, "for [IsAutomSemigroup]", [IsAutomSemigroup],
function (G)
  local g;
  if IsTrivial(G) then
    Info(InfoAutomGrp, 3, "Size(G): 1, G is trivial");
    return 1;
  fi;

  for g in Iterator(G) do od;

  return Size(G);
#  TryNextMethod();
end);

# ###############################################################################
# ##
# #M  \= (<G>, <H>)
# ##
# InstallMethod(\=, "for [IsAutomGroup, IsAutomGroup]",
#               IsIdenticalObj, [IsAutomGroup, IsAutomGroup],
# function(G, H)
#   local fgens1, fgens2, fam;
#
#   if HasIsGroupOfAutomFamily(G) and HasIsGroupOfAutomFamily(H) then
#     if IsGroupOfAutomFamily(G) <> IsGroupOfAutomFamily(H) then
#       Info(InfoAutomGrp, 3, "G = H: false, exactly one is GroupOfAutomFamily");
#       return false;
#     fi;
#     if IsGroupOfAutomFamily(G) then
#       Info(InfoAutomGrp, 3, "G = H: true, both are GroupOfAutomFamily");
#       return true;
#     fi;
#   fi;
#
#   fgens1 := List(GeneratorsOfGroup(G), g -> Word(g));
#   fgens2 := List(GeneratorsOfGroup(H), g -> Word(g));
#   fam := UnderlyingAutomFamily(G);
#
#   if fam!.rws <> fail then
#     fgens1 := AsSet(AG_ReducedForm(fam!.rws, fgens1));
#     fgens2 := AsSet(AG_ReducedForm(fam!.rws, fgens2));
#   fi;
#
#   if GroupWithGenerators(fgens1) = GroupWithGenerators(fgens2) then
#     Info(InfoAutomGrp, 3, "G = H: true, by subgroups of free group");
#     return true;
#   fi;
#
#   TryNextMethod();
# end);

# ###############################################################################
# ##
# #M  IsSubset (<G>, <H>)
# ##
# InstallMethod(IsSubset, "for [IsAutomGroup, IsAutomGroup]",
#               IsIdenticalObj, [IsAutomGroup, IsAutomGroup],
# function(G, H)
#   local h, fam, fgens1, fgens2;
#
#   if HasIsGroupOfAutomFamily(G) and IsGroupOfAutomFamily(G) then
#     Info(InfoAutomGrp, 3, "IsSubgroup(G, H): true");
#     Info(InfoAutomGrp, 3, "  G is GroupOfAutomFamily");
#     return true;
#   fi;
#
#   fgens1 := List(GeneratorsOfGroup(G), g -> Word(g));
#   fgens2 := List(GeneratorsOfGroup(H), g -> Word(g));
#   fam := UnderlyingAutomFamily(G);
#
#   if fam!.rws <> fail then
#     fgens1 := AsSet(AG_ReducedForm(fam!.rws, fgens1));
#     fgens2 := AsSet(AG_ReducedForm(fam!.rws, fgens2));
#   fi;
#
#   if IsSubgroup(GroupWithGenerators(fgens1), GroupWithGenerators(fgens2)) then
#     Info(InfoAutomGrp, 3, "IsSubgroup(G, H): true");
#     Info(InfoAutomGrp, 3, "  by subgroups of free group");
#     return true;
#   fi;
#
#   TryNextMethod();
# end);

###############################################################################
##
#M  <g> in <G>
##
InstallMethod(\in, "for [IsAutom, IsAutomGroup]",
              [IsAutom, IsAutomSemigroup],
function(g, G)
  local fam, fgens, w;

  if HasIsGroupOfAutomFamily(G) and IsGroupOfAutomFamily(G) then
    return true;
  fi;

  fgens := List(GeneratorsOfSemigroup(G), g -> Word(g));
  w := Word(g);

  fam := UnderlyingAutomFamily(G);

  if fam!.rws <> fail then
    fgens := AsSet(AG_ReducedForm(fam!.rws, fgens));
    w := AG_ReducedForm(fam!.rws, w);
  fi;

  if w in SemigroupByGenerators(fgens) then
    Info(InfoAutomGrp, 3, "g in G: true");
    Info(InfoAutomGrp, 3, "  by elements of free group");
    Info(InfoAutomGrp, 3, "  g = ", g, "; G = ", G);
    return true;
  fi;

  TryNextMethod();
end);

###############################################################################
##
#M  Random(<G>)
##
InstallMethodWithRandomSource(Random, "for a random source and [IsAutomSemigroup]",
              [IsRandomSource, IsAutomSemigroup],
function(rs, G)
  local w, monoid, F, gens, pi;

  if IsAutomatonSemigroup(G) then
    monoid := UnderlyingFreeMonoid(G);

    if IsTrivial(monoid) then
      w := One(monoid);
    else
      while true do
        w := Random(rs, monoid);
        if not IsOne(w) then
          break;
        fi;
      od;
    fi;
    return Autom(w, UnderlyingAutomFamily(G));
  else
    gens := GeneratorsOfSemigroup(G);
    F := FreeGroup(Length(gens));
    pi := GroupHomomorphismByImagesNC(F,                      UnderlyingFreeGroup(G),
                                      GeneratorsOfGroup(F),   List(gens, Word)        );
    return Autom( Random( rs, SemigroupByGenerators( GeneratorsOfGroup(F)))^pi, UnderlyingAutomFamily(G));
  fi;
end);

###############################################################################
##
#M  UnderlyingFreeMonoid( <G> )
##
InstallMethod(UnderlyingFreeMonoid, "for [IsAutomSemigroup]",
              [IsAutomSemigroup],
function(G)
  return UnderlyingFreeMonoid(UnderlyingAutomFamily(G));
end);

###############################################################################
##
#M  UnderlyingFreeGroup( <G> )
##
InstallMethod(UnderlyingFreeGroup, "for [IsAutomSemigroup]",
              [IsAutomSemigroup],
function(G)
  return UnderlyingFreeGroup(UnderlyingAutomFamily(G));
end);

###############################################################################
##
#M  UnderlyingFreeGenerators( <G> )
##
InstallMethod(UnderlyingFreeGenerators, "for [IsAutomSemigroup]",
              [IsAutomSemigroup],
function(G)
  return List(GeneratorsOfSemigroup(G), g -> Word(g));
end);

# ###############################################################################
# ##
# #M  UnderlyingFreeSubgroup(<G>)
# ##
# InstallMethod(UnderlyingFreeSubgroup, "for [IsAutomGroup]",
#               [IsAutomGroup],
# function(G)
#   local f;
#   if HasIsGroupOfAutomFamily(G) and IsGroupOfAutomFamily(G) then
#     return UnderlyingFreeGroup(G);
#   fi;
#   f := Subgroup(UnderlyingFreeGroup(G), UnderlyingFreeGenerators(G));
#   if f = UnderlyingFreeGroup(G) then
#     SetIsGroupOfAutomFamily(G, true);
#   fi;
#   return f;
# end);

###############################################################################
##
#M  UnderlyingFreeGroup( <G> )
##
InstallMethod(UnderlyingFreeGroup, "for [IsAutomSemigroup]",
              [IsAutomSemigroup],
function(G)
  return UnderlyingAutomFamily(G)!.freegroup;
end);

###############################################################################
##
#M  UnderlyingFreeGenerators( <G> )
##
InstallMethod(UnderlyingFreeGenerators, "for [IsAutomSemigroup]",
              [IsAutomSemigroup],
function(G)
  return List(GeneratorsOfSemigroup(G), g -> Word(g));
end);

InstallMethod(SphericalIndex, "for IsAutomSemigroup",
              [IsAutomSemigroup],
function(G)
  return SphericalIndex(GeneratorsOfSemigroup(G)[1]);
end);
InstallMethod(DegreeOfTree, "for IsAutomSemigroup",
              [IsAutomSemigroup],
function(G)
  return UnderlyingAutomFamily(G)!.deg;
end);
InstallMethod(TopDegreeOfTree, "for IsAutomSemigroup",
              [IsAutomSemigroup],
function(G)
  return UnderlyingAutomFamily(G)!.deg;
end);

###############################################################################
##
#M  UnderlyingAutomaton(<G>)
##
InstallMethod(UnderlyingAutomaton, "for [IsAutomSemigroup]",
              [IsAutomSemigroup],
function(G)
  local fam, numstates, automatonlist, d, i, j;
  fam := UnderlyingAutomFamily(G);
  numstates := fam!.numstates;
# if you do the next 2 operations in one "List", it will remove unbinded spaces
  automatonlist := List(fam!.automatonlist);
  Apply(automatonlist, x -> ShallowCopy(x));
  d := fam!.deg;

# in case we have 1 in the list we move it to the numstates+1 postion
  if Length(automatonlist)=2*numstates+1 then
    for i in [1..numstates] do
      for j in [1..d] do
        if automatonlist[i][j]=2*numstates+1 then automatonlist[i][j] := numstates+1; fi;
      od;
    od;
    automatonlist[numstates+1] := List([1..d], x -> numstates+1);
    Add(automatonlist[numstates+1], automatonlist[2*numstates+1][d+1]);
    numstates := numstates+1;
  fi;
  return MealyAutomaton(automatonlist{[1..numstates]});
end);

###############################################################################
##
#M  IsAutomatonSemigroup(<G>)
##
InstallMethod(IsAutomatonSemigroup, "for [IsAutomSemigroup]",
              [IsAutomSemigroup],
function(G)
  if not HasIsAutomatonSemigroup(G) then return false; fi;
end);

###############################################################################
##
#M  SemigroupOfAutomFamily(<G>)
##
InstallMethod(SemigroupOfAutomFamily, "for [IsAutomSemigroup]",
                   [IsAutomSemigroup],
function(G)
  return SemigroupOfAutomFamily(UnderlyingAutomFamily(G));
end);

#E

max_stars_repo_name: gap-packages/automgrp

max_issues_repo_name: gap-packages/automgrp

max_forks_repo_name: gap-packages/automgrp


text:
#############################################################################
##
#W  union.gi                                                    Karel Dekimpe
#W                                                               Bettina Eick
##

#############################################################################
##
#F GenericDeterminantMat( mat )
##
GenericDeterminantMat := function( mat )
    local d, det, sig, i, sub;

    # set up
    d := Length( mat );
    if ForAny( mat, x -> Length(x) <> d ) then return fail; fi;

    # the trivial cases
    if d = 1 then return mat[1][1]; fi;
    if d = 2 then return mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1]; fi;

    # otherwise use first row and recursion
    det := 0;
    sig := 1;
    for i in [1..d] do
        sub := Concatenation( [1..i-1], [i+1..d] );
        sub := mat{[2..d]}{sub};
        det := det + sig * mat[1][i] * GenericDeterminantMat( sub );
        sig := - sig;
    od;

    return det;
end;

#############################################################################
##
#F NullspaceIntMod( base, vec, p )  . . . . . . . . . . . . b * vec = 0 mod p
##
## computes those elements b in <base> with b * vec = 0 mod p. Returns a
## triangulized basis with respect to <base>.
##
NullspaceIntMod := function( base, vec, p )
    local imgs, d, null;

    # get images
    imgs := List( base, x -> x * vec );
    imgs := List( imgs, x -> x mod p );
    d    := Length( imgs );
    if ForAll( imgs, x -> x = 0 ) then return IdentityMat( d ); fi;

    # compute kernel of imgs vector - must have rank d-1
    Add( imgs, p );
    null := NullspaceIntMat( TransposedMat( [imgs] ) );
    null := List( null, x -> x{[1..d]} );
    null := NormalFormIntMat( null, 0 ).normal;

    # return this images nullspace
    return Filtered( null, x -> PositionNonZero(x) <= d );
end;

#############################################################################
##
#F FindMaximals( sub ) . . . . . . . .subgroups which are maximal within sub
##
FindMaximals := function( sub )
    local new, i, tmp;
    new := [];
    for i in [1..Length(sub)] do
        if Size( sub[i] ) > 1 then
            if Length(new) = 0 then
                Add( new, sub[i] );
            elif not ForAny( new, x -> IsSubset( x, sub[i] ) ) then
                tmp := Filtered( new, x -> not IsSubset( sub[i], x ) );
                Add( new, sub[i] );
            fi;
        fi;
    od;
    return new;
end;

if not IsBound( SizeOfUnion ) then SizeOfUnion := false; fi;
#############################################################################
##
#F SizeOfUnionRec( sub ) -- recursive version
##
SizeOfUnionRec := function( list )
    local s, i, int, t, n;
    if Length( list ) = 0 then return 1; fi;
    s := Size( list[1] );
    n := Length( list );
    for i in [2..n] do
        int := List( list{[1..i-1]}, x -> Intersection( x, list[i]));
        t := SizeOfUnion( int );
        s := s + Size( list[i] ) - t;
    od;
    return s;
end;

#############################################################################
##
#F SizeOfUnionTriv -- trivial version
##
SizeOfUnionTriv := function( list )
    return Length( Union( List( list, Elements ) ) );
end;

#############################################################################
##
#F SizeOfUnion -- main function
##
SizeOfUnion := function( sub )
    local list;
    list := FindMaximals( sub );
    if Length(list) = 0 then return 1; fi;
    if Length(list) = 1 then return Size(list[1]); fi;
    if Sum(List(list, Size)) < 2000 then
        return SizeOfUnionTriv(list);
    fi;
    return SizeOfUnionRec(list);
end;

#############################################################################
##
#F SizeOfUnionMod( subs, e ) . . . . . . . . .size of the union of subs mod e
##
## <subs> is a list of bases containing (eZ)^d. Compute the size of the union
## of <subs> in (Z/eZ)^d.
##
## This function needs to be profiled.
##
SizeOfUnionMod := function( subs, e )
    local d, F, V, b, news;

    # the trivial case
    if Length( subs ) = 0 then return 1; fi;
    d := Length( subs[1] );

    if IsPrimeInt( e ) then
        F := GF(e);
        V := F^d;
        b := BasisVectors( Basis( V ) );
        news := List( subs, x -> x * b );
        news := List( news, x -> Subspace( V, x ) );
        if ForAny( news, x -> Size(x) = e^d ) then return e^d; fi;
        # return SizeOfUnion( news );
        return Length( Union( List( news, Elements ) ) );
    fi;

    V := AbelianGroup( List( [1..d], x -> e ) );
    b := GeneratorsOfGroup(V);
    news := List( subs, x -> List( x, y -> MappedVector( y, b ) ) );
    news := List( news, x -> Subgroup( V, x ) );
    if ForAny( news, x -> Size(x) = e^d ) then return e^d; fi;

    # return SizeOfUnion( news );
    return Length( Union( List( news, Elements ) ) );
end;

max_stars_repo_name: alex-konovalov/aclib

max_issues_repo_name: alex-konovalov/aclib

max_forks_repo_name: alex-konovalov/aclib


text:
Class(OpenMP_POWERUnparser, paradigms.smp.OpenMP_UnparseMixin, packages.powerisa.power9.p9macro.POWER9Unparser);
Class(OpenMP_POWERUnparser_ParFor, paradigms.smp.OpenMP_UnparseMixin_ParFor, packages.powerisa.power9.p9macro.POWER9Unparser);

max_stars_repo_name: franzfranchetti/spiral-package-fftx

max_issues_repo_name: franzfranchetti/spiral-package-fftx

max_forks_repo_name: franzfranchetti/spiral-package-fftx


text:
# Copyright (c) 2018-2021, Carnegie Mellon University
# See LICENSE for details

Declare(PropagateUnification);

list_dup := (n, a) -> Cond(IsList(a), Checked(Length(a)=n, a), Replicate(n, a));

getXTypeOpts := (t, opts) -> When(IsBound(t.a.t_in),  List(t.a.t_in,  e -> TPtr(e)), list_dup(t.arity()[2], opts.XType));
getYTypeOpts := (t, opts) -> When(IsBound(t.a.t_out), List(t.a.t_out, e -> TPtr(e)), list_dup(t.arity()[1], opts.YType));

# TypeOpts(<sums>, <opts>) returns options with 'XType' and 'YType' fields set acording
# to 't_in' and 't_out' attributes (when available).
# Playing ping-pong between SumsUnification and TypeOpts allows to perform search for OL and
# multi data type formulas on existing DP search infrastructure.

TypeOpts := (sums, opts) -> Cond( not IsBound(opts.XType) or not IsBound(opts.YType), opts,
    CopyFields(opts, rec(
        XType := StripList(getXTypeOpts(sums, opts)),
        YType := StripList(getYTypeOpts(sums, opts))
)));

# SumsUnification sets data type 't_in' and 't_out' attributes on sums tree..
# It propagates 't_in' types of sums root object (if available) or opts.XType.
SumsUnification := function(sums, opts)
    if IsBound(opts.doSumsUnification) and opts.doSumsUnification then
        sums := PropagateUnification(sums, List(getXTypeOpts(sums, opts), e -> e.t));
    fi;
    return sums;
end;

Class(PropagateUnification, HierarchicalVisitor, rec(

     # makes sure the object is not a polymorphic function/SPL like L,
     # which would cause trouble because they are handled below via
     # FuncClass, and thus .a.t_in and t_out are not set
     _isSPL  := x -> IsSPL(x) and not IsFunction(x),

     _unify  := (t1, t2) -> Cond( IsList(t1) and IsList(t2) and Length(t1)=Length(t2),
                                   List([1..Length(t1)], i -> UnifyPair(t1[i], t2[i])),
                               IsList(t1) and IsType(t2),
                                   List([1..Length(t1)], i -> UnifyPair(t1[i], t2)),
                               IsType(t1) and IsType(t2),
                                   UnifyPair(t1, t2),
                               Error("PropagateUnification: type dimensions mismatch")),
     _to_cplx := (lst) -> List(lst, e -> T_Complex(e)),
     _to_real := (lst) -> List(lst, e -> e.realType()),

     TCast := meth(self, o, inputs)
         o.params[3] := UnifyTypes(inputs :: [o.params[3]]);
         return o.withA( t_in  => [o.params[3]],
                         t_out => [o.params[2]]);
     end,

     TCvt := (self, o, inputs) >>
         o.withA( t_in  => [o.isa_from().t.base_t()],
                  t_out => [o.isa_to().t.base_t()  ] ),

     Cvt := ~.TCvt,

     Compose := (self, o, inputs) >> let(
         ch    := DropLast(ScanR( o.rChildren(), (l, e) -> When(IsList(l), self(e, l), self(e, l.a.t_out)), inputs ), 1),
         o.from_rChildren(ch).withA( t_in  => Last(ch).a.t_in,
                                     t_out => ch[   1].a.t_out)),

     RC := ( self, o, inputs) >> let(
         ch := List(o.rChildren(), e -> When( self._isSPL(e), self(e, self._to_cplx(inputs)), e)),
         o.from_rChildren(ch).withA(t_in => self._to_real(ch[1].a.t_in), t_out => self._to_real(ch[1].a.t_out))),
     GRC := ~.RC,
     TRC := ~.RC,

     BaseMat := (self, o, inputs) >> let(
         dmn := List(o.dmn(), e -> e.t.base_t()),
         rng := List(o.rng(), e -> e.t.base_t()),
         inp := List(Zip2(dmn, inputs), UnifyTypes),
         # we have no idea how to propagate data types for objects which arity is other than [1, 1],
         # let's hope object defined it's rng() and dmn() data types or uses same data type everywhere
         # or data type propagated horizontally in the case of same arity on the left and on the right.
         rnginp := Replicate( Length(rng), UnifyTypes(inp) ),
         out := List(Zip2(rng, Cond( Length(dmn)=Length(rng), inp, rnginp )), UnifyTypes),
         o.withA( t_in => inp, t_out => out )),

     ClassSPL := (self, o, inputs) >> let(
         ch := List(o.rChildren(), e -> When( self._isSPL(e), self(e, inputs), e )),
         spl_ch := Filtered(ch, self._isSPL),
         When( spl_ch <> [],
             o.from_rChildren(ch).withA(
                 t_in  => FoldL1(List(spl_ch, x->x.a.t_in ), self._unify),
                 t_out => FoldL1(List(spl_ch, x->x.a.t_out), self._unify)),
             self.BaseMat(o, inputs))),

     Cross := meth (self, o, inputs)
         local ch, i, j, d;
         ch := o.rChildren();
         j  := 1;
         for i in [1..Length(ch)] do
             d     := Length(ch[i].dmn());
             ch[i] := self(ch[i], Flat(inputs{[j..j+d-1]}));
             j     := j + d;
         od;
         return o.from_rChildren(ch).withA(
                     t_in  => Flat(List(ch, e -> e.a.t_in)),
                     t_out => Flat(List(ch, e -> e.a.t_out)));
     end,

     SMPBarrier := (self, o, inputs) >> let( ch := self(o.child(1), inputs),
         ObjId(o)(o.nthreads, o.tid, ch).withA( t_in => ch.a.t_in, t_out => ch.a.t_out)),

     VContainer := ~.ClassSPL,
     SumsBase   := ~.ClassSPL,
     VPerm      := ~.BaseMat,

     TaggedNonTerminal := meth( self, o, inputs)
         local tag, tL, tR, tM;
         tag := FirstDef(o.getTags(), t -> Global.autolib.IsACVec(t), false);
         if tag=false then
             return self.ClassSPL(o, inputs);
         else
             tL := List(Flat([tag.isaL()]), e->e.t.base_t());
             tR := List(Flat([tag.isaR()]), e->e.t.base_t());
             tM := List(list_dup(o.arity()[2], tag.isaM()), e->e.t.base_t());
             if IsComplexT(inputs[1]) then
                 tL := List(tL, e -> T_Complex(e));
                 tR := List(tR, e -> T_Complex(e));
                 tM := List(tM, e -> T_Complex(e));
             fi;
             return self.ClassSPL(o, tM).withA(t_in => tR, t_out => tL);
         fi;
     end,

     GTBase := ~.TaggedNonTerminal,

     FuncClass     := Ignore,
     FuncClassOper := Ignore,
));


text:
dimX := 0;
dimY := 1;
dimZ := 2;
dimW := 3;

dimXYZ := [dimX, dimY, dimZ];

max_stars_repo_name: spiral-software/spiral-package-mpi

max_issues_repo_name: spiral-software/spiral-package-mpi

max_forks_repo_name: spiral-software/spiral-package-mpi


text:
# Built-in, using sets (which are also lists)
a := [ 1, 2, 3, 1, [ 4 ], 5, 5, [4], 6 ];
# [ 1, 2, 3, 1, [ 4 ], 5, 5, [ 4 ], 6 ]
b := Set(a);
# [ 1, 2, 3, 5, 6, [ 4 ] ]
IsSet(b);
# true
IsList(b);
# true

ext: gap

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text:
##  Copyright (c) 2018-2021, Carnegie Mellon University
##  See LICENSE for details

Class(RulesDropSymDecl, RuleSet);
RewriteRules(RulesDropSymDecl, rec(
    drop_decl := Rule(@@(1, decl, (e,cx)->ForAny(cx.opts.symbol, i->i in e.vars)),
        (e,cx)-> decl(Filtered(@@(1).val.vars, i->not i in cx.opts.symbol), @@(1).val.cmd))
));


text:
#
#
#

Read("~/Workspace/Chevalley.gap/init.gi");

Read(Filename(home_dir,"lib/rsys.gd"));
Read(Filename(home_dir,"lib/rsys.gi"));

Read(Filename(home_dir,"lib/chvadj.gd"));
Read(Filename(home_dir,"lib/chvadj.gi"));

Read(Filename(home_dir,"lib/nilchv.gd"));
Read(Filename(home_dir,"lib/nilchv.gi"));

Read(Filename(home_dir,"lib/algU.gd"));
Read(Filename(home_dir,"lib/algU.gi"));

#
# UipotentChv needs SolveRelations
#
Read(Filename(home_dir,"lib/poly.gd"));
Read(Filename(home_dir,"lib/poly.gi"));

#
# Unialg needed
#
Read(Filename(home_dir,"lib/unichv.gd"));
Read(Filename(home_dir,"lib/witt.gd"));
Read(Filename(home_dir,"lib/unialg.gd"));

Read(Filename(home_dir,"lib/unichv.gi"));


text:
SelectionSort := function(v)
   local i, j, k, n, m;
   n := Size(v);
   for i in [1 .. n] do
      k := i;
      m := v[i];
      for j in [i + 1 .. n] do
         if v[j] < m then
            k := j;
            m := v[j];
         fi;
      od;
      v[k] := v[i];
      v[i] := m;
   od;
end;

v := List([1 .. 100], n -> Random([1 .. 100]));
SelectionSort(v);
v;


text:
#############################################################################
##
##  Functions.gi        JConvex package
##                      Martin Bies
##
##  Copyright 2021      University of Pennsylvania
##
##  A Gap package to do convex geometry by Polymake
##
## Chapter: Functionality
##
#############################################################################

#############################################################################
## Availability of required software
#############################################################################

##
InstallMethod( PolymakeAvailable, [  ],
  function( )
    # Check if Polymake and the gap-julia interface are available
    if IsPackageMarkedForLoading( "JuliaInterface", ">= 0.5.2" ) then
        if IsBoundGlobal("_Polymake_jl") and IsJuliaObject(ValueGlobal("_Polymake_jl")) then
            return true;
        fi;
    fi;
    return false;
end );

##
InstallMethod( CddInterfaceAvailable, [  ],
  function( )
    return IsPackageMarkedForLoading( "CddInterface", ">= 2020.06.24" );
end );

##
InstallMethod( NormalizInterfaceAvailable, [  ],
  function( )
    return IsPackageMarkedForLoading( "NormalizInterface", ">= 1.2.0" );
end );


text:
# Copyright (c) 2018-2021, Carnegie Mellon University
# See LICENSE for details

# ==========================================================================
# Conjugate(<spl>, <perm>)
#    perm^-1 * spl * perm
# ==========================================================================
Class(Conjugate, BaseOperation, rec(
    new := (self, spl, conj) >> Checked(IsSPL(spl), IsSPL(conj),
	When(Dimensions(conj) = [1,1], spl,
	     SPL(WithBases( self,
		     rec( _children := [spl, conj],
			 dimensions := spl.dimensions ))))),
    #-----------------------------------------------------------------------
    isPermutation := False,
    #-----------------------------------------------------------------------
    dims := self >> self.dimensions,
    #-----------------------------------------------------------------------
    toAMat := self >>
        ConjugateAMat(AMatSPL(self._children[1]), AMatSPL(self._children[2])),
    #-----------------------------------------------------------------------
    print := (self, i, is) >>
        Print("(", self.child(1).print(i, is), ") ^ (", self.child(2).print(i+is, is), ")",
              self.printA()),
    #-----------------------------------------------------------------------
    transpose := self >>
       Inherit(self, rec(_children := [ TransposedSPL(self._children[1]),
					self._children[2] ],
	                dimensions := Reversed(self.dimensions))),
    #-----------------------------------------------------------------------
    arithmeticCost := (self, costMul, costAddMul) >>
        self._children[1].arithmeticCost(costMul, costAddMul)
));


text:
package main

import (
	"net/http"

	"github.com/unrolled/secure" // or "gopkg.in/unrolled/secure.v1"
)

var myHandler = http.HandlerFunc(func(w http.ResponseWriter, r *http.Request) {
	w.Write([]byte("hello world"))
})

func main() {
	secureMiddleware := secure.New(secure.Options{
		AllowedHosts:          []string{"example.com", "ssl.example.com"},
		HostsProxyHeaders:     []string{"X-Forwarded-Host"},
		SSLRedirect:           true,
		SSLHost:               "ssl.example.com",
		SSLProxyHeaders:       map[string]string{"X-Forwarded-Proto": "https"},
		STSSeconds:            315360000,
		STSIncludeSubdomains:  true,
		STSPreload:            true,
		FrameDeny:             true,
		ContentTypeNosniff:    true,
		BrowserXssFilter:      true,
		ContentSecurityPolicy: "script-src $NONCE",
		PublicKey:             `pin-sha256="base64+primary=="; pin-sha256="base64+backup=="; max-age=5184000; includeSubdomains; report-uri="https://www.example.com/hpkp-report"`,
		IsDevelopment: false,
	})

	app := secureMiddleware.Handler(myHandler)
	http.ListenAndServe("127.0.0.1:3000", app)
}

max_stars_repo_name: visuo-dev/yaegi

max_issues_repo_name: visuo-dev/yaegi

max_forks_repo_name: visuo-dev/yaegi


text:
# Copyright (c) 2018-2021, Carnegie Mellon University
# See LICENSE for details

#F   tSPL Twiddle factors
Class(TTwiddle, Tagged_tSPL_Container, rec(
    abbrevs :=  [
        (mn, n) -> [mn, n, 1],
        (mn, n, k) -> [mn, n, k]
     ],
    dims := self >> let(n:=self.params[1], [n,n]),
    terminate := self >> Diag(Tw1(self.params[1], self.params[2], self.params[3])),
    transpose := self >> Copy(self),
    isReal := False,
    doNotMeasure := true,
    noCodelet := true,
    normalizedArithCost := self >> 6*Rows(self)
));

Class(TTwiddle_Stockham, Tagged_tSPL_Container, rec(
    abbrevs :=  [
        (n, blocksize, rdx, j) -> [n, blocksize, rdx, j]
     ],
    dims := self >> let(n:=self.params[1], [n,n]),
    #params := self >> let(params := [self.params[1], self.params[2], self.params[3], self.params[4]]),
    terminate := self >> Tw1(self.params[1], self.params[1]/self.params[3], 1),
#F 	D tensor I has to be done in the TTwiddle_St rule, not here
#F    terminate := self >> diagTensor(Tw1(self.params[1], self.params[2], self.params[3]), fConst(TComplex, self.params[2], 1)),
#F    terminate := self >> fCompose(dOmega(self.params[1],1), fId(self.params[1]/self.params[2])),
    transpose := self >> Copy(self),
    isReal := False,
    doNotMeasure := true,
    noCodelet := true,
    domain := self >> TComplex,
    normalizedArithCost := self >> 6*Rows(self)
));

NewRulesFor(TTwiddle_Stockham, rec(
    TTwiddle_St := rec(
        minSize := false,
        maxSize := false,
        applicable := (self, nt) >> ((IsBool(self.minSize) and not self.minSize) or nt.params[1] >= self.minSize) and
            ((IsBool(self.maxSize) and not self.maxSize) or nt.params[1] <= self.maxSize),
        forTransposition := false,
        apply := function(t, C, Nonterms)
            local access, blocksize, i, i1, j, n, d, k, rdx, stride, t;

			n := t.params[1];
			blocksize := t.params[2];
			rdx := t.params[3];
			j := t.params[4];
			t := LogInt(n, rdx);
			i := Ind(n);
			i1 := Ind(n/(rdx^j));
			stride := rdx^j;
			access := Lambda(i1, i1 * stride);

#F	Older versions, not using an access function, but returning the full matrix
#F			TDiag(fPrecompute(diagTensor(fConst(TComplex, rdx^j, 1), Stockham_radix.gen(rdx, LogInt(N, rdx)-j-1))))
#F			omega := fCompose(FData(fCompose(dOmega(n,1),diagTensor(dLin(div(n,blocksize), 1, 0, TInt), dLin(blocksize, 1, 0, TInt))).tolist()), Lambda(i, i));
#F			omega := fCompose(fPrecompute(fCompose(dOmega(n,1),diagTensor(dLin(div(n,blocksize), 1, 0, TInt), dLin(blocksize, 1, 0, TInt))).tolist()), Lambda(i, i));

			omega := fCompose(fPrecompute(Tw1(n, n/rdx,1)), access);
			omega := diagTensor(omega, fConst(TComplex, rdx^j, 1));

        	return Diag(omega);

		end
    ),
));

NewRulesFor(TTwiddle, rec(
    TTwiddle_Tw1 := rec(
        minSize := false,
        maxSize := false,
        applicable := (self, nt) >> ((IsBool(self.minSize) and not self.minSize) or nt.params[1] >= self.minSize) and
            ((IsBool(self.maxSize) and not self.maxSize) or nt.params[1] <= self.maxSize),
        forTransposition := false,
        apply := (t, C, Nonterms) -> Diag(fPrecompute(Tw1(t.params[1], t.params[2], t.params[3])))
    ),

    TTwiddle_TwoTables := rec(
        minSize := false,
        maxSize := false,
        applicable := (self, nt) >> ((IsBool(self.minSize) and not self.minSize) or nt.params[1] >= self.minSize) and
            ((IsBool(self.maxSize) and not self.maxSize) or nt.params[1] <= self.maxSize),
        forTransposition := false,
        apply := function(t, C, Nonterms)
            local mn, n, _k, rdx, dp, i, i1, i2, j, k, f, omega_1, omega_2;

            #t in this case is the TTwiddle(mn,n,k) function, so the t.params[1]  = first parameter of TTwiddle() = mn
            mn := t.params[1];
            n := t.params[2];
            _k := t.params[3];
            dp := DivisorPairs(mn);
            rdx := Maximum(dp[Int(Length(dp)/2+1)]);

            i := Ind(mn);
            i1 := Ind(mn);
            i2 := Ind(mn);
            j := Ind(mn/rdx);
            k := Ind(rdx);

            f := Lambda(i, idiv(i, n) * imod(i, n));	# -> guarantees contiguous access

#   NOTE: fPrecompute does not work, so I have to force table generation through FData()
#            omega_1 := fCompose(fPrecompute(dOmega(mn/rdx, _k)), fComputeOnline(Lambda(i1, idiv(f.at(i1), rdx))));
#            omega_2 := fCompose(fPrecompute(dOmega(mn, _k)), fComputeOnline(Lambda(i2, imod(f.at(i2), rdx))));
#            omega_1 := fCompose(fPrecompute(fCompose(dOmega(mn/rdx, _k), fId(mn/rdx))), fComputeOnline(Lambda(i1, idiv(f.at(i1), rdx))));
#            omega_2 := fCompose(fPrecompute(fCompose(dOmega(mn, _k), fId(rdx))), fComputeOnline(Lambda(i2, imod(f.at(i2), rdx))));

#   NOTE: the domain of dOmega cannot be set to an integer, so I need to compose with fId()...
            omega_1 := fCompose(FData(fCompose(dOmega(mn/rdx, _k), fId(mn/rdx)).tolist()), Lambda(i1, idiv(f.at(i1), rdx)));
            omega_2 := fCompose(FData(fCompose(dOmega(mn, _k), fId(rdx)).tolist()), Lambda(i2, imod(f.at(i2), rdx)));

            return Diag(omega_1) * Diag(omega_2);

# old version with FData without fCompose
#            omega_1 := FData(Lambda(j, omega(mn/rdx, j*imod(_k, mn/rdx))).tolist());
#            omega_2 := FData(Lambda(k, omega(mn, k*_k)).tolist());
#            return Diag(Lambda(i1, omega_1.at(idiv(f.at(i1), rdx)))) * Diag(Lambda(i2, omega_2.at(imod(f.at(i2), rdx))));
        end
    ),
));


text:
IsKaprekar := function(n)
	local a, b, p, q;
	if n = 1 then
		return true;
	fi;
	q := n*n;
	p := 10;
	while p < q do
		a := RemInt(q, p);
		b := QuoInt(q, p);
		if a > 0 and a + b = n then
			return true;
		fi;
		p := p*10;
	od;
	return false;
end;

Filtered([1 .. 10000], IsKaprekar);
# [ 1, 9, 45, 55, 99, 297, 703, 999, 2223, 2728, 4879, 4950, 5050, 5292, 7272,
#   7777, 9999 ]

Size(last);
# 17

Filtered([1 .. 1000000], IsKaprekar);
# [ 1, 9, 45, 55, 99, 297, 703, 999, 2223, 2728, 4879, 4950, 5050, 5292, 7272,
#   7777, 9999, 17344, 22222, 38962, 77778, 82656, 95121, 99999, 142857,
#   148149, 181819, 187110, 208495, 318682, 329967, 351352, 356643, 390313,
#   461539, 466830, 499500, 500500, 533170, 538461, 609687, 627615, 643357,
#   648648, 670033, 681318, 791505, 812890, 818181, 851851, 857143, 961038,
#   994708, 999999 ]

Size(last);
# 54

IsKaprekarAndHow := function(n, base)
	local a, b, p, q;
	if n = 1 then
		return true;
	fi;
	q := n*n;
	p := base;
	while p < q do
		a := RemInt(q, p);
		b := QuoInt(q, p);
		if a > 0 and a + b = n then
			return [a, b];
		fi;
		p := p*base;
	od;
	return false;
end;

IntegerToBaseRep := function(n, base)
	local s, digit;
	if base > 36 then
		return fail;
	elif n = 0 then
		return "0";
	else
		s := "";
		digit := "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
		while n <> 0 do
			Add(s, digit[RemInt(n, base) + 1]);
			n := QuoInt(n, base);
		od;
		return Reversed(s);
	fi;
end;

PrintIfKaprekar := function(n, base)
	local v;
	v := IsKaprekarAndHow(n, base);
	if IsList(v) then
		Print(n, "(10) or in base ", base, ", ",
			IntegerToBaseRep(n, base), "^2 = ",
			IntegerToBaseRep(n^2, base), " and ",
			IntegerToBaseRep(v[2], base), " + ",
			IntegerToBaseRep(v[1], base), " = ",
			IntegerToBaseRep(n, base), "\n");
	fi;
	return fail;
end;

# In base 17...
Perform([1 .. 1000000], n -> PrintIfKaprekar(n, 17));
# 16(10) or in base 17, G^2 = F1 and F + 1 = G
# 64(10) or in base 17, 3D^2 = E2G and E + 2G = 3D
# 225(10) or in base 17, D4^2 = A52G and A5 + 2G = D4
# 288(10) or in base 17, GG^2 = GF01 and GF + 1 = GG
# 1536(10) or in base 17, 556^2 = 1B43B2 and 1B4 + 3B2 = 556
# 3377(10) or in base 17, BBB^2 = 8093B2 and 809 + 3B2 = BBB
# 4912(10) or in base 17, GGG^2 = GGF001 and GGF + 1 = GGG
# 7425(10) or in base 17, 18BD^2 = 24E166G and 24E + 166G = 18BD
# 9280(10) or in base 17, 1F1F^2 = 39B1B94 and 39B + 1B94 = 1F1F
# 16705(10) or in base 17, 36DB^2 = B992C42 and B99 + 2C42 = 36DB
# 20736(10) or in base 17, 43CD^2 = 10DE32FG and 10DE + 32FG = 43CD
# 30016(10) or in base 17, 61EB^2 = 23593F92 and 2359 + 3F92 = 61EB
# 36801(10) or in base 17, 785D^2 = 351E433G and 351E + 433G = 785D
# 37440(10) or in base 17, 7A96^2 = 37144382 and 3714 + 4382 = 7A96
# 46081(10) or in base 17, 967B^2 = 52G94382 and 52G9 + 4382 = 967B
# 46720(10) or in base 17, 98B4^2 = 5575433G and 5575 + 433G = 98B4
# 53505(10) or in base 17, AF26^2 = 6GA43F92 and 6GA4 + 3F92 = AF26
# 62785(10) or in base 17, CD44^2 = 9A5532FG and 9A55 + 32FG = CD44
# 66816(10) or in base 17, DA36^2 = AEG42C42 and AEG4 + 2C42 = DA36
# 74241(10) or in base 17, F1F2^2 = D75F1B94 and D75F + 1B94 = F1F2
# 76096(10) or in base 17, F854^2 = E1F5166G and E1F5 + 166G = F854
# 83520(10) or in base 17, GGGG^2 = GGGF0001 and GGGF + 1 = GGGG
# 266224(10) or in base 17, 33334^2 = A2C52A07G and A2C5 + 2A07G = 33334


text:
#
# ---------- Constructors ----------
#

InstallMethod(Unipotent,
              "Unipotent element in simple adjoint type",
              [IsChevalleyAdj,IsList,IsList,IsInt],
              function(sys,coeffs,ordering,order)
              local object;

              if Filtered(coeffs,i->not i[2] in Integers)=[] then coeffs:=List(coeffs,i->[i[1],One(ring(sys))*i[2]]);
              elif Filtered(coeffs,i->not i[2] in ring(sys))<>[] then Error("The coefficients ar not in indicated ring."); fi;

              object:=Objectify(NewType(NewFamily("UnipotentFamily"),
                                        IsAttributeStoringRep and
                                        IsUnipotent and
                                        IsMultiplicativeElementWithInverse and
                                        IsMultiplicativeElementWithOne),
                           rec());

              SetchevalleyAdj(object,sys);
              SetOrdering(object,ordering);

              Setcoefficients(object,Canonic(sys,coeffs,ordering));
#             SetInverse(object,InverseOp(object)); This is set at the first call Inverse(obj);

              if order = -1 and Characteristic(sys)>0 then
                  SetOrder(object,OrderOp(object));
              elif order > -1 then
                  SetOrder(object,order);
              fi;

              SetName(object,Concatenation("<unipotent element for ",type(sys),String(rank(sys)),
                                      " in characteristic ",String(Characteristic(sys)),">"));
              return object;
end);

InstallMethod(Unipotent,
              "Unipotent element in simple adjoint type",
              [IsChevalleyAdj,IsList,IsList],
              function(sys,coeffs,ordering)

              return Unipotent(sys,coeffs,ordering,-2);
end);

InstallMethod(Unipotent,
              "Unipotent element in simple adjoint type",
              [IsChevalleyAdj,IsList,IsInt],
              function(sys,coeffs,order)
              local ordering,tmp,i,slot;

              tmp:=ShallowCopy(coeffs);
              Sort(tmp,function(a,b) return a[1]<b[1]; end);

              ordering:=[1..Length(positiveRoots(sys))];
              tmp:=List(tmp,i->ordering[i[1]]);

              for i in [1..Length(coeffs)] do
                  ordering[coeffs[i][1]]:=tmp[i];
              od;

              return Unipotent(sys,Filtered(coeffs,i->i[2]<>Zero(ring(sys))),ordering,order);
end);

InstallMethod(Unipotent,
              "Unipotent element in simple adjoint type",
              [IsChevalleyAdj,IsList],
              function(sys,coeffs)

              return Unipotent(sys,coeffs,-2);#-1);
end);

InstallMethod(Unipotent,
              "Unipotent element in simple adjoint type",
              [IsNilpotentChv],
              function(e)

              return Unipotent(chevalleyAdj(e),coefficients(e));
end);

#
# ---------- Canonic form of element ----------
#

InstallMethod(Canonic,
              "Canonic form of coefficients of unipotent element in given ordering",
              [IsChevalleyAdj,IsList,IsList],
              function(sys,coeffs,ordering)
              local lista,
              pr,Cijrs,
              flag,i,temp,suma,param,cc,k;

              lista:=List(coeffs,i->ShallowCopy(i));
              pr:=positiveRoots(sys);

              Cijrs:=C(sys);
              flag:=true;
              while flag do
                  flag:=false;
                  i:=Length(lista);
                  while 0 < i do
                            if lista[i][2] = Zero(ring(sys)) then
                                Remove(lista,i);
                                flag:=true;
                            elif 1 < i and lista[i][1] = lista[i-1][1] then
                                lista[i-1][2]:=lista[i-1][2]+lista[i][2];
                                Remove(lista,i);
                                flag:=true;
                            elif 1 < i and ordering[lista[i][1]] < ordering[lista[i-1][1]] then
                                temp:=lista[i]; lista[i]:=lista[i-1]; lista[i-1]:=temp;

                                #aici r si s sunt pe pozitia i-1 si repsectiv i
                                k:=1;
                                for cc in Cijrs[lista[i-1][1]][lista[i][1]] do
                                    suma:=cc[1]*pr[lista[i-1][1]]+cc[2]*pr[lista[i][1]];
                                    param:=cc[3]*((-1)^cc[1])*(lista[i-1][2]^cc[1])*(lista[i][2]^cc[2]);
                                    Add(lista,[Position(pr,suma),param],i+k);
                                    k:=k+1;
                                od;
                                flag:=true;
                            fi;
                            i:=i-1;
                  od;
              od;

              return Immutable(lista);
end);

#
# ---------- Arithmetic Operations ----------
#

InstallMethod(InverseMutable,
              "Inverse of unipotent element",
              [IsUnipotent],
              function(u)
              local i,lista,len,reverse_order;

              lista:=List(coefficients(u),i->StructuralCopy(i));
              len:=Length(lista);
              for i in [1..len] do
                  Add(lista,[lista[len-i+1][1],-lista[len-i+1][2]]);
                  Remove(lista,len-i+1);
              od;

              reverse_order:=function(o)
                  local maxim;
                  maxim:=Maximum(o);
                  return List(o,i->maxim-i+1);
              end;

              return Unipotent(chevalleyAdj(u),lista,reverse_order(Ordering(u)),-2);#0);
end);

InstallMethod(OneMutable,
              "Inverse of unipotent element",
              [IsUnipotent],
              function(u)

              return One(chevalleyAdj(u),[],0);
end);

InstallMethod(\*,
              "Multiplication for unipotent elements a*b with the ordering of b",
              [IsUnipotent,IsUnipotent],
              function(u1,u2)
              local sys1,sys2,
              generic,result,coeffs1,coeffs2,
              APR,avars,pr_len,
              i;

              sys1:=chevalleyAdj(u1);
              sys2:=chevalleyAdj(u2);
              if type(sys1)<>type(sys2) or
                 rank(sys1)<>rank(sys2) or
                 ring(sys1)<>ring(sys2) then Error("Not in the same family."); fi;

              return Unipotent(sys2,Concatenation(coefficients(u1),coefficients(u2)),Ordering(u2));
end);

InstallMethod(\*,
              "Action of unipotent element u on nilpotent e: u*e",
              [IsUnipotent,IsNilpotentChv],
              function(u,e)
              local result,
              sys1,sys2,Mrsi,pr,
              uas,r,s,i,pos;

              sys1:=chevalleyAdj(u);
              sys2:=chevalleyAdj(e);
              if type(sys1)<>type(sys2) or
                 rank(sys1)<>rank(sys2) or
                 ring(sys1)<>ring(sys2) then Error("Not in the same family.");
              fi;

              pr:=positiveRoots(sys1);
              Mrsi:=M(sys1);

              uas:=Reversed(coefficients(u));
              result:=List(coefficients(e),i->ShallowCopy(i));

              # for this see Carter pg.64
              for r in uas do
                  for s in result do
                      for i in [1..Length(Mrsi[r[1]][s[1]])] do
                          pos:=Position(pr,i*pr[r[1]]+pr[s[1]]);
                          result:=Concatenation(result,[[pos,Mrsi[r[1]][s[1]][i]*r[2]^i*s[2]]]);
                      od;
                  od;
              od;

              return NilpotentChv(sys2,result);
end);

InstallMethod(OrderOp,
              "Order of unipotent element",
              [IsUnipotent],
              function(u)
              local result,
              p,u_to_p,putere,i;

              if coefficients(u)=[] then return 1; fi;

              p:=Characteristic(chevalleyAdj(u));

              u_to_p:=u;
              for i in [1..p-1] do u_to_p:=u_to_p*u; od;

              result:=1;
              putere:=u_to_p;
              while coefficients(putere)<>[] do
                  putere:=putere*u_to_p;
                  result:=result+1;
              od;

              return p*result;
end);

InstallMethod(Comm,
              "Commutator for unipotent elements a^-1b^-1ab with the ordering of b",
              [IsUnipotent,IsUnipotent],
              function(u1,u2)
              local sys1,sys2;

              sys1:=chevalleyAdj(u1);
              sys2:=chevalleyAdj(u2);
              if type(sys1)<>type(sys2) or
                 rank(sys1)<>rank(sys2) or
                 ring(sys1)<>ring(sys2) then Error("Not in the same family."); fi;
              return Unipotent(sys2,
                                  Concatenation(coefficients(u1^-1),coefficients(u2^-1),
                                                     coefficients(u1),coefficients(u2)),
                                  Ordering(u2));
end);

InstallMethod(Conj,
              "Conjugation for unipotent elements b^-1ab with the ordering of b",
              [IsUnipotent,IsUnipotent],
              function(u1,u2)
              local sys1,sys2;

              sys1:=chevalleyAdj(u1);
              sys2:=chevalleyAdj(u2);
              if type(sys1)<>type(sys2) or
                 rank(sys1)<>rank(sys2) or
                 ring(sys1)<>ring(sys2) then Error("Not in the same family."); fi;
              return Unipotent(sys2,
                                  Concatenation(coefficients(u2^-1),coefficients(u1),coefficients(u2)),
                                  Ordering(u2));
end);

InstallMethod(Subword,
              "Subword of a unipotent element by list of roots",
              [IsUnipotent,IsList],
              function(u,l)
              local coeffs;
              coeffs:=Filtered(coefficients(u),i->i[1] in l);

              return Unipotent(chevalleyAdj(u),coeffs,Ordering(u));
end);

InstallMethod(InsertWord,
              "Shuffles a word into a unipotent element u according to the ordering of u",
              [IsUnipotent,IsList],
              function(u,word)
              local coeffs,
              sys,check,ord;

              sys:=chevalleyAdj(u);
              check:=ForAll(word,i-> i[2] in ring(sys));
              if not check then
                  Error("The word cannot be shuffled into u: not over the same ring!\n");
              fi;

              coeffs:=coefficients(u);
              coeffs:=Concatenation(coeffs,word);
              ord:=Ordering(u);
              Sort(coeffs,function(a,b) return ord[a[1]]<ord[b[1]]; end);
              coeffs:=Canonic(sys,coeffs,ord);

#              Setcoefficients(u,coeffs); Doesn't work this way

              return Unipotent(sys,coeffs,Ordering(u));
end);

#
# ---------- Lie Algebra Action ----------
#

InstallMethod(AduAlphat,
              "Adu_{alpha}(t) on the whole Lie algebra",
              [IsChevalleyAdj,IsPosInt,IsRingElement],
              function(sys,root,param)
              local result,
              pr,pr_len,roots,
              A_rs,M_rsi,
              L,CB,Rank,B,cochars,h,
              tmp,zero,unu,s,i,poz;

              pr:=positiveRoots(sys);
              pr_len:=Length(pr);
              roots:=Concatenation(pr,-pr);

              L:=lieAlgebra(sys);
              CB:=ChevalleyBasis(L);
              Rank:=rank(sys);
              B:=Basis(L);
              cochars:=Cocharacters(L);
              cochars:=List(cochars,i->Coefficients(B,i));

              zero:=Zero(ring(sys));
              unu:=One(ring(sys));

              A_rs:=A(sys);
              M_rsi:=M(sys);

              result:=[];
              # s in positive roots
              for s in [1..pr_len] do
                  tmp:=List([1..Dimension(L)],i->zero);

                  if s = root then
                      tmp[s]:=unu;
                  else
                      tmp[s]:=unu;
                      for i in [1..Length(M_rsi[root][s])] do
                          poz:=Position(roots,i*roots[root]+roots[s]);
                          tmp[poz]:=M_rsi[root][s][i]*param^i;
                      od;
                  fi;
                  Add(result,tmp);
              od;

              # s in negative roots
              for s in [1..pr_len]+pr_len do
                  tmp:=List([1..Dimension(L)],i->zero);

                  if s = root+pr_len then
                      # e_{-r}
                      tmp[s]:=unu;
                      # -t^2*e_{r}
                      tmp[root]:=-param^2;
                      # t*h_{r}
                      tmp:=tmp+param*cochars[root];
                  else
                      tmp[s]:=unu;
                      for i in [1..Length(M_rsi[root][s])] do
                          poz:=Position(roots,i*roots[root]+roots[s]);
                          tmp[poz]:=M_rsi[root][s][i]*param^i;
                      od;
                  fi;
                  Add(result,tmp);
              od;

              # h_s
              for s in [1..Rank] do
                  tmp:=List([1..Dimension(L)],i->zero);
                  # h_s
                  tmp:=tmp+cochars[s];
                  # -A_{sr}*t*e_r
                  tmp[root]:=-A_rs[s][root]*param;
                  Add(result,tmp);
              od;

              return result;#TransposedMat(result);
end);

InstallMethod(Adu,
              "Adu on the whole Lie algebra",
              [IsUnipotent],
              function(u)
              return Product(coefficients(u),i->AduAlphat(chevalleyAdj(u),i[1],i[2]));
end);

InstallMethod(AduJordanBlocksOp,
              "Computes Jordan Blocks of Adu using SageMath",
              [IsUnipotent],
              function(u)
              local sage_script,current_dir,tmp_file,result,
              ad,dim,i,j;

              sage_script:=Concatenation("~/sage/sage-4.7.2-linux-32bit-ubuntu_10.04_lts-i686-Linux/sage ",
                                      "~/workspace/ChevalleyAdj/sage/tmp/jordan.sage");
              current_dir:=Directory("~/workspace/ChevalleyAdj/sage/tmp/");
              tmp_file:=Filename(current_dir,"matrix.sage");
              result:=Filename(current_dir,"blocks.gap");

              #while IsPolynomialRing(ring(chevalleyAdj(u))) do u:=Descend(u,CoefficientsRing(ring(chevalleyAdj(u)))); od;
              if IsAlgebraicU(chevalleyAdj(u)) then u:=Descend(u,chevalleyAdj(chevalleyAdj(u))); fi;
              ad:=Adu(u);
              ad:=List(ad,i->List(i,j->Int(j)));
              PrintTo(tmp_file,"p=",String(Characteristic(chevalleyAdj(u))),"\nmat=");
              AppendTo(tmp_file,"[");
              dim:=Length(ad);
              for i in [1..dim] do
                  AppendTo(tmp_file,"[");
                  for j in [1..dim] do
                      AppendTo(tmp_file,String(ad[i][j]));
                          if j<>dim then AppendTo(tmp_file,","); fi;
                  od;
                  AppendTo(tmp_file,"]");
                  if i<>dim then AppendTo(tmp_file,",\n"); fi;
              od;
              AppendTo(tmp_file,"]");

              Exec(sage_script);

              result:=ReadAsFunction(result);
              result:=result();

              return result;
end);

InstallMethod(PositiveAdu,
              "Ad(u) restricted to Lie(U) where U is spanned by positive root groups",
              [IsUnipotent],
              function(u)
              local result,
              sys,
              baza,b,pr,pr_len;

              sys:=chevalleyAdj(u);
              pr:=positiveRoots(sys);
              pr_len:=Length(pr);
              baza:=List([1..pr_len],i->NilpotentChv(sys,[[i,1]]));

              result:=[];
              for b in baza do
                  result:=Concatenation(result,[u*b]);
              od;

              return List(result,i->LieAlgebraCoeffs(i){[1..pr_len]});
end);

InstallMethod(RestrictedAdu,
              "Ad(u) restricted to a subspace of Lie(U) spanned given (positive) roots",
              [IsUnipotent,IsList],
              function(u,roots)
              local sys,adu,pr,
              croots,check;

              if not IsInt(roots[1]) then TryNextMethod(); fi;

              sys:=chevalleyAdj(u);
              adu:=PositiveAdu(u);
              pr:=positiveRoots(sys);
              croots:=Filtered([1..Length(pr)],i->not i in roots);

              check:=Filtered(Concatenation(List(croots,i->adu[i]{roots})),j->j<>Zero(ring(sys)));
              if check<>[] then Error("Ade cannot be restricted like this."); fi;
              check:=Filtered(Concatenation(List(roots,i->adu[i]{croots})),j->j<>Zero(ring(sys)));
              if check<>[] then Error("Ade cannot be restricted like this."); fi;

              return List(roots,i->adu[i]{roots});
end);

InstallMethod(RestrictedAdu,
              "Ad(u) restricted to a subspace of Lie(U) spanned given (positive) roots",
              [IsUnipotent,IsList],
              1,
              function(u,vectors)
              local result,vectors_len,
              adu,img,
              CPR,cvar_names,cvars,cvectors,sysCPR,uu,
              tmp,generic,relations,relations_len,sol,i;

              vectors_len:=Length(vectors);

              cvar_names:=List([1..vectors_len],i->Concatenation("c_",String(i)));
              CPR:=PolynomialRing(ring(chevalleyAdj(u)),cvar_names);
              cvars:=IndeterminatesOfPolynomialRing(CPR);
              sysCPR:=ChevalleyAdj(chevalleyAdj(u),CPR);

              cvectors:=List(vectors,i->Ascend(i,sysCPR));
              uu:=Ascend(u,sysCPR);
              img:=List(cvectors,i->uu*i);

              result:=[];
              for i in [1..vectors_len] do
                  tmp:=ShallowCopy(cvars);
                  generic:=Sum(List([1..vectors_len],i->tmp[i]*cvectors[i]));
                  relations:=List(coefficients(img[i]-generic),i->i[2]);
                  relations:=Filtered(relations,i->i<>Zero(CPR));
              #Error("!");
                  tmp:=SolveRelations(relations,cvars,tmp,CPR);
                  if tmp[2]<>[] then
                      #Error("!");
                      Print("RestrictedAdu: did not solve all relations! ",relations,"\n");
                  fi;
                  tmp:=List(tmp[1],i->Value(i,cvars,cvars));
                  #tmp:=Sum(List([1..vectors_len],i->tmp[i]*vectors[i])); # I want the matrix
                  Add(result,tmp);
              od;

              return result;
end);

InstallMethod(PositiveFixedspace,
              "Ker of PositiveAdu",
              [IsUnipotent],
              function(u)
              local result,
              adu,sys,unu,pr,pr_len,ns;

              sys:=chevalleyAdj(u);
              pr:=positiveRoots(sys);
              pr_len:=Length(pr);

              adu:=PositiveAdu(u);
              unu:=DiagonalMat(List([1..pr_len],i->One(ring(sys))));
              ns:=NullspaceMat(adu-unu);

              result:=List(ns,v->List([1..pr_len],j->[j,v[j]]));
              result:=List(result,i->NilpotentChv(sys,i));

              return result;
end);

InstallMethod(RestrictedFixedspace,
              "Ker of RestrictedAdu",
              [IsUnipotent,IsList],
              function(u,roots)
              local result,
              radu,sys,unu,roots_len,ns;

              sys:=chevalleyAdj(u);
              roots_len:=Length(roots);

              radu:=RestrictedAdu(u,roots);
              unu:=DiagonalMat(List([1..roots_len],i->One(ring(sys))));
              ns:=NullspaceMat(radu-unu);

              result:=List(ns,v->List([1..roots_len],j->[roots[j],v[j]]));
              result:=List(result,i->NilpotentChv(sys,i));

              return result;
end);


text:
# Copyright (c) 2018-2021, Carnegie Mellon University
# See LICENSE for details

#P Permutations
#P ------------
#P
#P Under different circumstances different objects are called permutations,
#P even in the context of linear algrebra.
#P
#P Following list describes these objects and their representation in SPIRAL:
#P
#P  1. Regular permutations
#P     (a) in cycle notation : (1,2)(3,4)
#P     (b) as lists          : [2,1,4,3]
#P     (c) as index space mapping function:  i -> i+1 mod N
#P
#P  2. Parametrized permutation classes
#P     [For example stride permutations L(size,str) | size mod str = 0]
#P     Represented as a construction function, which returns an index mapping
#P     function - (c) above.
#P
#P  Conversions: ListPerm     (b) <- (a)   [gap.perm]
#P               PermList     (a) <- (b)   [gap.perm]
#P               PermFunc     (a) <- (c)   [spiral.spl.perm]
#P               ListPermFunc (b) <- (c)   [spiral.spl.perm]
#P
#P -------------------

#F PermFunc(<func>, <size>) . . . . . convert perm. function into explicit GAP perm.
#F
#F  PermFunc converts a 0-based permutation function used in SPLs, into an explicit
#F  GAP permutation. Recall that GAP permutations are 1-based, and are in cycle
#F  representation.
#F
PermFunc := (func, size) ->
    PermList( List([0..size-1], func) + 1 );

#F PermFunc(<func>, <size>) . . . . . convert perm. function into explicit list perm.
#F
#F  PermFunc converts a 0-based permutation function used in SPLs, into an explicit
#F  1-based list permutation. List permutations can be converted to GAP permutations
#F  using PermList. Alternatively PermFunc can be used, wihch returns GAP permutation.
#F
ListPermFunc := (func, size) ->
    List([0..size-1], func) + 1;

# ==========================================================================
# FuncClass
#
# Base class for symbolic functions
# ==========================================================================
Class(FuncClass, BaseMat, Function, rec(
    #-----------------------------------------------------------------------
    # Must be implemented in subclasses
    #-----------------------------------------------------------------------
    lambda := self >> Error("not implemented"),
    domain := self >> Error("not implemented"),
    range := self >> Error("not implemented"),

    #-----------------------------------------------------------------------
    _perm := true,
    isReal := True,
    isPermutation := self >> false,
    perm := self >> Checked(self.isPermutation(), PermFunc(x->self.lambda().at(x), self.range())),

    #-----------------------------------------------------------------------
    print   := (self,i,is) >> Print(
	self.name, "(", PrintCS(self.params), ")", When(self.transposed, ".transpose()")),
    #-----------------------------------------------------------------------
    equals := (self, other) >> ObjId(other) = ObjId(self) and self.params=other.params,
    dims := self >> [self.range(), self.domain()],

    advdims   := self >> [ [[ self.range() ]], [[ self.domain() ]] ],

    # size along each dimension, for a multidimensional range
    # for compatibility with ClassSPL, we wrap this into another list
    # (list of outputs, since ClassSPL can have >1 output)
    advrange  := self >> [[ self.range() ]],

    # size along each dimension, for a multidimensional domain
    # for compatibility with ClassSPL, we wrap this into another list
    # (list of outputs, since ClassSPL can have >1 output)
    advdomain := self >> [[ self.domain() ]],

    # dimensionality of range (1-d, 2-d, etc)
    advrangeDim := self >> Length(self.advrange()[1]),
    # dimensionality of domain (1-d, 2-d, etc)
    advdomainDim := self >> Length(self.advdomain()[1]),
    #-----------------------------------------------------------------------
    toAMat := self >> Gath(self).toAMat(),
    #-----------------------------------------------------------------------
    arithmeticCost := (self, costMul, costAddMul) >> costMul(0) - costMul(0),
    #-----------------------------------------------------------------------
    transpose := self >> CopyFields(self, rec(transposed := not self.transposed )),
    #-----------------------------------------------------------------------
    conjTranspose := self >> self.transpose(),
    #-----------------------------------------------------------------------
    normalizedArithCost := self >> 0,
    #-----------------------------------------------------------------------
    free := self >> Union(List(self.params, FreeVars)),
    #-----------------------------------------------------------------------
    # Rewrite rules support
    #
    from_rChildren := (self, rch) >> CopyFields(ApplyFunc(ObjId(self), rch),
	rec(transposed:=self.transposed)),

    # NOTE: self.transposed not exposed
    rChildren := self >> self.params,

    rSetChild := meth(self, n, newChild)
        self.params[n] := newChild;
        # self.canonizeParams(); ??
	self.dimensions := self.dims();
    end,
    # ----------------------------------------------------------------------
    checkParams := meth(self, params)
        local nargs, nump;
	nargs := NumArgs(self.def);
	nump := Length(params);
	if nargs <> -1 and nargs <> nump then
            Error(self.name, " needs ", NumArgs(self.def), " parameters: ",
		ParamsFunc(self.def), "\n");
	fi;
    end,
    #-----------------------------------------------------------------------
    canonizeParams := meth(self, params)
        local A, nump;
	nump := Length(params);
	if IsBound(self.abbrevs) then
            for A in self.abbrevs do
                if NumArgs(A) = -1 or NumArgs(A) = nump then
		    return ApplyFunc(A, params);
		fi;
	    od;
            return params;
	else
            return params;
	fi;
    end,
    #-----------------------------------------------------------------------
    __call__ := meth(arg)
        local result, self, params, lkup, h;
        self := arg[1];
        params := arg{[2..Length(arg)]};
        params := self.canonizeParams(params);
        self.checkParams(params);

        h := self.hash;
        if h<>false then
            lkup := h.objLookup(self, params);
            if lkup[1] <> false then return lkup[1]; fi;
        fi;

        result := SPL(WithBases(self, rec(params := params, transposed := false)));
        result := Inherit(result, ApplyFunc(result.def, params));
	result.dimensions := result.dims();

        if h<>false then return h.objAdd(result, lkup[2]);
        else return result;
        fi;
    end,

# obsolete
#    checkInverse := self >> Checked(self.isPermutation(),
#        PermFunc(self.direct, self.range()) * PermFunc(self.inverse, self.domain())),
));

Class(PermClass, FuncClass, rec(
    isPermutation := self >> true,
    toAMat := self >> Gath(self).toAMat(),
    equals := (self, other) >> ObjId(other) = ObjId(self) and self.params=other.params
));


text:
#
# Read("~/Workspace/groupsSB/epi/ext.gi");
#
# highest_weight:=[1,1];# A2: corresponding highest weight module is the Lie algebra
#

L:= SimpleLieAlgebra(type,rank,Rationals);
cb0:=CanonicalBasis(L);
V:= HighestWeightModule( L, highest_weight );
extA:=List([1..Dimension(V)],i->ExteriorPowerOfAlgebraModule( V, i ));
extA_basis:=List(extA,i->Basis(i));

ext_e:=function(k,e)
	local W,extb,result,v;
    W:=extA[k];
    extb:=extA_basis[k];

	result:=[];
	for v in extb do
		Append(result,[Coefficients(extb,e^v)]);
	od;
	result:=TransposedMat(result);
	return result;
end;

ade0:=function(e)
	local result,v;
	result:=[];
	for v in cb0 do
		Append(result,[Coefficients(cb0,e*v)]);
	od;
	result:=TransposedMat(result);
	return result;
end;

ext_root_group:=function(k,index,t)
	local ee,tmp,result,i;
	ee:=ext_e(k,cb0[index]);
	tmp:=ee;
	result:=tmp^0;
	i:=1;
	while Length(Set(Concatenation(tmp)))<>1 do
		result:=result+t^i*tmp/Factorial(i);
		i:=i+1;
		tmp:=tmp*ee;
	od;
	return result;
end;

ext_root_nil_degree:=function(k,index)
	local ee,tmp,i;
	ee:=ext_e(k,cb0[index]);
	tmp:=ee;
	i:=1;
	while Length(Set(Concatenation(tmp)))<>1 do
		i:=i+1;
		tmp:=tmp*ee;
	od;
	return i;
end;

ext_lin_comb:=function(k,coeffs)
	local extAb,result;
	extAb:=extA_basis[k];
	result:=List([1..Length(extAb)],i->coeffs[i]*extAb[i]);
	return Sum(result);
end;

max_stars_repo_name: iuliansimion/groupsSB

max_issues_repo_name: iuliansimion/groupsSB

max_forks_repo_name: iuliansimion/groupsSB