text:
#
# Directory of where my stuff is
#
home_dir:=Directory("~/Workspace/Chevalley.gap/");
lib_dir:=Directory("~/Workspace/Chevalley.gap/lib");
test_dir:=Directory("~/Workspace/Chevalley.gap/test");
data_dir:=Directory("~/Workspace/Chevalley.gap/lib/data");
g2_dir:=Directory("~/Workspace/Chevalley.gap/lib/cases/g2");
g2_dir_char2:=Directory("~/Workspace/Chevalley.gap/lib/cases/g2/char2");
g2_dir_char3:=Directory("~/Workspace/Chevalley.gap/lib/cases/g2/char3");
f4_dir:=Directory("~/Workspace/Chevalley.gap/lib/cases/f4");
f4_dir_char2:=Directory("~/Workspace/Chevalley.gap/lib/cases/f4/char2");
f4_dir_char3:=Directory("~/Workspace/Chevalley.gap/lib/cases/f4/char3");
e6_dir:=Directory("~/Workspace/Chevalley.gap/lib/cases/e6");
e6_dir_char2:=Directory("~/Workspace/Chevalley.gap/lib/cases/e6/char2");
e6_dir_char3:=Directory("~/Workspace/Chevalley.gap/lib/cases/e6/char3");
e7_dir:=Directory("~/Workspace/Chevalley.gap/lib/cases/e7");
e7_dir_char2:=Directory("~/Workspace/Chevalley.gap/lib/cases/e7/char2");
e7_dir_char3:=Directory("~/Workspace/Chevalley.gap/lib/cases/e7/char3");
e8_dir:=Directory("~/Workspace/Chevalley.gap/lib/cases/e8");
e8_dir_char2:=Directory("~/Workspace/Chevalley.gap/lib/cases/e8/char2");
e8_dir_char3:=Directory("~/Workspace/Chevalley.gap/lib/cases/e8/char3");
e8_dir_char5:=Directory("~/Workspace/Chevalley.gap/lib/cases/e8/char5");
components_dir:=Directory("~/Workspace/Chevalley.gap/lib/components");
#
# Read("~/Workspace/Chevalley.gap/init.gap"); Read(Filename(home_dir,"load.gap"));
#

ext: gi


text:
#############################################################################
####
##
#A  anupq.gi                    ANUPQ package                  Eamonn O'Brien
#A                                                             & Frank Celler
##
#Y  Copyright 1992-1994,  Lehrstuhl D fuer Mathematik,  RWTH Aachen,  Germany
#Y  Copyright 1992-1994,  School of Mathematical Sciences, ANU,     Australia
##

#############################################################################
##
#F  ANUPQDirectoryTemporary( <dir> ) . . . . .  redefine ANUPQ temp directory
##
##  calls the UNIX command `mkdir' to create <dir>, which must be  a  string,
##  and if successful a directory  object  for  <dir>  is  both  assigned  to
##  `ANUPQData.tmpdir' and returned. The field  `ANUPQData.outfile'  is  also
##  set to be a file in `ANUPQData.tmpdir', and on exit from {\GAP} <dir>  is
##  removed.
##
InstallGlobalFunction(ANUPQDirectoryTemporary, function(dir)
local created;

  # check arguments
  if not IsString(dir) then
    Error(
      "usage: ANUPQDirectoryTemporary( <dir> ) : <dir> must be a string.\n");
  fi;

  # create temporary directory
  CreateDir(dir);
  if not IsDirectoryPath(dir) then
    return fail;
  fi;

  Add( GAPInfo.DirectoriesTemporary, dir );
  ANUPQData.tmpdir  := Directory(dir);
  ANUPQData.outfile := Filename(ANUPQData.tmpdir, "PQ_OUTPUT");
  return ANUPQData.tmpdir;
end);

#############################################################################
##
#F  ANUPQerrorPq( <param> ) . . . . . . . . . . . . . . . . . report an error
##
InstallGlobalFunction( ANUPQerrorPq, function( param )
    Error(
    "Valid Options:\n",
    "    \"ClassBound\", <bound>\n",
    "    \"Prime\", <prime>\n",
    "    \"Exponent\", <exponent>\n",
    "    \"Metabelian\"\n",
    "    \"OutputLevel\", <level>\n",
    "    \"Verbose\"\n",
    "    \"SetupFile\", <file>\n",
    "    \"PqWorkspace\", <workspace>\n",
    "Illegal Parameter: \"", param, "\"" );
end );

#############################################################################
##
#F  ANUPQextractPqArgs( <args> )  . . . . . . . . . . . . . extract arguments
##
InstallGlobalFunction( ANUPQextractPqArgs, function( args )
    local   CR,  i,  act,  match;

    # allow to give only a prefix
    match := function( g, w )
    	return 1 < Length(g)
               and Length(g) <= Length(w)
               and w{[1..Length(g)]} = g;
    end;

    # extract arguments
    CR := rec();
    i  := 2;
    while i <= Length(args)  do
        act := args[i];
        if not IsString( act ) then ANUPQerrorPq( act ); fi;

    	# "ClassBound", <class>
        if match( act, "ClassBound" ) then
            i := i + 1;
            CR.ClassBound := args[i];

    	# "Prime", <prime>
        elif match( act, "Prime" )  then
            i := i + 1;
            CR.Prime := args[i];

    	# "Exponent", <exp>
        elif match( act, "Exponent" )  then
            i := i + 1;
            CR.Exponent := args[i];

        # "Metabelian"
        elif match( act, "Metabelian" ) then
            CR.Metabelian := true;

    	# "Output", <level>
        elif match( act, "OutputLevel" )  then
            i := i + 1;
            CR.OutputLevel := args[i];
    	    CR.Verbose     := true;

    	# "SetupFile", <file>
        elif match( act, "SetupFile" )  then
    	    i := i + 1;
            CR.SetupFile := args[i];

    	# "PqWorkspace", <workspace>
        elif match( act, "PqWorkspace" )  then
    	    i := i + 1;
            CR.PqWorkspace := args[i];

    	# "Verbose"
        elif match( act, "Verbose" ) then
            CR.Verbose := true;

    	# signal an error
    	else
            ANUPQerrorPq( act );

    	fi;
    	i := i + 1;
    od;
    return CR;

end );

#############################################################################
##
#V  ANUPQGlobalVariables
##
InstallValue( ANUPQGlobalVariables,
              [ "F",          #  a free group
                "MapImages"   #  images of the generators in G
                ] );

#############################################################################
##
#F  ANUPQReadOutput . . . . read pq output without affecting global variables
##
InstallGlobalFunction( ANUPQReadOutput, function( file )
    local globalvars, var, result;

    globalvars := [ "ANUPQmagic", "ANUPQautos", "ANUPQgroups" ];

    for var in globalvars do
        HideGlobalVariables( var );
    od;

    Read( file );

    result := rec();

    for var in globalvars do
        if IsBoundGlobal( var ) then
            result.(var) := ValueGlobal( var );
        else
            result.(var) := fail;
        fi;
    od;

    for var in globalvars do
        UnhideGlobalVariables( var );
    od;

    return result;
end );

#############################################################################
##
#F  PqEpimorphism( <arg> : <options> ) . . . . .  epimorphism onto p-quotient
##
InstallGlobalFunction( PqEpimorphism, function( arg )
    return PQ_EPI_OR_PCOVER(arg : PqEpiOrPCover := "pQepi");
end );

#############################################################################
##
#F  Pq( <arg> : <options> ) . . . . . . . . . . . . . . . . . . .  p-quotient
##
InstallGlobalFunction( Pq, function( arg )
    return PQ_EPI_OR_PCOVER(arg : PqEpiOrPCover := "pQuotient");
end );

#############################################################################
##
#F  PqPCover( <arg> : <options> ) . . . . . .  p-covering group of p-quotient
##
InstallGlobalFunction( PqPCover, function( arg )
    return PQ_EPI_OR_PCOVER(arg : PqEpiOrPCover := "pCover");
end );

#############################################################################
##
#F  PQ_GROUP_FROM_PCP(<datarec>,<out>) . extract gp from pq pcp file into GAP
##
InstallGlobalFunction( PQ_GROUP_FROM_PCP, function( datarec, out )
    local gens;
    HideGlobalVariables( "F", "MapImages" );
    Read( datarec.outfname );
    if out = "pCover" then
      datarec.pCover := ValueGlobal( "F" );
      IsPGroup( datarec.pCover );
    else
      if IsBound(datarec.pcgs) then
        gens := datarec.pcgs;
      else
        gens := GeneratorsOfGroup( datarec.group );
      fi;
      datarec.pQepi := GroupHomomorphismByImagesNC(
                           datarec.group,
                           ValueGlobal( "F" ),
                           gens,
                           ValueGlobal( "MapImages" )
                           );
      SetIsSurjective( datarec.pQepi, true );
      datarec.pQuotient := Image( datarec.pQepi );
      IsPGroup( datarec.pQuotient );
    fi;
    UnhideGlobalVariables( "F", "MapImages" );
end );

#############################################################################
##
#F  TRIVIAL_PQ_GROUP(<datarec>, <out>) . . . extract gp when trivial into GAP
##
InstallGlobalFunction( TRIVIAL_PQ_GROUP, function( datarec, out )
local Q;
    Q := TrivialGroup( IsPcGroup );
    if out = "pCover" then
      datarec.pCover := Q;
      IsPGroup( datarec.pCover );
    else
      datarec.pQepi := GroupHomomorphismByFunction(
                           datarec.group, Q, g -> One(Q) );
      SetIsSurjective( datarec.pQepi, true );
      datarec.pQuotient := Image( datarec.pQepi );
      IsPGroup( datarec.pQuotient );
    fi;
end );

#############################################################################
##
#F  PQ_EPI_OR_PCOVER(<args>:<options>) .  p-quotient, its epi. or its p-cover
##
InstallGlobalFunction( PQ_EPI_OR_PCOVER, function( args )
    local   out, datarec, AtClass, trivial;

    out := ValueOption("PqEpiOrPCover");
    datarec := ANUPQ_ARG_CHK("Pq", args);
    datarec.filter := ["Output file in", "Group presentation"];
    VALUE_PQ_OPTION("Identities", [], datarec);
    if datarec.calltype = "GAP3compatible" then
        # ANUPQ_ARG_CHK calls PQ_EPI_OR_PCOVER itself in this case
        # (so datarec.(out) has already been computed)
        return datarec.(out);
    fi;
    trivial := IsEmpty( datarec.group!.GeneratorsOfMagmaWithInverses );
    if trivial then
        ; #the `pq' binary spits out nonsense if given a trivial pres'n
    elif datarec.calltype = "interactive" and
         ( IsBound(datarec.pQuotient) or IsBound(datarec.pCover) ) then
        AtClass := function()
          return IsBound(datarec.complete) and datarec.complete or
                 IsBound(datarec.class) and datarec.class = datarec.ClassBound;
        end;

        if IsBound(datarec.pcoverclass) and
           datarec.pcoverclass = datarec.class and not AtClass() then
            # ``reduce'' the p-cover to a p-class
            PQ_FINISH_NEXT_CLASS( datarec );
        fi;
        while not AtClass() do
            PQ_NEXT_CLASS( datarec );
        od;
        # the following is not executed if the while-loop is
        # executed at least once
        if IsBound( datarec.(out) ) then
            return datarec.(out); # it had already been computed
        fi;
    else
        PQ_PC_PRESENTATION(datarec, "pQ");
        if datarec.class < Minimum(63, datarec.ClassBound) then
            datarec.complete := true;
        fi;
    fi;

    trivial := trivial or IsEmpty(datarec.ngens) or datarec.ngens[1] = 0;
    if not trivial then
        if out = "pCover" then
          PQ_P_COVER( datarec );
        fi;

        PushOptions( rec(nonuser := true) );
        PQ_WRITE_PC_PRESENTATION(datarec, datarec.outfname);
        PopOptions();
    fi;

    if datarec.calltype = "non-interactive" then
        PQ_COMPLETE_NONINTERACTIVE_FUNC_CALL(datarec);
        if IsBound( datarec.setupfile ) then
          if trivial then
            return fail;
          fi;
          return true;
        fi;
    fi;

    if trivial then
        TRIVIAL_PQ_GROUP( datarec, out );
    else
        # read group and images from file
        PQ_GROUP_FROM_PCP( datarec, out );
    fi;
    return datarec.(out);
end );

#############################################################################
##
#F  PqRecoverDefinitions( <G> ) . . . . . . . . . . . . . . . . . definitions
##
##  This function finds a definition for each generator of the p-group <G>.
##  These definitions need not be the same as the ones used by pq.  But
##  they serve the purpose of defining each generator as a commutator or
##  power of earlier ones.  This is useful for extending an automorphism that
##  is given on a set of minimal generators of <G>.
##
InstallGlobalFunction( PqRecoverDefinitions, function( G )
    local   col,  gens,  definitions,  h,  g,  rhs,  gen;

    col  := ElementsFamily( FamilyObj( G ) )!.rewritingSystem;
    gens := GeneratorsOfRws( col );

    definitions := [];

    for h in [1..NumberGeneratorsOfRws( col )] do
        rhs := GetPowerNC( col, h );
        if Length( rhs ) = 1 then
            gen := Position( gens, rhs );
            if not IsBound( definitions[gen] ) then
                definitions[gen] := h;
            fi;
        fi;

        for g in [1..h-1] do
            rhs := GetConjugateNC( col, h, g );
            if Length( rhs ) = 2 then
                gen := SubSyllables( rhs, 2, 2 );
                gen := Position( gens, gen );
                if not IsBound( definitions[gen] ) then
                    definitions[gen] := [h, g];
                fi;
            fi;
        od;
    od;
    return definitions;
end );

