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

Declare(CodeBlock, BRAMPermStreamOne, TPrmMulti);

Class(TPrmMulti, Tagged_tSPL_Container, rec(
    abbrevs :=  [ (tspl_list, which) -> [tspl_list, which] ],
    dims := self >> Cond(IsBound(self.params[1][1].domain), [self.params[1][1].domain(), self.params[1][1].range()], self.params[1][1].dims()),

    terminate := meth(self)
       local prms, res, len, p, i;
       prms := self.params[1];
       len := Length(prms);
       res := prms[len];
       for i in Reversed([1..len-1]) do
          p := Cond(IsBound(prms[i].terminate), prms[i].terminate(), Prm(prms[i]));
          res := COND(eq(self.params[2], (i-1)), p, res);
       od;

       return res.terminate();
    end,

    transpose := self >> TPrmMulti(List(self.params[1], i->i.transpose()), self.params[2]).withTags(self.getTags()),

    isReal := True,
    doNotMeasure := true,
    noCodelet := true,
    normalizedArithCost := self >> 0
));

#F PermutationModMatrix ( <perm>, <n>, <m> )
#F   constructs for the permutation <perm> on <n> points an
#F   <m> x <m> matrix M. The (i,j) entry of M contains the number
#F   of pairs (k, k^<perm>) with k mod m = i and k^{perm> mod m = j.
#F   <m> has to divide <n>.
#F   We view <perm> as a permutation of the points {0,.., n-1} even
#F   though in gap it is represented as permutation on {1,..,n}
#F   (i.e shifted by 1).
#F

PermutationModMatrix := function ( p, n, m)
  local M, i;

  # error checking
  if not IsPerm(p) then
    Error("<p> is not a permutation");
  fi;
  if not ( IsInt(n) and IsInt(m) ) then
    Error("<n> and <m> must be integers");
  fi;
  if not n mod m = 0 then
    Error("<m> must divide <n>");
  fi;
  if not ( p = () or LargestMovedPointPerm(p) <= n ) then
    Error("<n> is not a valid degree for <perm>");
  fi;

  # initialize m x m matrix M
  M := List([1..m], r -> List([1..m], c -> 0));

  # construct M
  for i in [1..n] do
    M[(i-1) mod m + 1][(i^p-1) mod m + 1] :=
      M[(i-1) mod m + 1][(i^p-1) mod m + 1] + 1;
  od;

  return M;
end;

## ============================================================
##
##  PermMatrixRowVals(p) takes permutation matrix p and returns a list
##    of the '1' location in each row of the matrix.
##
##  Thus, the value a in list element b means that the input in
##    location b will be permuted to location a in the output.
##
## ============================================================
PermMatrixRowVals := meth(x)

   return List(Prm(FPerm(x)).func.lambda().tolist(), i->i.ev());
end;

## ============================================================
##
##  PermMatrixRowVals(p) takes permutation matrix p and returns a list
##    of the '1' location in each row of the matrix.
##
##  Thus, the value a in list element b means that the input in
##    location b will be permuted to location a in the output.
##
## ============================================================
PermMatrixRowValsOld := 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;

## ============================================================
##
##  min_list(l) returns the smallest value in l, a list of lists.
##
## ============================================================
min_list := meth(l)
   local t, i;
   t := [];
   for i in [1 .. Length(l)] do
      t[i] := Minimum(l[i]);
   od;
   return Minimum(t);
end;

## ============================================================
##
##  max_list(l) returns the largest value in l, a list of lists.
##
## ============================================================
max_list := meth(l)
   local t, i;
   t := [];
   for i in [1 .. Length(l)] do
      t[i] := Maximum(l[i]);
   od;
   return Maximum(t);
end;

## ============================================================
##
##  RandomPerm(n) returns a random permutation on n points.
##
##  RandomPermSPL(n) is not used becuase 1/3 of the time it
##  returns a stride permutation, and 1/3 of the time it
##  returns J(n).
##
## ============================================================
RandomPerm := meth(n)
   return Perm(Random(SymmetricGroup(n)), n);
end;

## ============================================================
##
##   IsPermMult(p) returns 1 if matrix p is a permutation of a
##   constant multiple of a permutation, and 0 otherwise.
##
## ============================================================

IsPermMult := meth(p)
   local n, i, j, total;
   n := Length(p);
   for i in [1..n] do
      total := 0;
      for j in [1..n] do
         if (p[i][j] <> 0) then
            total := total+1;
	    if (total > 1) then
		return 0;
	    fi;
         fi;
      od;
   od;
   return 1;
end;

## ============================================================
##   OnesVector(n)
##   Returns a vector of 1s of length n
## ============================================================
OnesVector := meth(n)
   local res, i;
   res := [];
   for i in [1..n] do
      res[i] := 1;
   od;
   return res;
end;

## ============================================================
##  Given the mux addresses, i.e., the multiplexor settings
##  that say "given an output location, where does it get its
##  input from on each cycle," determine the inverse, i.e.,
##  "given an input location, determine where its output will
##  go for each cycle.
##
##  This function is intended as a helper function, used in
##  StreamPermGen.createCode().
## ============================================================
PermHelperFindPortAddresses := meth(muxaddr)
    local ports, i, row, j, rowres;
    ports := [];
    for i in [1..Length(muxaddr)] do
       rowres := [];
       row := muxaddr[i];
       for j in [1..Length(row)] do
          rowres[row[j]+1] := j-1;
       od;
       Append(ports, [rowres]);
    od;
    return ports;
end;

_optsInColumn := meth(col, rows_left)
    local res, i, size;
    size := Length(col);
    res := [];
    for i in [1..size] do
        if ((rows_left[i] = 1) and (col[i] > 0)) then
            Append(res, [i]);
        fi;
    od;
    return res;
end;

# Used in debugging.
#_findPermBacktracks := 0;

try_column := meth(whichCol, board, rows_left)
    local size, col, options, i, spec_rows_left, tres;

    size := Length(board[1]);

    col := TransposedMat(board)[whichCol];

    options := _optsInColumn(col, rows_left);

    for i in options do
        if (whichCol = size) then
            return [i];
        fi;

        rows_left[i] := 0;
        tres := try_column(whichCol+1, board, rows_left);
        if (tres <> [-1]) then
            Append(tres, [i]); # backwards
            return tres;
        fi;
        rows_left[i] := 1;
    od;

    # Debugging
    #_findPermBacktracks := _findPermBacktracks+1;
    return [-1];

end;

_randomPermGood := meth(board)
    local rows_left, permVec, size, permMat, col, row, choices, failed, choice, i;
    size      := Length(board[1]);
    rows_left := List([1..size], i->1);
    permVec := [];

    failed := false;
    col := 1;
    while ((col <= size) and (not failed)) do

        # Collect list of choices in this column.
        choices := [];
        for row in [1..size] do
            if ((rows_left[row] <> 0) and (board[row][col] > 0)) then
                Append(choices, [row]);
            fi;
        od;

        if (Length(choices) = 0) then
            failed := true;
        else
            choice := Random(choices);
            rows_left[choice] := 0;
            Append(permVec, [choice]);
        fi;

        col := col + 1;
    od;

    if (failed) then
        # If we fail, then we can assume that the square contains at least one 0 entry.  So,
        # we can send back the all 1 matrix, and be guaranteed that P - [1] < 0 for some value.
        return(List([1..size], r -> List([1..size], c -> 1)));
    else
        permMat   := [];
        for i in [1..size] do
            permMat[i] := BasisVec(size, permVec[i]-1);
        od;

        permMat := TransposedMat(permMat);

        return permMat;
    fi;
end;

_findPerm := meth(board)
    local rows_left, permVec, size, permMat, i;
    size      := Length(board[1]);
    rows_left := List([1..size], i->1);
    permVec   := Reversed(try_column(1, board, rows_left));
    permMat   := [];

    for i in [1..size] do
        permMat[i] := BasisVec(size, permVec[i]-1);
    od;

    permMat := TransposedMat(permMat);

    return permMat;
end;

_minisat_path := "";

_findPermUsingSat := meth(board)
    local size, r, c, zeros, dir, i, zerofile, outfile, res, permmat;

    size := Length(board[1]);
    zeros := [];

    for r in [1..size] do
       for c in [1..size] do
           if (board[r][c] = 0) then
               Append(zeros, [ (r-1)*size + c ]);
           fi;
       od;
    od;

    dir := Concat("/tmp/spiral/", String(GetPid()));
    MakeDir(dir);
    zerofile := Concat(dir, "/zeros");
    outfile := Concat(dir, "/out");
    PrintTo(zerofile, "");

    for i in [1..Length(zeros)] do
        AppendTo(zerofile, zeros[i], "\n");
    od;

    if (_minisat_path = "") then
        Error("Error: path to SAT solver not set: paradigms.stream._minisat_path is not bound.");
    fi;

#    IntExec(Concat("/Users/pam/minisat/core/minisat ", String(size), " ", zerofile, " ", outfile));
    IntExec(Concat(_minisat_path, " ", String(size), " ", zerofile, " ", outfile));

    res := ReadVal(outfile);

    permmat := [];
    for i in [1..size] do
       permmat[i] := BasisVec(size, res[i]);
    od;

    return permmat;
end;

_SumList := meth(m)
    local res, i;
    res := 0;
    for i in [1..Length(m)] do
        res := res + m[i];
    od;
    return res;
end;

_whichMeth := 1;

Declare(routePerm);

## ============================================================
##     listofperms = FactorSemiMagicSquare(square)
##
##  Given a semi-magic square of size n, find a (non-unique)
##  decomposition into a sum of permutations.  Returns the
##  decomposition as a list of permutations (with each perm
##  being represented in the "row value" form).
##
##  Each semi-magic square can be decomposed into at most
##  (n-1)^2 + 1 different permutations.  See:
##      Leep et al. Marriage, Magic, and Solitaire. The American
##      Mathematical Monthly (1999) vol. 106 (5) pp. 419-429.
##
##  This method uses randomly generated permutations.  A better
##  method would construct them based upon the square.  A more
##  realistic approach would be a hybrid of the two.
##
##  This function is intended as a helper function, used in
##  StreamPermGen.createCode().
## ============================================================
FactorSemiMagicSquare := meth(square)
   local length, multplr, rndprm, j, res, rows_left, cols_left, startTime, endTime, tries, which, ct;
   length := Length(square);
   res := [];

   # Debugging
   startTime := TimeInSecs();
   tries := 0;

   ct := 0;

   which := stream._whichMeth;

   ##  While the semi-magic square is not all zero...
   while (max_list(square) > 0) do

	  ## What is left?  Is it a multiple of a permutation matrix?
	  ## (This step not needed for correctness, but is a performance
	  ## optimization.)
	  if (IsPermMult(square) = 1) then

	      ## If what is left is a constant times a perm matrix,
	      ## set that perm as the random step, and complete the
	      ## final step.
	      rndprm := square / max_list(square);
	      multplr := 0;
	      while (min_list(square - multplr*rndprm) >= 0) do
	         multplr := multplr+1;
	      od;
	      multplr := multplr-1;
	      square := square - multplr*rndprm;
	      for j in [1..multplr] do