#############################################################################
##
#F  PqAutomorphism( <epi>, <autoimages> ) . . . . . . . . . . . . definitions
##
##  Take an automorphism of the preimage and produce the induced automorphism
##  of the image of the epimorphism.
##
InstallGlobalFunction( PqAutomorphism, function( epi, autoimages )
    local   G,  p,  gens,  definitions,  d,  epimages,  i,  pos,  def,
            phi;

    G      := Image( epi );
    p      := PrimePGroup( G );
    gens   := GeneratorsOfGroup( G );

    autoimages := List( autoimages, im->Image( epi, im ) );

    ##  Get a definition for each generator.
    definitions := PqRecoverDefinitions( G );
    d := Number( [1..Length(definitions)],
                 i->not IsBound( definitions[i] ) );

    ##  Find the images for the defining generators of G under the
    ##  automorphism.  We have to be careful, as some of the generators for
    ##  the source might be redundant as generators of G.
    epimages := List( GeneratorsOfGroup(Source(epi)), g->Image(epi,g) );
    for i in [1..d] do
        ##  Find G.i ...
        pos := Position( epimages, G.(i) );
        if pos = fail then
            Error( "generators ", i, "not image of a generators" );
        fi;
        ##  ... and set its image.
        definitions[i] := autoimages[pos];
    od;

    ##  Replace each definition by its image under the automorphism.
    for i in [d+1..Length(definitions)] do
        def := definitions[i];
        if IsInt( def ) then
            definitions[i] := definitions[ def ]^p;
        else
            definitions[i] := Comm( definitions[ def[1] ],
                                    definitions[ def[2] ] );
        fi;
    od;

    phi := GroupHomomorphismByImages( G, G, gens, definitions );
    SetIsBijective( phi, true );

    return phi;
end );

#############################################################################
##
#F  PqLeftNormComm( <words> ) . . . . . . . . . . . . .  left norm commutator
##
##  returns for a list <words> of words in the generators of a group the left
##  norm commutator of <words>, e.g.~if <w1>, <w2>, <w3>  are  words  in  the
##  generators of some free or fp group then  `PqLeftNormComm(  [<w1>,  <w2>,
##  <w3>] );' is equivalent to `Comm( Comm( <w1>, <w2> ), <w3> );'. Actually,
##  the only restrictions on <words> are that <words> must constitute a  list
##  of group elements of the  same  group  (so  a  list  of  permutations  is
##  allowed, for example) and that <words> must contain at least *two* words.
##
InstallGlobalFunction( PqLeftNormComm, function( words )
local fam, comm, word;
  if not IsList(words) or 2 > Length(words) or
     not ForAll(words, IsMultiplicativeElementWithInverse) then
    Error( "<words> should be a list of at least 2 group elements\n" );
  else
    fam := FamilyObj(words[1]);
    if not ForAll(words, w -> IsIdenticalObj(FamilyObj(w), fam)) then
      Error( "<words> should belong to the same group\n" );
    fi;
  fi;
  comm := words[1];
  for word in words{[2 .. Length(words)]} do
    comm := Comm(comm, word);
  od;
  return comm;
end );

#############################################################################
##
#F  PqGAPRelators( <group>, <rels> ) . . . . . . . . pq relators as GAP words
##
##  returns a list of words that {\GAP} understands, given a list  <rels>  of
##  strings in the string representations of the generators of the  fp  group
##  <group> prepared as a list of relators for the `pq' program.
##
##  *Note:*
##  The `pq' program does not  use  `/'  to  indicate  multiplication  by  an
##  inverse and uses square brackets to represent (left normed)  commutators.
##  Also, even though the `pq' program accepts  relations,  all  elements  of
##  <rels> *must* be in relator form, i.e.~a relation of form `<w1>  =  <w2>'
##  must be written as `<w1>*(<w2>)^-1'.
##
##  Here is an example:
##
##  \beginexample
##  gap> F := FreeGroup("a", "b");
##  gap> PqGAPRelators(F, [ "a*b^2", "[a,b]^2*a", "([a,b,a,b,b]*a*b)^2*a" ]);
##  [ a*b^2, a^-1*b^-1*a*b*a^-1*b^-1*a*b*a, b^-1*a^-1*b^-1*a^-1*b*a*b^-1*a*b*a^
##      -1*b*a^-1*b^-1*a*b*a*b^-1*a^-1*b^-1*a^-1*b*a*b^-1*a*b^-1*a^-1*b*a^-1*b^
##      -1*a*b*a*b*a^-1*b*a*b^-1*a*b*a^-1*b*a^-1*b^-1*a*b*a*b^-1*a^-1*b^-1*a^
##      -1*b*a*b^-1*a*b^-1*a^-1*b*a^-1*b^-1*a*b*a*b^2*a*b*a ]
##  \endexample
##
InstallGlobalFunction( PqGAPRelators, function( group, rels )
local gens, relgens, diff, g;
  if not( IsFpGroup(group) ) then
    Error("<group> must be an fp group\n");
  fi;
  gens := List( FreeGeneratorsOfFpGroup(group), String );
  if not ForAll(rels, rel -> Position(rel, '/') = fail) then
    Error( "pq binary does not understand `/' in relators\n" );
  fi;
  relgens := Set( Concatenation(
                      List( rels, rel -> Filtered(
                                             SplitString(rel, "", "*[]()^, "),
                                             str -> Int(str) = fail) ) ) );
  diff := Difference(relgens, gens);
  if not IsEmpty(diff) then
    Error( "generators: ", diff,
           "\nare not among the generators of the group supplied\n" );
  fi;
  CallFuncList(HideGlobalVariables, gens);
  for g in FreeGeneratorsOfFpGroup(group) do
    ASS_GVAR(String(g), g);
  od;
  rels := List( rels, rel -> EvalString(
                                 ReplacedString(
                                     ReplacedString(rel, "]", "])"),
                                     "[", "PqLeftNormComm(["
                                     ) ) );
  CallFuncList(UnhideGlobalVariables, gens);
  return rels;
end );

#############################################################################
##
#F  PqParseWord( <F>, <word> ) . . . . . . . . . . . . parse word through GAP
#F  PqParseWord( <n>, <word> )
##
##  parse <word> through {\GAP}, where <word> is a string representing a word
##  in the generators of <F> (the first form  of  `PqParseWord')  or  <n>  pc
##  generators `x1,...,x<n>'. `PqParseWord' is provided as a  rough-and-ready
##  check of <word> for syntax errors. A syntax error will cause the entering
##  of a `break'-loop,  in  which  the  error  message  may  or  may  not  be
##  meaningful (depending on whether the syntax  error  gets  caught  at  the
##  {\GAP} or kernel level).
##
##  *Note:*
##  The reason the generators *must* be `x1,...,x<n>' in the second  form  of
##  `PqParseWord' is that these are the pc generator names used by  the  `pq'
##  program (as distinct from the generator names for the group  provided  by
##  the user to a function like `Pq' that invokes the `pq' program).
##
InstallGlobalFunction( PqParseWord, function( n, word )
local ParseOnBreak, ParseOnBreakMessage, NormalOnBreak, NormalOnBreakMessage,
      parts, gens;

  if IsGroup(n) or
     Position(word, '[') <> fail or Position(word, '(') <> fail then
    #pass word through GAP's parser to see if it's ok

    NormalOnBreak := OnBreak;
    ParseOnBreak := function()
      Where(0);
      OnBreak := NormalOnBreak;
    end;
    OnBreak := ParseOnBreak;

    if IsFunction(OnBreakMessage) then
      NormalOnBreakMessage := OnBreakMessage;
      ParseOnBreakMessage := function()
        Print( " syntax error in: ", word, "\n" );
        Print( " you can type: 'quit;' to quit to outer loop.\n" );
        OnBreakMessage := NormalOnBreakMessage;
      end;
      OnBreakMessage := ParseOnBreakMessage;
    fi;

    if IsGroup(n) then
      PqGAPRelators(n, [ word ]);
    else
      PqGAPRelators(FreeGroup(n, "x"), [ word ]);
    fi;

    OnBreak := NormalOnBreak;
    if IsFunction(OnBreakMessage) then
      OnBreakMessage := NormalOnBreakMessage;
    fi;

  else
    parts := List( SplitString(word, "*"), part -> SplitString(part, "^") );
    if ForAny( parts, part -> 2 < Length(part) or
                              2 = Length(part) and not IsInt( Int(part[2]) ) )
       then
      Error( "detected invalid exponent in argument <word>: ", word, "\n");
    fi;
    if ForAny( parts, part -> IsEmpty( part[1] ) or part[1][1] <> 'x' ) then
      Error( "generators in argument <word> must all be of form:\n",
             "`x<i>' for some integer <i>\n" );
    fi;
    gens := List( parts, part -> Int( part[1]{[2 .. Length(part[1])]} ) );
    if not ForAll(gens, gen -> IsPosInt(gen) and gen <= n) then
      Error( "generators in argument <word> must be in the range: ",
             "x1,...,x", n, "\n" );
    fi;
  fi;
  return true;
end );

#############################################################################
##
#F  PQ_EVALUATE( <string> ) . . . . . . . . . evaluate a string emulating GAP
##
##  For each substring of the string <string> that is a statement (i.e.  ends
##  in a `;'), `PQ_EVALUATE( <string> )' evaluates it in the same way  {\GAP}
##  would. If the substring is further followed by  a  `;'  (i.e.  there  was
##  `;;'), this is an indication that the statement would produce no  output;
##  otherwise the output that the user would normally see if  she  typed  the
##  statement interactively is displayed.
##
InstallGlobalFunction(PQ_EVALUATE, function(string)
local from, pos, statement, parts, var;
  from := 0;
  pos := Position(string, ';', from);
  while pos <> fail do
    statement := string{[from + 1..pos]};
    statement := ReplacedString(statement," last "," ANUPQData.example.last ");
    if pos < Length(string) and string[pos + 1] = ';' then
      Read( InputTextString(statement) );
      from := pos + 1;
    else
      parts := SplitString(statement, "", " \n");
      if 1 < Length(parts) and parts[2] = ":=" then
        Read( InputTextString(statement) );
        Read( InputTextString(
                  Concatenation( "View(", parts[1], "); Print(\"\\n\");" ) ) );
        ANUPQData.example.last := parts[1];
      else
        var := EvalString(statement);
        View( var );
        Print( "\n" );
        ANUPQData.example.last := var;
      fi;
      from := pos;
    fi;
    pos := Position(string, ';', from);
  od;
end );