#! fails sometimes	         Append(res, [routePerm(PermMatrixRowValsOld(rndprm))]);
	         Append(res, [PermMatrixRowValsOld(rndprm)]);
	      od;

          ct := ct + 1;

      else

          ##  Generate a random permutation.
          # set 'which' as: paradigms.stream.whichMeth := 0;

#          rndprm := _randomPermGood(square);
          tries := tries+1;

#            if (which = 0) then
#                rndprm := MatSPL(RandomPerm(length));
#                #tries := tries + 1;
#            else
#                if (which = 1) then
#                    rndprm := _findPerm(square);
#                else
                    rndprm := _findPermUsingSat(square);
#                fi;
# 	       fi;
      fi;

      ##  How many times can this permutation be
      ##  subtracted from the square?
      multplr := 0;
      while (min_list(square - multplr*rndprm) >= 0) do
         multplr := multplr+1;
      od;
      multplr := multplr-1;

      ## If this permutation has been subtracted from the square,
      ## record the decision, and check if we are done.
      if (multplr > 0) then

          ## Subtract the permutation the correct number of times.
          square := square - multplr*rndprm;

     	  ## Record the decision
	      for j in [1..multplr] do
	          Append(res, [PermMatrixRowValsOld(rndprm)]);
	      od;

          endTime := TimeInSecs();
#          Print(which, ", ", (endTime - startTime), ", ", tries, "\n");
          startTime := endTime;
          tries := 0;
          ct := ct + 1;

      fi;
   od;

   return res;
end;

## ============================================================
##  row = PermHelperFindRow(permrows, x (readdport),
##                           y (writeport), streaming-width)
##
##  Given a permutation as a set of row values, and a selected
##  read and write port (x, y), find which permutation word can
##  be scheduled that will lead to a read from port x to a
##  write in port y.
##
##  As the algorithm progresses, words that have been scheduled
##  are removed from permrows and replaced with -1.
##
##  This function is intended as a helper function, used in
##  StreamPermGen.createCode().
## ============================================================
PermHelperFindRow := meth(p, x, y, w)
   local a, b;

   for b in [1 .. Length(p)] do
      a := p[b];

      # The a <> -1 isn't needed for correctness but it makes the
      # code more clear.
      if ((a <> -1) and (Mod(a,w) = x) and (Mod(b-1,w) = y)) then
         return b;
      fi;
   od;

   Error("Cannot find a word to schedule.");
end;

## ============================================================
##  [rd_addr, wr_addr] = PermHelperFindAddresses(permrows,
##                          muxaddr, w);
##
##  Given a permutation as a set of row values and a cycle-by-
##  cycle read/write port mapping (muxaddr), and streaming width
##  w, determine the assoicated memory read and write addresses.
##
##  This function is intended as a helper function, used in
##  StreamPermGen.createCode().
## ============================================================
PermHelperFindAddresses := meth(permrows, muxaddr, w)
   local rd_addr, wr_addr, i, x, y, b, a, rdcyc, wrcyc, wrcyc2, ports;

   ports := PermHelperFindPortAddresses(muxaddr);
   rd_addr := []; wr_addr := [];

   for i in [1..Length(permrows)/w] do
      rdcyc := []; wrcyc := []; wrcyc2 := [];

      for x in [1..w] do
         y := ports[i][x]+1;
	 b := PermHelperFindRow(permrows, x-1, y-1, w);
	 a := permrows[b];
	 Append(rdcyc, [floor(a/w).v]);
	 Append(wrcyc, [floor((b-1)/w).v]);
	 Append(wrcyc2, []);
	 permrows[b] := -1;
      od;
      Append(rd_addr, [rdcyc]);

      # Here, wrcyc tells the write address for each *read port* at this
      # time.  I need to permute these values so it tells the write address
      # for each *write port* at this time.
      # So, I just need to assign wrcyc to wrcyc2 based upon this mapping.
      for x in [1..w] do
         wrcyc2[x] := wrcyc[muxaddr[i][x]+1];
      od;

      Append(wr_addr, [wrcyc2]);
   od;
   return [rd_addr, wr_addr];
end;

## ==================================================================
##  BRAMPermGeneral([rdaddr], [swnetctrl], [wraddr], [perm], it)
##    where rdaddr is the set of read addresses, swnetctrl is the
##    control bits for the \Omega^{-1}\Omega network, wraddr is the
##    set of write addresses, and perm is the permutation as an SPL
##    object.  it is the iteration variable that specifies which perm
##    to perform.
##
##  This represents a streaming multi-permutation structure that has
##  been factored using our "general" method.
##
##  This is not ideal for permutations that are linear on the bits.
##
##  Do not attempt to build this object by hand. Instead use
##  StreamPerm(perm, streamsize) to build a generic streaming
##  perm.  Then, StreamPerm.createCode() will construct the
##  BRAMPermGeneral given the permutation and streaming width.
## ==================================================================
Class(BRAMPermGeneral, BaseMat, SumsBase, rec(
	abbrevs := [(r,m,w,p,it)-> [r,m,w,p,it]],

	new := (self, r, m, w, p, it) >> SPL(WithBases(self, rec(
		    _children  := [r, m, w, p, it],
		    dimensions := [ Length(r[1])*Length(r[1][1]), Length(r[1])*Length(r[1][1])],
		    streamSize := Length(r[1][1]),
	     ))),

	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, "(", self.dims()[1],
	    ", ", self.streamSize, ", ", self.child(1), ", ", self.child(2),
	    ", ", self.child(3), ", ", self.child(5), ")"),

    toAMat := meth(self)
        local prms, len, it, i, res;
        prms := self.child(4);
        len  := Length(prms);
        it   := self.child(5);
        res  := prms[len];

        for i in Reversed([1..len-1]) do
            res := COND(eq(it, i-1), prms[i], res);
        od;

        return res.toAMat();
    end,

	dims   := self >> [Length(self._children[1][1])*Length(self._children[1][1][1]),
	                   Length(self._children[1][1])*Length(self._children[1][1][1])],

	sums   := self >> self
));

PadToTwoPowSize := meth(muxaddr)
   local n, n2p, i, res, dif, line;
   n := Length(muxaddr[1]);

   if (2^Log2Int(n) = n) then
       return muxaddr;
   fi;

   n2p := 2^(Log2Int(n)+1);
   dif := n2p-n;

   res := [];

   for i in [1..Length(muxaddr)] do
       line := Concatenation(muxaddr[i], [n..n2p-1]);
       Append(res, [line]);
   od;

   return res;
end;

RandomVectorVector := meth(domain, inner, outer)
   local res, i, j;
   res := [];
   for i in [1..outer] do
      res[i] := [];
      for j in [1..inner] do
         res[i][j] := RandomList([0..domain-1]);
      od;
   od;
   return res;
end;

## ======================================================
##  A container for a streaming multi-permutation structure
##
##  StreamPerm([<SPL>], <m>, <w>, <i>)
##
##  where the permutations are I(<M>) x <SPL>, <SPL> is
##  the permutation as an SPL object, <w> is the
##  streaming width, and <i> selects between the <SPL>.
## ======================================================