#############################################################################
##
#F  PqExample() . . . . . . . . . . execute a pq example or display the index
#F  PqExample( <example>[, PqStart][, Display] )
#F  PqExample( <example>[, PqStart][, <filename>] )
##
##  With no arguments,  or  with  single  argument  `"index"',  or  a  string
##  <example> that is not the name of a file in the `examples' directory,  an
##  index of available examples is displayed.
##
##  With just the one argument <example> that is the name of a  file  in  the
##  `examples' directory, the example contained in that file is  executed  in
##  its simplest form. Some examples accept options  which  you  may  use  to
##  modify some of the options used in the commands of the example.  To  find
##  out which options an example accepts,  use  one  of  the  mechanisms  for
##  displaying the example described below.
##
##  Some examples have both non-interactive and interactive forms; those that
##  are non-interactive only have a name ending  in  `-ni';  those  that  are
##  interactive only have a name ending in `-i'; examples with  names  ending
##  in  `.g'  also  have  only  one  form;  all  other  examples  have   both
##  non-interactive and interactive forms and for these giving  `PqStart'  as
##  second argument invokes `PqStart' initially  and  makes  the  appropriate
##  adjustments  so  that  the  example  is  executed  or   displayed   using
##  interactive functions.
##
##  If `PqExample' is called with last (second or third)  argument  `Display'
##  then the example  is  displayed  without  being  executed.  If  the  last
##  argument is a non-empty  string  <filename>  then  the  example  is  also
##  displayed without being executed but is also written to a file with  that
##  name. Passing an empty string as last argument has  the  same  effect  as
##  passing `Display'.
##
##  *Note:*
##  The  variables  used  in  `PqExample'  are  local  to  the   running   of
##  `PqExample', so there's no  danger  of  having  some  of  your  variables
##  over-written. However, they are not  completely  lost  either.  They  are
##  saved to a record `ANUPQData.examples.vars', i.e.~if `F'  is  a  variable
##  used in the example then you will be able to access it after  `PqExample'
##  has finished as `ANUPQData.examples.vars.F'.
##
InstallGlobalFunction(PqExample, function(arg)
local name, file, instream, line, input, doPqStart, vars, var, printonly,
      filename, DoAltAction, GetNextLine, PrintLine, action, datarec, optname,
      linewidth, sizescreen, CheckForCompoundKeywords, hasFunctionExpr, parts,
      iscompoundStatement, compoundDepth;

  sizescreen := SizeScreen();
  if sizescreen[1] < 80 then
    SizeScreen([80, sizescreen[2]]);
    linewidth := 80;
  else
    linewidth := sizescreen[1];
  fi;

  if IsEmpty(arg) then
    name := "index";
  else
    name := arg[1];
  fi;

  if name = "README" then
    file := fail;
  else
    file := Filename(DirectoriesPackageLibrary( "anupq", "examples"), name);
  fi;
  if file = fail then
    Info(InfoANUPQ + InfoWarning, 1,
         "Sorry! There is no ANUPQ example with name `", name, "'",
         " ... displaying index.");
    name := "index";
    file := Filename(DirectoriesPackageLibrary( "anupq", "examples"), name);
  fi;

  if name <> "index" then
    doPqStart := false;
    if Length(arg) > 1 then
      # At this point the name of the variable <printonly> doesn't make
      # sense; however, if the value assigned to <printonly> is `Display'
      # or an empty string then we ``print only'' and if it is a non-empty
      # string then it is assumed to be a filename and we `LogTo' that filename.
      printonly := arg[Minimum(3, Length(arg))];
      if arg[2] = PqStart then
        if 2 < Length(name) and
           name{[Length(name) - 1 .. Length(name)]} in ["-i", "ni", ".g"] then
          Error( "example does not have a (different) interactive form\n" );
        fi;
        doPqStart := true;
      fi;
    else
      printonly := false;
    fi;

    DoAltAction := function()
      if doPqStart then
        if action[2] = "do" then
          # uncomment line
          line := line{[2..Length(line)]};
        else
          # replace a variable with a proc id
          line := ReplacedString( line, action[5], action[3] );
        fi;
      fi;
    end;

    if printonly = Display or IsString(printonly) then
      GetNextLine := function()
        local from, to;
        line := ReadLine(instream);
        if line = fail then
          return;
        elif IsBound(action) then
          action := SplitString(action, "", "# <>\n");
          DoAltAction();
          Unbind(action);
        elif 3 < Length(line) and line{[1..4]} = "#alt" then
          # only "#alt" actions recognised
          action := line;
        elif IsMatchingSublist(line, "#comment:") then
          line := ReplacedString(line, " supplying", "");
          from := Position(line, ' ');
          to   := Position(line, '<', from);
          Info(InfoANUPQ, 1,
               "In the next command, you may", line{[from .. to - 1]});
          from := to + 1;
          to   := Position(line, '>') - 1;
          Info(InfoANUPQ, 1, "supplying to `PqExample' the option: `",
                             line{[from .. to]}, "'");
        fi;
      end;

      if IsString(printonly) and printonly <> "" then
        filename := printonly;
        LogTo( filename ); #Make sure it's empty and writable
      fi;
      PrintLine := function()
        if IsMatchingSublist(line, "##") then
          line := line{[2..Length(line)]};
        elif line[1] = '#' then
          return;
        fi;
        Print( ReplacedString(line, ";;", ";") );
      end;
      printonly := true; #now the name of the variable makes sense
    else
      printonly := false;
      ANUPQData.example := rec(options := rec());
      datarec := ANUPQData.example.options;

      CheckForCompoundKeywords := function()
        local compoundkeywords;
        compoundkeywords := Filtered( SplitString(line, "", "( ;\n"),
                                      w -> w in ["do", "od", "if", "fi",
                                                 "repeat", "until",
                                                 "function", "end"] );
        hasFunctionExpr := "function" in compoundkeywords;
        compoundDepth := compoundDepth
                         + Number(compoundkeywords,
                                  w -> w in ["do", "if", "repeat", "function"])
                         - Number(compoundkeywords,
                                  w -> w in ["od", "fi", "until",  "end"]);
        return not IsEmpty(compoundkeywords);
      end;

      GetNextLine := function()
        local from, to, bhsinput;
        repeat
          line := ReadLine(instream);
          if line = fail then return; fi;
        until not IsMatchingSublist(line, "#comment:");
        if IsBound(action) then
          action := SplitString(action, "", "# <>\n");
          if action[1] = "alt:" then
            DoAltAction();
          else
            # action[2] = name of a possible option passed to `PqExample'
            # action[4] = string to be replaced in <line> with the value
            #             of the option if ok and set
            optname := action[2];
            if IsDigitChar(optname[ Length(optname) ]) then
              optname := optname{[1..Length(optname) - 1]};
            fi;
            datarec.(action[2]) := ValueOption(action[2]);
            if datarec.(action[2]) = fail then
              Unbind( datarec.(action[2]) );
            else
              if not ANUPQoptionChecks.(optname)( datarec.(action[2]) ) then
                Info(InfoANUPQ, 1, "\"", action[2], "\" value must be a ",
                                   ANUPQoptionTypes.(optname),
                                   ": option ignored.");
                Unbind( datarec.(action[2]) );
              else
                if action[1] = "add" then
                  line[1] := ' ';
                fi;
                if IsString( datarec.(action[2]) ) then
                  line := ReplacedString( line, action[4],
                                          Flat(['"',datarec.(action[2]),'"']) );
                else
                  line := ReplacedString( line, action[4],
                                          String( datarec.(action[2]) ) );
                fi;
              fi;
            fi;
          fi;
          Unbind(action);
        elif IsMatchingSublist(line, "##") then
          ; # do nothing
        elif 3 < Length(line) and line{[1..4]} in ["#sub", "#add", "#alt"] then
          action := line;
        elif line[1] = '#' then
          # execute instructions behind the scenes
          bhsinput := "";
          repeat
            Append( bhsinput,
                    ReplacedString(line{[2..Length(line)]},
                                   "datarec",
                                   "ANUPQData.example.options") );
            line := ReadLine(instream);
          until line[1] <> '#' or
                (3 < Length(line) and line{[1..4]} in ["#sub", "#add", "#com"]);
          Read( InputTextString(bhsinput) );
        fi;
      end;

      PrintLine := function()
        if IsMatchingSublist(line, "##") then
          line := line{[2..Length(line)]};
        elif line[1] = '#' then
          return;
        fi;
        if input = "" then
          Print("gap> ");
        else
          Print(">    ");
        fi;
        Print( ReplacedString(line, ";;", ";") );
      end;
    fi;
  fi;

  instream := InputTextFile(file);
  if name <> "index" then
    FLUSH_PQ_STREAM_UNTIL( instream, 10, 1, ReadLine,
                           line -> IsMatchingSublist(line, "#Example") );
    line := FLUSH_PQ_STREAM_UNTIL( instream, 1, 10, ReadLine,
                                   line -> IsMatchingSublist(line, "#vars:") );
    if Length(line) + 21 < linewidth then
      Info(InfoANUPQ, 1, line{[Position(line, ' ')+1..Position(line, ';')-1]},
                         " are local to `PqExample'");
    else
      #this assumes one has been careful to ensure the `#vars:' line is not
      #longer than 72 characters.
      Info(InfoANUPQ, 1, line{[Position(line, ' ')+1..Position(line, ';')-1]},
                         " are");
      Info(InfoANUPQ, 1, "local to `PqExample'");
    fi;
    vars := SplitString(line, "", " ,;\n");
    vars := vars{[2 .. Length(vars)]};
    if not printonly then
      CallFuncList(HideGlobalVariables, vars);
    fi;
    line := FLUSH_PQ_STREAM_UNTIL(instream, 1, 10, ReadLine,
                                  line -> IsMatchingSublist(line, "#options:"));
    input := "";
    GetNextLine();
    while line <> fail do
      PrintLine();
      if line[1] <> '#' then
        if not printonly then
          if input = "" then
            compoundDepth := 0;
            iscompoundStatement := CheckForCompoundKeywords();
          elif iscompoundStatement and compoundDepth > 0 then
            CheckForCompoundKeywords();
          fi;
          if line <> "\n" then
            Append(input, line);
            if iscompoundStatement then
              if compoundDepth = 0 and Position(input, ';') <> fail then
                Read( InputTextString(input) );
                if hasFunctionExpr then
                  parts := SplitString(input, "", ":= \n");
                  Read( InputTextString(
                            Concatenation(
                                "View(", parts[1], "); Print(\"\\n\");" ) ) );
                  ANUPQData.example.last := parts[1];
                fi;
                iscompoundStatement := false;
                input := "";
              fi;
            elif Position(input, ';') <> fail then
              PQ_EVALUATE(input);
              input := "";
            fi;
          fi;
        fi;
      fi;
      GetNextLine();
    od;
    if printonly then
      if IsBound(filename) then
        LogTo();
      fi;
    else
      ANUPQData.example.vars := rec();
      for var in Filtered(vars, ISBOUND_GLOBAL) do
        ANUPQData.example.vars.(var) := VALUE_GLOBAL(var);
      od;
      Info(InfoANUPQ, 1, "Variables used in `PqExample' are saved ",
                         "in `ANUPQData.example.vars'.");
      CallFuncList(UnhideGlobalVariables, vars);
    fi;
  else
    FLUSH_PQ_STREAM_UNTIL(instream, 1, 10, ReadLine, line -> line = fail);
  fi;
  CloseStream(instream);
  if linewidth <> sizescreen[1] then
    SizeScreen( sizescreen ); # restore what was there before
  fi;
end);

#############################################################################
##
#F  AllPqExamples() . . . . . . . . . .  list the names of all ANUPQ examples
##
InstallGlobalFunction( AllPqExamples, function()
  local dir,  files;

  dir   := DirectoriesPackageLibrary( "anupq", "examples" )[1];
  files := DirectoryContents( Filename( dir, "" ));
  # Remove certain files
  files := Difference( files, [".", "..", "index", "README", "CVS",
                                         "5gp-PG-e5-i", "7gp-a-x-Rel-i"] );
  # Remove files ending with a tilde
  files := Filtered( files, file -> file[ Length(file) ] <> '~' );
  return files;
end );

#############################################################################
##
#F  GrepPqExamples( <string> ) . . . . . . . grep ANUPQ examples for a string
##
##  runs the UNIX command `grep <string>'  over  the  {\ANUPQ}  examples  and
##  returns the list of examples for which  there  is  a  match.  The  actual
##  matches are `Info'-ed at `InfoANUPQ' level 2.
##
InstallGlobalFunction( GrepPqExamples, function( string )
  local dir,  str,  grep,  out,  opts,  lines,  matches,  line;

  dir := DirectoriesPackageLibrary( "anupq", "examples" )[1];
  grep := Filename( DirectoriesSystemPrograms(), "grep" );
  str := "";
  out := OutputTextString( str, true );
  opts := Concatenation( [ string ], AllPqExamples() );
  Process( dir, grep, InputTextNone(), out, opts );
  CloseStream( out );
  lines := SplitString( str, "",  "\n" );
  matches := [];
  for line in lines do
    Info(InfoANUPQ, 2, line);
    Add( matches, SplitString(line, "", ":")[1] );
  od;
  return Set(matches);
end );

#E  anupq.gi  . . . . . . . . . . . . . . . . . . . . . . . . . . . ends here

max_stars_repo_name: gap-system/anupq

max_issues_repo_name: gap-system/anupq

max_forks_repo_name: gap-system/anupq


text:
Pascal := function(n)
	local i, v;
	v := [1];
	for i in [1 .. n] do
		Display(v);
		v := Concatenation([0], v) + Concatenation(v, [0]);
	od;
end;

Pascal(9);
# [ 1 ]
# [ 1, 1 ]
# [ 1, 2, 1 ]
# [ 1, 3, 3, 1 ]
# [ 1, 4, 6, 4, 1 ]
# [ 1, 5, 10, 10, 5, 1 ]
# [ 1, 6, 15, 20, 15, 6, 1 ]
# [ 1, 7, 21, 35, 35, 21, 7, 1 ]
# [ 1, 8, 28, 56, 70, 56, 28, 8, 1 ]

ext: gap

max_stars_repo_name: mullikine/RosettaCodeData

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

#F Replaces RemoteData inside a BB with RemoteDataUnrolled
UnrollRemoteData := function(e, cx)
   local extvars, loopvar, i, unrolledLoopVars, BBLoopVars, unrolledRange, newobj, newvar, newFofs, newofs, newfdataofs, f;

   # Compute list of loopvars that will be unrolled (ones inside BB)
   # These are all loops under the BB
   i := Length(cx.parents);
   unrolledLoopVars := [];
   unrolledRange := 1;
   while (ObjId(cx.parents[i])<>MultiBufISum and ObjId(cx.parents[i]) <> MultiBufISumFinal) do
      i := i-1;
      if ObjId(cx.parents[i])=ISum then
         unrolledLoopVars := unrolledLoopVars :: [cx.parents[i].var];
      fi;
   od;

   #Error("BP");

  # Take intersection because twiddle might not involve all loops under BB.
  BBLoopVars := Intersection(e.ofs.free(), unrolledLoopVars);
  unrolledRange := ListProduct(List(unrolledLoopVars, i->i.range));
  unrolledRange := When(unrolledRange=0, 1, unrolledRange);

  #NOTE: why was this an issue?
  #if BBLoopVars = [] then
  #   PrintLine("------------------> Going ahead with no Isum inside BB");
  #fi;

  # a) return RemoteDataUnrolled
  # b) with FDataOfs_mbuf that has a new range

  #Error("BP");

  newobj := e.child(1);    # This is the VRCDiag or similar obj

  # Change the FDataOfs_mbuf to include a different length and offset

  # Offset for FDataOfs_mbuf is RemoteData.ofs, with freevars that are not inside the BB set to zero.
  newFofs := Copy(e.ofs);

  extvars := newFofs.free();
  SubtractSet(extvars, BBLoopVars);

  for loopvar in extvars do
     SubstVars(newFofs, rec( (loopvar.id) := 0));
  od;

  newobj := SubstTopDown(newobj, FDataOfs_mbuf, f->FDataOfs_mbuf(f.var, f.len*unrolledRange, newFofs));

  e.var.t.size    := e.var.t.size    * unrolledRange;
  e.altbuf.t.size := e.altbuf.t.size * unrolledRange;

  # Must manually set unrolled loopvars to zero here because this will go to
  # init, and won't be taken care of by unrolling
  newofs := e.ofs;
  for loopvar in BBLoopVars do
     SubstVars(newofs, rec( (loopvar.id) := V(0) ));
  od;

  return(RemoteDataUnrolled(e.var, e.altbuf, e.value, e.ofs, newobj));
end;

Class(MBufCodegen, DistCodegen, rec(

    RemoteDataInit := meth(self, o, y, x, opts)
    #-----------------------------------------------------------------------------
        local r;
        r := Collect(o, RemoteDataUnrolled);
        # NOTE: handle multiple RemoteDataUnrolleds here.

        if Length(r) > 1 then
           Error("ERROR: Multi RemoteDataInit case NOT implemented!");
        fi;
        r := r[1];
        return(chain(
          self( RemoteDataNoBody(r.var, r.value, r.ofs, r.child(1)), y, x, opts ),
          self(o.child(1), y, x, opts)
        ));
    end,

    # This is called only from the BB body. So codegen only o.child(1)
    RemoteDataUnrolled := (self, o, y, x, opts) >> self( o.child(1), y, x, opts ),
    #-----------------------------------------------------------------------------

    # This is called in situations where there is no unrolled code.
    RemoteDataNoBody := meth(self, o, y, x, opts)
    #-----------------------------------------------------------------------------

      local pkSizeInBytes, dmaCommand;

      dmaCommand := var("GATHMEM_DIAG");
      pkSizeInBytes := o.child(1).element.len * opts.vector.vlen * (opts.vector.isa.bits/8);

      if pkSizeInBytes > 16384 then
         if pkSizeInBytes = 32768 then
            dmaCommand := var("GATHMEM_DIAG_32K");
         elif pkSizeInBytes = 65536 then
            dmaCommand := var("GATHMEM_DIAG_64K");
         else
           Error("Too large a packet size for getting twiddles");
         fi;
      fi;
     return(chain(
       call(dmaCommand,
       o.altbuf.id,    # Destination base address
       0,              # Destination offset
       "spe_info."::o.value.id,     # Source base address
       o.ofs*opts.vector.vlen,       # Source offset
       o.child(1).element.len * opts.vector.vlen # Packet size in elements
       )

       # Only if

       #call(var("spu_writech"),
       #"MFC_WrTagMask",
       #"1 << 3"),
       #
       #call(var("spu_mfcstat"),
       #"MFC_TAG_UPDATE_ALL"
       #)
     ));
    end,

    # # These should've gotten converted to RemoteDataUnrolled and RemoteDataInit
    # RemoteData := meth(self, o, y, x, opts)
    #     Error("# These should've gotten converted to RemoteDataUnrolled and RemoteDataInit");
    # end,

    # These should've gotten converted to RemoteDataUnrolled and RemoteDataInit
    # But during DP, these might not always be contained in a MultiBufISum. In
    # this case, we should ideally add a new MultiBufISum. Here's a hack instead:
    #RemoteData := (self, o, y, x, opts) >> self(o.child(1), y, x, opts),
    RemoteData := (self, o, y, x, opts) >> skip(),

#-----------------------------------------------------------------------------
    Formula := meth(self, o, y, x, opts)
#-----------------------------------------------------------------------------
#NOTE: Formula needs major cleanup Tue 31 Mar 2009 12:58:40 AM EDT
        local icode, datas, datas_ppe, prog, params, sub, initsub, initsub_ppe, io, e, d;

        [x, y] := self.initXY(x, y, opts);

        o :=  Process_fPrecompute(o, opts);

        o := o.child(1);
        params := Set(Concatenation(Collect(o, param), Filtered(Collect(o, var), IsParallelLoopIndex)));

        SubstTopDown(o, @@(1, VRCDiag,
            (e,cx)->(ObjId(e.element)=FDataOfs
                     and e.element.var.t.size * opts.vector.vlen >= opts.maxSPETwiddles
                     #and IsBound(cx.MultiBufISum)
                     #and IsBound(cx.BB)
                     )
            ),
          e->let(v := @@(1).val,
                 f     := v.element,
                 dmbuf     := var.fresh_t("Tmbuf",  TArray(f.var.t.t, f.len)),
                 dmbuf_alt := var.fresh_t("TmbufAlt", TArray(f.var.t.t, f.len)),
                 RemoteData(dmbuf, dmbuf_alt, f.var, f.ofs, VRCDiag(FDataOfs_mbuf(dmbuf, f.len, 0), v.v))
              )
        );

        # Replace FDataOfs that won't fit inside an SPE with FDataOfs_mbuf
        #SubstTopDown(o, @@(1, FDataOfs,
        #    (e,cx)->(e.var.t.size*opts.vector.vlen >= opts.maxSPETwiddles
        #             #and IsBound(cx.MultiBufISum)
        #             #and IsBound(cx.BB)
        #             )
        #    ),
        #e->let(f:=@@(1).val, FDataOfs_mbuf(f.rChildren()[1], f.rChildren()[2], f.rChildren()[3]))
        #);

        datas := Collect(o, FDataOfs);

        datas_ppe := Collect(o, RemoteData);

        #datas_ppe := [];
        #for d in datas do
        #   if d.var.t.size * opts.vector.vlen >= opts.maxSPETwiddles then
        #      datas_ppe := datas_ppe :: [d];
        #      datas := RemoveList(datas, d);
        #   fi;
        #od;

        # Replace RemoteData with RemoteDataUnroll where necessary
        #NOTE: Checking for context is not enough: must check for length
        #SubstTopDown(o, @@(1, RemoteData, (e,cx)->IsBound(cx.BB)), UnrollRemoteData);

        #SubstTopDown(o, @@(1, RemoteData, (e,cx)->(IsBound(cx.MultiBufISum) or IsBound(cx.MultiBufISumFinal))),
        #    UnrollRemoteData);

        SubstTopDown(o, @@(1, RemoteData, (e,cx)->IsBound(cx.MultiBufISum) or IsBound(cx.MultiBufISumFinal)),
            UnrollRemoteData);

        # Enclose BBs within a RemoteDataInit if needed (if the BB contains one or more
        # RemoteDataUnrolled's)
        #SubstTopDown(o, @@(1, MultiBufISum, (e,cx)->(
        #    Length(Collect(e, RemoteDataUnrolled))>=1
        #    and (not IsBound(cx.RemoteDataInit) or Length(cx.RemoteDataInit) = 0)
        #  )),
        #  f->RemoteDataInit(f)
        #);

        #return(o);

        o := BlockSumsOpts(o, opts);
        icode := self(o, y, x, opts);
        icode := ESReduce(icode, opts);
        icode := RemoveAssignAcc(icode);
        icode := BlockUnroll(icode, opts);
        # icode := PowerOpt(icode);
        icode := DeclareHidden(icode);
        if IsBound(opts.isFixedPoint) and opts.isFixedPoint then
            icode := FixedPointCode(icode, opts.bits, opts.fracbits);
        fi;

        io := When(x=y, [x], [y, x]);
        sub := Cond(IsBound(opts.subName), opts.subName, "transform");
        initsub := Cond(IsBound(opts.subName), Concat("init_", opts.subName), "init");
        initsub_ppe := Cond(IsBound(opts.subName), Concat("init_ppe_", opts.subName), "init_ppe");
        icode := func(TVoid, sub, Concatenation(io, params), icode);

        #Error("BP");

        if IsBound(opts.generateInitFunc) and opts.generateInitFunc then
            prog := program(
                #decl(List(datas, x->x.var)::List(datas_ppe, x->x.var),
                decl(List(datas, x->x.var),
                    chain(
                        func(    TVoid, initsub,     params, chain(List(datas,     x -> SReduce(x.var.init, opts)))),
                        func_ppe(TVoid, initsub_ppe, params, decl(List(datas_ppe, x->x.value), chain(List(datas_ppe, x -> SReduce(x.value.init, opts)))) ),
                        icode
                    )));
        else
            prog := program( func(TVoid, initsub, params, chain()), icode);
        fi;

        # FF: I really don't know why suddenly AVX_8x32f requires me to do that !!
        if IsBound(opts.vector.isa.fixProblems) then prog := opts.vector.isa.fixProblems(prog, opts); fi;

        prog.dimensions := o.dims();

        # Perform software pipelining
        # prog := MarkPreds(prog);
        # prog := MarkDefUse(prog);
        # SubstTopDown(prog, @(1, loop, e->(Length(Collect(e, loop))=1 and Length(e.range) >= 4)),
        #     e->let(l := @(1).val, loop_sw(l.rChildren()[1], l.rChildren()[2], l.rChildren()[3]) )
        # );
        # SubstTopDown(prog, loop_sw, e->SoftwarePipeline(e));

        return prog;
    end,

    GathMem := meth(self, o, y, x, opts)
        #Gath: loop(i, o.func.domain(), assign(nth(y,i), nth(x, func.at(i))));
        #DO_DMA(Xalt, spe_info.X, IDEAL_DMA_SIZE_BYTES, MFC_GET_CMD); // Xalt = spe_info.X
        local func, pkSize, i, numPkts, pkSizeInBytes, dmaCommand;
        func        := o.func.lambda();
        pkSize      := o.pkSize;
        numPkts     := o.func.domain();
        #i          := 0;
        i           := Ind(numPkts);

        dmaCommand := var("GATHMEM_GET");
        pkSizeInBytes := o.pkSize * (opts.vector.isa.bits/8);

        #NOTE: 16384 is the largest packet we can send on the Cell
        # If we're not using tables, then we can go higher using backend macros
        if (pkSizeInBytes > 16384) then
          # NOTE: How to handle cases where pkSize is greater than 16k?
          if (pkSizeInBytes = 32768) then
             dmaCommand := var("GATHMEM_GET_32K");
          elif (pkSizeInBytes = 65536) then
             dmaCommand := var("GATHMEM_GET_64K");
          elif (pkSizeInBytes = 131072) then
             dmaCommand := var("GATHMEM_GET_128K");
          else
          Error("Selected packet size - ", o.pkSize, " - is higher than the Cell's max allowed (using macros).");
          fi;
       fi;

        # NOTE: remove this after testing
        if (numPkts <= 16 ) then
            #16 is the size of the SPU's DMA Command queue. Exceeding this
            #means taking a huge performance loss because in effect, DMAs won't
            #be done in the background.
            # NOTE: this (16) is hardcoded for the implementation than for the ISA
            return(loop(i, numPkts, call(
                   dmaCommand,
                   #y,
                   var(Concatenation(y.id, "alt")),
                   i*pkSize,
                   var(Concatenation("spe_info.", y.id)),
                   func.at(i)*pkSize,
                   pkSize
                )));
        fi;

        if (numPkts > 2048) then
          Error("Selected attemping to DMA more than 2048 packets at one time.  The Cell architecture has no facility for this. (Doing multiple sets is possible, but in most cases, a better algorithm is needed).");
        fi;

        #NOTE: 16384 is the largest packet we can send on the Cell using DMA tables
        if (pkSizeInBytes > 16384) then
          Error("Selected packet size - ", o.pkSize, " - is higher than the Cell's max allowed.");
        fi;

        return(chain(
            # Build DMA list
            loop(i, numPkts, chain(
              call(var("GATH_LIST_SETSIZE"), pkSize, i.id),
              call(var("GATH_LIST_SETADDR"), func.at(i)*pkSize, i.id, var(Concatenation("spe_info.", y.id)))
            )),

            # Execute DMA list
            call(var("GATHMEM_LIST"),
                var(Concatenation(y.id, "alt")),
                var(Concatenation("spe_info.", y.id)),
                numPkts
            )
        ));

    end,

    ScatMem := meth(self, o, y, x, opts)
        #Scat: return loop(i, o.func.domain(), assign(nth(y,func.at(i)), nth(x, i)));
        local func, pkSize, i, numPkts, pkSizeInBytes, dmaCommand;
        func          := o.func.lambda();
        pkSize        := o.pkSize;
        numPkts     := o.func.domain();
        #i          := 0;
        i           := Ind(numPkts);

        pkSizeInBytes := o.pkSize * (opts.vector.isa.bits/8);

        #NOTE: 16384 is the largest packet we can send on the Cell
        # If we're not using tables, then we can go higher using backend macros
        dmaCommand := var("SCATMEM_PUT");
        pkSizeInBytes := o.pkSize * (opts.vector.isa.bits/8);

        #NOTE: 16384 is the largest packet we can send on the Cell
        # If we're not using tables, then we can go higher using backend macros
        if (pkSizeInBytes > 16384) then
          # NOTE: How to handle cases where pkSize is greater than 16k?

          if (pkSizeInBytes = 32768) then
             dmaCommand := var("SCATMEM_PUT_32K");
          elif (pkSizeInBytes = 65536) then
             dmaCommand := var("SCATMEM_PUT_64K");
          elif (pkSizeInBytes = 131072) then
             dmaCommand := var("SCATMEM_PUT_128K");
          else
          Error("Selected packet size - ", o.pkSize, " - is higher than the Cell's max allowed (using macros).");
          fi;
        fi;

        #DO_DMA(Yalt, spe_info.Y, IDEAL_DMA_SIZE_BYTES, MFC_PUT_CMD); // spe_info.Y = Yalt

        if (numPkts <= 16 ) then
            return(loop(i, numPkts, call(
                dmaCommand,
                #x,
                var(Concatenation(x.id, "alt")),
                i*pkSize,
                var(Concatenation("spe_info.", x.id)),
                func.at(i)*pkSize,
                pkSize
            )));
        fi;

        if (numPkts > 2048) then
          Error("Selected attemping to DMA more than 2048 packets at one time. The Cell architecture has no facility for this");
        fi;

        #NOTE: 16384 is the largest packet we can send on the Cell using DMA tables
        if (pkSizeInBytes > 16384) then
          Error("Selected packet size - ", o.pkSize, " - is higher than the Cell's max allowed for DMA lists.");
        fi;

        return(chain(
            # Build DMA list
            loop(i, numPkts, chain(
              call(var("SCAT_LIST_SETSIZE"), pkSize, i.id),
              call(var("SCAT_LIST_SETADDR"), func.at(i)*pkSize, i.id, var(Concatenation("spe_info.", x.id)))
            )),

            # Execute DMA list
            call(var("SCATMEM_LIST"),
                var(Concatenation(x.id, "alt")),
                var(Concatenation("spe_info.", x.id)),
                numPkts
            )
        ));

    end,

     R2Sum := (self, o, y, x, opts) >>
            loop(o.var1, o.domain1, loop(o.var2, o.domain2, self(o.child(1), y, x, opts))),

    # Note: to perform multibuffering, x,y are reversed in some of the following calls.
    MultiBufISum := meth(self, o, y, x, opts)
       local r, mbufvars, mloop, loopbody;

        r := Collect(o, RemoteDataUnrolled);

        # NOTE: handle multiple RemoteDataUnrolleds here.
        if Length(r) > 1 then Error("ERROR: Multi RemoteDataInit case NOT implemented!"); fi;

        mbufvars := [];
        if Length(r) > 0 then
           r := r[1];
           mbufvars := [r.var]::[r.altbuf];
        fi;

        mloop := multibuffer_loop;
        if IsBound(opts.nombuf) and opts.nombuf=true then mloop := mem_loop; fi;

#loopvar, range, y, x, gathmem, twiddles, bufs, cmd, scatmem)
        return(decl(mbufvars,
          mloop(o.var, o.domain, y, x,
               self(o.gathmem, x, y, opts),
               When(r=[], [], self( RemoteDataNoBody(r.var, r.altbuf, r.value, r.ofs, r.child(1)), y, x, opts )),
               mbufvars,
               self(o.child(1), y, x, opts),
               self(o.scatmem, x, y, opts)
          )
        ));
    end,

    MultiBufISumFinal := (self, o, y, x, opts) >>
        self(MultiBufISum(o.var, o.domain, o.scatmem, o.child(1), o.gathmem), y, x, opts),

    #HACK: this won't work if twiddles need to be streamed.
    MemISum := (self, o, y, x, opts) >>
       mem_loop(o.var, o.domain, y, x,
            self(o.gathmem, x, y, opts),
            [],
            [],
            self(o.child(1), y, x, opts),
            self(o.scatmem, x, y, opts)
       ),

    MemISumFinal := (self, o, y, x, opts) >>
       mem_loop(o.var, o.domain, y, x,
            self(o.gathmem, x, y, opts),
            self(o.child(1), y, x, opts),
            self(o.scatmem, x, y, opts)
       ),

    MultiBufDistSum := (self, o, y, x, opts) >>
        self(MultiBufISum(o.var, o.domain, o.scatmem, o.child(1), o.gathmem), y, x, opts)
));


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