Class(StreamPerm, BaseMat, SumsBase, rec(
	abbrevs := [(s, m, w, i) -> [s, m, w, i]],

	new := (self, prm, m, ss, i) >> SPL(WithBases(self, rec(
	     _children := [prm, m, ss, i],
	     streamSize := ss,
         par := m,
         it := i,
	     dimensions := Cond(IsBound(prm[1].domain), [m * prm[1].domain(), m * prm[1].range()], m * prm[1].dims()),
	     ))),

	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 := meth(self)
        local lin, prm, w, P, length, permobj, permrows, modmatrix, muxaddr, muxaddr2, ports, addrs,
              rdaddr, wraddr, network, res, p3, p1, r, k, M, N, P, x, i, j, l, K, size, count, ptrans,
              z, G, q1, q3, q3b, p1b, G, dm, plist, bitlist, prms, Pt, thisPerm, thisPrm, prmCode, it,
              prmTens, dm_no_tens, prm_no_tens;

        prm := List(self.child(1), p->Tensor(I(self.child(2)), p));
        prm_no_tens := self.child(1);
        w := self.child(3);
        it := self.child(4);

        P := [];

        # If the permutation is a composition of multiple perms, it is more efficient to
        # compute the bit representation of each individually, and then combine them.
        for i in [1..Length(self.child(1))] do
            thisPerm := self.child(1)[i];
            thisPrm := prm[i];

	    if (ObjId(thisPrm) = L) then
		thisPrm := TL(thisPerm.params[1], thisPerm.params[2]);
	    fi;

	    if (ObjId(thisPrm) = Tensor and ObjId(thisPrm.child(1)) = I and ObjId(thisPrm.child(2)) = L) then
		thisPrm := TL(thisPrm.child(2).params[1], thisPrm.child(2).params[2], thisPrm.child(1).params[1], 1);
	    fi;

            if (ObjId(thisPerm) = Compose) then
                plist := thisPerm.children();

		# If we just have an L, turn it into a TL so we can use its .permBits() function
		plist := List(plist, i -> Cond(ObjId(i) = L, TL(i.params[1], i.params[2]), i));

		# If we have I x L, turn it into TL
		plist := List(plist, i -> Cond(ObjId(i) = Tensor and ObjId(i.child(1)) = I and ObjId(i.child(2)) = L, TL(i.child(2).params[1], i.child(2).params[2], i.child(1).params[1], 1), i));

                bitlist := List(plist, i ->
		    Cond(IsBound(i.permBits),
			i.permBits(),
			Cond(ObjId(i) = Tensor and ObjId(i.child(1)) = I and IsBound(i.child(2).permBits),
				DirectSumMat(MatSPL(I(Log2Int(i.child(1).params[1])))*GF(2).one, i.child(2).permBits()),
				PermMatrixToBits(MatSPL(i))

			)
		    )
		);

                if (-1 in bitlist) then
                    Pt := -1;
                else
                    Pt := Product(bitlist);
                    if (self.child(2) > 1) then
                        Pt := DirectSumMat(MatSPL(I(Log2Int(self.child(2))))*GF(2).one, Pt);
                    fi;
                fi;
            else
		Pt := Cond(IsBound(thisPrm.permBits),
                    # If P.permBits defined, use it
		    thisPrm.permBits(),
	            # If P = Tensor(I, Q) and Q.permBits defined, use it.
		    Cond(ObjId(thisPrm) = Tensor and ObjId(thisPrm.child(1)) = I and IsBound(thisPrm.child(2).permBits),
			DirectSumMat(MatSPL(I(Log2Int(thisPrm.child(1).params[1])))*GF(2).one, thisPrm.child(2).permBits()),

			# Otherise, use PermMatrixToBits function (slow)
			PermMatrixToBits(MatSPL(thisPrm))
		    )
		);
            fi;

            Append(P, [Pt]);
        od;

        # so, Make a list of Ps, one for each perm.  then set lin if all are <> -1.
        lin := Cond(ForAny(P, i->i=-1) or 2^Log2Int(w) <> w or Length(prm) > 1, 0, 1);

        # Find the overall dimensions and the dimension of the "non-tensor" part
        dm := Cond(IsBound(prm[1].domain), prm[1].domain(), prm[1].dims()[1]);

	# PM: This is not very robust.
        dm_no_tens := Cond(IsBound(self.child(1)[1].domain), self.child(1)[1].domain(), self.child(1)[1].dims()[1]);

        # If the permutations are on w words, then we can generate code.
        if (self.child(3) = dm_no_tens) then

            ## If this perm is multiple permutation functions composed, separate them to make the
            ## code generation work as expected.

            prmCode := List(prm_no_tens, j ->
                Cond(ObjId(j) = Compose,
                    Compose(List(j.children(), i->Prm(FPerm(i)))),
                    Prm(FPerm(j))
                    )
                );

            res := prmCode[Length(prmCode)];

            for i in Reversed([1..Length(prmCode)-1]) do
            res := COND(eq(it, (i-1)), prmCode[i], res);
            od;

            if (self.child(2) > 1) then
                return STensor(CodeBlock(res).createCode(), self.child(2), self.child(3));
            fi;

            return CodeBlock(res).createCode();

        fi;

        # If our perm is not linear on the bits, or if we have multiple permutations, or
        # if we have manually overridden this to avoid the patented permutations:
        if ((lin = 0) or (stream.avoidPermPatent = 1)) then

            # We want to ignore the tensor product until the end
            prm := self.child(1);
            prmTens := List(self.child(1), p->Tensor(I(self.child(2)), p));

            it := self.child(4);

            if (self.child(3) = 1) then
                return BRAMPermStreamOne(
#                           List(prmTens, i->PermMatrixRowValsOld(MatSPL(i))),
                           List(prmTens, i->PermMatrixRowVals(i)),
                           prmTens,
                           it);
            fi;

	        length    := self._children[1][1].dims()[1];

	        permobj   := List(self._children[1], i -> PermSPL(i));
	        permrows  := List(permobj, i -> PermMatrixRowValsOld(MatPerm(i, length)));
	        modmatrix := List(permobj, i -> PermutationModMatrix(i, length, self.child(3)));
	        muxaddr   := List(modmatrix, i -> FactorSemiMagicSquare(i));

            muxaddr2  := List(muxaddr, i -> PadToTwoPowSize(i));

            network   := List(muxaddr2, j -> List([1..Length(j)], i -> routePerm(j[i])));

	        addrs     := List([1..Length(permobj)], i -> PermHelperFindAddresses(permrows[i], muxaddr[i], self.child(3)));
	        rdaddr    := List(addrs, i -> i[1]);
	        wraddr    := List(addrs, i -> i[2]);

            if (self.child(2) > 1) then
                return STensor(BRAMPermGeneral(rdaddr, network, wraddr, self._children[1], it), self.child(2), w);
            else
                return BRAMPermGeneral(rdaddr, network, wraddr, self._children[1], it);
            fi;

        else
            # If we are here, we only have one permutation
            P := P[1];

	        k := Log2Int(self.child(3));
	        size := Log2Int(self.dimensions[1]);

            # If our streaming width is 1, then we are done.  P = I*P.
	        if (self.child(3) = 1) then
                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 # This is algorithm 5.2 of J. ACM 2009 paper.

                     # Is P a permutation matrix?
                     if (CheckPerm(P) = 1) then
	                  j := NonZeroRows(p3);
	                  i := ZeroRows(p1);
	             else
		          p1b := TransposedMat(Copy(p1));
		          G := MyTriangulizeMat(p1b);
		          G := TransposedMat(G);
		          G := G * GF(2).one;
		          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 .. Length(q3)] do
 		               Append(q3b, [q3[l]{[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.child(3));
	        else
	            return STensor(BRAMPerm(N,M,self.child(3)), 2^count, self.child(3));
	        fi;
        fi;
    end,

    print := (self,i,is) >> Print(self.name, "(", self._children[1], ", ", self.child(2), ", ", self.child(3), ", ", self._children[4], ")"),

    toAMat := meth(self)
        local prms, len, it, i, res;
        prms := self._children[1];
        len  := Length(prms);
        it   := self._children[4];
        res  := Tensor(I(self.child(2)), prms[len]);

        for i in Reversed([1..len-1]) do
            res := COND(eq(it, (i-1)), Tensor(I(self.child(2)), prms[i]), res);
        od;

        return res.toAMat();
    end,

	dims := self >> self.child(2) * self._children[1][1].dims(),
	sums := self >> self

));

## ======================================================
##  A container for a streaming RC(permutation)
##
##  RCStreamPerm([<SPL>], <m>, <w>, <i>)
##
##  where the permutations is I(<M>) x [<SPL>] x I(2), [<SPL>] is
##  list of permutations as SPL objects, <w> is the
##  streaming width, and <i> selects which of the perms.
## ======================================================
Class(RCStreamPerm, BaseMat, SumsBase, rec(
	abbrevs := [(s, m, w, i) -> [s, m, w, i]],

	new := (self, prm, m, ss, i) >> SPL(WithBases(self, rec(
	     _children := [prm, m, ss, i],
	     streamSize := ss,
         par := m,
         it := i,
	     dimensions := Cond(IsBound(prm[1].domain), [2 * m * prm[1].domain(), 2 * m * prm[1].range()], 2 * m * prm[1].dims()),
	     ))),

	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 := meth(self)
        local r;
        r := StreamPerm(self.child(1), self.par, self.streamSize, self.it).createCode();
        if (ObjId(r) = STensor) then
            return STensor(RC(r.child(1)), r.p, r.bs*2);
        else
             return RC(r);
        fi;
    end,

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

	toAMat := self >> RC(StreamPerm(self._children[1], self.par, self.streamSize, self.it)).toAMat(),

	dims := self >> 2 * self.par * self._children[1][1].dims(),
	sums := self >> self

));

## ======================================================
##  A container for a general streaming permutation (not
##  necessarily linear on the bits)
##
##  StreamPermGen(<SPL>, <q>)
##
##  where <SPL> is the permutation as an SPL object, and
##  <q> is the streaming width.
## ======================================================
## NB: Deprecated.  Use StreamPerm, which covers all cases
# Class(StreamPermGen, BaseMat, SumsBase, rec(
# 	abbrevs := [(s, q) -> [s, q]],

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

# 	     dimensions := Cond(IsBound(prm.domain), [prm.domain(), prm.range()], 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 length, permobj, permrows, modmatrix, muxaddr, muxaddr2, ports,
#              addrs, rdaddr, wraddr, network, dms, prms;

#        dms := Cond(IsBound(self.child(1).domain), [self.child(1).domain(), self.child(1).range()], self.child(1).dims());

#        if (self.streamSize = 1) then
#            return BRAMPermStreamOne(
#                     [PermMatrixRowVals(MatSPL(self._children[1]))],
#                     [self._children[1]], 0);
#        fi;

#        if (self.streamSize = dms[1]) then

#            ## If this perm is multiple permutation functions composed, separate them to make the
#            ## code generation work as expected.
#            if (ObjId(self.child(1)) = Compose) then
#                prms := Compose(List(self.child(1).children(), i->Prm(FPerm(i))));
#            else
#                prms := Prm(FPerm(self.child(1)));
#            fi;

#            return CodeBlock(prms).createCode();

#        fi;

# 	   length    := self._children[1].dims()[1];
# 	   permobj   := PermSPL(self._children[1]);
# 	   permrows  := PermMatrixRowVals(MatPerm(permobj, length));
# 	   modmatrix := PermutationModMatrix(permobj, length, self.streamSize);
# 	   muxaddr   := FactorSemiMagicSquare(modmatrix);

#        muxaddr2  := PadToTwoPowSize(muxaddr);

# #       Print("muxaddr = ", muxaddr, "\nmuxaddr2 = ", muxaddr2);

#        network   := List([1..Length(muxaddr2)], i -> routePerm(muxaddr2[i]));

# 	   addrs     := PermHelperFindAddresses(permrows, muxaddr,
# 	                   self.streamSize);
# 	   rdaddr    := addrs[1];
# 	   wraddr    := addrs[2];

# #!	   rdaddr := RandomVectorVector(length/self.streamSize, self.streamSize, length/self.streamSize);
# #!	   wraddr := RandomVectorVector(length/self.streamSize, self.streamSize, length/self.streamSize);
# #!	   muxaddr := RandomVectorVector(self.streamSize, self.streamSize, length/self.streamSize);

# # Previously, we used 'muxaddr.'
# #	   return BRAMPermGeneral(rdaddr, muxaddr, wraddr, self._children[1]);
# # Now I am switching to the omega^-1 omega network control bits
#        return BRAMPermGeneral(rdaddr, network, wraddr, self._children[1]);

#         end,

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

# 	toAMat := self >> self._children[1].toAMat(),
# 	dims := self >> self._children[1].dims(),
# 	sums := self >> self

# ));

## ======================================================
##  A container for a streaming multi-permutation structure
##
##  StreamPermGen(<l>, <i>, <q>)
##
##  where <l> is a list of permutations as SPL objects,
##  <i> is a variable that selects which permutation to
##  perform, and <q> is the streaming width.
## ======================================================
## Deprecated.  Use StreamPerm, which supports multiple perms (or will soon!)

##  !!!!!!!!!!! IMPORTANT !!!!!!!!
#     Do not delete this code until I have successfully rolled its multiple perm
#     support into StreamPermGen.
##  !!!!!!!!!!! IMPORTANT !!!!!!!!

# Class(StreamPermGenMult, BaseMat, SumsBase, rec(
# 	abbrevs := [(s, i, q) -> [s, i, q]],

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

# 	     dimensions := Cond(IsBound(prm[1].domain), [prm[1].domain(), prm[1].range()], prm[1].dims()),

# 	     ))),

# 	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 := meth(self)
# 	   local length, permobj, permrows, modmatrix, muxaddr, muxaddr2, ports,
#              addrs, rdaddr, wraddr, network, dms, prms, it, prmcode, i, res;

#        dms := Cond(IsBound(self.child(1)[1].domain), [self.child(1)[1].domain(), self.child(1)[1].range()], self.child(1)[1].dims());

#        it   := self.child(2);
#        prms := self.child(1);

#        if (self.streamSize = 1) then
#            return BRAMPermStreamOne(List(prms, i->PermMatrixRowVals(MatSPL(i))),
#                       prms,
#                       it);
#        fi;

#        if (self.streamSize = dms[1]) then

#            ## If this perm is multiple permutation functions composed, separate them to make the
#            ## code generation work as expected.

#            prmcode := List(self.child(1), j ->
#                         Cond(ObjId(j)=Compose,
#                            Compose(List(self.child(1).children(), i->Prm(FPerm(i)))),
#                            Prm(FPerm(j))
#                         )
#            );

#            res := prmcode[Length(prmcode)];

#            for i in Reversed([1..Length(prmcode)-1]) do
#               res := COND(eq(it, i), prmcode[i], res);
#            od;

#            return CodeBlock(res).createCode();

#        fi;

# 	   length    := self._children[1].dims()[1];
# 	   permobj   := PermSPL(self._children[1]);
# 	   permrows  := PermMatrixRowVals(MatPerm(permobj, length));
# 	   modmatrix := PermutationModMatrix(permobj, length, self.streamSize);
# 	   muxaddr   := FactorSemiMagicSquare(modmatrix);

#        muxaddr2  := PadToTwoPowSize(muxaddr);

#        network   := List([1..Length(muxaddr2)], i -> routePerm(muxaddr2[i]));

# 	   addrs     := PermHelperFindAddresses(permrows, muxaddr,
# 	                   self.streamSize);
# 	   rdaddr    := addrs[1];
# 	   wraddr    := addrs[2];

# # Previously, we used 'muxaddr.'
# #	   return BRAMPermGeneral(rdaddr, muxaddr, wraddr, self._children[1]);
# # Now I am switching to the omega^-1 omega network control bits
#        return BRAMPermGeneral(rdaddr, network, wraddr, self._children[1]);

#         end,

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

#     toAMat := meth(self)
#         local prms, len, it, i, res;

#          prms := self._children[1];
#          len  := Length(prms);
#          it   := self._children[2];
#          res  := prms[len];

#          for i in Reversed([1..len-1]) do
#              res := COND(eq(it, i), prms[i], res);
#          od;

#          return res.toAMat();
#     end,

# #    toAMat := self >> COND(eq(self._children[2], 0), self._children[1][1], self._children[1][2]).toAMat(),
# 	dims := self >> self._children[1][1].dims(),
# 	sums := self >> self

# ));

## ==================================================================
##   BRAMPermStreamOne([rd-addr, ...], [<SPL>, ... ], i)
##     A container for a streaming multi-permutation structure where width = 1.
##     This merits its own class because the hardware implementation
##     is greatly simplified.
##
##     The permutations are given as a list of vectors of read addresses, and
##     as a list of SPL objects.
##
##     i is the variable that selects which permutation to perform
## ==================================================================
Class(BRAMPermStreamOne, BaseMat, SumsBase, rec(
	abbrevs := [(r, p, i)-> [r, p, i]],

	new := (self, r, p, i) >> SPL(WithBases(self, rec(
		    _children  := [r, p, i],
		    dimensions := [ Length(r[1]), Length(r[1])]
	     ))),

	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, "(", self.child(1), ", ", self.child(3), ")"),

#	toAMat := self >> self.child(2).toAMat(),

    toAMat := meth(self)
        local prms, len, it, i, res;

         prms := self._children[2];
         len  := Length(prms);
         it   := self._children[3];
         res  := prms[len];

         for i in Reversed([1..len-1]) do
             res := COND(eq(it, i-1), prms[i], res);
         od;

         return res.toAMat();
    end,

	dims   := self >> [Length(self._children[1][1]),
	                   Length(self._children[1][1])],

	sums   := self >> self
));

Declare(InitStreamUnrollHw, SumStreamStrategy, StreamStrategy);

# CountPermSwitches(t)
# Returns number of switches needed for permutation t.
# We assume t is a TPrm().withTags([AStream(w)]);
CountPermSwitches := function (t)
  local i, M, M2, M2_t, N, N2, N2_t, n, k, r, s, s2, s3, s_m, s_n, opts;

  opts := InitStreamUnrollHw();
  r := RandomRuleTree(t, opts);
  s := SumsRuleTreeStrategy(r, SumStreamStrategy, opts);
  s2 := ApplyStrategy(s, StreamStrategy, UntilDone, opts);
  s3 := s2.createCode();

  k := Log2Int(s3.streamSize);
  M := s3.child(2);
  n := Length(M);

  # Extract M2.  First just get the bottom k rows.
  M2_t := Sublist(M, [n-k+1 .. n]);

  M2 := [];
  # Now get the first n-k columns
  for i in [1..k] do
     Append(M2, [Sublist(M2_t[i], [1..k])]);
  od;

  s_m := Rank(M2);

  N := s3.child(1);
  # Extract N2.  First just get the bottom k rows.
  N2_t := Sublist(N, [n-k+1 .. n]);

  N2 := [];
  # Now get the first n-k columns
  for i in [1..k] do
     Append(N2, [Sublist(N2_t[i], [1..k])]);
  od;

  s_n := Rank(N2);

  return (s_m + s_n) * s3.streamSize/2;
end;


text:
#
# Read("~/Workspace/Chevalley.gap/init.gi"); Read(Filename(test_dir,"_allTests.gi"));
#

LogTo(Filename(test_dir,"_allTests.log"));

#
# Polynomial stuff
#
Print("[Polynomial stuff]\n");
Read(Filename(test_dir,"poly.test.init.gi"));
Read(Filename(test_dir,"poly.test.gi"));

#
# Root system
#
Print("[Root system]\n");
Read(Filename(test_dir,"rsys.test.init.gi"));
Read(Filename(test_dir,"rsys.test.gi"));

#
# Parabolic system
#
Print("[Parabolic system]\n");
Read(Filename(test_dir,"psys.test.init.gi"));
Read(Filename(test_dir,"psys.test.gi"));

#
# Chevalley group stuff
#
Print("[Chevalley group stuff]\n");
Read(Filename(home_dir,"lib/chvadj.gd"));
Read(Filename(home_dir,"lib/chvadj.gi"));
Read(Filename(test_dir,"chvadj.test.gi"));

#
# Nilpotent elements
#
Print("[Nilpotent elements]\n");
Read(Filename(home_dir,"lib/nilchv.gd"));
Read(Filename(home_dir,"lib/nilchv.gi"));
Read(Filename(test_dir,"nilchv.test.gi"));

#
# Algebraic unipotent Sylow stuff
#
Print("[Algebraic unipotent Sylow stuff]\n");
Read(Filename(home_dir,"lib/algU.gd"));
Read(Filename(home_dir,"lib/algU.gi"));
Read(Filename(test_dir,"algU.test.gi"));

#
# Unipotent elements
#
Print("[Unipotent elements]\n");
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"));
Read(Filename(test_dir,"unichv.test.gi"));

#
# Arithmetics modulo p
#
Print("[Arithmetics modulo p]\n");
Read(Filename(home_dir,"lib/unimod.gd"));
Read(Filename(home_dir,"lib/unimod.gi"));
#Read(Filename(test_dir,"unimod.test.gi"));

#
# Witt groups
#
Print("[Witt groups]\n");
Read(Filename(home_dir,"lib/witt.gi"));
Read(Filename(test_dir,"witt.test.gi"));

#
# Unipotent elements over polynomials
#
Print("[Unipotent elements over polynomials]\n");
Read(Filename(home_dir,"lib/unialg.gi"));
Read(Filename(test_dir,"unialg.test.gi"));

ext: gi


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

EPS:=1E-14;
RoundEPS:=(v)->Cond((v>-EPS) and (v<EPS),0,v);
IsVal:=(x,v)->((x>(v-EPS)) and (x<(v+EPS)));
IsZero:=(x)->IsVal(x,0);
IsOne:=(x)->IsVal(x,1);
IsMinusOne:=(x)->IsVal(x,-1);

# Printing to C

Declare(CGen);

ISum.cprint := meth(self,i,is)
    Print("for( ", self.var.cprint(), " = ", self.domain, ") { \n");
    self.expr.cprint(i+is,i);
    Print(Blanks(i), "}\n");
end;

Compose.cprint := meth(self,i,is)
    local c;
    for c in self.children() do
	    Print(Blanks(i+is));
        c.cprint(i+is, is);
        Print(";\n");
    od;
    Print("\n");
end;

# -------------------
# TYPES
# -------------------

re.cprint := self >> Print(self.args[1].cprint(), ".r");
im.cprint := self >> Print(self.args[1].cprint(), ".i");

TComplex.ctype := "struct{double r,i;}";
TDouble.ctype  := "double";
TInt.ctype     := "int";

Typ.cdecl := (self, var) >> When(IsList(var),
    Print(self.ctype, " ", PrintCS(List(var, v->v.id))),
    Print(self.ctype, " ", var.id));

TArray.cdecl  := meth(self,var)
    local dims, elemtype;
    if IsList(var) then
	    DoForAll(var, v->Print(v.t.cdecl(v), "; "));
    elif not IsArray(self.t) then
        Print(self.t.ctype, Print(" ", var.id, "[", self.size, "]"));
    else
        dims := [];
        elemtype := self;
        while IsArray(elemtype) do
            Add(dims, elemtype.size);
            elemtype := elemtype.t;
        od;
        Print(elemtype.ctype, " ", var.id);
        DoForAll(dims, d->Print("[",d,"]"));
    fi;
end;

# -------------------
# EXPRESSIONS
# -------------------

_cprintcs := function ( lst )
    local  i;
    if Length(lst) = 0 then
	    return;
	fi;
    lst[1].cprint();
    for i  in [ 2 .. Length(lst) ]  do
        Print(", ", lst[i].cprint());
    od;
end;

ccprint := (v,t) -> v;

Command.cprint := (self,i,is) >> self.print(i,is);
Exp.cprint := self >> self.print();
Loc.cprint := self >> self.print();

Value.idx_cprint := self >> Print(self.v);

Value.cprint := self >> Cond(
    self.t = TComplex, let(c:=Complex(self.v),
        Print("{", ReComplex(c), ",", ImComplex(c), "}")),
    self.t = TDouble and IsCyc(self.v),
        Print("(", ccprint(ReComplex(Complex(self.v)), self.t), ")"),
    IsArray(self.t),
        Print("{",
		(self.v), "}"),
    self.v < 0,
        Print("(", ccprint(self.v, self.t), ")"),
    Print(ccprint(self.v, self.t)));

infix_print := function(lst, sep)
    local first, c;
    first := true;
    for c in lst do
        if first then
		    first := false;
        else
		    Print(sep);
        fi;
        c.cprint();
    od;
    return "";
end;

pinfix_print := function(lst, sep)
    local first, c;
    first := true;
    for c in lst do
        if first then
		    first := false;
        else
		    Print(sep);
        fi;
        if not IsBound(c.rChildren) or c.rChildren()=[] then
		    c.cprint();
        else
		    Print("(",c.cprint(), ")");
		fi;
    od;
    return "";
end;

add.idx_cprint := self >> Print( "(", infix_print(self.args, " + "), ")");
sub.idx_cprint := self >> Print( "(", infix_print(self.args, " - "), ")");
mul.idx_cprint := self >> Print( "", infix_print(self.args, "*"), "");
neg.idx_cprint := self >> Print("-(", self.args[1].cprint(), ")");

ExpCommand.cprint := (self,i,is) >>
    Print(Blanks(i), self.name, "(", infix_print(self.args, ", "), ");\n");

nth.cprint := self >> Print(self.loc.cprint(), "[", self.idx.cprint(), "]");
Exp.cprint := self >> Print(self.name, "(", infix_print(self.args, ", "), ")");
add.cprint := self >> Print( "(", infix_print(self.args, " + "), ")");
sub.cprint := self >> Print( "(", infix_print(self.args, " - "), ")");
mul.cprint := self >> Print( "", infix_print(self.args, "*"), "");
neg.cprint := self >> Print("-(", self.args[1].cprint(), ")");

div.cprint := self >> Print( "", pinfix_print(self.args, "/"), "");
idiv.cprint := self >> Print( "(", pinfix_print(self.args, "/"), ")");
imod.cprint := self >> Print( "(", pinfix_print(self.args, "%"), ")");
no_mod.cprint := self >> self.args[1].cprint();

cond.cprint := self >> Checked(Length(self.args)=3,
    Print("((", self.args[1].cprint(), ") ? (",
            self.args[2].cprint(), ") : (",
        self.args[3].cprint(), "))"));

leq.cprint := self >> Chain(
        Print("(1"),
        FoldL1(self.args, (a,b) ->
        Chain(Print(" && ", "(", a.cprint(), " <= ", b.cprint(), ")"), b)),
        Print(")"));

# -------------------
# COMMANDS
# -------------------

assign.cprint := meth(self, i, is)
    Print(Blanks(i), self.loc.cprint(), " = ", self.exp.cprint(), ";\n");
end;

chain.cprint := meth(self,i,is)
    local c, vars;
    Print(Blanks(i), "{", Blanks(is-1));
    vars := Set(
    List(Collect(self, [assign, @.target(var), @(1)]),
		e->e.loc));
    if false and Length(vars) > 0 then
	    Print("double ");
	    PrintCS(vars);
	    Print(";\n");
    else
	    Print("\n");
    fi;
    for c in self.cmds do
	    c.cprint(i+is, is);
    od;
    Print(Blanks(i), "}\n");
end;

data.cprint := meth(self,i,is)
    Print(Blanks(i), "{", Blanks(is-1), When(IsArray(self.value.t), "static ", ""),
        self.value.t.cdecl(self.var), " = ", self.value.cprint(), ";");
    Print("\n");
    self.cmd.cprint(i+is, is);
    Print(Blanks(i), "}\n");
end;

decl.cprint := meth(self,i,is)
    local v, first, types;
    Print(Blanks(i), "{", Blanks(is-1));
    if Length(self.vars) > 1 and Length(Set(List(self.vars, v->v.t)))=1 then
        # all types are the same
        self.vars[1].t.cdecl(self.vars);
        Print(";\n");
    else
        first := self.vars[1];
        if IsArray(first.t) then
		    Print("static ");
		fi;
        first.t.cdecl(first);
		Print(";\n");
        for v in Drop(self.vars, 1) do
            Print(Blanks(i+is));
            if IsArray(v.t) then
			    Print("static ");
			fi;
            v.t.cdecl(v);
            Print(";\n");
        od;
    fi;
    self.cmd.cprint(i+is, is);
    Print(Blanks(i), "}\n");
end;

loop.cprint_c99 := meth(self,i,is)
       local v, lo, hi;
       Constraint(IsRange(self.range));
       v := self.var;
	   lo := self.range[1];
	   hi := Last(self.range);
       Print(Blanks(i));
       Print("for(int ", v, " = ", lo, "; ", v, " <= ", hi, "; ", v, "++) \n");
       self.cmd.cprint(i+is, is);
end;

loop.cprint := meth(self,i,is)
       local v, lo, hi;
       Constraint(IsRange(self.range));
       v := self.var;
	   lo := self.range[1];
	   hi := Last(self.range);
       Print(Blanks(i));
       Print("{int ",v,"; for(", v, " = ", lo, "; ", v, " <= ", hi, "; ", v, "++) \n");
       self.cmd.cprint(i+is, is);
       Print("}");
end;

IF.cprint := meth(self,i,is)
       local v, lo, hi;
       Print(Blanks(i));
       Print("if(", self.cond.cprint(), ") {\n");
       self.then_cmd.cprint(i+is, is);
       Print(Blanks(i), "} else {\n");
       self.else_cmd.cprint(i+is, is);
       Print(Blanks(i), "}\n");
end;

CGenScalarANSI_C := function(subname, code)
   Constraint(IsCommand(code));
   Print(
    "extern void init_",subname," ();\n",
    "extern void ",subname," (",TDouble.ctype , " *, ",TDouble.ctype , " *);\n"
   );
   Print("void ", subname, "(", TDouble.ctype, " *Y, ", TDouble.ctype, " *X) {\n");
   code.cprint(4,2);
   Print("}\n");

   Print("void init_", subname, "() { }\n");
end;

CGenScalarIntel_C := function(subname, code)
   Constraint(IsCommand(code));
   Print(
    "extern void init_",subname," ();\n",
    "extern void ",subname," (",TDouble.ctype , " *, ",TDouble.ctype , " *);\n"
   );
   Print("void ", subname, "(", TDouble.ctype, " *restrict Y, ", TDouble.ctype, " *restrict X) {\n",
   "    __assume_aligned(X, 16)\n",
   "    __assume_aligned(Y, 16)\n");
   code.cprint(4,4);
   Print("}\n");

   Print("void init_", subname, "() { }\n");
end;

CGenScalarVectorC := function(subname, code)
   Constraint(IsCommand(code));
   Print(
    "extern void init_",subname," ();\n",
    "extern void ",subname," (",TDouble.ctype , " *, ",TDouble.ctype , " *);\n"
   );
   Print("void ", subname, "(__declspec(alignedvalue(16)) ", TDouble.ctype, " *restrict Y, __declspec(alignedvalue(16)) ", TDouble.ctype, " *restrict X) {\n");
   code.cprint(4,4);
   Print("}\n");

   Print("void init_", subname, "() { }\n");
end;

CGenScalarCPP := function(subname, code)
   Constraint(IsCommand(code));
   Print(
    "#ifdef __cplusplus\n",
    "extern \"C\" {\n",
    "#endif\n",
    "extern void init_",subname," ();\n",
    "extern void ",subname," (",TDouble.ctype , " *, ",TDouble.ctype , " *);\n",
    "#ifdef __cplusplus\n",
    "}\n",
    "#endif\n"
   );
   Print("void ", subname, "(", TDouble.ctype, " *Y, ", TDouble.ctype, " *X) {\n");
   code.cprint(4,4);
   Print("}\n");

   Print("void init_", subname, "() { }\n");
end;

CGen:=CGenScalarANSI_C;

# CGenDMP is set to a dummy-value here.
# Will be set to DMPUnparser.gen as soon as this is defined
# CGenDMP:=CGenScalarANSI_C;

CRuntimeTest := function ( testFunc, cmd, opts )
    local prog, result, cgenFunc;
    opts := Copy(opts);
    Constraint(IsCommand(cmd));
    Constraint(IsRec(opts));
    if not IsBound(opts.dataType) then
	    opts.dataType := "real";
	fi;
    opts := DeriveSPLOptions(cmd, opts);
    prog := ProgSPL(cmd, opts);
    prog.type := ProgInputType.TargetSource;

    cgenFunc := When(IsBound(opts.unparser), (f, cmd) -> opts.unparser.gen(f, cmd, opts), CGen);
    PrintTo(prog.file, cgenFunc(prog.sub_name, cmd));
    result := testFunc(cmd, opts, prog);

    if IsBool(result) and result=false then
        Error("Compiled code execution failed, offending program kept in '", prog.file, "'");
    else
        return result;
    fi;
end;

Declare(CodeRuleTreeOpts);
CVerifyRT := function(rt, opts)
    local c, m1,m2,nt;
    nt:=rt.node;
    c := CodeRuleTreeOpts(rt, opts);
    m1 := MatSPL(nt);
    m2 := CMatrix(c, opts);
    return Maximum(Flat(MapMat(m1-m2, i->AbsComplex(ComplexAny(i)))));
end;

CVerifyMMM:=(c,opts,t)->opts.profile.verify(c,opts,t);


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

ImportAll(paradigms.scratchpad);

doNormalWHT := function(arg)
   local N, t, rt, srt, code, a;

   N := When(Length(arg) >=1, arg[1], 2);

   t := WHT(N);
   rt := RandomRuleTree(t, SpiralDefaults);
   srt := SumsRuleTree(rt, SpiralDefaults);
   code := CodeSums(srt, SpiralDefaults);
   a := CMatrix(code, SpiralDefaults);

   return a;
end;

doScratchWHT := function(arg)
   local N, ls, opts, t, rt, srt, code, b;
   N := When(Length(arg) >=1, arg[1], 2);
   ls := When(Length(arg) >=2, 2^arg[2],2);

   opts := ScratchX86CMContext.getOpts(ls,1,1,N);
   t := WHT(opts.size).withTags(opts.tags);
   rt := RandomRuleTree(t, opts);
   srt := SumsRuleTree(rt, opts);
   code := CodeSums(srt, opts);
   b := CMatrix(code, opts);

   return b;
end;

doWHT := function(arg)
    local N, N1, i, j, a, b;

    SpiralDefaults.includes := ["\"scratchc.h\""];
    SpiralDefaults.profile.makeopts.CFLAGS := "-O2 -Wall -fomit-frame-pointer -msse4.1 -std=gnu99 -static";

    N := When(Length(arg) >= 1, arg[1], 2);

    for i in [2..N] do
        a := doNormalWHT(i);
        N1 := i - 1;
        for j in [1..N1] do
            b := doScratchWHT(i,j);
            PrintLine("Size of WHT:", 2^i);
            PrintLine("Scratch buffer size:", 2^j);
            PrintLine(a=b);
            PrintLine("------------------------------------------");
        od;
    od;
end;


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

#F SPLAMat( <amat> )
#F   returns an spl corresponding to <amat>.
#F

SPLAMat := function ( A )
  if not IsAMat(A) then
    Error("<A> must be an amat");
  elif A.char <> 0 then
    Error("<A> has char <> 0");
  fi;

  # catch identity case
  if IsIdentityMat(A) then
    return Checked(A.dimensions[1]=A.dimensions[2], I(A.dimensions[1]));
  fi;

  # dispatch on types
  if A.type = "perm" then

    return Perm(A.element, A.dimensions[1]);

  elif A.type = "mon" then

    # catch diagonal case
    if A.element.perm = () then
      return Diag(A.element.diag);
    fi;
    return Perm(A.element.perm, A.dimensions[1]) * Diag(A.element.diag);

  elif A.type = "mat" then

    # catch DFT case
    if IsBound(A.isDFT) and A.isDFT = true then
      return spiral.transforms.DFT(A.dimensions[1]);
    fi;

    # catch RDFT case
    if  IsBound(A.isRDFT) and A.isRDFT = true then
      if A.RDFTinfo[1] = 1 then
        if A.RDFTinfo[2] = "not transposed" then
          return spiral.transforms.SRDFT(A.dimensions[1]);
        else # transposed
          return spiral.transforms.SRDFT(A.dimensions[1]).transpose();
        fi;
      else # A.RDFTinfo[1] = 3
        if A.RDFTinfo[2] = "not transposed" then
          return spiral.transforms.SRDFT3(A.dimensions[1]);
        else # transposed
          return spiral.transforms.SRDFT3(A.dimensions[1]).transpose();
        fi;
      fi;
    fi;

    # catch rotation case
    if IsBound(A.isRotation) and A.isRotation = true then
      return Rot(A.angle);
    fi;

    return Mat(A.element);

  elif A.type = "scalarMultiple" then

    if IsIdentityMat(A.element) then
      return
        Diag(
          List([1..A.dimensions[1]], i -> A.scalar)
        );
    else
      return
        Compose(
          Diag(
            List([1..A.dimensions[2]], i -> A.scalar)
          ),
          SPLAMat(A.element)
        );
    fi;

  elif A.type = "product" then

    return Compose(List(A.factors, SPLAMat));

  elif A.type = "power" then

    # no negative exponents
    if A.exponent < 0 then
      Error("<A>.exponent must not be negative");
    fi;
    return Compose(Replicate(A.exponent, SPLAMat(A.element)));

  elif A.type = "conjugate" then

    return Conjugate(SPLAMat(A.element), SPLAMat(A.conjugation));

  elif A.type = "directSum" then

    return DirectSum(List(A.summands, SPLAMat));

  elif A.type = "tensorProduct" then

    return Tensor(List(A.factors, SPLAMat));

  else

    Error("unrecognized type of amat <A>");

  fi;
end;


text:
InstallMethod(ChevalleyAdj,
              "Simple algebraic of adjoint type",
              [IsString,IsPosInt,IsRing],
              function(Type,Rank,Ring)
              local object,test;
              object:=Objectify(NewType(NewFamily("ChevalleyAdjFamily"),
                                        IsAttributeStoringRep and
                                        IsChevalleyAdj),
                                rec());

              Settype(object,Type);
              Setrank(object,Rank);

              SetlieAlgebra(object,SimpleLieAlgebraTypeA_G(Type,Rank,Ring));
              SetrootSystem(object,RootSystem(SimpleLieAlgebraTypeA_G(Type,Rank,Integers)));
              SetpositiveRoots(object,positiveRoots(rootSystem(object)));
              SetallRoots(object,Concatenation(
                  positiveRoots(rootSystem(object)),
                  -positiveRoots(rootSystem(object))));
              SetweylGroup(object,WeylGroup(rootSystem(object)));

              Setring(object,Ring);

              SetCharacteristic(object,Characteristic(Ring));

              SetA(object,A_rs(object));
              SetN(object,N_rs(object));
              SetH(object,Eta(object));
              SetM(object,M_rsi(object));
              SetC(object,C_ijrs(object));

              SetName(object,Concatenation("<simple adjoint ",Type,String(Rank),
                                           " in characteristic ",String(Characteristic(Ring)),">"));

              return object;
end);

InstallMethod(ChevalleyAdj,
              "Simple algebraic of adjoint type",
              [IsChevalleyAdj,IsRing],
              function(sys,Ring)
              local object;
              object:=Objectify(NewType(NewFamily("ChevalleyAdjFamily"),
                                        IsAttributeStoringRep and
                                        IsChevalleyAdj),
                                rec());

              Settype(object,type(sys));
              Setrank(object,rank(sys));

              SetlieAlgebra(object,SimpleLieAlgebraTypeA_G(type(object),rank(object),Ring));
              SetrootSystem(object,rootSystem(sys));
              SetpositiveRoots(object,positiveRoots(sys));
              SetallRoots(object,allRoots(sys));
              SetweylGroup(object,weylGroup(sys));

              Setring(object,Ring);

              SetCharacteristic(object,Characteristic(Ring));

              SetA(object,A(sys));
              SetH(object,H(sys));
              SetN(object,N(sys));
              SetM(object,M(sys));
              SetC(object,C(sys));

              SetName(object,Concatenation("<simple adjoint ",type(object),String(rank(object)),
                                           " in characteristic ",String(Characteristic(object)),">"));

              return object;
end);

InstallMethod(A_rs,
              "Extended Cartan Matrix",
              [IsChevalleyAdj],
              function(sys)
              local A,
              root_len,L,B,cochars,
              r,s;

              root_len:=2*Length(positiveRoots(sys));
              L:=UnderlyingLieAlgebra(rootSystem(sys));
              cochars:=Cocharacters(L);
              cochars:=Concatenation(cochars,-cochars);
              B:=Basis(L);

              A:=NullMat(root_len,root_len,Integers);
              for r in [1..root_len] do
                  for s in [1..root_len] do
                      A[r][s]:=Sum(Coefficients(Basis(L),cochars[r]*B[s]));
                  od;
              od;

              return A;
end);

eta:=function(r,s,roots,NN)
    local p,q,sus,jos;

    p:=0;
    while Position(roots,-(p+1)*roots[r]+roots[s])<>fail do p:=p+1; od;
    q:=0;
    while Position(roots,(q+1)*roots[r]+roots[s])<>fail do q:=q+1; od;

    sus:=List([1..p],i->SignInt(NN[r][Position(roots,(i-1-p)*roots[r]+roots[s])]));
    if sus=[] then sus:=1; else sus:=Product(sus); fi;
    if sus=0 then sus:=1; fi;

    jos:=List([1..q],i->SignInt(NN[r][Position(roots,(i-1-p)*roots[r]+roots[s])]));
    if jos=[] then jos:=1; else jos:=Product(jos); fi;
    if jos=0 then jos:=1; fi;

    return (-1)^p*sus/jos;
end;

InstallMethod(Eta,
              "The matrix \eta with the signs (see Carter page 95)",
              [IsChevalleyAdj],
              function(sys)
              local H,
              NN,roots,roots_len,r,s;

              roots:=allRoots(sys);
              roots_len:=Length(roots);

              NN:=N_rs(sys);

              H:=NullMat(roots_len,roots_len,Integers);
#              for r in [1..roots_len/2] do
#                  H[r][r]:=-1;
#                  H[r][Position(roots,-roots[r])]:=-1;
#                  H[Position(roots,-roots[r])][r]:=-1;
#                  H[Position(roots,-roots[r])][Position(roots,-roots[r])]:=-1;
#              od;

              for r in [1..roots_len] do
                  for s in [1..roots_len] do
                      if r<>s and r<>Position(roots,-roots[s]) then H[r][s]:=eta(r,s,roots,NN);
                      else H[r][s]:=-1; fi;
                  od;
              od;

              return H;
end);

#InstallMethod(Eta,
#              "The matrix \eta with the signs (see Carter page 95)",
#              [IsChevalleyAdj],
#              function(sys)
#              local H,
#              A_rs,pr,roots,pr_len,root_refs,
#              r,s,
#              wr_on_s,flag;
#
#              A_rs:=A(sys);
#              pr:=positiveRoots(sys);
#              roots:=allRoots(sys);
#              pr_len:=Length(pr);
#              root_refs:=RootReflections(weylGroup(sys));
#
#              H:=NullMat(2*pr_len,2*pr_len,Integers);
#              for r in [1..pr_len] do
#                  H[r][r]:=-1;
#                  H[r][Position(roots,-roots[r])]:=-1;
#                  H[Position(roots,-roots[r])][r]:=-1;
#                  H[Position(roots,-roots[r])][Position(roots,-roots[r])]:=-1;
#              od;
#
#              for s in [1..2*pr_len] do
#                  flag:=false;
#                  for r in [1..2*pr_len] do
#                      wr_on_s:=Position(roots,roots[s]*root_refs[((r-1) mod pr_len) +1]);
#                      if H[r][s]<>0 and H[r][wr_on_s]=0 then
#                          H[r][wr_on_s]:=H[r][s]*(-1)^A_rs[r][s];
#                          H[r][Position(roots,-roots[wr_on_s])]:=H[r][wr_on_s];
#                          flag:=true;
#                      elif H[r][s]=0 and H[r][wr_on_s]<>0 then
#                          H[r][s]:=H[r][wr_on_s]*(-1)^A_rs[r][s];
#                          H[r][Position(roots,-roots[s])]:=H[r][s];
#                          flag:=true;
#                      fi;
#                  od;
#                  if not flag then
#                      for r in [1..2*pr_len] do
#                          wr_on_s:=Position(roots,roots[s]*root_refs[((r-1) mod pr_len)+1]);
#                          if H[r][s]=0 then
#                              H[r][s]:=1;
#                              H[r][Position(roots,-roots[s])]:=1;
#                              H[r][wr_on_s]:=(-1)^A_rs[r][s];
#                              H[r][Position(roots,-roots[wr_on_s])]:=H[r][wr_on_s];
#                          fi;
#                      od;
#                  fi;
#              od;
#
#              return H;
#end);

InstallMethod(N_rs,
              "Lie Algebra structure constants (see Carter page 52)",
              [IsChevalleyAdj],
              function(sys)
              local N,
              L,B,pr_len,
              r,s;

              L:=UnderlyingLieAlgebra(rootSystem(sys));
              B:=Basis(L);
              pr_len:=Length(positiveRoots(sys));

              N:=NullMat(2*pr_len,2*pr_len,Integers);
              for r in [1..pr_len*2] do
                  for s in [1..pr_len*2] do
                      if r=s+pr_len or s=r+pr_len then N[r][s]:=fail;
                      else N[r][s]:=Sum(Coefficients(B,B[r]*B[s])); fi;
                  od;
              od;

              return N;
end);

InstallMethod(M_rsi,
              "The M_rsi (see Certer page 61)",
              [IsChevalleyAdj],
              function(sys)
              local M,
              pr,pr_len,roots,n,
              i,j,k,tmp,poz;

              pr:=positiveRoots(sys);
              roots:=allRoots(sys);
              pr_len:=Length(pr);
              n:=N(sys);

              M:=List([1..pr_len*2],i->List([1..pr_len*2],j->[]));
              for i in [1..pr_len*2] do
                  for j in [1..pr_len*2] do
                      if i=j+pr_len or j=i+pr_len then M[i][j]:=fail;
                      else
                          tmp:=[n[i][j]];
                          for k in [2..4] do
                              poz:=Position(roots,(k-1)*roots[i]+roots[j]);
                              if poz<>fail then Add(tmp,n[i][poz]); fi;
                          od;
                          M[i][j]:=List([1..Length(tmp)],i->1/Factorial(i)*Product(tmp{[1..i]}));
                          M[i][j]:=Filtered(M[i][j],k->k<>0);
                      fi;
                  od;
              od;

              return M;
end);

InstallMethod(C_ijrs,
              "Algebraic group structure constants (see Carter page 76)",
              [IsChevalleyAdj],
              function(sys)
              local C,
              pr,pr_len,roots,m,
              r,s,poz;

              pr:=positiveRoots(sys);
              roots:=allRoots(sys);
              pr_len:=Length(pr);
              m:=M(sys);

              C:=List([1..pr_len*2],i->List([1..pr_len*2],j->[]));
              for r in [1..pr_len*2] do
                  for s in [1..pr_len*2] do
                    if r=s+pr_len or s=r+pr_len then C[r][s]:=fail;
                    else
                      #C_i1rs=M_rsi
                      Add(C[r][s],List([1..Length(m[r][s])],i->[i,1,m[r][s][i]]));
                      #C_1jrs=(-1)^j*M_srj
                      Add(C[r][s],List([2..Length(m[s][r])],j->[1,j,m[s][r][j]*(-1)^j]));#(1,1) apare mai sus
                      #C_32rs=1/3*M_(r+s)r2
                      poz:=Position(roots,3*roots[r]+2*roots[s]);
                      if poz <> fail then
                          poz:=Position(roots,roots[r]+roots[s]);
                          Add(C[r][s],[[3,2,(1/3)*m[poz][r][2]]]);
                      fi;
                      #C_23rs=1/3*M_(r+s)r2
                      poz:=Position(roots,2*roots[r]+3*roots[s]);
                      if poz <> fail then
                          poz:=Position(roots,roots[r]+roots[s]);
                          Add(C[r][s],[[2,3,-(2/3)*m[poz][s][2]]]);
                      fi;
                      C[r][s]:=Concatenation(C[r][s]);
                      Sort(C[r][s],function(a,b) return a[1]+a[2]<b[1]+b[2]; end);
                    fi;
                  od;
              od;

              return C;
end);

InstallMethod(PermutationSign,
              "Uses the \eta's to give the sign of the permutation on an e_alpha",
              [IsChevalleyAdj,IsList,IsPosInt],
              function(sys,perm,root_index)
              local result,
              eta,ref_mats,roots,
              i;

              eta:=H(sys);
              ref_mats:=RootReflections(weylGroup(sys));
              roots:=allRoots(sys);

              result:=1;
              for i in Reversed(perm) do  # !!! ATENTIE AICI .. conteaza daca e actiune stanga sau dreapta
              #for i in perm do
                  result:=result*eta[i][root_index];
                  root_index:=Position(roots,roots[root_index]*ref_mats[i]);
#                  root_index:=Position(roots,ref_mats[i]*roots[root_index]);
              od;

              return result;
end);

InstallMethod(PermutationSign,
              "Uses the \eta's to give the sign of the permutation on an e_alpha",
              [IsChevalleyAdj,IsList,IsList],
              function(sys,perm,root)
              return PermutationSign(sys,perm,Position(positiveRoots(sys),root));
end);

InstallMethod(InversePermutationSign,
              "Uses the \eta's to give the sign of the permutation on an e_alpha",
              [IsChevalleyAdj,IsList,IsPosInt],
              function(sys,perm,root_index)
              local result,
              eta,ref_mats,roots,
              i;

              eta:=H(sys);
              ref_mats:=RootReflections(weylGroup(sys));
              roots:=allRoots(sys);

              result:=1;
              #for i in Reversed(perm) do  # !!! ATENTIE AICI .. conteaza daca e actiune stanga sau dreapta
              for i in perm do
                  #result:=result*eta[i][root_index];
                  root_index:=Position(roots,roots[root_index]*ref_mats[i]);
                  result:=result*eta[i][root_index];
#                  root_index:=Position(roots,ref_mats[i]*roots[root_index]);
              od;

              return result;
end);

InstallMethod(InversePermutationSign,
              "Uses the \eta's to give the sign of the permutation on an e_alpha",
              [IsChevalleyAdj,IsList,IsList],
              function(sys,perm,root)
              return InversePermutationSign(sys,perm,Position(positiveRoots(sys),root));
end);


text:
Load(graph);
Import(graph);
conf := SparseDefaultConf;
opts := conf.getOpts(V(1));
opts.vec := true;
i := Ind(opts.symbol[1]);
front := var.fresh_t("frontier",TArray(TSemiring_Boolean(TInt), opts.symbol[1]));
fdataofs := FDataOfs(front, opts.symbol[1], 0);
md := MakeDiag(fdataofs);
eye := I(i);
sid := SUM(eye, md);
matmul := MatMul(opts.symbol[1],opts.symbol[1]);
spmv := matmul * sid;
spmv2 := Rewrite(spmv, RulesIDiag, opts);
push := Rewrite(spmv2, RulesMatMul, opts);
fdos := FDataOfs(front, opts.symbol[1], i);
rv := RowVec(fdos);
pull := IterVStack(i, opts.symbol[1], rv);
f := var.fresh_t("nnz", TInt);
push_pull := COND(lt(div(f, opts.symbol[1]), V(0.01)), PushBFSKernel(push), PullBFSKernel(pull));
rt := RandomRuleTree(last, opts);
srt := SumsRuleTree(rt, opts);
srt2 := Rewrite(Copy(srt), RulesSigSPMV2, opts);
push_pullfinal := Rewrite(Copy(srt2), RulesScatRow, opts);
cs := CodeSums(BFSKernel(last), opts);
PrintCode("dobfs", cs, opts);


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

unpacklo := (l1,l2,n,k) -> Flat(List([1..n/(2*k)], i-> [
        List([1..k], j->l1[(i-1)*k+j]),
        List([1..k], j->l2[(i-1)*k+j])
]));

unpackhi := (l1,l2,n,k) -> Flat(List([1..n/(2*k)], i -> [
        List([1..k], j->l1[n/2+(i-1)*k+j]),
        List([1..k], j->l2[n/2+(i-1)*k+j])
]));

# saturation in pack_semantic is not taken into account
pack_semantic := (l1, l2, n) -> Flat(List( [1..n], i -> [l1[i], l2[i]]));

sparams := (l,n) -> List([1..l], i->[1..n]);

shuffle := (in1, in2, p, n, k) -> Flat([
    List([1..n/(2*k)],     i->List([1..k], j->in1[(p[i]-1)*k+j])),
    List([n/(2*k)+1..n/k], i->List([1..k], j->in2[(p[i]-1)*k+j]))
]);

inverse_ushuffle := (inp, p, n) -> List([1..n], i->inp[p[i]]);

shufflehi := (in1, p, n, k) -> Concat(
    Sublist(in1, [1..n/2]),
    let(l := Sublist(in1, [n/2+1..n]), shuffle(l, l, p, n/2, k)));

shufflelo := (in1, p, n, k) -> Concat(
    let(l := Sublist(in1, [1..n/2]), shuffle(l, l, p, n/2, k)),
    Sublist(in1, [n/2+1..n]));

iclshuffle := p -> Print("_MM_SHUFFLE", When(Length(p)=2, "2", ""), "(", PrintCS(Reversed(p-1)), ")");

iclprintop := self >> Print(self.icl, "(", _vcprintcs(self.args), ")");

iperm4 := self >> Filtered(Cartesian(self.params()), i->i[1]<>i[2] and i[3] <> i[4]);

vtakehi := (v) -> Checked(IsValue(v) and IsVecT(v.t) and v.t.size>1,  TakeLast(v.v, Length(v.v)/2));
vtakelo := (v) -> Checked(IsValue(v) and IsVecT(v.t) and v.t.size>1,  Take(v.v, Length(v.v)/2));


text:
method vForm0 vForm0.mLoad <alias=vForm0_mLoad>( )
{
//	this->vForm.mCreateWin()
	ustr ustmp
	uint comp
	comp as this
	with comp
	{
		.AutoLang=1
		.Border=$fbrdDialog
		.Bottom=0
		.Caption=ustmp.fromutf8("Gentee Manager")
		.Enabled=1
		.FormStyle=$fsChild
		.Height=400
		.HelpTopic=ustmp.fromutf8("")
		.Hint=ustmp.fromutf8("")
		.HorzAlign=$alhLeft
		.IconName=ustmp.fromutf8("")
		.Left=0
		.Name="Form0"
		.Right=0
		.StartPos=$spScreenCenter
		.Style=""
		.Tag=0
		.Top=0
		.TopMost=0
		.VertAlign=$alvTop
		.Visible=1
		.Width=595
		.WindowState=$wsNormal
		.OnCreate.Set( this, Form0_appinit )

		uint comp
		comp as this.Tab0
		comp.Owner = this
		with comp
		{
			.AutoLang=1
			.Bottom=52
			.Enabled=1
			.Height=305
			.HelpTopic=ustmp.fromutf8("")
			.Hint=ustmp.fromutf8("")
			.HorzAlign=$alhLeftRight
			.Left=15
			.Name="Tab0"
			.Right=12
			.Style=""
			.TabOrder=0
			.Tag=0
			.Top=10
			.VertAlign=$alvTopBottom
			.Visible=1
			.Width=560

			uint comp
			comp as this.TabItem0
			comp.Owner = this.Tab0
			with comp
			{
				.AutoLang=1
				.Bottom=0
				.Caption=ustmp.fromutf8("Application")
				.Enabled=1
				.Height=276
				.HelpTopic=ustmp.fromutf8("")
				.Hint=ustmp.fromutf8("")
				.HorzAlign=$alhLeft
				.Left=4
				.Name="TabItem0"
				.Right=0
				.Style=""
				.Tag=0
				.Top=25
				.VertAlign=$alvTop
				.Visible=1
				.Width=552

				uint comp
				comp as this.Btn0
				comp.Owner = this.TabItem0
				with comp
				{
					.AutoLang=1
					.Bottom=0
					.BtnStyle=$bsClassic
					.Caption=ustmp.fromutf8("Gentee IDE (Debugger)")
					.Checked=0
					.Enabled=1
					.Height=35
					.HelpTopic=ustmp.fromutf8("")
					.Hint=ustmp.fromutf8("")
					.HorzAlign=$alhLeft
					.Left=25
					.Name="Btn0"
					.Right=0
					.Style=""
					.TabOrder=0
					.Tag=0
					.Top=30
					.VertAlign=$alvTop
					.Visible=1
					.Width=225
					.OnClick.Set( this, Form0_runide )
				}
				comp as this.Btn1
				comp.Owner = this.TabItem0
				with comp
				{
					.AutoLang=1
					.Bottom=0
					.BtnStyle=$bsClassic
					.Caption=ustmp.fromutf8("Gentee Studio")
					.Checked=0
					.Enabled=1
					.Height=35
					.HelpTopic=ustmp.fromutf8("")
					.Hint=ustmp.fromutf8("")
					.HorzAlign=$alhLeft
					.Left=25
					.Name="Btn1"
					.Right=0
					.Style=""
					.TabOrder=1
					.Tag=0
					.Top=80
					.VertAlign=$alvTop
					.Visible=1
					.Width=225
					.OnClick.Set( this, Form0_runstudio )
				}
				comp as this.Btn2
				comp.Owner = this.TabItem0
				with comp
				{
					.AutoLang=1
					.Bottom=0
					.BtnStyle=$bsClassic
					.Caption=ustmp.fromutf8("VisEdit (Demo)")
					.Checked=0
					.Enabled=1
					.Height=35
					.HelpTopic=ustmp.fromutf8("")
					.Hint=ustmp.fromutf8("")
					.HorzAlign=$alhLeft
					.Left=25
					.Name="Btn2"
					.Right=0
					.Style=""
					.TabOrder=2
					.Tag=0
					.Top=180
					.VertAlign=$alvTop
					.Visible=1
					.Width=225
					.OnClick.Set( this, Form0_runvis )
				}
				comp as this.Panel0
				comp.Owner = this.TabItem0
				with comp
				{
					.AutoLang=1
					.Border=$brdGroupBox
					.Bottom=0
					.Caption=ustmp.fromutf8("Documentation")
					.Enabled=1
					.Height=160
					.HelpTopic=ustmp.fromutf8("")
					.Hint=ustmp.fromutf8("")
					.HorzAlign=$alhLeft
					.Left=290
					.Name="Panel0"
					.Right=0
					.Style=""
					.TabOrder=3
					.Tag=0
					.Top=20
					.VertAlign=$alvTop
					.Visible=1
					.Width=250

					uint comp
					comp as this.Btn5
					comp.Owner = this.Panel0
					with comp
					{
						.AutoLang=1
						.Bottom=0
						.BtnStyle=$bsClassic
						.Caption=ustmp.fromutf8("Help v2.5")
						.Checked=0
						.Enabled=1
						.Height=35
						.HelpTopic=ustmp.fromutf8("")
						.Hint=ustmp.fromutf8("")
						.HorzAlign=$alhLeft
						.Left=20
						.Name="Btn5"
						.Right=0
						.Style=""
						.TabOrder=0
						.Tag=0
						.Top=95
						.VertAlign=$alvTop
						.Visible=1
						.Width=205
						.OnClick.Set( this, Form0_openchm25 )
					}
					comp as this.Btn3
					comp.Owner = this.Panel0
					with comp
					{
						.AutoLang=1
						.Bottom=0
						.BtnStyle=$bsClassic
						.Caption=ustmp.fromutf8("Help v3")
						.Checked=0
						.Enabled=1
						.Height=35
						.HelpTopic=ustmp.fromutf8("")
						.Hint=ustmp.fromutf8("")
						.HorzAlign=$alhLeft
						.Left=20
						.Name="Btn3"
						.Right=0
						.Style=""
						.TabOrder=1
						.Tag=0
						.Top=35
						.VertAlign=$alvTop
						.Visible=1
						.Width=205
						.OnClick.Set( this, Form0_openchm3 )
					}
				}
				comp as this.Btn8
				comp.Owner = this.TabItem0
				with comp
				{
					.AutoLang=1
					.Bottom=0
					.BtnStyle=$bsClassic
					.Caption=ustmp.fromutf8("GeViewer")
					.Checked=0
					.Enabled=1
					.Height=35
					.HelpTopic=ustmp.fromutf8("")
					.Hint=ustmp.fromutf8("")
					.HorzAlign=$alhLeft
					.Left=25
					.Name="Btn8"
					.Right=0
					.Style=""
					.TabOrder=4
					.Tag=0
					.Top=130
					.VertAlign=$alvTop
					.Visible=1
					.Width=225
					.OnClick.Set( this, Form0_runviewer )
				}
			}
			comp as this.TabItem1
			comp.Owner = this.Tab0
			with comp
			{
				.AutoLang=1
				.Bottom=0
				.Caption=ustmp.fromutf8("Associate Ext.")
				.Enabled=1
				.Height=276
				.HelpTopic=ustmp.fromutf8("")
				.Hint=ustmp.fromutf8("")
				.HorzAlign=$alhLeft
				.Left=4
				.Name="TabItem1"
				.Right=0
				.Style=""
				.Tag=0
				.Top=25
				.VertAlign=$alvTop
				.Visible=0
				.Width=552

				uint comp
				comp as this.Label0
				comp.Owner = this.TabItem1
				with comp
				{
					.AutoLang=1
					.Bottom=0
					.Caption=ustmp.fromutf8("Current Gentee compiler:")
					.Enabled=1
					.Height=25
					.HelpTopic=ustmp.fromutf8("")
					.Hint=ustmp.fromutf8("")
					.HorzAlign=$alhLeft
					.Left=20
					.Name="Label0"
					.Right=0
					.Style=""
					.Tag=0
					.TextHorzAlign=$talhLeft
					.TextVertAlign=$talvTop
					.Top=20
					.VertAlign=$alvTop
					.Visible=1
					.Width=680
					.WordWrap=0
					.AutoSize=0
				}
				comp as this.vRegcur
				comp.Owner = this.TabItem1
				with comp
				{
					.AutoLang=1
					.Bottom=0
					.Enabled=0
					.Height=30
					.HelpTopic=ustmp.fromutf8("")
					.Hint=ustmp.fromutf8("")
					.HorzAlign=$alhLeft
					.Left=20
					.MaxLen=32768
					.Multiline=0
					.Name="vRegcur"
					.Password=0
					.ReadOnly=0
					.Right=0
					.ScrollBars=$sbNone
					.Style=""
					.TabOrder=1
					.Tag=0
					.Text=ustmp.fromutf8("")
					.Top=45
					.VertAlign=$alvTop
					.Visible=1
					.Width=480
					.WordWrap=0
					.Border=1
				}
				comp as this.vRegnew
				comp.Owner = this.TabItem1
				with comp
				{
					.AutoLang=1
					.Bottom=0
					.Enabled=1
					.Height=30
					.HelpTopic=ustmp.fromutf8("")
					.Hint=ustmp.fromutf8("")
					.HorzAlign=$alhLeft
					.Left=20
					.MaxLen=32768
					.Multiline=0
					.Name="vRegnew"
					.Password=0
					.ReadOnly=0
					.Right=0
					.ScrollBars=$sbNone
					.Style=""
					.TabOrder=2
					.Tag=0
					.Text=ustmp.fromutf8("")
					.Top=95
					.VertAlign=$alvTop
					.Visible=1
					.Width=480
					.WordWrap=0
					.Border=1
				}
				comp as this.Btn6
				comp.Owner = this.TabItem1
				with comp
				{
					.AutoLang=1
					.Bottom=0
					.BtnStyle=$bsClassic
					.Caption=ustmp.fromutf8("Change .G command line")
					.Checked=0
					.Enabled=1
					.Height=30
					.HelpTopic=ustmp.fromutf8("")
					.Hint=ustmp.fromutf8("")
					.HorzAlign=$alhLeft
					.Left=20
					.Name="Btn6"
					.Right=0
					.Style=""
					.TabOrder=3
					.Tag=0
					.Top=145
					.VertAlign=$alvTop
					.Visible=1
					.Width=225
					.OnClick.Set( this, Form0_changegline )
				}
				comp as this.Btn7
				comp.Owner = this.TabItem1
				with comp
				{
					.AutoLang=1
					.Bottom=0
					.BtnStyle=$bsClassic
					.Caption=ustmp.fromutf8("Set default command line")
					.Checked=0
					.Enabled=1
					.Height=35
					.HelpTopic=ustmp.fromutf8("")
					.Hint=ustmp.fromutf8("")
					.HorzAlign=$alhLeft
					.Left=20
					.Name="Btn7"
					.Right=0
					.Style=""
					.TabOrder=4
					.Tag=0
					.Top=210
					.VertAlign=$alvTop
					.Visible=1
					.Width=225
					.OnClick.Set( this, Form0_setdefault )
				}
			}
			.CurIndex=0
		}
		comp as this.Btn4
		comp.Owner = this
		with comp
		{
			.AutoLang=1
			.Bottom=0
			.BtnStyle=$bsClassic
			.Caption=ustmp.fromutf8("Exit")
			.Checked=0
			.Enabled=1
			.Height=30
			.HelpTopic=ustmp.fromutf8("")
			.Hint=ustmp.fromutf8("")
			.HorzAlign=$alhLeft
			.Left=225
			.Name="Btn4"
			.Right=0
			.Style=""
			.TabOrder=1
			.Tag=0
			.Top=325
			.VertAlign=$alvTop
			.Visible=1
			.Width=130
			.OnClick.Set( this, Form0_exit )
		}
	}

	return this
}

method vForm0 vForm0.init( )
{
   this.pTypeId = vForm0
   return this
}
func init_vForm0 <entry>()
{

   regcomp( vForm0, "vForm0", vForm, $vForm_last,
      %{ %{$mLoad,     vForm0_mLoad}},
      0->collection )

}

max_stars_repo_name: gentee/c-gentee

max_issues_repo_name: gentee/c-gentee

max_forks_repo_name: gentee/c-gentee


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

####################################
##
##  Attributes of polyopes
##
####################################

InstallMethod( ExternalPolymakePolytope,
               "for polyopes",
               [ IsPolytope ],
    function( poly )

    if IsBound( poly!.input_points ) then
        return Polymake_PolytopeByGenerators( poly!.input_points );
    fi;

    # otherwise our fallback is poly from inequalities
    return Polymake_PolytopeFromInequalities( poly!.input_ineqs );

end );

InstallMethod( VerticesOfPolytope,
               "for polyopes",
               [ IsPolytope ],
  function( poly )

    return Set( Polymake_V_Rep( ExternalPolymakePolytope( poly ) )!.vertices );

end );

InstallMethod( DefiningInequalities,
               "for polyopes",
               [ IsPolytope ],
  function( poly )
    local ineqs, eqs;

    ineqs := Polymake_Inequalities( ExternalPolymakePolytope( poly ) );
    eqs := Polymake_Equalities( ExternalPolymakePolytope( poly ) );
    return Set( Concatenation( eqs, (-1) * eqs, ineqs ) );

end );

InstallMethod( FacetInequalities,
               " for external polyopes",
               [ IsExternalPolytopeRep ],
  function( poly )

    return Polymake_Inequalities( ExternalPolymakePolytope( poly ) );

end );

InstallMethod( IsBounded,
               " for external polyopes.",
               [ IsPolytope ],
    function( poly )

    return Polymake_IsBounded( ExternalPolymakePolytope( poly ) );

end );

InstallMethod( LatticePoints,
               [ IsPolytope ],
    function( poly )

    return Polymake_LatticePoints( ExternalPolymakePolytope( poly ) );

end );

InstallMethod( Dimension,
               [ IsPolytope ],
    function( poly )

    return Polymake_Dimension( ExternalPolymakePolytope( poly ) );

end );

InstallMethod( IntersectionOfPolytopes,
               [ IsPolytope, IsPolytope ],
  function( poly1, poly2 )
    local poly, ext_polytope;

    if not Rank( ContainingGrid( poly1 ) ) = Rank( ContainingGrid( poly2 ) ) then
        Error( "polytopes are not of the same dimension" );
    fi;

    # compute the intersection
    ext_polytope:= Polymake_Intersection( ExternalPolymakePolytope( poly1 ), ExternalPolymakePolytope( poly2 ) );

    # find the generating vertices of this new polytope and construct the polytope
    # currently, this ignores the lineality space Polymake_Linealities( ext_polytope ) completely!
    poly := Polytope( Polymake_Vertices( ext_polytope ) );

    # set properties
    SetExternalPolymakePolytope( poly, ext_polytope );
    SetContainingGrid( poly, ContainingGrid( poly1 ) );
    SetAmbientSpaceDimension( poly, AmbientSpaceDimension( poly1 ) );

    # return the result
    return poly;

end );