InstallMethod( ExternalPolymakeCone,
               [ IsCone ],
   function( cone )

    if IsBound( cone!.input_rays ) and Length( cone!.input_rays ) = 1 and IsZero( cone!.input_rays ) then
        return Polymake_ConeByGenerators( cone!.input_rays[ 1 ] );
    fi;

    if IsBound( cone!.input_rays ) then
        return Polymake_ConeByGenerators( cone!.input_rays );
    fi;

    if IsBound( cone!.input_equalities ) then
        return Polymake_ConeFromInequalities( cone!.input_inequalities, cone!.input_equalities );
    fi;

    # otherwise our fallback is cone from inequalities
    return Polymake_ConeFromInequalities( cone!.input_inequalities );

end );

InstallMethod( RayGenerators,
               [ IsCone ],
    function( cone )
    local rays, lin;

    rays := Polymake_Rays( ExternalPolymakeCone( cone ) );
    lin := Polymake_Lineality( ExternalPolymakeCone( cone ) );
    return Set( Concatenation( rays, lin, (-1) * lin ) );

end );

InstallMethod( IsPointed,
               [ IsCone ],
    function( cone )

    return Polymake_IsPointed( ExternalPolymakeCone( cone ) );

end );

InstallMethod( Dimension,
               [ IsCone ],
    function( cone )

    return Polymake_Dimension( ExternalPolymakeCone( cone ) );

end );

InstallMethod( DefiningInequalities,
               [ IsCone ],
  function( cone )
    local ineqs, eqs;

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

end );

InstallMethod( LinealitySpaceGenerators,
               [ IsCone ],
  function( cone )

    return Set( Polymake_Lineality( ExternalPolymakeCone( cone ) ) );

end );

InstallMethod( RaysInFacets,
               [ IsCone ],
    function( cone )

    return Polymake_RaysInFacets( ExternalPolymakeCone( cone ) );

end );

InstallMethod( RaysInFaces,
               " for cones",
               [ IsCone ],

  function( cone )

    return Set( Polymake_RaysInFaces( ExternalPolymakeCone( cone ) ) );

end );

InstallMethod( IntersectionOfCones,
               "for homalg cones",
               [ IsCone, IsCone ],
  function( cone1, cone2 )
    local ext_cone, ineqs, eqs, cone;

    if not Rank( ContainingGrid( cone1 ) ) = Rank( ContainingGrid( cone2 ) ) then
        Error( "cones are not from the same grid" );
    fi;

    # compute the intersection cone in Polymake
    ext_cone := Polymake_Intersection( ExternalPolymakeCone( cone1 ), ExternalPolymakeCone( cone2 ) );
    ineqs := ext_cone!.inequalities;
    eqs := ext_cone!.equalities;
    cone := ConeByEqualitiesAndInequalities( eqs, ineqs );
    SetContainingGrid( cone, ContainingGrid( cone1 ) );

    # and return it
    return cone;

end );


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

Declare(DST3, DST7);

#F DST1(<n>) - Discrete Sine Transform, Type I, non-terminal
#F Definition: (n x n)-matrix [ sin((k+1)*(l+1)*pi/(n+1)) | k,l = 0...n-1 ]
#F Note:       The natural size for a DST1 is 2^k - 1
#F Note:       DST4 is symmetric
#F Example:    DST1(7)
#F
Class(DST1, DTTBase, rec(
    terminate := self >> When(self.params[1] > 1, Mat(DST_Iunscaled(self.params[1])), I(1)),
    transpose := self >> Copy(self),
    SmallRandom := () -> Random([2,3,5,7,11,15,23,31])
));

#F DST2(<n>) - Discrete Sine Transform, Type II, non-terminal
#F Definition: (n x n)-matrix [ sin(k*(l+1/2)*pi/n) | k,l = 0...n-1 ]
#F Note:       DST2 is the transpose of DST3
#F Example:    DST2(8)
Class(DST2, DTTBase, rec(
    terminate := self >> Mat(DST_IIunscaled(self.params[1])),
    transpose := self >> DST3(self.params[1]),
    SmallRandom := () ->  Random([2,3,4,6,8,9,12,16,18,24,27,32])
));

#F DST3(<n>) - Discrete Sine Transform, Type III, non-terminal
#F Definition: (n x n)-matrix [ sin((k-1/2)*l*pi/n) | k,l = 1...n ]
#F Note:       DST3 is the transpose of DST2
#F Example:    DST3(8)
#F Scaled variant (not supported) is:
#F                [ a_l*sin(k*(l-1/2)*pi/n) | k,l = 1...n ]
#F            with  a_j = 1/sqrt(2) for j = n and = 1 else.
Class(DST3, DTTBase, rec(
    terminate := self >> Mat(DST_IIIunscaled(self.params[1])),
    transpose := self >> DST2(self.params[1]),
    SmallRandom := () ->  Random([2,3,4,6,8,9,12,16,18,24,27,32])
));

#F DST4(<n>) - Discrete Sine Transform, Type IV, non-terminal
#F Definition: (n x n)-matrix [ sin((k-1/2)*(l-1/2)*pi/n) | k,l = 1...n ]
#F Note:       DST4 is symmetric
#F Example:    DST4(8)
Class(DST4, DTTBase, rec(
    terminate := self >> Mat(DST_IVunscaled(self.params[1])),
    transpose := self >> Copy(self),
    SmallRandom := () ->  Random([2,3,4,6,8,9,12,16,18,24,27,32])
));

#F DST5(<n>) - Discrete Sine Transform, Type V, non-terminal
#F Definition: (n x n)-matrix [ sin((k+1)*(l+1)*pi/(n+1/2)) | k,l = 0...n-1 ]
#F Note:       DST5 is symmetric
#F Example:    DST5(8)
Class(DST5, DTTBase, rec(
    terminate := self >> Mat(DST_Vunscaled(self.params[1])),
    transpose := self >> Copy(self),
));

#F DST6(<n>) - Discrete Sine Transform, Type VI, non-terminal
#F Definition: (n x n)-matrix [ sin((k+1)*(l+1/2)*pi/(n+1/2)) | k,l = 0...n-1 ]
#F Note:       The transpose of DST6 is DST7
#F Example:    DST6(8)
Class(DST6, DTTBase, rec(
    terminate := self >> Mat(DST_VIunscaled(self.params[1])),
    transpose := self >> DST7(self.params[1])
));

#F DST7(<n>) - Discrete Sine Transform, Type VII, non-terminal
#F Definition: (n x n)-matrix [ sin((k+1/2)*(l+1)*pi/(n+1/2)) | k,l = 0...n-1 ]
#F Note:       The transpose of DST7 is DST6
#F Example:    DST7(8)
Class(DST7, DTTBase, rec(
    terminate := self >> Mat(DST_VIIunscaled(self.params[1])),
    transpose := self >> DST6(self.params[1])
));

#F DST8(<n>) - Discrete Sine Transform, Type VIII, non-terminal
#F Definition: (n x n)-matrix [ sin((k+1/2)*(l+1/2)*pi/(n-1/2)) | k,l = 0...n-1 ]
#F Note:       DST8 is symmetric
#F Example:    DST8(8)
Class(DST8, DTTBase, rec(
    terminate := self >> Mat(DST_VIIIunscaled(self.params[1])),
    transpose := self >> Copy(self),
));


text:
# This file is a script for GAP and tests a list of polynomials in hexadecimal
# for irreducibility over F_2

# create x over F_2 = GF(2)
x := Indeterminate(GF(2), "x");

# test if polynomial is irreducible, i.e. the number of factors is one
IrredPoly := function (poly)
    return (Length(Factors(poly)) = 1);
end;;

# create a polynomial in x from the hexadecimal representation of the
# coefficients
Hex2Poly := function (s)
    return ValuePol(CoefficientsQadic(IntHexString(s), 2), x);
end;;

# list of candidates, in hex
candidates := [ "3DA3358B4DC173" ];

# create real polynomials
L := List(candidates, Hex2Poly);

# filter and display the list of irreducible polynomials contained in L
Display(Filtered(L, x -> (IrredPoly(x))));

max_stars_repo_name: swaitw/restic

max_issues_repo_name: LCYLYM/restic

max_forks_repo_name: LCYLYM/restic


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

TempArrayOL:=function(y,x,child)
  local y1,x1,mytype,accu;
  return StripList(List(child.dmn(),l->TempVec(l)));
end;

#Class(OLCodegen,VectorCodegen, rec(
Class(OLCodegen, RecCodegenMixin, SMPCodegenMixin, VectorCodegen, rec(

    LinkIO := (self, o, y, x, opts) >> self(o.params[1], y, x, opts),

    Compose := meth(self, o,y,x,opts)
        local ch, numch, vecs, allow, i, j, indices, output, chcomp, veclen, rveclen, dmnoffset, rngoffset, outvar, cc;

    ch := Filtered(o.children(), i-> not i.name in ["DMPGath","DMPScat"]);
    numch := Length(ch);
    vecs := [y];
    allow := (x<>y);

    for i in [1..numch-1] do
        if allow and ObjId(ch[i])=Inplace then
            vecs[i+1] := vecs[i];
        else
            vecs[i+1] := TempArrayOL(y,x,ch[i]);
        fi;
    od;
    vecs[numch+1] := x;
    for i in Reversed([1..numch]) do
            if allow and ObjId(ch[i])=Inplace
        then vecs[i] := vecs[i+1]; fi;
    od;

        # everything was inplace, make it go from x -> y as expected
    if vecs[1] = vecs[numch+1] then vecs[1] := y; fi;
    if vecs[1] = vecs[numch+1] then vecs[numch+1] := x; fi;

    #Error("Codegn-Compose -> Wire");
    #Start from the most right and handle Wire

    i := numch;
    cc:=ch[i];
    if ObjId(cc)=BB then
	cc:=cc._children[1];
    fi;

    while (i>1 and (ObjId(cc)=Wire or ObjId(cc) = Cross)) do

        if ObjId(cc) = Wire then
        #Variable length
            if not IsList(vecs[i+1]) then
                veclen := 1;
                indices := Compacted([1..veclen]);
                output := List(cc.params[2]*indices,x->vecs[i+1]);
            else
                veclen := Length(vecs[i+1]);
                indices := Compacted([1..veclen]);
                output := List(cc.params[2]*indices,x->vecs[i+1][x]);
	fi;
	vecs[i] := output;
    else
    # Wires in Cross
        dmnoffset := 0;
        rngoffset := 0;
        for chcomp in cc._children do
	    if ObjId(cc)=BB then
	       cc:=cc._children[1];
	    fi;

        if (ObjId(chcomp) = Wire) then
            if not IsList(chcomp.dims()[2]) then
                veclen := 1;
            else
                veclen := Length(chcomp.dims()[2]);
            fi;
            if not IsList(chcomp.dims()[1]) then
                rveclen := 1;
            else
                rveclen := Length(chcomp.dims()[1]);
            fi;
            indices := Compacted([1..veclen]);
            output := List(chcomp.params[2]*indices,x->vecs[i+1][x+dmnoffset]);
            j:=1;
            for outvar in output do
                vecs[i][rngoffset + j] := outvar;
                j:=j+1;
            od;
            dmnoffset:=dmnoffset+veclen;
            rngoffset:=rngoffset+rveclen;
        fi;
    od;

    fi;
        i:=i-1;
        cc:=ch[i];
        if ObjId(cc)=BB then
	  cc:=cc._children[1];
	fi;
    od;

#    for i in [1..numch-1] do
#      if ObjId(ch[i])=Wire then
#          indices := Compacted([1..Length(vecs[i])]);
#          output := List(ch[i].params[2]*indices,x->vecs[i][x]);
#          vecs[i+1]:=output;
#      fi;
#    od;

    ## HACK For the A Cross I that were I is not computed
    ## Could be (advantageously) replaced by a Wire(Id) when
    ## the new system will be around

    # Disabled by Hao because it is broken in some way#
    for i in [2..numch] do
       if ObjId(ch[i])=Cross then
           for j in [1..Length(ch[i]._children)] do
              if (ObjId(ch[i]._children[j])=Prm and
                          ObjId(ch[i]._children[j].func)=fId and
                          ObjId(vecs[i][j].t)=TArray)or
                      (ObjId(ch[i]._children[j])=BB and
                          ObjId(ch[i]._children[j]._children[1])=Prm and
                          ObjId(ch[i]._children[j]._children[1].func)=fId and
                          ObjId(vecs[i][j].t)=TArray) or
                      (ObjId(ch[i]._children[j])=BB and
                          ObjId(ch[i]._children[j]._children[1])=I and
                          ObjId(vecs[i][j].t)=TArray) then
	      skip();
              #vecs[i][j]:=vecs[i+1][j];
          fi;
           od;
       fi;
    od;

    [vecs, ch] := [Reversed(vecs), Reversed(ch)];
    return decl( Difference(Flat(vecs{[2..Length(vecs)-1]}), Flat([x,y])),
        chain( List([1..numch], i -> When(vecs[i+1]=vecs[i],
            self(ch[i], vecs[i],   vecs[i], CopyFields(opts, rec(_inplace:=true))),
            self(ch[i], vecs[i+1], vecs[i], opts)))));
    end,

    ScatInit := (self, o, y, x, opts) >> let(ii := Ind(),
    condition := When(Length(o.cond)=0,V(1),FoldL1(List(o.cond,x->eq(x,V(0))),logic_and)),
    chain(
    IF(condition,
        loop(ii, o.func.domain(), assign(nth(y,o.func.at(ii)), When(ObjId(y)=var,y.t.t.zero(),y.t.zero()))),
        skip()),
    self(o._children,y,x,opts))
    ),

    Wire := (self, o, y, x, opts) >> skip(),

    ScatAcc := (self, o, y, x, opts) >>
        self._acc(self(Scat(o.func),y,x,opts),y),

    ICScatAcc := meth(self,o,y,x,opts)
      local i, func;
      i := Ind(); func := o.func.lambda();
      return loop(i, o.func.domain(), chain(assign_acc(nth(y,func.at(i)), nth(x, i)),
                        assign_acc(nth(y,add(func.at(i),1)),nth(x,add(i,1)))));
    end,

     # NOTE: handle unaligned case - IsUnalignedPtrT(y)
     VScatAcc := (self, o, y, x, opts) >>
         self._acc(self(VScat(o.func,o.v),y,x,opts),y),

     VTensor_OL := meth(self, o, y, x, opts)
       local CastToVect;
       CastToVect:=x->tcast(TPtr(TVect(x.t.t, o.vlen)), x);
       return self(o.child(1), StripList(List(Flat([y]),t->CastToVect(t))), StripList(List(Flat([x]),t->CastToVect(t))), opts);
     end,

##########
### to FIX: this Multiplication implicitely takes 2 inputs
    Multiplication:= meth(self, o, y, x, opts)
    local iterator;

    iterator:=Ind();
    return loop(iterator, [ 0 .. o.element[2]-1 ],
            assign(nth(StripList(y), iterator), mul(nth(x[1], iterator),nth(x[2], iterator))));
    end,

    ICMultiplication:=meth(self,o,y,x,opts)
    local len,iterator;

    iterator:=Ind();
    len := When(o.element[2] = 2,1,o.element[2]);
    len := When(Length(o.element) = 3,len/2,len);
    return loop(iterator, [ 0 .. (len)- 1 ],
            chain(
                assign(nth(StripList(y),mul(iterator,2)),sub(mul(nth(x[1], mul(iterator,2)),nth(x[2],mul(iterator,2))),
                                 mul(nth(x[1], add(mul(iterator,2),1)),nth(x[2], add(mul(iterator,2),1))))),
                assign(nth(StripList(y), add(mul(iterator,2),1)),add(mul(nth(x[1], mul(iterator,2)),nth(x[2], add(mul(iterator,2),1))),
                                             mul(nth(x[1], add(mul(iterator,2),1)),nth(x[2],mul(iterator,2)))))));

    end,

    Addition:= meth(self, o, y, x, opts)
        local iterator;
        iterator:=Ind();
        return loop(iterator, [ 0 .. _unwrap(o.element[2])-1 ],
            assign(nth(StripList(y), iterator), add(nth(x[1], iterator),nth(x[2], iterator))));
    end,

   Subtraction:= meth(self, o, y, x, opts)
        local iterator;
        iterator:=Ind();
        return loop(iterator, [ 0 .. o.element[2]-1 ],
            assign(nth(StripList(y), iterator), sub(nth(x[1], iterator),nth(x[2], iterator))));
    end,

    Codelet := meth(self, o, y, x, opts)
        local code,inputlist,outputlist,outputTypes,inputTypes;

    o := o.child(1);
        o := SubstBottomUp(o,BB,e->e.rChildren()[1]);
        o := OLQuickAndDirtyHackForCodelet(o);

        o:=ApplyStrategy(o,opts.formulaStrategies.postProcess, UntilDone, opts);
        o := OLRulesBufferFinalize(o);
        outputTypes := opts.OutputTypes;
        inputTypes := opts.InputTypes;
        inputlist :=[];
        outputlist :=[];
        for i in [1..DimLength(o.dims()[1])] do
          Append(outputlist, [TPtr(outputTypes[i])]);
        od;
        for i in [1..DimLength(o.dims()[2])] do
          Append(inputlist, [TPtr(inputTypes[i])]);
        od;
        Unification(o, StripList(inputlist));

        if List(outputlist,l->l.t)<>List(o.rng(),l->l.t) then
          Error("Unification of final output failed");
        fi;

    ## Generating code : main body
#        o := BlockSums(opts.libgen.basesUnrolling, o);
    code := SReduce(self(o, y, x, opts), opts);
        code := ESReduce(code, opts);
    code := RemoveAssignAcc(code);
    code := BlockUnroll(code, opts);
    code := DeclareHidden(code);
    return code;
    end,

   Formula := meth(self, o, y, x, opts)
        local code, init_code, datas, prog, params, sub, inputTypes, outputTypes,
    initsub,i,inputlist,outputlist,smp, num_threads, buffers, bufalloc, dalloc, dvars, map,ignore, codelet_codes, codelet_recs, chash;

    o := o.child(1);
        o := OLRulesBufferFinalize(o);
    o := OLRulesCode(o);

    params := Set(Collect(o, param));
    smp := Collect(o, SMPSum);
    num_threads := When(smp=[], 1, smp[1].p);
    if not ForAll(smp, x->x.p=num_threads) then Error("Non-uniform num_threads in SMPSum's"); fi;
        outputTypes := opts.OutputTypes;
        inputTypes :=  opts.InputTypes;
    inputlist :=[];
    outputlist :=[];
    for i in [1..DimLength(o.dims()[1])] do
      Append(outputlist, [TPtr(outputTypes[i])]);
    od;
    for i in [1..DimLength(o.dims()[2])] do
      Append(inputlist, [TPtr(inputTypes[i])]);
    od;
    Unification(o, StripList(inputlist));

    if List(outputlist,l->l.t)<>List(o.rng(),l->l.t) then
      Error("Unification of final output failed");
    fi;

    inputlist :=[];
    outputlist :=[];
    for i in [1..DimLength(o.dims()[1])] do
      Append(outputlist, [var(ConcatenationString("Y",String(i)), TPtr(outputTypes[i]))]);
    od;
    for i in [1..DimLength(o.dims()[2])] do
      Append(inputlist, [var(ConcatenationString("X",String(i)), TPtr(inputTypes[i]))]);
    od;

        #We're dropping y and x right away because previous stuff is backward compatibility
    y:=StripList(outputlist);
        x:= StripList(inputlist);

        ## HACK for the unroll of ROIs
        if IsBound(opts.libgen) then
            chash:=spiral.libgen.CreateCodeletHashTable();
            for i in Flat(opts.libgen.codeletTab.entries) do
                if Collect(i.data.sums,Multiplication)=[] then
                    i.data.unrolling:=300;         #this is the unrolling threshold for non kernels
                fi;
                HashAdd(chash,i.key, i.data);
            od;
            opts.libgen.codeletTab:=chash;
        elif Collect(o,Multiplication)=[] then
            o := SubstBottomUp(o,BB,e->e.rChildren()[1]);
            SubstTopDownNR(o, @(1).cond(IsMarkedBlock),
                function(e)
                    local f;
                    f:=@(1).val;
                    f.isBlock:=false;
                    return f;
                end);
            o := compiler.BlockSumsOpts(o,CopyFields(opts,rec(globalUnrolling:=32)));
        fi;

    ## Generating code : codelets
    codelet_recs := spiral.libgen.CompileCodelets(o, opts);
    codelet_codes := List(codelet_recs,
            function(clrec)
              local inputs, outputs, myf;
              outputs:=StripList(List([1..Length(Flat([clrec.sums.dims()[1]]))],x->var(Concat("Y",When(x>1,String(x),"")), TPtr(outputTypes[x]))));
              inputs:=StripList(List([1..Length(Flat([clrec.sums.dims()[2]]))],x->var(Concat("X",When(x>1,String(x),"")), TPtr(inputTypes[x]))));
              myf:=func(TVoid, clrec.name, Concatenation(Flat([outputs, inputs]), clrec.params), clrec.code);
              myf.inline:=true;
              return myf;
            end
        );

    map := tab();
    ## Generating code : main body
    datas := Collect(o, FDataOfs);
        code := self(o, y, x, opts);
    code := SReduce(code, opts);
        code := ESReduce(code, opts);
    code := RemoveAssignAcc(code);
    code := BlockUnroll(code, opts);
#   code := SimplifySingularArray(code);
        # code := PowerOpt(code);
    code := DeclareHidden(code);
    if IsBound(opts.isFixedPoint) and opts.isFixedPoint then
        code := FixedPointCode(code, opts.bits, opts.fracbits);
    fi;
    [buffers, bufalloc, code] := spiral.libgen.AllocBuffersSMP("buf", code, map, opts);

    ## Generating code : initialization
    [dvars, dalloc, ignore] := spiral.libgen.AllocBuffersSMP("dat", decl(datas, skip()), map, opts);

    # NOTE: is there a better way?
    sub := opts.subName;
    initsub := opts.subInitName;
    code := func(TVoid, sub, Concatenation(params,
                When(IsBound(opts.subParams), opts.subParams, []), [y,x]), code);

        init_code := func(TVoid, initsub, [],
            chain(bufalloc, dalloc, List(datas, x -> SReduce(x.var.init, opts))));

        prog := SubstVars(chain(init_code, code), map);

        if Length(Filtered(opts.subParams,x->x.id="num_threads"))>0 then
            prog := data(var("NUM_THREADS", TInt), V(num_threads),prog);
        fi;

        prog := program(
            codelet_codes,
            decl(Concatenation(dvars, buffers, List(datas, x->x.var)),
        prog));

        prog:= GenerateBench(o, prog,opts);

    return FlattenCode(prog);
    end,
));

DefaultCodegen.ExpDiag := (self, o, y, x, opts) >> let(
   i   := Ind(),     d := o.d,
   elt := o.element, s := Length(elt.vars),
   xx  := List([0..s-1], j -> nth(x, i*s + j)),
   loop(i, d,
       assign(nth(y,i), elt.at(xx)))
);


text:
local reps,tori,labels,I,II,ords,blocks;

reps:=[[[1,1],[2,1],[3,1],[4,1],[5,1],[6,1],[7,1],[8,1]],#E_8
[[1,1],[2,1],[10,1],[11,1],[5,1],[6,1],[7,1],[8,1]],#E_8(a_1)
[[1,1],[2,1],[3,1],[10,1],[12,1],[13,1],[14,1],[8,1]],#E_8(a_2)
[[7,1],[8,1],[9,1],[10,1],[11,1],[12,1],[13,1],[19,1]],#E_8(a_3)
[[9,1],[10,1],[11,1],[12,1],[13,1],[14,1],[15,1],[19,1]],#E_8(a_4)
[[1,1],[2,1],[3,1],[4,1],[5,1],[6,1],[7,1]],#E_7
[[7,1],[8,1],[9,1],[17,1],[18,1],[19,1],[20,1],[21,1]],#E_8(b_4)
[[9,1],[14,1],[15,1],[17,1],[18,1],[19,1],[20,1],[26,1]],#E_8(a_5)
[[1,1],[3,1],[10,1],[11,1],[5,1],[6,1],[7,1]],#E_7(a_1)
[[8,1],[16,1],[17,1],[18,1],[20,1],[21,1],[24,1],[33,1]],#E_8(b_5)
[[2,1],[3,1],[4,1],[5,1],[6,1],[7,1],[8,1]],#D_7
[[4,1],[10,1],[12,1],[14,1],[15,1],[26,1],[27,1],[30,1]],#E_8(a_6)
[[1,1],[2,1],[3,1],[10,1],[12,1],[13,1],[14,1]],#E_7(a_2)
[[1,1],[2,1],[3,1],[4,1],[5,1],[6,1],[8,1]],#E_6+A_1
#[[2,1],[3,1],[10,1],[12,1],[6,1],[7,1],[8,1],[83,1]],#(D_7(a_1))^{(2)}
[[2,1],[3,1],[10,1],[12,1],[6,1],[7,1],[8,1]],#D_7(a_1)
[[15,1],[22,1],[23,1],[24,1],[25,1],[26,1],[27,1],[28,1]],#E_8(b_6)
[[1,1],[4,1],[7,1],[12,1],[13,1],[17,1],[19,1]],#E_7(a_3)
[[1,1],[2,1],[11,1],[12,1],[5,1],[6,1],[8,1]],#E_6(a_1)+A_1
#[[1,1],[3,1],[4,1],[5,1],[6,1],[7,1],[8,1],[67,1]],#(A_7)^{(3)}
[[1,1],[3,1],[4,1],[5,1],[6,1],[7,1],[8,1]],#A_7
[[2,1],[3,1],[10,1],[12,1],[13,1],[14,1],[8,1]],#D_7(a_2)
[[1,1],[2,1],[3,1],[4,1],[5,1],[6,1]],#E_6
[[2,1],[3,1],[4,1],[5,1],[6,1],[7,1]],#D_6
#[[1,1],[2,1],[3,1],[4,1],[5,1],[7,1],[8,1],[51,1]],#(D_5+A_2)^{(2)}
[[1,1],[2,1],[3,1],[4,1],[5,1],[7,1],[8,1]],#D_5+A_2
[[1,1],[2,1],[11,1],[12,1],[5,1],[6,1]],#E_6(a_1)
[[1,1],[4,1],[7,1],[11,1],[25,1],[26,1],[27,1]],#E_7(a_4)
[[1,1],[2,1],[4,1],[5,1],[6,1],[7,1],[8,1]],#A_6+A_1
[[2,1],[3,1],[10,1],[12,1],[6,1],[7,1]],#D_6(a_1)
[[1,1],[3,1],[4,1],[5,1],[6,1],[7,1]],#A_6
[[30,1],[33,1],[35,1],[36,1],[40,1],[45,1],[53,1],[58,1]],#E_8(a_7)
[[1,1],[2,1],[3,1],[4,1],[5,1],[7,1]],#D_5+A_1
[[4,1],[11,1],[14,1],[16,1],[24,1],[25,1],[26,1]],#E_7(a_5)
[[4,1],[6,1],[9,1],[11,1],[18,1],[19,1],[8,1]],#E_6(a_3)+A_1
[[2,1],[3,1],[10,1],[12,1],[13,1],[14,1]],#D_6(a_2)
[[1,1],[10,1],[11,1],[2,1],[5,1],[7,1],[8,1]],#D_5(a_1)+A_2
[[1,1],[3,1],[4,1],[5,1],[6,1],[8,1]],#A_5+A_1
[[1,1],[3,1],[4,1],[2,1],[6,1],[7,1],[8,1]],#A_4+A_3
[[1,1],[2,1],[3,1],[4,1],[5,1]],#D_5
[[4,1],[6,1],[9,1],[11,1],[18,1],[19,1]],#E_6(a_3)
#[[2,1],[3,1],[4,1],[5,1],[7,1],[8,1],[63,1]],#(D_4+A_2)^{(2)}
[[2,1],[3,1],[4,1],[5,1],[7,1],[8,1]],#D_4+A_2
[[1,1],[3,1],[2,1],[5,1],[6,1],[7,1],[8,1]],#A_4+A_2+A_1
[[1,1],[10,1],[11,1],[2,1],[5,1],[7,1]],#D_5(a_1)+A_1
[[1,1],[3,1],[4,1],[5,1],[6,1]],#A_5
[[1,1],[3,1],[4,1],[2,1],[6,1],[7,1]],#A_4+A_2
[[1,1],[3,1],[4,1],[2,1],[6,1],[8,1]],#A_4+2A_1
[[1,1],[10,1],[11,1],[2,1],[5,1]],#D_5(a_1)
[[1,1],[3,1],[4,1],[6,1],[7,1],[8,1]],#2A_3
[[1,1],[3,1],[4,1],[5,1],[7,1]],#A_4+A_1
[[2,1],[3,1],[10,1],[12,1],[7,1],[8,1]],#D_4(a_1)+A_2
[[2,1],[3,1],[4,1],[5,1],[7,1]],#D_4+A_1
[[1,1],[3,1],[2,1],[5,1],[6,1],[7,1]],#A_3+A_2+A_1
[[1,1],[3,1],[4,1],[5,1]],#A_4
#[[1,1],[3,1],[4,1],[6,1],[7,1],[44,1]],#(A_3+A_2)^{(2)}
[[1,1],[3,1],[4,1],[6,1],[7,1]],#A_3+A_2
[[2,1],[3,1],[10,1],[12,1],[7,1]],#D_4(a_1)+A_1
[[1,1],[3,1],[4,1],[6,1],[8,1]],#A_3+2A_1
[[1,1],[3,1],[2,1],[5,1],[6,1],[8,1]],#2A_2+2A_1
[[2,1],[3,1],[4,1],[5,1]],#D_4
[[2,1],[3,1],[10,1],[12,1]],#D_4(a_1)
[[1,1],[3,1],[4,1],[6,1]],#A_3+A_1
[[1,1],[3,1],[2,1],[5,1],[6,1]],#2A_2+A_1
[[1,1],[3,1],[5,1],[6,1]],#2A_2
[[1,1],[3,1],[2,1],[5,1],[7,1]],#A_2+3A_1
[[1,1],[3,1],[4,1]],#A_3
[[1,1],[3,1],[5,1],[7,1]],#A_2+2A_1
[[1,1],[3,1],[5,1]],#A_2+A_1
[[1,1],[2,1],[5,1],[7,1]],#4A_1
[[1,1],[3,1]],#A_2
[[1,1],[2,1],[5,1]],#3A_1
[[1,1],[2,1]],#2A_1
[[1,1]]#A_1
];

labels:=["E_8",
"E_8(a_1)",
"E_8(a_2)",
"E_8(a_3)",
"E_8(a_4)",
"E_7",
"E_8(b_4)",
"E_8(a_5)",
"E_7(a_1)",
"E_8(b_5)",
"D_7",
"E_8(a_6)",
"E_7(a_2)",
"E_6+A_1",
#"(D_7(a_1))^{(2)}",
"D_7(a_1)",
"E_8(b_6)",
"E_7(a_3)",
"E_6(a_1)+A_1",
#"(A_7)^{(3)}",
"A_7",
"D_7(a_2)",
"E_6",
"D_6",
#"(D_5+A_2)^{(2)}",
"D_5+A_2",
"E_6(a_1)",
"E_7(a_4)",
"A_6+A_1",
"D_6(a_1)",
"A_6",
"E_8(a_7)",
"D_5+A_1",
"E_7(a_5)",
"E_6(a_3)+A_1",
"D_6(a_2)",
"D_5(a_1)+A_2",
"A_5+A_1",
"A_4+A_3",
"D_5",
"E_6(a_3)",
#"(D_4+A_2)^{(2)}",
"D_4+A_2",
"A_4+A_2+A_1",
"D_5(a_1)+A_1",
"A_5",
"A_4+A_2",
"A_4+2A_1",
"D_5(a_1)",
"2A_3",
"A_4+A_1",
"D_4(a_1)+A_2",
"D_4+A_1",
"A_3+A_2+A_1",
"A_4",
#"(A_3+A_2)^{(2)}",
"A_3+A_2",
"D_4(a_1)+A_1",
"A_3+2A_1",
"2A_2+2A_1",
"D_4",
"D_4(a_1)",
"A_3+A_1",
"2A_2+A_1",
"2A_2",
"A_2+3A_1",
"A_3",
"A_2+2A_1",
"A_2+A_1",
"4A_1",
"A_2",
"3A_1",
"2A_1",
"A_1"
];

I:=[[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[]];#,[],[],[],[],[]];

II:=[[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[]];#,[],[],[],[],[]];

tori:=[[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[]];#,[],[],[],[],[]];

blocks:=[];

#ords:=[];
ords:=[ 125, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 5, 25, 25, 25, 5, 25, 25, 5, 5, 25, 5, 5, 5, 25, 5, 5, 5, 5, 5, 5, 25, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5 ];

return [reps,tori,labels,I,II,ords,blocks];


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

Class(CConv, TaggedNonTerminal, rec(
    abbrevs := [
        (n,h) -> Checked(IsPosIntSym(n),
            [_unwrap(n), h, false]),
        (n,h,isFreqData) -> Checked(IsPosIntSym(n),
            [_unwrap(n), h, isFreqData])
        ],
    dims      := self >> [self.params[1], self.params[1]],
    terminate := self >> let(size := self.params[1],
        taps := When(self.params[7], MatSPL(DFT(self.params[1], -1)) * When(IsList(self.params[2]), self.params[2], self.params[2].list), self.params[2]),
        spiral.transforms.filtering.Circulant(size, taps, -size).terminate()),
    isReal    := self >> false,
    normalizedArithCost := self >> let(n := self.params[1], IntDouble(5 * n * d_log(n) / d_log(2))),
    TType := TComplex,
    hashAs := self >> ApplyFunc(ObjId(self), [self.params[1], fUnk(self.params[2].range(), self.params[2].domain())]::Drop(self.params, 2))
));

NewRulesFor(CConv, rec(
    CConv_DFT_IDFT := rec(
       forTransposition := false,
       applicable :=  (self, nt) >> not nt.hasTags(),
       children := nt -> [[ DFT(nt.params[1], 1),
                            DFT(nt.params[1], -1) ]],
       apply := (nt, C, cnt) -> Cond(ObjId(nt.params[2]) = var, # a variable is frequency domain data -> FData
                                    C[1] * Diag(FData(nt.params[2])) * C[2], # if data is provided it needs to be frequency domain
                                    ObjId(nt.params[2]) = Lambda,  # data is inside a larger array
                                    C[1] * Diag(nt.params[2]) * C[2],
                                    nt.params[7],   # if params[7] is true then data is frequency domain data
                                    C[1] * Diag(When(IsList(nt.params[2]), FList(TComplex, nt.params[2]), nt.params[2])) * C[2],
                                    let(cxfftdiag := 1/nt.params[1]*ComplexFFT(List(nt.params[2].tolist(), i->ComplexAny(_unwrap(i)))),
                                        C[1] * Diag(FData(cxfftdiag)) * C[2]))
    )
));

NewRulesFor(MDRConv, rec(
    MDPRConv_Base := rec(
        info := "MDRConv -> IMDPRDFT * RCDiag(sym) * MDPRDFT",
        applicable     := nt -> nt.params[3], # need symbol to be in frequency domain

        children       := nt -> [[ MDPRDFT(nt.params[1], -1).withTags(nt.getTags()),
                                   IMDPRDFT(nt.params[1], 1).withTags(nt.getTags())
                                 ]],

        # nonterminal, children, children non-terminals
        apply          := (nt, C, Nonterms) -> C[2] *
                            RCDiag(FDataOfs(nt.params[2], 2*Product(DropLast(nt.params[1], 1))* (Last(nt.params[1])/2+1), 0)) *
                            C[1]
    ),

    MDRConv_3D_2trip_zyx_freqdata := rec(
       forTransposition := false,
       applicable :=  (self, nt) >> not nt.hasTags() and Length(nt.params[1]) = 3 and nt.params[3],
       children := nt -> let(nlist := nt.params[1],
                            diag := nt.params[2],
                            nfreq := nlist[1]/2+1,
                            i := Ind(nfreq*nlist[2]),
                            hfunc := Cond(ObjId(diag) = Lambda,
                                let(j := Ind(nlist[3]),
                                    # Lambda(j, cxpack(diag.at(2*(j + i*nlist[3])), diag.at(2*(j + i*nlist[3])+1)))
                                    pos := i +j*nfreq*nlist[2],
                                    Lambda(j, cxpack(diag.at(2*pos), diag.at(2*pos+1)))
                                ),
                                IsVar(diag),
                                let(j := Ind(nlist[3]),
                                    pos := i +j*nfreq*nlist[2],
                                    Lambda(j, cxpack(nth(diag, 2*pos), nth(diag, 2*pos+1))).setRange(TComplex)
                                ),
                                ObjId(diag) = fUnk,
                                fUnk(TComplex, nlist[3]),
                                let(list := nt.params[2].list,  # here we assume FList(TReal, [...])
                                    clist := List([1..Length(list)/2], i->Cplx(list[2*i-1], list[2*i])),
                                    fc := FList(TComplex, clist),
                                    gf := fTensor(fBase(i), fId(nlist[3])),
                                    fCompose(fc, gf)
                                )
                            ),
                            [[ PRDFT(nlist[1], -1),  # stage 1: PRDFT z
                               DFT(nlist[2], -1),    # stage 2: DFT y
                               CConv(nlist[3], hfunc), # stage 3+4+5: complex conv in x
                               DFT(nlist[2], 1), # stage 6: iDFT in y
                               IPRDFT(nlist[1], 1),   # stage 7: iPRDFT in z
                               InfoNt(i)
                            ]]),

       apply := (nt, C, cnt) -> let(prdft1d := C[1],
                                    pdft1d := C[2],
                                    iopconv := C[3],
                                    ipdft1d := C[4],
                                    iprdft1d := C[5],
                                    i := cnt[6].params[1],
                                    nlist := nt.params[1],
                                    n1 := nlist[1],
                                    nfreq := nlist[1]/2+1,
                                    n2 := nlist[2],
                                    n3 := nlist[3],
                                    stage1 := L(2*nfreq*n3*n2, n3) * Tensor(I(n2), Tensor(L(2*nfreq, 2) * prdft1d, I(n3))) * Tensor(L(n2*n1, n2), I(n3)),
                                    stage2 := Tensor(I(n3), Tensor(RC(pdft1d), I(nfreq))),
                                    pp := Tensor(L(n3*n2*nfreq, n2*nfreq), I(2)) * Tensor(I(n3), L(2*nfreq*n2, nfreq)),
                                    ppi := Tensor(I(n3), L(2*nfreq*n2, 2*n2)) * Tensor(L(n3*n2*nfreq, n3), I(2)),
                                    stage543 := Grp(ppi * IDirSum(i, RC(iopconv)) * pp),
                                    stage76 := Tensor(L(n2*n1, n1), I(n3)) * Grp(Tensor((Tensor(I(n2), iprdft1d * L(2*nfreq, nfreq)) *
                                       Tensor(RC(ipdft1d), I(nfreq))), I(n3)) * L(2*nfreq*n3*n2, 2*nfreq*n2)),
                                    conv3dr := stage76 * stage543 * stage2 * stage1,
                                    conv3dr
                            )
    )
));

NewRulesFor(MDRConvR, rec(
    MDPRConvR_Base := rec(
        info := "MDRConv -> IMDPRDFT * Diag(sym) * MDPRDFT",
        applicable     := nt -> nt.params[3], # need symbol to be in frequency domain

        children       := nt -> [[ MDPRDFT(nt.params[1], -1).withTags(nt.getTags()),
                                   IMDPRDFT(nt.params[1], 1).withTags(nt.getTags())
                                 ]],

        # nonterminal, children, children non-terminals
        apply          := (nt, C, Nonterms) -> C[2] *
                            Diag(diagTensor(FDataOfs(nt.params[2], Product(DropLast(nt.params[1], 1))* (Last(nt.params[1])/2+1), 0), fConst(TReal, 2, 1))) *
                            C[1]
    )
));

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

NewRulesFor(TRaderMid, rec(
    Pad_vec := rec(
        applicable := (self, t) >> t.isTag(1, AVecReg) or t.isTag(1, AVecRegCx),
        forTransposition := false,
        apply := (t, C, Nonterms) -> let(
            v := t.firstTag().v,
            ds := Rows(t)-1,
            When(IsInt(ds/v),
                DelayedDirectSum(VScat_sv(fId(1), v , 1), I(ds))
                * VecRaderMid(t.params[1], t.params[2], t.params[3], v)
                * DelayedDirectSum(VGath_sv(fId(1), v , 1), I(ds)),
                DelayedDirectSum(VScat_sv(fId(1), v , 1), VScat_sv(fId(ds), v, 1))
                * VecRaderMid(t.params[1], t.params[2], t.params[3], v)
                * DelayedDirectSum(VGath_sv(fId(1), v , 1), VGath_sv(fId(ds), v, 1))
            )
        )
    )
));


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

setLeft := function(tags)
   local i, retval;
   retval := Copy(tags);
   for i in retval do
      i.isLeftChild := true;
      i.isRightChild := false;
 od;
 return(retval);
end;

setRight := function(tags)
   local i, retval;
   retval := Copy(tags);
   for i in retval do
      i.isLeftChild := false;
      i.isRightChild := true;
 od;
 return(retval);
end;

NewRulesFor(WHT, rec(
    WHT_GT := rec(
    switch := false,
    maxSize := false,
    minSize := false,
    codeletSize := false,
    inplace := false,

    applicable := (self, t) >> let(n := Rows(t),
        n > 2 and
        (self.maxSize=false or n <= self.maxSize) and
        (self.minSize=false or n >= self.minSize) and
        not IsPrime(n)),

        children := (self, t) >> Map2(Filtered(DivisorPairs(Rows(t)), d->When(IsInt(self.codeletSize), d[1]<=self.codeletSize, true)),
        (m,n) -> [
        GT(WHT(LogInt(m, 2)), XChain([0, 1]), XChain([0, 1]), [n]).withTags(t.getTags()),
        GT(WHT(LogInt(n, 2)), XChain([1, 0]), XChain([1, 0]), [m]).withTags(t.getTags()),
        #GT(WHT(LogInt(m, 2)), XChain([0, 1]), XChain([0, 1]), [n]).withTags( let(ts := t.getTags(),  setLeft(ts)) ),
        #GT(WHT(LogInt(n, 2)), XChain([1, 0]), XChain([1, 0]), [m]).withTags( let(ts := t.getTags(), setRight(ts)) ),
        ]),

    apply := (self, t, C, Nonterms) >>  C[1] * C[2]
    )
));

NewRulesFor(DFT, rec(
    DFT_GT_CT := rec(
    switch := false,
    maxSize := false,
    minSize := 2,
    codeletSize := false,
    inplace := false,

    a := rec(
        precompute := true
    ),

    applicable := (self, t) >> let(n := Rows(t),
        n > 2 and
        (self.maxSize=false or n <= self.maxSize) and
        (self.minSize=false or n >= self.minSize) and
        not IsPrime(n)),

        #GT(DFT(m, t.params[2] mod m), XChain([0, 1]), XChain([0, 1]), [n]).withTags(  let(a := t.getTags(),  ),
        #GT(DFT(n, t.params[2] mod n), XChain([0, 1]), XChain([1, 0]), [m]).withTags(  t.getTags())

        #GT(DFT(m, t.params[2] mod m), XChain([0, 1]), XChain([0, 1]), [n]).withTags(t.getTags()),
        #GT(DFT(n, t.params[2] mod n), XChain([0, 1]), XChain([1, 0]), [m]).withTags(t.getTags())

    children := (self, t) >> Map2(Filtered(DivisorPairs(Rows(t)), d->When(IsInt(self.codeletSize), d[1]<=self.codeletSize, true)),
        (m,n) -> [
        GT(DFT(m, t.params[2] mod m), XChain([0, 1]), XChain([0, 1]), [n]).withTags(t.getTags()),
        GT(DFT(n, t.params[2] mod n), XChain([0, 1]), XChain([1, 0]), [m]).withTags(t.getTags()),
        #GT(DFT(m, t.params[2] mod m), XChain([0, 1]), XChain([0, 1]), [n]).withTags( let(ts := t.getTags(),  setLeft(ts)) ),
        #GT(DFT(n, t.params[2] mod n), XChain([0, 1]), XChain([1, 0]), [m]).withTags( let(ts := t.getTags(), setRight(ts)) ),
         ]),

    apply := (self, t, C, Nonterms) >>  let(
            inplace := When(self.inplace, Inplace, Grp),
            n := Rows(Nonterms[2].params[1]),
            rot := t.params[2],
            compute := When(self.a.precompute, fPrecompute, fComputeOnline),

            inplace(Buf(C[1] * Diag(compute(Tw1(Rows(t), n, rot))))) * C[2])
    )
));

# GT(DFT(), ...) rules
#
NewRulesFor(GT, rec(
    GT_DFT_Base2 := rec(
    applicable := (self, t) >> let(rank := Length(t.params[4]),
        rank = 0 and PatternMatch(t, [GT, DFT, @, @, @, @, @], empty_cx()) and Rows(t.params[1])=2),
    apply := (t, C, Nonterms) -> F(2)
    ),

    GT_DFT_CT := rec(
        maxSize       := false,
        minSize       := false,
        minRank       := 0,
        maxRank       := 1,
        codeletSize := 32,
        inplace := true,

        a := rec(
            precompute := true
        ),

        applicable := (self, t) >> let(
            rank := Length(t.params[4]),
            dft := t.params[1],

            rank >= self.minRank
            and rank <= self.maxRank
            and When(rank>0, t.getTags()=[], true)
            and (self.maxSize=false or Rows(dft) <= self.maxSize)
            and (self.minSize=false or Rows(dft) >= self.minSize)
            and PatternMatch(t, [GT, DFT, XChain, XChain, @, @, @], empty_cx())
            and DFT_GT_CT.applicable(dft)
        ),

        children := (self, t) >> let(
            dft := t.params[1],
            g := t.params[2], s := t.params[3],
            rot := dft.params[2],
            loop_dims := t.params[4],
            nloops := Length(loop_dims),
            tags := t.getTags(),

            Map2( Filtered(DivisorPairs(Rows(dft)), d->d[1]<=self.codeletSize),
                (m,n) -> [
                    GT(DFT(m, rot mod m),
                        s.composeWith(XChain([0, 1])), s.composeWith(XChain([0, 1])),
                        Concatenation([n], loop_dims)
                    ).withTags(tags),

                    GT(DFT(n, rot mod n),
                        g.composeWith(XChain([0, 1])), s.composeWith(XChain([1, 0])),
                        Concatenation([m], loop_dims)
                    ).withTags(tags)
                ]
            )
        ),

        apply := (self, t, C, Nonterms) >> let(
            loop_dims := t.params[4],
            s := t.params[3],
            N := Rows(t.params[1]),
            k := t.params[1].params[2],
            n := Rows(Nonterms[2].params[1]),

            inplace := When(self.inplace, Inplace, Grp),
            compute := When(self.a.precompute, fPrecompute, fComputeOnline),

            inplace(Buf(C[1] * Diag(compute(s.toDiag(loop_dims, Tw1(N,n,k)))))) * C[2]
        )
    )
));


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

NewRulesFor(TICompose, rec(
    TICompose_DropTag := rec(
        applicable := t -> true,
        children := t -> [[ t.params[3] ]],
        apply := (t, C, Nonterms) -> ICompose(t.params[1], t.params[2], C[1])
    )
));

#F InitStreamSw(<radix>)
#F
#F Example:
#F   opts := InitStreamSw(2);
#F   c := CodeRuleTree(opts.transforms.DFT(64), opts);
#F   me   := CMatrix(c, opts);
#F   them := MatSPL(opts.transforms.DFT(64));
#F   InfinityNormMat(me-them);
#F
InitStreamSw := radix -> CopyFields(
    CplxSpiralDefaults,
    IntelC99Mixin, # uncomment this for C99 output (default = SSE2 intrinsics)
    rec(
    breakdownRules := rec(
        DFTDR :=  [CopyFields(DFTDR_tSPL_Pease, rec(precompute:=true)) ],
        DFT   :=  [DFT_Base, DFT_CT ],
        TTensorI := [IxA_base, IxA_L_base, L_IxA_base, AxI_base],
        TCompose := [TCompose_tag],
        TICompose := [TICompose_DropTag],
        TDiag    := [TDiag_base]
    ),

    transforms := rec(
        DFT := n -> DFTDR(n, radix, [AStream(radix)])
    ),

#    unparser := CMacroUnparserProg,

    globalUnrolling := radix
));

#F TestStreamSw(<n>, <radix>, <use_transpose>)  - software test of Pease DFT used for streaming hardware
#F
#F use_transpose = true:  DRDFT
#F use_transpose = false: DFTDR
#F
TestStreamSw := function(n, radix, use_transpose)
    local opts, t, rt, c, err, me, them;
    opts := InitStreamSw(radix);

    PrintErr("Compiling...\n");
    t := When(use_transpose, opts.transforms.DFT(n).transpose(), opts.transforms.DFT(n));
    rt := RandomRuleTree(t, opts);
    if rt = false
        then Error("No ruletrees."); fi;
    c := CodeRuleTree(rt, opts);

    PrintErr("Running and computing the code matrix...\n");
    me := CMatrix(c, opts);

    PrintErr("Comparing to the definition...\n");
    them := MatSPL(t);
    err := InfinityNormMat(me-them);

    if err < 1e-6 then PrintErr(t, " -- ", GreenStr("OK\n"));
    else PrintErr(t, " -- ", RedStr("FAIL\n")); fi;
end;

TestStreamSwPrint := function(n, radix, use_transpose)
    local opts, t, rt, c, err, me, them;
    opts := InitStreamSw(radix);

    t := When(use_transpose, opts.transforms.DFT(n).transpose(), opts.transforms.DFT(n));
    rt := RandomRuleTree(t, opts);
    if rt = false
        then Error("No ruletrees."); fi;
    c := CodeRuleTree(rt, opts);

    return PrintCode("function_name", c, opts);
#    return rt;
end;