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

#############################################################################
#F
#F  sched.g                    Bag of Schedulers                Ernest Chan
#F
##

#############################################################################
##
##  Constants
##
PRINTDAG := false;
RED := 0;           # INPUT NODES
BLACK := 1;         # TEMP NODES
BLUE := 2;          # OUTPUT NODES
GREEN := 3;         # WEIRD! no pred and no succ

############################################################################
##
##  Building the DAG
##
##  Node structure
##    cmd:     the original statement
##    pred:    a list of predecessors
##    succ:    a list of successors
##    id:      list of command outputs
##

############################################################################
##
#F  AddNode( <dag>, <cmd> ) . . makes a node for a variable/input/output
##
AddNode := function(dag, cmd)
  local a;
  a := rec(
    id := Set(cmd.op_out() :: cmd.op_inout()),     # either var or nth
    index := Length(dag),
    pred := [],
    succ := [],
    cmd := cmd,
  );

  Add(dag, a);

  return a;
end;

############################################################################
##
#F  GetNode( <dag> , <id> ) . . . . . . if node with the specified <id> is
#F      in <dag>, then it is returned,  otherwise returns false
##
GetNode := function(dag, var)
  local i;

  i := Length(dag);
  while i >= 1 do
    if var in dag[i].id then
      return dag[i];
    fi;
    i := i-1;
  od;

  return false;
end;

############################################################################
##
#F  AddEdges( <dag> , <cmd> ) . . . parses out the variables from
#F      <exp>, and addes an edge from each of the variables to <node> in <dag>
##
AddEdges := function(dag, cmd)
    local v, n, ops_in, node;

    ops_in := Set(cmd.op_in() :: cmd.op_inout());

    node := AddNode(dag, cmd);

    for v in ops_in do
        n := GetNode(dag, v);
        if n <> false then
            Add(n.succ, node);
            Add(node.pred, n);
        fi;
    od;
end;

##############################################################################
##
#F  BuildDag( <chain> ) . . . . builds the Dag by examining
#F      each command statement in <chain>
##
BuildDag := function(chain_object)
    local i, j, k, dag, newnode, old_id, newnode_name, f, objid, assignstmt, node;
    if (ObjId(chain_object) = chain) then
        dag := [];
        for i in [1..Length(chain_object.cmds)] do
            # for each assign statement, make appropriate changes to the Dag have to
            #   1) add nodes if not in the graph
            #   2) add edges
            AddEdges(dag, chain_object.cmds[i]);
        od;
        return dag;
    else
        Print("in function BuildDag: argument not a chain structure");
    fi;
end;

##############################################################################
##
#F  PrintDag ( <dag> , <bool> ) . . . . if true, prints the dag with its orders,
#F      otherwise prints the dag with the names of the nodes
##
PrintDag := function(dag, order)
  local count, _PrintDag, __PrintDag;

  Print("digraph DAG {\n");

  count := 0;

  _PrintDag := function(dag, order)
    local i,j,templist;

    if Length(dag)<1 then return; fi;

    if IsList(dag[1]) then
      for i in [1..Length(dag)] do
        Print("  subgraph cluster");
        Print(count);
        count := count + 1;
        Print(" {\n");
        _PrintDag(dag[i], order);
        Print("  }\n");
      od;
    else
      for i in [1..Length(dag)] do
        Print("  ");
        if order then
          Print(dag[i].order);
        else
          Print("\"", dag[i].id, "\"");
        fi;
      od;
      Print("\n");
    fi;
  end;

  _PrintDag(dag, order);

  __PrintDag := function(dag,order)
    local i,j;

    if Length(dag)<1 then return; fi;

    if IsList(dag[1]) then
      for i in [1..Length(dag)] do
        __PrintDag(dag[i], order);
      od;
    else
      for i in [1..Length(dag)] do
        for j in [1..Length(dag[i].succ)] do
          Print("  ");
          if order then
            Print(dag[i].order);
          else
            Print("\"", dag[i].id, "\"");
          fi;
          Print(" -> ");
          if order then
            Print(dag[i].succ[j].order);
          else
            Print("\"", dag[i].succ[j].id, "\"");
          fi;
          Print("\n");
        od;
      od;
    fi;
  end;

  __PrintDag(dag, order);

  Print("}\n");
end;

##############################################################################
##
#F  OutputDag( <dag> , <string> ) . . . . . . . outputs a postscript file for
#F      the dag with filename <string>.ps
##
OutputDag := function(dag, file)
  if PRINTDAG then
    PrintTo("temp", PrintDag(dag, false));
    SYS_EXEC(Concat("dot -Tps -o",file," temp"));
    SYS_EXEC("rm temp");
  fi;
end;

##############################################################################
##
#F  InitColor ( <dag> )  . . . . . . . . initialize color of nodes in a dag
##
InitColor := function(dag)
  local i, j, pred_exists, succ_exists;

  for i in [1..Length(dag)] do
    pred_exists := Filtered(dag[i].pred, e->e in dag) <> [];
    succ_exists := Filtered(dag[i].succ, e->e in dag) <> [];

    if pred_exists then
      if succ_exists then
        dag[i].color := BLACK;
      else
        dag[i].color := BLUE;
      fi;
    else
      if succ_exists then
        dag[i].color := RED;
      else
        Error("Not a DAG");
      fi;
    fi;
  od;

end;

###############################################################################
##
#F  ColorPhase ( <dag> , <color> ) . . advances one step in the coloring phase
##
ColorPhase := function(dag, color)
  local i, j, pos, all_has_color, colorlist;

  colorlist := [];

  if color=RED then
    for i in [1..Length(dag)] do
      if dag[i].color = BLACK then
        all_has_color := true;
        for j in dag[i].pred do
          if j in dag then
            if j.color<>RED then
              all_has_color := false;
            fi;
          fi;
        od;
        if all_has_color then
          Add(colorlist, i);
        fi;
      fi;
    od;
  elif color=BLUE then
    for i in [1..Length(dag)] do
      if dag[i].color = BLACK then
        all_has_color := true;
        for j in dag[i].succ do
          if j in dag then
            if j.color<>BLUE then
              all_has_color := false;
            fi;
          fi;
        od;
        if all_has_color then
          dag[i].color := BLUE;
        fi;
      fi;
    od;
  fi;

  for i in [1..Length(colorlist)] do
    dag[colorlist[i]].color := color;
  od;

end;

###############################################################################
##
#F  SplitDag( <dag> ) . . . . . splits a dag into two (roughly equal) halves
##
SplitDag := function(dag)
  local getComponent, set, newset, temp, answer, i;

  getComponent := function(dag, set, node)
    local i, j, pred, succ, tempset;

    pred := Filtered(node.pred, e->e in dag);
    succ := Filtered(node.succ, e->e in dag);

    if not node in set then
      Add(set, node);
      node.taken := true;
    fi;
    for i in pred do
      if not i in set then
        Add(set, i);
        i.taken := true;
        set := getComponent(dag, set, i);
      fi;
    od;
    for i in succ do
      if not i in set then
        Add(set, i);
        i.taken := true;
        set := getComponent(dag, set, i);
      fi;
    od;
    return set;
  end;

  answer := [];
  for i in dag do
    i.taken := false;
  od;

  while not Length(Filtered(dag, x -> not x.taken)) = 0  do
    set := Filtered(dag, x -> not x.taken);
    temp := getComponent(dag, [], set[1]);
    Add(answer, temp);
  od;

  if Length(answer) = 1 then
    return rec(graph := answer[1], modified := false);
  else
    return rec(graph := answer, modified := true);
  fi;
end;

###############################################################################
##
#F  BreakDeps ( <dag> )  . . . . . . . . .  breaking one-to-many dependencies
##

BreakDeps := function (dag)
  local subset, split, graph, i;
  subset := List([1..Length(dag)], i -> rec(
    node_index := i,
    pred := Filtered(List(dag[i].pred, e->Position(dag, e)), e->e <> false),
    succ := Filtered(List(dag[i].succ, e->Position(dag, e)), e->e <> false)));
  split := Filtered(subset, r -> let( max := Maximum0(List(r.succ, i->Length(subset[i].succ))),
             Length(r.pred)=0 and Length(r.succ) > 0 and max > 0 and Length(r.succ) > 2*max));
  split := List(split, e->e.node_index);
  if Length(split) > 0 then
    graph := [[],[]];
    for i in [1..Length(dag)] do
      if i in split then
        Add(graph[1], dag[i]);
      else
        Add(graph[2], dag[i]);
      fi;
    od;
    return rec(graph := graph, modified := true);
  else
    return rec(graph := dag, modified := false);
  fi;
end;

###############################################################################
##
#F  PartitionDag ( <dag> )  . . . . . . . . . . . . returns a partitioned dag
##
PartitionDag := function(dag)
  local hasBlackNode, bBlackNodes, i, redlist, bluelist, serial;

  serial := BreakDeps(dag);
  if serial.modified then
    return List(serial.graph, x->PartitionDag(x));
  else
    serial := SplitDag(dag);
    if serial.modified then
      return List(serial.graph, x->PartitionDag(x));
    else

    # Parallel Split
    if Length(dag)>=2 then

      InitColor(dag);
      bBlackNodes := false;
      while (Length(Filtered(dag, x -> x.color = BLACK))>0) do
        bBlackNodes := true;
        ColorPhase(dag, RED);
        ColorPhase(dag, BLUE);
      od;

      if bBlackNodes then
        redlist := Filtered(dag, x->x.color=RED);
        bluelist := Filtered(dag, x->x.color=BLUE);
        return [PartitionDag(redlist), PartitionDag(bluelist)];
      fi;
    fi;
    return dag;
    fi;
  fi;
end;

##########################################################################
##
#F  ScheduleDag( <partitioned_dag> ) . . . . schedules a paritioned dag
##
ScheduleDag := function(pdag)
  local f, newlist, bigorder;

  newlist := [];
  bigorder := 0;

  f := function(pdag)
    local i,j,k,m;
    if Length(pdag)>0 then
      if IsRec(pdag[1]) then          # use bfs for the order of the nodes
        for j in pdag do j.taken := false; j.order:=0; od;
        while Length(Filtered(pdag, x->not x.taken))>0 do
          k := Filtered(pdag, x -> not x.taken);
          # want all nodes not taken, with preds either outside pdag or taken
          for j in k do
            m := true;
            for i in j.pred do
              if (i in pdag) and (not i.taken) then
                m := false;
              fi;
            od;
            if m then
               j.taken := true;
               Add(newlist, j.cmd);
               bigorder := bigorder + 1;
               j.order := bigorder;
            fi;
          od;
        od;
      elif IsList(pdag[1]) then
        for i in [1..Length(pdag)] do
           f(pdag[i]);
        od;
      fi;
    fi;
  end;

  f(pdag);
  return newlist;

end;

###############################################################################
##
#F  RandomScheduleDag( <dag> ) . . . . . . schedules a dag in a random manner
##
RandomScheduleDag := function(dag)
  local f, newlist, bigorder, list_of_recs;

  newlist := [];
  bigorder := 0;

  f := function(dag)
    local i,j,k,m;
    if Length(dag)>0 then
      if IsRec(dag[1]) then
        for j in dag do j.taken := false; j.order:=0; od;
        while Length(Filtered(dag, x->not x.taken))>0 do
          k := Filtered(dag, x -> not x.taken);
          list_of_recs := [];
          # want all nodes not taken, with preds either outside dag or taken
          for j in k do
            m := true;
            for i in j.pred do
              if (not i.taken) then
                m := false;
              fi;
            od;
            if m then
              Add(list_of_recs, j);
            fi;
          od;
          j := Random(list_of_recs);
          j.taken := true;
          Add(newlist, j.cmd);
          bigorder := bigorder + 1;
          j.order := bigorder;
        od;
      fi;
    fi;
  end;

  f(dag);
  return chain(Filtered(newlist, x -> x <> 0));

end;

###############################################################################
##
#F  GetCompleteSchedules( <dag>, <int_limit> ) . . . . . get complete set of
#F      schedules of a given dag, with the number of limit of the number of
#F      schedules
##
GetCompleteSchedules := function(dag, limit)
  local f, newlist, bigorder, li, dag_length, j;

  newlist := [];
  dag_length := Length(dag);

  f := function(dag, li, bigorder)
    local i,j,k,m, sel, newdag, li2;

    if Length(newlist) >= limit then
      return;
    fi;

    if Length(li) = dag_length then
      m := Filtered(li, x -> x<>0);
      if Length(Filtered(newlist, x->x=m))<=0 then
        Add(newlist, m);
      fi;
    elif Length(dag)>0 then
      if IsRec(dag[1]) then          # use bfs for the order of the nodes
        newdag := Copy(dag);
        # want all nodes not taken, with preds either outside dag or taken
        for j in [1..Length(dag)] do
          if not dag[j].taken  then
            m := true;
            for i in dag[j].pred do
              if (not i.taken) then
                m := false;
              fi;
            od;
            if m then
              dag[j].taken := true;
              li2 := Copy(li);
              if dag[j].cmd <> 0 then
                Add(li2, j);
              else
                Add(li2, 0);
              fi;
              dag[j].order := bigorder;
              newdag := Copy(dag);
              f(newdag, li2, bigorder + 1);
              dag[j].order := 0;
              dag[j].taken := false;
            fi;
          fi;
        od;
      fi;
    fi;
  end;

  for j in dag do j.taken := false; j.order:=0; od;
  f(dag, [], 0);

  newlist := List(newlist, x -> chain(List(x, y -> dag[y].cmd)));

  return newlist;
end;

###############################################################################
##
#F  FFTWScheduleAssignments( <chain> ) . . . schedule assignments
#F      using FFTW's scheduler
##
FFTWScheduleAssignments := function(chain_obj)
  local dag, pdag, schedule;

  dag := BuildDag(chain_obj);
  pdag := PartitionDag(dag); #SplitDagByDataType(dag);

  schedule := chain(ScheduleDag(pdag));

  if PRINTDAG then
    PrintTo("SA_dag", PrintDag(dag, false));
    SYS_EXEC("dot -Tps -odag.ps SA_dag");
    SYS_EXEC("rm SA_dag");

    PrintTo("SA_dag", PrintDag(dag, true));
    SYS_EXEC("dot -Tps -oorder.ps SA_dag");
    SYS_EXEC("rm SA_dag");

    PrintTo("SA_pdag", PrintDag(pdag, false));
    SYS_EXEC("dot -Tps -opdag.ps SA_pdag");
    SYS_EXEC("rm SA_pdag");

    PrintTo("SA_oldschedule", chain_obj);
    PrintTo("SA_newschedule", schedule);

    PrintTo("SA_order", PrintDag(pdag, true));
    SYS_EXEC("dot -Tps -oporder.ps SA_order");
    SYS_EXEC("rm SA_order");
  fi;

  return schedule;
end;

###############################################################################
##
#F  RandomScheduleAssignments(<chain>, <int>) . . . generates a random schedule
##
RandomScheduleAssignments := function(chain_obj)
  local dag, schedule;

  dag := BuildDag(chain_obj);

  schedule := RandomScheduleDag(dag);

  if PRINTDAG then
    PrintTo("SA_dag", PrintDag(dag, false));
    SYS_EXEC("dot -Tps -ordag.ps SA_dag");
    SYS_EXEC("rm SA_dag");

    PrintTo("SA_oldschedule", chain_obj);
    PrintTo("SA_newschedule", schedule);

    PrintTo("SA_dag", PrintDag(dag, true));
    SYS_EXEC("dot -Tps -ororder.ps SA_dag");
    SYS_EXEC("rm SA_dag");

  fi;

  return schedule;
end;

###############################################################################
##
#F  CompleteScheduleAssignments (<chain>, <int_limit>) . . .
#F      generates all the schedules (with limit as the number of schedules
#F      generated)
##
CompleteScheduleAssignments := function(chain_obj, number)
  local dag;

  dag := BuildDag(chain_obj);
  return GetCompleteSchedules(dag, number);

end;

###############################################################################
##
#F  EdwardScheduleAssignments( <chain>, <int_reg>) . . . using
#F      Edward's scheduling algorithm, generate a schedule
##
EdwardScheduleAssignments := function(chain_obj, n)
  local dag,i,j,k,l,m,mincost,minindex,readylist, final, registers, getcost, costlist;

  dag := BuildDag(chain_obj);

  # init readylist and taken
  readylist := Filtered(dag, x -> x.pred=[]);
  for i in dag do
    i.taken := false;
  od;

  final := [];
  registers := List([1..n], x->1);

  getcost := function(node)
    local i,j;
    i:=0;
    for j in node.pred do
      if not (i in registers) then
        i:=i+2;
      else
        i:=i-1;
      fi;
    od;
    return i;
  end;

  while (Length(readylist)>0) do
     # for each node in readylist, determine which to pick
     #   among the nodes in the readylist, pick one with
     #   least cost
     #   1) determine cost
     #   2) determine min
     #   3) take min
     #
     #   N.B. maintain registers
     costlist := List(readylist, x->getcost(x));
     minindex := 1;
     mincost := costlist[1];
#     Print(costlist, "\n");
#     Print("reg = ", Length(Filtered(registers, x->x<>1)), "\n");

     for i in [1..Length(costlist)] do
       if mincost > costlist[i] then
         minindex := i;
         mincost := costlist[i];
       fi;
     od;

     k := readylist[i];
#     Print(k.original, "\n");

     # update registers
     if k in registers then
       i := Position(registers, k);
       j := [k];
       for l in [1..i-1] do
         Add(j, registers[l]);
       od;
       for l in [i+1..n] do
         Add(j, registers[l]);
       od;
       registers := j;
     else
       j := [k];
       for l in [1..n-1] do
         Add(j, registers[l]);
       od;
       registers := j;
     fi;

     readylist := Filtered(readylist, x -> not x=k);
     Add(final, k.cmd);
     k.taken := true;

     for i in k.succ do
       if not i.taken then
         l := true;
         for j in i.pred do
           if not j.taken then
             l := false;
           fi;
         od;
         if l then
           Add(readylist, i);
         fi;
       fi;
     od;

  od;

  return chain(Filtered(final, x->x<>0));

end;

###############################################################################
##
#F  ImFFTWScheduleAssignments( <chain> ) . . . using an algorithm
#F      similar to Edwards scheduler, generate a schedule
##
ImFFTWScheduleAssignments := function(chain_obj)
  local dag, first_level, second_level, output_nodes,i,j,k, color, color_no, li;

  dag := BuildDag(chain_obj);

  # label second level nodes
  for i in dag do
    i.level := 3;        # OTHER LEVEL
  od;
  first_level := Filtered(dag, x -> x.pred=[]);
  second_level := [];
  for i in first_level do
    i.level:= 1;
    for j in i.succ do
      j.level:= 2;
      if not (j in second_level) then
        Add(second_level, j);
      fi;
    od;
  od;

  # get outputnodes;
  output_nodes := Filtered(dag, x->x.succ = []);

  # color output nodes
  for i in dag do
    i.color := 0;  # means uninitialized
  od;
  color := function(node, c)
    local i;
    if not node.color = 0 then
      return;
    elif node.level=3 then
      node.color := c;
      for i in node.pred do
        color(i, c);
      od;
      for i in node.succ do
        color(i, c);
      od;
    elif node.level = 2 then
      if node.color = 0 then
        node.color := c;
      fi;
    fi;
  end;

  color_no := 0;
  for i in output_nodes do
    if i.color=0 then
      color_no := color_no + 1;
      color(i, color_no);
    fi;
  od;

  # for each color, schedule graph
  li := [];
  for i in [1..color_no] do
    j := Filtered(dag, x -> x.color=i);
    j := ScheduleDag(j);
    for k in [1..Length(j)] do
      Add(li, j[k]);
    od;
  od;

  # output statements;
  return chain(li);
end;

# tochain := x -> Cond(IsChain(x), Copy(x), chain(Copy(x)));

# bsplit := function(code)
#     local dims;
#     dims := code.dimensions;
#     code := When(IsBound(code.cmd), code.cmd, code);
#     code := FlattenCode(BinSplit(tochain(code)));
#     code := Compile.declareVars(code);
#     code.dimensions := dims;
#     return code;
# end;

# sched := function(code)
#     local dims;
#     dims := code.dimensions;
#     code := When(IsBound(code.cmd), code.cmd, code);
#     code := FFTWScheduleAssignments(tochain(code));
#     code := Compile.declareVars(code);
#     code.dimensions := dims;
#     return code;
# end;

# sched2 := function(code)
#     local dims;
#     dims := code.dimensions;
#     code := When(IsBound(code.cmd), code.cmd, code);
#     code := EdwardScheduleAssignments(tochain(code), 8);
#     code := Compile.declareVars(code);
#     code.dimensions := dims;
#     return code;
# end;

# randsched := function(code)
#     local dims;
#     dims := code.dimensions;
#     code := When(IsBound(code.cmd), code.cmd, code);
#     code := RandomScheduleAssignments(tochain(code));
#     code := Compile.declareVars(code);
#     code.dimensions := dims;
#     return code;
# end;


text:
#
# francy: Interactive Discrete Mathematics in GAP
#

#############################################################################
##
#M  Menu( <title> ) . . . . . . . . . . . . . . . a new menu entry
##
InstallMethod(Menu,
  "a title string",
  true,
  [IsString,
   IsCallback],
  0,
function(title, callback)
  return Objectify(MenuObjectType, rec(
    id       := GenerateID(),
    title    := title,
    callback := callback,
    menus    := rec()
  ));
end);

InstallOtherMethod(Menu,
  "a title string",
  true,
  [IsString],
  0,
function(title)
  return Menu(title, NoopCallback());
end);

#############################################################################
##
#M  Add( <menu>, <menu> ) . . . . . add menu to canvas
##
InstallOtherMethod(Add,
  "a menu, a menu",
  true,
  [IsMenu,
   IsMenu],
  0,
function(menu, object)
  menu!.menus!.(object!.id) := object;
  return menu;
end);

InstallOtherMethod(Add,
  "a canvas, a list of francy objects",
  true,
  [IsMenu,
   IsList],
  0,
function(menu, objects)
  local object;
  for object in objects do
    Add(menu, object);
  od;
  return menu;
end);

#############################################################################
##
#M  Remove( <menu>, <menu> ) . . . . . remove menu from menu
##
InstallOtherMethod(Remove,
  "a menu, a menu",
  true,
  [IsMenu,
   IsMenu],
  0,
function(menu, object)
  Unbind(menu!.menus!.(object!.id));
  return menu;
end);

InstallOtherMethod(Remove,
  "a menu, a list of francy objects",
  true,
  [IsMenu,
   IsList],
  0,
function(menu, objects)
  local object;
  for object in objects do
    Remove(menu, object);
  od;
  return menu;
end);

ext: gi

max_stars_repo_name: LaGuer/francy

max_issues_repo_name: LaGuer/francy

max_forks_repo_name: LaGuer/francy


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

#F IsNonTerminal ( <obj> )
#F   returns true, if <obj> is a non-terminal
#F
IsNonTerminal := N -> IsRec(N) and IsBound(N.isNonTerminal) and N.isNonTerminal = true;

Class(NonTerminalOps, SPLOps, rec(
));

Declare(NonTerminal);

# ==========================================================================
# NonTerminal
# ==========================================================================
Class(NonTerminal, BaseMat, rec(
    #F create getter methods
    #F    ex: NonTerminal.WithParams("param1_getter_name", "param2_getter_name")
    WithParams := arg >> CopyFields(arg[1],
	RecList( ConcatList( [2..Length(arg)],
	    i -> [arg[i], DetachFunc(Subst(self >> self.params[$(i-1)]))] ))),

    isNonTerminal := true, # for nonterm.g and IsNonTerminal() function
    _short_print := false,

    #--------- Transformation rules support --------------------------------
    from_rChildren := (self, rch) >> let(
        t := ApplyFunc(ObjId(self), DropLast(rch, 1)),
        When(Last(rch), t.transpose(), t)
    ),

    rChildren := self >>
        Concatenation(self.params, [self.transposed]),

    rSetChild := meth(self, n, newChild)
        if n <= 0 or n > Length(self.params) + 1
            then Error("<n> must be in [1..", Length(self.params)+1, "]"); fi;
        if n <= Length(self.params) then
            self.params[n] := newChild;
        else
            self.transposed := newChild;
        fi;
        # self.canonizeParams(); ??
        self.dimensions := self.dims();
    end,

    #-----------------------------------------------------------------------
    canonizeParams := meth(self)
        local A, nump;
        nump := Length(self.params);
    if nump = 0 then return; fi;
        for A in self.abbrevs do
            if NumArgs(A) = -1 or NumArgs(A) = nump then
                self.params := ApplyFunc(A, self.params);
                return;
            fi;
        od;
        Error("Nonterminal ", self.name, " can take ", List(self.abbrevs, NumArgs),
              "  arguments, but not ", nump);
    end,

    new := meth(arg)
        local result, self, params;
        self := arg[1];
        params := arg{[2..Length(arg)]};
        result := SPL(WithBases(self, rec(params := params, transposed := false )));
        result.canonizeParams();
        result.dimensions := result.dims();
        return result;
    end,
    #-----------------------------------------------------------------------
    __call__ := ~.new,
    #-----------------------------------------------------------------------
    isTerminal := False,
    #-----------------------------------------------------------------------
    isPermutation := False,
    #-----------------------------------------------------------------------
    transpose := self >>
        Inherit(self, rec(transposed := not self.transposed,
                  dimensions := Reversed(self.dimensions))),
    setTransposed := (self, v) >> When( self.transposed <> v,
        self.transpose(), self),
    #-----------------------------------------------------------------------
    # .conjTranspose() is needed for RC(T).transpose() == RC(T.conjTranspose())
    # The inert form will only work is non-terminal is not used inside RC,
    # or used inside RC, but not transposed
    conjTranspose := self >> InertConjTranspose(self),
    isInertConjTranspose := self >> self.conjTranspose = NonTerminal.conjTranspose,

    #-----------------------------------------------------------------------
    print := (self, i, is) >>
	Cond(
            not IsBound(self.params), Print(self.name),
	    Print(self._print(self.params, i, is),
		  When(self.transposed, ".transpose()"))),
    #-----------------------------------------------------------------------
    toAMat := self >> self.terminate().toAMat(),
    #-----------------------------------------------------------------------
    export := arg >> Error("Can't export a non-terminal"),
    #-----------------------------------------------------------------------
    arithmeticCost := arg >> Error("Can't compute arithmetic cost of a non-terminal"),

    transposeSymmetric := True,
    isSymmetric := False,

    area := self >> Product(self.dims())
));


text:
#
# walrus: Computational Methods for Finitely Generated Monoids and Groups
#
## Locations

# NewLocation : IsPregroupRelator -> Int -> IsPregroupLocation
InstallGlobalFunction(NewLocation,
function(R,i)
    return Objectify(PregroupLocationType, [R, i, [], R[i-1], R[i]]);
end);

InstallMethod(Relator, "for a location",
              [ IsPregroupLocationRep ],
              l -> l![1]);

InstallMethod(Position, "for a location",
              [ IsPregroupLocationRep ],
              l -> l![2]);

InstallMethod(InLetter, "for a location",
              [ IsPregroupLocationRep ],
              l -> l![4] );

InstallMethod(OutLetter, "for a location",
              [ IsPregroupLocationRep ],
              l -> l![5] );

#X Get rid of this
InstallMethod(__ID, "for a location",
              [ IsPregroupLocationRep ],
              l -> l![6]);

# Return list of Places that have this location
InstallMethod(Places, "for a location",
              [ IsPregroupLocationRep ],
function(l)
    return l![3];
end);

InstallMethod(NextLocation, "for a location",
              [ IsPregroupLocationRep ],
function(l)
    local loc;

    # All locations on the relator that l is on
    loc := Locations(l![1]);
    if l![2] < Length(loc) then
        return loc[l![2] + 1];
    else
        return loc[1];
    fi;
end);

InstallMethod(PrevLocation, "for a location",
              [ IsPregroupLocationRep ],
function(l)
    local loc;

    # All locations on the relator that l is on
    loc := Locations(l![1]);
    if l![2] > 1 then
        return loc[l![2] - 1];
    else
        return loc[Length(loc)];
    fi;
end);

# Two locations are the same if they are on the same
# relator at the same index
InstallMethod(\=, "for a location and a location",
              [ IsPregroupLocationRep, IsPregroupLocationRep],
function(l,r)
    return (l![1] = r![1])
           and (l![2] = r![2]);
end);

InstallMethod(PregroupPresentationOf, "for a location",
              [ IsPregroupLocationRep ],
              l -> PregroupPresentationOf(l![1]));

InstallMethod(ViewString, "for a location",
              [ IsPregroupLocationRep ],
function(l)
    return STRINGIFY(ViewString(l![1]), "("
                     , l![2], ","
                     , ViewString(InLetter(l)), ","
                     , ViewString(OutLetter(l)), ")"
                    );
end);

max_stars_repo_name: RussWoodroofe/walrus

max_issues_repo_name: RussWoodroofe/walrus

max_forks_repo_name: RussWoodroofe/walrus


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

#s:=SPLRuleTree(ExpandSPL(PolyBDFT3(16))[1]);

#s:=SPLRuleTree(ExpandSPL(BSkewDFT3(16))[1]);
# Even1 -- Odd0
s := K(16,4) *
     DirectSum(
	 BSkewDFT3(4, 1/16),
	 BSkewDFT3(4, 3/16),
	 BSkewDFT3(4, 5/16),
	 BSkewDFT3(4, 7/16)
     ) *
     L(16,8) *
     Tensor(I(2), PolyBDFT3(8, 1/4));
     L(16,2);

left := K(8,2) *
        DirectSum(
	  SkewDTT(DCT3(2), 1/8),
	  SkewDTT(DCT3(2), 3/8),
	  SkewDTT(DCT3(2), 5/8),
	  SkewDTT(DCT3(2), 7/8)
        );

 right:=MatSPL(left)^-1 * MatSPL(DCT3(8));

spiral> ps(MatSPL(L(8,2))*right*MatSPL(L(8,4)));

left := K(8,4) *
        DirectSum(
	  SkewDTT(DCT3(4), 1/4),
	  SkewDTT(DCT3(4), 3/4)
        );

 right:=MatSPL(left)^-1 * MatSPL(DCT3(8));

spiral> ps(MatSPL(L(8,4))*right*MatSPL(L(8,2)));

# X * BRDFT * IJ = BRDFT,      X = BRDFT * IJ * BRDFT^-1
# Y * BRDFT * IJ = BRDFT * J,  Y = BRDFT*J*IJ*BRDFT^-1
#
c3polyalgo := function(N,m,a)
    local k, left, right, zz, mid;
    k := N/2/m;
    zz := RulesDFTSymmetry(L(N, m) * GathExtend(N,Odd0));
    left := K(N/2,m) * DirectSum(List([0..k-1], i -> SkewDTT(DCT3(m), fr(k,i,a).ev())));
    right := DirectSum(
	SkewDTT(DCT3(k), a),
	List([1..(m-2)/2], i -> BRDFT3(2*k,a/2)*DirectSum(I(1), -J(2*k-1))),
	PolyDTT(SkewDTT(DCT4(k),a))
    ) *
    Compose(Drop(zz.children(), 1));

    mid := MatSPL(left)^-1 * MatSPL(SkewDTT(DCT3(N/2),a)) * MatSPL(right)^-1;
    #PrintMat(mid);
    VisualizeMat(mid, ", ");
    return mid;
end;

#N := 32; m := 4; k := N/2/m;

#[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]
mid := c3polyalgo(32,8);
p:=Prm(FList([0,2,3,6,7,10,11,14, 1,4,5,8,9,12,13, 15]).setRange(15));

ps(mid * MatSPL(p)^-1);


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

Class(TerminateRaderMid, RuleSet);
RewriteRules(TerminateRaderMid, rec(
    terminateRaderMid := Rule(@(1, VecRaderMid), e-> e.term())
));

Class(StretchRaderMid, RuleSet);
RewriteRules(StretchRaderMid, rec(

    stretchCT_A := ARule(Compose,
        [ [ DelayedDirectSum, @(2, VScat_sv, e->Cols(e)=getV(e)), [ @(3, Compose), ..., @(4, IxVGath_pc) ] ],
          @(5, VecRaderMid),
          [ DelayedDirectSum, @(7, VGath_sv, e->Rows(e)=getV(e)), [ @(8, Compose),
                @(9, IxVScat_pc, s-> let(g:=@(4).val, s.k=g.k and s.N=g.N and s.ofs=g.ofs and s.v=g.v )), ... ] ]
        ],
        e -> [ DelayedDirectSum(@(2).val, Compose(DropLast(@(3).val.children(), 1))),
               let(g:=@(4).val , @(5).val.stretchCT(g.N)),
               DelayedDirectSum(@(7).val, Compose(Drop(@(8).val.children(), 1))) ]),

    stretchCT := ARule(Compose,
        [ [ DelayedDirectSum, [ @(1, Compose), ..., @(2, VScat_sv, e->Cols(e)=getV(e)) ], [ @(3, Compose), ..., @(4, IxVGath_pc) ] ],
          @(5, VecRaderMid),
          [ DelayedDirectSum, [ @(6, Compose), @(7, VGath_sv, e->Rows(e)=getV(e)), ... ], [ @(8, Compose),
                @(9, IxVScat_pc, s-> let(g:=@(4).val, s.k=g.k and s.N=g.N and s.ofs=g.ofs and s.v=g.v )), ... ] ]
        ],
        e -> [ DelayedDirectSum(@(1).val, Compose(DropLast(@(3).val.children(), 1))),
               let(g:=@(4).val , @(5).val.stretchCT(g.N)),
               DelayedDirectSum(@(6).val, Compose(Drop(@(8).val.children(), 1))) ]),

#--
    stretchCT2_A := ARule(Compose,
        [ [ DelayedDirectSum, @(2, VScat_sv, e->Cols(e)=getV(e)), [ @(3, Compose), ..., @(4, IxVGath_pc), [VScat_sv, fId] ] ],
          @(5, VecRaderMid),
          [ DelayedDirectSum, @(7, VGath_sv, e->Rows(e)=getV(e)), [ @(8, Compose),
                [VGath_sv, fId], @(9, IxVScat_pc, s-> let(g:=@(4).val, s.k=g.k and s.N=g.N and s.ofs=g.ofs and s.v=g.v )), ... ] ]
        ],
        e -> [ DelayedDirectSum(@(2).val, Compose(DropLast(@(3).val.children(), 2))),
               let(g:=@(4).val , @(5).val.stretchCT(g.N)),
               DelayedDirectSum(@(7).val, Compose(Drop(@(8).val.children(), 2))) ]),

    stretchCT2 := ARule(Compose,
        [ [ DelayedDirectSum, [ @(1, Compose), ..., @(2, VScat_sv, e->Cols(e)=getV(e)) ], [ @(3, Compose), ..., @(4, IxVGath_pc), [VScat_sv, fId] ] ],
          @(5, VecRaderMid),
          [ DelayedDirectSum, [ @(6, Compose), @(7, VGath_sv, e->Rows(e)=getV(e)), ... ], [ @(8, Compose),
                [VGath_sv, fId], @(9, IxVScat_pc, s-> let(g:=@(4).val, s.k=g.k and s.N=g.N and s.ofs=g.ofs and s.v=g.v )), ... ] ]
        ],
        e -> [ DelayedDirectSum(@(1).val, Compose(DropLast(@(3).val.children(), 2))),
               let(g:=@(4).val , @(5).val.stretchCT(g.N)),
               DelayedDirectSum(@(6).val, Compose(Drop(@(8).val.children(), 2))) ]),

#--
    stretchPFA_A := ARule(Compose,
        [ [ DelayedDirectSum, @(2, VScat_sv, e->Cols(e)=getV(e)), [ @(3, Compose), ..., @(4, VStretchGath), [VScat_sv, fId] ] ],
          @(5, VecRaderMid),
          [ DelayedDirectSum, @(7, VGath_sv, e->Rows(e)=getV(e)), [ @(8, Compose),
                [VGath_sv, fId], @(9, VStretchScat, s-> let(g:=@(4).val, s.part=g.part and s.func.range()=g.func.range() and s.v=g.v )),  ... ] ]
        ],
        e -> [ DelayedDirectSum(@(2).val, Compose(DropLast(@(3).val.children(), 2))),
               let(g:=@(4).val , @(5).val.stretchPFA(g.part, g.func.range(), @(4).val.func)),
               DelayedDirectSum(@(7).val, Compose(Drop(@(8).val.children(), 2))) ]),

    stretchPFA := ARule(Compose,
        [ [ DelayedDirectSum, [ @(1, Compose), ..., @(2, VScat_sv, e->Cols(e)=getV(e)) ], [ @(3, Compose), ..., @(4, VStretchGath), [VScat_sv, fId] ] ],
          @(5, VecRaderMid),
          [ DelayedDirectSum, [ @(6, Compose), @(7, VGath_sv, e->Rows(e)=getV(e)), ... ], [ @(8, Compose),
                [VGath_sv, fId], @(9, VStretchScat, s-> let(g:=@(4).val, s.part=g.part and s.func.range()=g.func.range() and s.v=g.v )),  ... ] ]
        ],
        e -> [ DelayedDirectSum(@(1).val, Compose(DropLast(@(3).val.children(), 2))),
               let(g:=@(4).val , @(5).val.stretchPFA(g.part, g.func.range(), @(4).val.func)),
               DelayedDirectSum(@(6).val, Compose(Drop(@(8).val.children(), 2))) ]),
#--
    stretchPFA2_A := ARule(Compose,
        [ [ DelayedDirectSum, [ @(1, Compose), ..., @(2, VScat_sv, e->Cols(e)=getV(e)) ], [ @(3, Compose), ..., @(4, VStretchGath), ] ],
          @(5, VecRaderMid),
          [ DelayedDirectSum, [ @(6, Compose), @(7, VGath_sv, e->Rows(e)=getV(e)), ... ], [ @(8, Compose),
                @(9, VStretchScat, s-> let(g:=@(4).val, s.part=g.part and s.func.range()=g.func.range() and s.v=g.v )),  ... ] ]
        ],
        e -> [ DelayedDirectSum(@(1).val, Compose(DropLast(@(3).val.children(), 1))),
               let(g:=@(4).val , @(5).val.stretchPFA(g.part, g.func.range(), @(4).val.func)),
               DelayedDirectSum(@(6).val, Compose(Drop(@(8).val.children(), 1))) ]),

    stretchPFA2 := ARule(Compose,
        [ [ DelayedDirectSum, @(2, VScat_sv, e->Cols(e)=getV(e)), [ @(3, Compose), ..., @(4, VStretchGath), ] ],
          @(5, VecRaderMid),
          [ DelayedDirectSum, @(7, VGath_sv, e->Rows(e)=getV(e)), [ @(8, Compose),
                @(9, VStretchScat, s-> let(g:=@(4).val, s.part=g.part and s.func.range()=g.func.range() and s.v=g.v )),  ... ] ]
        ],
        e -> [ DelayedDirectSum(@(2).val, Compose(DropLast(@(3).val.children(), 1))),
               let(g:=@(4).val , @(5).val.stretchPFA(g.part, g.func.range(), @(4).val.func)),
               DelayedDirectSum(@(7).val, Compose(Drop(@(8).val.children(), 1))) ]),
#--
    stretchPFA3_A := ARule(Compose,
        [ [ DelayedDirectSum, [ @(1, Compose), ..., @(2, VScat_sv, e->Cols(e)=getV(e)) ], [ @(3, Compose), ..., @(4, DelayedPrm) ] ],
          @(5, VecRaderMid),
          [ DelayedDirectSum, [ @(6, Compose), @(7, VGath_sv, e->Rows(e)=getV(e)), ... ], [ @(8, Compose),
                @(9, DelayedPrm, s-> let(g:=@(4).val, 1=1)),  ... ] ]
        ],
        e -> [ DelayedDirectSum(@(1).val, Compose(DropLast(@(3).val.children(), 1))),
               let(g:=@(4).val , @(5).val.stretchPFA(1, Cols(g), @(4).val.func)),
               DelayedDirectSum(@(6).val, Compose(Drop(@(8).val.children(), 1))) ]),

    stretchPFA3 := ARule(Compose,
        [ [ DelayedDirectSum, @(2, VScat_sv, e->Cols(e)=getV(e)), [ @(3, Compose), ..., @(4, DelayedPrm) ] ],
          @(5, VecRaderMid),
          [ DelayedDirectSum, @(7, VGath_sv, e->Rows(e)=getV(e)), [ @(8, Compose),
                @(9, DelayedPrm, s-> let(g:=@(4).val, 1=1)),  ... ] ]
        ],
        e -> [ DelayedDirectSum(@(2).val, Compose(DropLast(@(3).val.children(), 1))),
               let(g:=@(4).val , @(5).val.stretchPFA(1, Cols(g), @(4).val.func)),
               DelayedDirectSum(@(7).val, Compose(Drop(@(8).val.children(), 1))) ])
));


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

# determine ldepth by object type
_ldepth := o -> let(
    x := When(ObjId(o)=tcast, o.args[2], o),

    # cond functions like a switch statement.
    Cond(
        IsBound(x.ldepth),
            x.ldepth,
        not IsVar(x),
            2^20,
        let(
            f := DefFrontier(x, 5),
            l := Length(f),
            last := Last(f),

            Cond(
                l = 1,
                    0,
                Length(last)=1 and IsLoopIndex(last[1]),
                    last[1].ldepth,
                2^20
            )
        )
    )
);

# added a high value for constants so that they sort at one end
# when the variables are sorted by ldepth
_ldepth2 := (o, params_first) -> let(
    x := When(ObjId(o)=tcast, o.args[2], o),

    # cond functions like a switch statement.
    Cond(
        IsBound(x.ldepth),
            x.ldepth,
        IsValue(x),
            (2^20)+1,
        ObjId(x) in [param, in_param], 100,
#            (2^20) + When(params_first, -2-(InternalHash(x.id) mod 1024),
#                                        +2+(InternalHash(x.id) mod 1024)),

        not IsVar(x),
            2^20,
        let(
            f := DefFrontier(x, 5),
            l := Length(f),
            last := Last(f),

            Cond(
                l = 1,
                    0,
                Length(last)=1 and IsLoopIndex(last[1]),
                    last[1].ldepth,
                2^21
            )
        )
    )
);

Class(BinSplit, rec(

    canSplit := (self, exp) >> ((IsExpCommand(exp) and not IsBound(exp.isExpComposite)) or
        (IsBound(exp.isExpComposite) and exp.isExpComposite
         and ForAny(exp.rChildren(), e->IsBound(e.isExpComposite))))
         and not (exp _is self._splitDisabled),

    expSplit := meth(self, exp)
        local c, v, composites, res, ch, args, i, prev, newexp;
        res := [];
        if self.canSplit(exp) then
            ch := exp.rChildren();
            composites := SplitBy([1..Length(ch)], i->IsBound(ch[i].isExpComposite))[1];
            for c in composites do
                v := var.fresh_t("a", ch[c].t);
                Append(res, self.assignSplit(assign(v, ch[c])));
                exp.rSetChild(c, v);
            od;
        fi;

        # check for non-binary add/mult
	# NOTE: use a list of assoc expr? use Exp property!!
        if (ObjId(exp) in [add, mul, xor, bin_xor, bin_and]) and Length(exp.args) > 2 then
            args := ShallowCopy(exp.args);
            for v in args do v.ldepth := _ldepth2(v, ObjId(exp)=mul); od;
            Sort(args, (a,b) -> a.ldepth < b.ldepth);

            prev := args[1];
            for i in [2..Length(args)-1] do
	        newexp := ObjId(exp)(prev, args[i]);
                v := var.fresh_t("b", newexp.t);
                Add(res, assign(v, newexp));
                prev := v;
            od;
            exp.args := [prev, Last(args)];
        fi;
	return [res, exp];
    end,

    assignSplit := meth(self, cmd)
        local loc, exp, cmds_loc, cmds_exp;
        [loc, exp] := [cmd.loc, cmd.exp];
	[cmds_exp, exp] := self.expSplit(exp);
	[cmds_loc, loc] := self.expSplit(loc);

        if ObjId(exp) in [div,mul,add,xor, stickyNeg] then
            # sets kdepth correctly for stuff like a+(b+c)
            loc.ldepth:=Minimum(List(exp.args,x->_ldepth2(x, false)));
            # Sets ldepth correctly for stuff like a2= 2*i3 => a2.ldepth=_ldepth(i3)
            if (loc.ldepth>=2^18) then
                loc.ldepth:=Minimum0(Filtered(List(exp.args,x->_ldepth2(x, false)), x->x<2^18));
            fi;
        fi;
        return cmds_exp :: cmds_loc :: [cmd];
    end,

    expCommandSplit := (self, cmd) >> self.expSplit(cmd),
#    expCommandSplit := meth(self, cmd)
#        local a, args, res, cmds;
#	res := [];
#	for a in cmd.args do
#	    [cmds, a] := self.expSplit(a);
#	    Append(res, cmds);
#	od;
#        return res :: [cmd];
#    end,

    # BinSplit(code[, opts])

    __call__ := meth(arg)
        local self, code, opts;
        [self, code, opts] := [arg[1], arg[2], Cond(Length(arg)>2, arg[3], rec())];
        if IsBound(opts.binsplitOnIndices) and not opts.binsplitOnIndices then
            self := CopyFields(self, rec( _splitDisabled := opts.simpIndicesInside ));
        else
            self := CopyFields(self, rec( _splitDisabled := [] ));
        fi;
        return self.apply(code);
    end,

    apply := meth(self, code)
        local chains, blk, i, cmd, aux, exp, split;
        chains := Collect(code, chain);
	for blk in chains do
	    for i in [1..Length(blk.cmds)] do
	        cmd := blk.cmds[i];
	        if IsAssign(cmd) then
                    # cmd could have parts aliased with different command
                    # since we do inplace modifications of cmd inside assignSplit, we copy it
                    split := self.assignSplit(Copy(cmd));
                    blk.cmds[i] := When(Length(split)=1, split[1], chain(split));
		elif IsExpCommand(cmd) then
                    split := self.expCommandSplit(Copy(cmd));
                    blk.cmds[i] := When(Length(split)=1, split[1], chain(split));
                fi;
            od;
        od;
        return FlattenCode(code);
    end,
));


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

_svctSizes := (N, p, v) -> Filtered([1..N], i->ForAll(Factors(i), j->j<=p) and IsInt(i/v^2));
_defaultSizes := (func, default) -> ((arg) -> When(Length(arg) = 1, func(When(IsList(arg[1]), arg[1], [arg[1]])), func(default)));

_parse := function(arglist)
    local isa, n, sizes, optrec;

    isa := "f64re";
    optrec := rec(interleavedComplex := true);

    if Length(arglist) = 0 then
        n := 1;
    elif Length(arglist) = 1 then
        if IsInt(arglist[1]) then
            n := arglist[1];
        elif IsList(arglist[1]) then
            n := arglist[1];
        else
            n := 0;
            isa := arglist[1];
        fi;
    elif Length(arglist) = 2 then
        n := arglist[1];
        isa := arglist[2];
    else
        n := arglist[1];
        isa := arglist[2];
        optrec := arglist[3];
    fi;

    sizes := When(IsList(n), n, List([1..n], i->2^i));

    return [sizes, isa, optrec];
end;

doDft := function(arg)
    local sizes, opts, dpr, isa, trafos, t, optrec, dpbench, defaults;

    [sizes, isa, optrec] := _parse(arg);

    defaults := rec(
        breakdownRules := rec(
            DFT := [ DFT_Base, DFT_PD, DFT_CT, DFT_GoodThomas, DFT_Rader ],
            WHT := [ WHT_Base, WHT_BinSplit ],
            MDDFT := [ MDDFT_Base, MDDFT_RowCol ],
            TRC := [ TRC_tag ]
        ),
        unparser := CUnparserProg
    );
    opts := CopyFields(When(optrec.generateComplexCode, CplxSpiralDefaults, SpiralDefaults), defaults, optrec);
    opts := InitDataType(opts, isa);

    if opts.generateComplexCode then
        t := i -> RC(DFT(i));
    else
        t := When(IsBound(opts.tags), i -> ComplexT(DFT(i), opts).withTags(opts.tags), i ->ComplexT(DFT(i), opts));
    fi;
    opts.benchTransforms := List(sizes, k -> t(k));

    dpr  := rec(verbosity := 0, timeBaseCases:=true);
    dpbench := spiral.libgen.DPBench(rec((Concat(isa, When(opts.interleavedComplex, "_ic", "_sc"))) := opts), dpr);
    if IsBound(optrec.verify) then dpbench.matrixVerify := optrec.verify; fi;

    return dpbench;
end;

# doDft auto-vectorization interface
#b := doDft(10, "f32re", _autovect(IntelC, rec(interleavedComplex := true, language := "c.icl.opt.core2")));
#b := doDft(16, "f32c", _autovect(IntelC, rec(interleavedComplex := true, language := "c.icl.opt.core2", generateComplexCode:=true)));

_autovect := (cc, opts) -> CopyFields(opts, rec(
    useRestrict := true,
    looppragma := cc.looppragma,
    restrict := cc.restrict,
    postalign := cc.postalign,
    arrayBufModifier := Concat("static ", cc.alignmentSpecifier()),
    arrayDataModifier := Concat("static ", cc.alignmentSpecifier()),
    compileStrategy:=NoCSE
));

benchScalar := () -> rec(
    float := rec(
        dft_ic := rec(
            medium := _defaultSizes(s->doDft(s, "f32re", rec(interleavedComplex := true)), 13)
        ),
        dft_sc := rec(
            medium := _defaultSizes(s->doDft(s, "f32re", rec(interleavedComplex := false)), 13)
        )
    ),
    double := rec(
        dft_ic := rec(
            medium := _defaultSizes(s->doDft(s, "f64re", rec(interleavedComplex := true)), 13)
        ),
        dft_sc := rec(
            medium := _defaultSizes(s->doDft(s, "f64re", rec(interleavedComplex := false)), 13)
        )
    )
);


text:
#############################################################################
##
#W    pcpgrp.gi                                                  Bettina Eick
##
##    This file contains the header functions to handle the 3- and
##    4-dimensional almost crystallographic pcp groups.
##

#############################################################################
##
#F IsAlmostCrystallographic( <G> )
##
InstallMethod( IsAlmostCrystallographic, "for pcp groups", true,
    [IsPcpGroup ], 0,
function( G )
    if HasAlmostCrystallographicInfo( G ) then return true; fi;
    return (IsInt(Index(G,FittingSubgroup(G))) and
            Size(NormalTorsionSubgroup(G)) = 1);
end );

#############################################################################
##
#F IsAlmostBieberbachGroup( <G> )
##
InstallMethod( IsAlmostBieberbachGroup, "for pcp groups", true,
    [IsPcpGroup ], 0,
function( G )
    return IsInt(Index(G,FittingSubgroup(G))) and IsTorsionFree(G);
end );

#############################################################################
##
#F IsolatorSubgroup( G, N )
##
InstallGlobalFunction( IsolatorSubgroup, function( G, N )
    local nat, F, T;
    if not IsNormal( G, N ) then return fail; fi;
    nat := NaturalHomomorphismByNormalSubgroup( G, N );
    F := Image( nat );
    T := TorsionSubgroup( F );
    if IsBool( T ) then return fail; fi;
    return PreImage( nat, T );
end);

#############################################################################
##
#F AlmostCrystallographicPcpDim3( type, param )
##
InstallGlobalFunction( AlmostCrystallographicPcpDim3,
function( type, param )
    local g, G, info;

    # type is integer or string
    if IsString( type ) then
        g := Int( type );
    elif IsInt( type ) then
        g := type;
    else
        g := false;
    fi;
    if not g in [1..17] then
        Error("type does not define a valid ac-group type");
    fi;

    # check parameters
    if IsBool( param ) then
        param := List( [1..ACDim3Param[g]], x -> Random(Integers) );
        if param[1] = 0 then param[1] := Random(Integers)^2+1; fi;
    elif IsInt( param ) then
        param := List( [1..ACDim3Param[g]], x -> param );
    elif Length( param ) <> ACDim3Param[g] then
        Error("parameter should be a list of length ", ACDim3Param[g] );
    fi;

    # get group
    if Length( param ) = 4 then
        G := ACPcpDim3Funcs[g](param[1], param[2], param[3], param[4]);
    elif Length(param) = 2 then
        G := ACPcpDim3Funcs[g](param[1], param[2] );
    elif Length(param) = 1 then
        G := ACPcpDim3Funcs[g](param[1]);
    fi;

    # set some information
    info := rec( dim := 3, type := g, param := param );
    SetAlmostCrystallographicInfo( G, info );
    SetIsAlmostCrystallographic( G, true );
    SetSize( G, infinity );
    return G;
end );

#############################################################################
##
#F AlmostCrystallographicPcpDim4( type, param )
##
InstallGlobalFunction( AlmostCrystallographicPcpDim4,
function( type, param )
    local g, G, info;

    # type is integer or string
    if IsString( type ) then
        g := Position( ACDim4Types, type );
    elif IsInt( type ) then
        g := type;
    else
        g := false;
    fi;
    if not g in [1..95] then
        Error("type does not define a valid ac-group type");
    fi;

    # check parameters
    if IsBool( param ) then
        param := List( [1..ACDim4Param[g]], x -> Random(Integers) );
        if param[1] = 0 then param[1] := Random(Integers)^2+1; fi;
    elif IsInt( param ) then
        param := List( [1..ACDim4Param[g]], x -> param );
    elif Length( param ) <> ACDim4Param[g] then
        Error("parameter should be a list of length ", ACDim4Param[g] );
    fi;

    # get group
    if Length(param) = 7 then
        G := ACPcpDim4Funcs[g]( param[1], param[2], param[3], param[4],
                                param[5], param[6], param[7]);
    elif Length(param) = 6 then
        G := ACPcpDim4Funcs[g]( param[1], param[2], param[3], param[4],
                                param[5], param[6]);
    elif Length(param) = 5 then
        G := ACPcpDim4Funcs[g]( param[1], param[2], param[3], param[4],
                                param[5]);
    elif Length(param) = 4 then
        G := ACPcpDim4Funcs[g]( param[1], param[2], param[3], param[4]);
    elif Length(param) = 3 then
        G := ACPcpDim4Funcs[g]( param[1], param[2], param[3]);
    fi;

    # set some information
    info := rec( dim := 4, type := g, param := param );
    SetAlmostCrystallographicInfo( G, info );
    SetIsAlmostCrystallographic( G, true );
    SetSize( G, infinity );
    return G;
end );

#############################################################################
##
#F AlmostCrystallographicPcpGroup( dim, type, param )
##
InstallGlobalFunction( AlmostCrystallographicPcpGroup,
function( dim, type, param )
    if dim = 3 then
        return AlmostCrystallographicPcpDim3( type, param );
    elif dim = 4 then
        return AlmostCrystallographicPcpDim4( type, param );
    else
        Error("dimension must be 3 or 4");
    fi;
end );

#############################################################################
##
#A  FittingSubgroup( < G > )
##
InstallMethod( FittingSubgroup, "for ac pcp groups", true,
      [IsPcpGroup and HasAlmostCrystallographicInfo], 0,
function( G )
    local pcp, rel, sub, F;
    pcp := Pcp( G );
    rel := RelativeOrdersOfPcp( pcp );
    sub := Filtered( [1..Length(pcp)], x -> rel[x] = 0 );
    F := Subgroup( G, pcp{sub} );
    SetIsNilpotentGroup( F, true );
    return F;
end );

#############################################################################
##
#F  NaturalHomomorphismOnHolonomyGroup( < G > )
#F  HolonomyGroup( < G > )
##
InstallMethod( NaturalHomomorphismOnHolonomyGroup, "for ac pcp groups", true,
    [IsPcpGroup and IsAlmostCrystallographic], 0,
function( G )
    return NaturalHomomorphismByNormalSubgroup(G,FittingSubgroup(G));
end );

InstallMethod( HolonomyGroup, "for ac pcp groups", true,
    [IsPcpGroup and IsAlmostCrystallographic], 0,
function( G )
    return Image( NaturalHomomorphismOnHolonomyGroup(G) );
end );

#############################################################################
##
#F IsomorphismPcpGroup( <G> )
##
InstallMethod( IsomorphismPcpGroup, "for ac groups", true,
    [IsMatrixGroup and HasAlmostCrystallographicInfo], 0,
function( G )
    local info, H, gensH, gensG, l, d, newsG, hom;

    # get the corresponding pcp group
    info := AlmostCrystallographicInfo( G );
    H := AlmostCrystallographicPcpGroup( info.dim, info.type, info.param );
    gensH := GeneratorsOfGroup(H);
    gensG := GeneratorsOfGroup(G);

    # sort the generators of G according to the generators of H
    l := Length( gensG );
    d := info.dim;
    newsG := Concatenation( Reversed( gensG{[d+1..l]} ), gensG{[1..d]} );
    hom := GroupHomomorphismByImagesNC( G, H, newsG, gensH );
    SetIsInjective( hom, true );
    SetIsSurjective( hom, true );
    return hom;
end );

max_stars_repo_name: alex-konovalov/aclib

max_issues_repo_name: alex-konovalov/aclib

max_forks_repo_name: alex-konovalov/aclib


text:
############################################################################
##
##  access.gi                    IRREDSOL                   Burkhard Höfling
##
##  Copyright © 2003–2016 Burkhard Höfling
##

############################################################################
##
#F  IndicesIrreducibleSolubleMatrixGroups(<n>, <q>, <d>)
##
##  see the IRREDSOL manual
##
InstallGlobalFunction(IndicesIrreducibleSolubleMatrixGroups,
    function( n, q, d )

        local data, inds, perm, i, k, l, max;

        if not IsPosInt(n) or not IsPosInt(d) or not IsPosInt(q)
                or not IsPPowerInt(q) then
            Error("n, q, and d must be positive integers, ",
                "q must be a prime power");
        fi;
        if n mod d <> 0 then
            return [];
        fi;

        LoadAbsolutelyIrreducibleSolubleGroupData(n/d, q^d);

        if d < n then
            data := IRREDSOL_DATA.GROUPS[n/d][q^d];
            inds := Filtered([1..Length(data)], i -> IsBound(data[i]));
                # as a last resort, redundant groups can be removed
        else
            inds := [1..Length(IRREDSOL_DATA.GROUPS_DIM1 [q^d])];
            if d > 1 then
                if IRREDSOL_DATA.GROUPS_DIM1 [q^d][1] = 1 then
                    RemoveSet(inds, 1); # rewriting the trivial group over a smaller field yields a reducible group
                fi;
            fi;
        fi;
        perm := IRREDSOL_DATA.GAL_PERM[n/d][q^d]^LogInt(q, SmallestRootInt(q));
        # permutation of a generator of Gal(GF(q^d)/GF(q)) on inds

        if perm <> () then
            i := PositionSorted(inds, SmallestMovedPoint(perm));
            max := LargestMovedPoint(perm);
            while i <= Length(inds) do
                k := inds[i];
                if k > max then
                    break;
                fi;
                l := k^perm;
                while l <> k do
                    RemoveSet(inds, l);
                    l := l^perm;
                od;
                i := i + 1;
            od;
        fi;
        MakeImmutable(inds);
        return inds;
    end);

############################################################################
##
#F  PermCanonicalIndexIrreducibleSolubleMatrixGroup(<n>, <q>, <d>, <k>
##
InstallGlobalFunction(PermCanonicalIndexIrreducibleSolubleMatrixGroup,
    function( n, q, d, k )

        local perm, pow, l, orb, min, powmin;

        LoadAbsolutelyIrreducibleSolubleGroupData(n/d, q^d);

        perm := IRREDSOL_DATA.GAL_PERM[n/d][q^d]^LogInt(q, SmallestRootInt(q));

        # permutation of a generator of Gal(GF(q^d)/GF(q)) on parameters

        powmin := 0;
        l := k^perm; # check whether k is least in orbit
        pow := 1;
        min := k;
        orb := [k];
        while l <> k do
            Add(orb, l);
            # we have l = k^(perm^pow)
            if l < k then
                powmin := pow;
                min := l;
            fi;
            l := l^perm;
            pow := pow + 1;
        od;
        return rec(perm := perm, pow := powmin, orb := orb, min := min);
    end);

############################################################################
##
#F  IrreducibleSolubleMatrixGroup(<n>, <q>, <d>, <k>)
##
##  see the IRREDSOL manual
##
InstallGlobalFunction(IrreducibleSolubleMatrixGroup,
    function( n, q, d, k )

        local perm, l, n0, q0, p, o, i, bas, mat, C, c, gddesc, desc, pres, gens, pcgs, grp, hom;

        if not IsPosInt(n) or not IsPosInt(d) or not IsPosInt(q)
                or not IsPPowerInt(q)  or not IsPosInt(k) or not n mod d = 0 then
            Error("n, q, d, and k must be positive integers, q must be a prime power ",
                "and d must divide n");
        fi;

        n0 := n;
        q0 := q;
        n := n0/d;
        q := q0^d;
        p := SmallestRootInt(q0);

        LoadAbsolutelyIrreducibleSolubleGroupData(n, q);

        if d > 1 then # rewrite as matrix group over subfield

            # switch to larger field

            # compute the permutation of a generator of Gal(GF(q)/GF(q0))
            # on the ccls of absolutely irreducible subgroups of GL(n,q)

            perm := IRREDSOL_DATA.GAL_PERM[n][q]^LogInt(q0, p);

            # check if k is a valid paramter, i. e., least in orbit

            l := k^perm;
            while l <> k do
                if l < k then
                    Error("inadmissible value for k");
                fi;
                l := l^perm;
            od;

            Info(InfoIrredsol, 3, "Constructing irreducible group with id ",
                [n0, q0, d, k]);

            bas := CanonicalBasis(AsVectorSpace(GF(q0), GF(q)));
                      # it is important to use CanonicalBasis here, in order to be sure
                      # that the result is the same when called multiple times
        fi;
        if n = 1 then
            if not IsBound(IRREDSOL_DATA.GROUPS_DIM1[q][k]) then
                Error("inadmissible value for k");
            fi;
            o := IRREDSOL_DATA.GROUPS_DIM1[q][k][1];

            mat := [[Z(q)^((q-1)/o)]];
            if d > 1 then
                mat := BlownUpMat(bas, mat);
            fi;
            grp := GroupWithGenerators([mat], IdentityMat(n0, GF(q0)));
            SetSize(grp, o);
            SetIsCyclic(grp, true);
            i := PositionSet(DivisorsInt(LogInt(q, p)), LogInt(q0, p));
            if d = 1 then
                SetMinimalBlockDimensionOfMatrixGroup(grp,  1);
            else
                SetMinimalBlockDimensionOfMatrixGroup(grp, IRREDSOL_DATA.GROUPS_DIM1[q][k][2][i]);
            fi;

            if o = 1 then
                hom := IdentityMapping(grp);
            else
                Assert(1, Length(IRREDSOL_DATA.GUARDIANS[1][q]) = 1);
                C := IRREDSOL_DATA.GUARDIANS[1][q][1];
                Assert(1, Length(MinimalGeneratingSet(C)) = 1);
                c := MinimalGeneratingSet(C)[1]^((q-1)/o);
                Assert(1, Order(c) = o);
                hom := GroupHomomorphismByImagesNC(SubgroupNC(C, [c]), grp, [c], [mat]);
                SetIsBijective(hom, true);
            fi;
        else
            if not IsBound(IRREDSOL_DATA.GROUPS[n][q][k]) then
                Error("inadmissible value for k");
            fi;

            # construct group and isomorphic pc group

            desc  := IRREDSOL_DATA.GROUPS[n][q][k];
            gddesc := IRREDSOL_DATA.GUARDIANS[n][q][desc[1]];
            pres := gddesc[3];

            pcgs := CanonicalPcgsByNumber(FamilyPcgs(Source(pres)), desc[2]);
            gens := List(pcgs, x -> ImageElm(pres, x));
            if d > 1 then
                gens := List(gens, x -> BlownUpMat(bas, x));
            fi;

            grp := GroupWithGenerators(gens, IdentityMat(n0, GF(q0)));
            SetSize( grp, Product(RelativeOrders(pcgs)) );

            hom := GroupGeneralMappingByImagesNC(GroupOfPcgs(pcgs), grp,
                pcgs, gens);
            SetIsGroupHomomorphism(hom, true);
            SetIsBijective(hom, true);

            # look up minimal block dimension
            if d = 1 then
                SetMinimalBlockDimensionOfMatrixGroup(grp, gddesc[4]);
            else
                i := PositionSet(DivisorsInt(LogInt(q, p)), LogInt(q0, p));
                SetMinimalBlockDimensionOfMatrixGroup(grp, IRREDSOL_DATA.GROUPS[n][q][k][3][i]);
            fi;
        fi;

        SetIdIrreducibleSolubleMatrixGroup(grp, [n0, q0, d, k]);
        SetFieldOfMatrixGroup(grp, GF(q0));
        SetDefaultFieldOfMatrixGroup(grp, GF(q0));
        SetTraceField(grp, GF(q0));
        SetConjugatingMatTraceField(grp, One(grp));
        SetRepresentationIsomorphism(grp, hom);
        SetIsPrimitiveMatrixGroup(grp, MinimalBlockDimensionOfMatrixGroup(grp) = n);
        SetIsIrreducibleMatrixGroup(grp, true);
        SetIsAbsolutelyIrreducibleMatrixGroup(grp, d = 1);
        SetIsSolvableGroup(grp, true);
        return grp;
    end);

############################################################################
##
#F  IndicesMaximalAbsolutelyIrreducibleSolubleMatrixGroups(<n>, <q>)
##
##  see the IRREDSOL manual
##
InstallGlobalFunction(IndicesMaximalAbsolutelyIrreducibleSolubleMatrixGroups,
    function( n, q )

        if not IsPosInt(n) or not IsPPowerInt(q)  then
            Error("n and q must be positive integers and q must be a prime power");
        fi;
        LoadAbsolutelyIrreducibleSolubleGroupData(n, q);
        if n = 1 then
            return Immutable([Length(IRREDSOL_DATA.GROUPS_DIM1 [q])]);
        else
            return IRREDSOL_DATA.MAX[n][q];
        fi;
    end);

############################################################################
##
#E
##

max_stars_repo_name: fingolfin/irredsol

max_issues_repo_name: fingolfin/irredsol

max_forks_repo_name: fingolfin/irredsol


text:
# Matrix exponentiation is built-in
A := [[0 , 1], [1, 1]];
PrintArray(A);
#   [ [  0,  1 ],
#     [  1,  1 ] ]
PrintArray(A^10);
#   [ [  34,  55 ],
#     [  55,  89 ] ]

ext: gap

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text:
# Copyright 2018-2019, Carnegie Mellon University
# See LICENSE for details

_HCOLProof_apply := (R, nt, children, nonterms) -> Checked(IsRule(R), IsSPL(nt),
    When(IsNewRule(R),
     R.apply(nt, children, nonterms),
     When(NumGenArgs(R.rule)=2,
          R.rule(nt.params, children),
          R.rule(nt.params, children, nonterms)))
);

HCOLProof_Codegen := function(ss, opts)
    local c;
    trace_log.beginStage("Sigma-SPL","icode", ss);
    c := opts.funcgen.generate(ss, opts);
    trace_log.endStage("Sigma-SPL","icode", c);
    return c;
end;

HCOLProof_CodeConversion := function(c, visitor, opts)
    local c2;
    # dead code?
    trace_log.beginStage("icode","icode", c);
    c2 := visitor(c, opts);
    trace_log.beginStage("icode","icode", c2);
    return c2;
end;

max_stars_repo_name: spiral-software/spiral-package-hcol

max_issues_repo_name: spiral-software/spiral-package-hcol

max_forks_repo_name: spiral-software/spiral-package-hcol


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

ImportAll(paradigms.vector);
ImportAll(paradigms.smp);

Class(RulesFuncSimpScratch,RulesFuncSimp);
RewriteRules(RulesFuncSimpScratch, rec(
    #L(mn,m) o (fBase(m,j) (X) f) -> f (X) fBase(j,m)
    LTensorFlip := ARule(fCompose,
       [ @(1,L), [ @(3,fTensor), @(2).cond(e->range(e) = @(1).val.params[2]), ...] ], e -> [ fTensor(Copy(Drop(@(3).val.children(), 1)), Copy(@(2).val)) ] ),
));

Class(RulesSumsScratch, RulesSums);
RewriteRules(RulesSumsScratch, rec(

    merge_LSKernel := ARule(Compose, [@(1, LSKernel), @(2, LSKernel)],
        e -> [LSKernel(Compose([@(1).val.child(1), @(2).val.child(1)]), @(1).val.mergeInfo(@(1).val.info, @(2).val.info))]),

    mark_DMAGS := Rule([Compose, @(1, Scat), @(2, LSKernel), @(3, Gath)], e -> DMAScat(@(1).val.func) * @(2).val * DMAGath(@(3).val.func) ),

    mark_SWPSum := Rule([DMAFence, [@(1, ISum), @(2)]], e -> let(
            var := @(1).val.var,
            dom := @(1).val.domain,
            comp := Compose(@(2).val._children),
            DMAFence(SWPSum(var, dom, comp)))),

    Scat_DMAFence := ARule( Compose,
       [ @(1, [Scat, Prm, Diag]),
         @(2, [DMAFence, ISum]) ],
        e -> [ CopyFields(@(2).val, rec(
             _children :=  List(@(2).val._children, c -> @(1).val * c),
             dimensions := [Rows(@(1).val), Cols(@(2).val)] )) ]),

    Gath_DMAFence := ARule( Compose,
       [ @(1, [DMAFence, ISum]),
         @(2, [Gath, Prm, Diag]) ],
        e -> [ CopyFields(@(1).val, rec(
                _children := List(@(1).val._children, c -> c * @(2).val),
                dimensions := [Rows(@(1).val), Cols(@(2).val)] )) ]),

    DMAFence_fold := Rule([@(1, DMAFence), @(2, DMAFence)],
        e -> @(2).val),

    Diag_fPrecompute := Rule([Diag, @(2, fPrecompute)], e -> Diag(@(2).val._children[1])),

    Diag_LSKernel := ARule( Compose,
       [ @(1, Diag), @(2, LSKernel) ],
        e -> [LSKernel(@(1).val * @(2).val.child(1), CopyFields(@(2).val.info, rec(opcount := @(2).val.info.opcount + 6 * Rows(@(1).val))))]),

    LSKernel_Diag := ARule( Compose,
        [ @(1, LSKernel), @(2, Diag) ],
        e -> [
            LSKernel(@(1).val.child(1) * @(2).val,
                CopyFields(@(1).val.info,
                    rec(opcount := @(1).val.info.opcount + 6 * Rows(@(2).val))
                )
            )
        ]
    )
));

Class(RulesRCScratch, RulesRC);
RewriteRules(RulesRCScratch, rec(
   RC_ISum := Rule([RC, @(1, ISum)],
        e -> let(s := @(1).val, ISum(s.var, s.domain, RC(s.child(1))))),

   RC_SWPSum := Rule([RC, @(1, SWPSum)],
        e -> let(s := @(1).val, SWPSum(s.var, s.domain, RC(s.child(1))))),

   RC_DMAFence := Rule([RC, @(1, DMAFence)],
        e -> let(s:=@(1).val, DMAFence(RC(s.child(1))))),

   RC_LSKernel := Rule([RC, @(1, LSKernel)],
        e -> let(s:=@(1).val, LSKernel(RC(s.child(1)), s.info))),

   RC_DMAGath := Rule([RC, @(1, DMAGath)], e -> DMAGath(fTensor(@(1).val.func, fId(2)))),

   RC_DMAScat := Rule([RC, @(1, DMAScat)], e -> DMAScat(fTensor(@(1).val.func, fId(2)))),
));

DMAGath.needInterleavedRight := False;
DMAGath.needInterleavedLeft := False;
DMAGath.cannotChangeDataFormat := True;
DMAGath.totallyCannotChangeDataFormat := True;

DMAScat.needInterleavedRight := False;
DMAScat.needInterleavedLeft := False;
DMAScat.cannotChangeDataFormat := True;
DMAScat.totallyCannotChangeDataFormat := True;


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

########################################################################
#   TRC rules

_VBlkMatSPL := function(mat, vtype)
   local i, n, m, nv, a, v, start, final, diff, m_zeros;
   [n, m] := DimensionsMat(mat);
   start := 1;
   a := [];
   v := vtype.size;
   nv := approx.CeilingRat(n/v);
   m_zeros := Replicate(m, vtype.t.zero().v);
   for i in [1 .. nv] do
       final := start + v - 1;
       if (final > n)  then
           diff := final - n;
           Add(a, List(TransposedMat(mat{[start .. n]} :: Replicate(diff, m_zeros)),
		       vec -> vtype.value(vec)));
       else
           Add(a, List(TransposedMat(mat{[start .. final]}),
		       vec -> vtype.value(vec)));
       fi;
       start := final + 1;
   od;
   return a;
end;

VectorizedMatSPL := (isa, spl) -> let(
    r   := Rows(spl),
    c   := Cols(spl),
    v   := isa.v,
    mat := Cond(IsMat(spl), spl, MatSPL(spl)),
    tt  := VBlk(_VBlkMatSPL(mat, isa.t), v) * VGath_dup(fId(c), v),
    When(Mod(r,v)=0, tt, VScat_pc(r,r,0,v)*tt)
);

NewRulesFor(TRC, rec(
#   base cases
    TRC_vect := rec(
        info := "TRC vect",
        forTransposition := false,
        applicable := t -> t.hasTag(AVecReg),
        children := t -> let(
            tags := t.getTags(),
            vecTag := t.getTag(AVecReg, 1),
            [[ t.params[1].withTags(tags).setWrap(VWrapTRC(vecTag.isa)) ]]
        ),
        apply := (t, c, nt) -> VRC(c[1], t.getTag(AVecReg, 1).v)
    ),

    TRC_VRCLR := rec(
        info := "TRC vect",
        forTransposition := false,
        applicable := t -> t.hasTags() and t.isTag(1, AVecReg),
        children := t -> [[ t.params[1].withTags(t.getTags()) ]],
        apply := (self, t, C, Nonterms) >> RC(C[1]),
        switch := false
    ),

    TRC_cplx := rec(
        forTransposition := false,
        applicable := t -> t.hasTag(AVecRegCx),
        children := t -> let(
            tags := t.getTags(),
            tag  := t.getTag(AVecRegCx, 1),
            nt   := t.params[1],
            [[ nt.withTags(tags).setWrap(VWrapTRCcplx(tag.isa)) ]]
        ),
        apply := (self, t, C, Nonterms) >> vRC(C[1]),
        switch := false
    ),

    # TRC_cplxvect will fire even when asked for real vectorization
    TRC_cplxvect := rec(
        forTransposition := false,
        applicable := t -> t.hasTag(AVecRegCx) or t.hasTag(AVecReg),
        children := t -> let(
            tags := t.getTags(),
            tag := Filtered([t.getTag(AVecReg, 1), t.getTag(AVecRegCx, 1)], i->not IsBool(i))[1],
            nt := t.params[1],
            ctag := When(ObjId(tag)=AVecRegCx, tag, AVecRegCx(tag.isa)),
            newTags := List(tags, i->When(i=tag, ctag, i)),
            [[ nt.withTags(newTags).setWrap(VWrapTRCcplx(tag.isa)) ]]
        ),
        apply := (self, t, C, Nonterms) >> vRC(C[1]),
        switch := false
    ),

    TRC_cplx_v2 := rec(
        forTransposition := false,
        applicable := t -> t.hasTag(AVecRegCx) and t.getTag(AVecRegCx, 1).isa.v = 2,
        children := t -> let(
            tags := t.getTags(),
            tag := t.getTag(AVecRegCx, 1),
            nt := t.params[1],
            newTags := Filtered(tags, i -> i<>tag),
            [[ nt.withTags(newTags).setWrap(VWrapTRCcplx(tag.isa)) ]]
        ),
        apply := (self, t, C, Nonterms) >> vRC(VTensor(C[1], 1)),
        switch := false
    ),

    TRC_By_Def := rec(
         forTransposition := false,
	 maxSize := 128,
         applicable := (self, t) >> Rows(t) <= self.maxSize and Cols(t)<=self.maxSize and
       	                            t.hasTags() and t.isTag(1, AVecReg) and t.free()=[],
         apply := (self, t, C, Nonterms) >> let(
	     isa := t.getTag(AVecReg, 1).isa,
	     VectorizedMatSPL(isa, RC(t.params[1]))
	 )
    )
));


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

fId.canTensorSplit := (self, div) >> self.domain() mod div = 0;
fId.tensorSplit    := (self, div) >> [fId(div), fId(self.size/div)];
fBase.canTensorSplit := (self, div) >> self.range() mod div = 0;
fBase.tensorSplit := (self, div) >> [ fBase(div,        idiv(self.params[2], self.range()/div)),
                                      fBase(self.range()/div, imod(self.params[2], self.range()/div)) ];

CanTensorSplit := (what, div) -> When(not IsRec(what) or not IsBound(what.canTensorSplit), false,
                                      what.canTensorSplit(div));

TensorSplit    := (what, div) -> what.tensorSplit(div);

try_split := function(i, flist, fdim, N)
    local ffrun, next, len, split;
    len := Length(flist);
    ffrun := fdim(flist[i]);
    while i < len and ffrun < N do
        next := flist[i+1];
	if ffrun * fdim(next) > N then
            if N mod ffrun <> 0 or not CanTensorSplit(next, N/ffrun) then
		return false;
        #Error("can't merge incompatible tensor chains (could not split)",
         #     "N=", N, " ffrun=", ffrun, " next=", next);
            else
		split := TensorSplit(next, N/ffrun);
		flist := Concat(flist{[1..i]}, split, flist{[i+2..len]});
		ffrun := ffrun * N/ffrun; # = N
            fi;
	else
            ffrun := ffrun * fdim(next);
	fi;
	i := i+1;
    od;
    #Print([i, flist], "\n");
    return [i, flist];
end;

full_merge_tensor_chains := function(target, ff, gg, compose, ftensor, gtensor, fid, gid, fdom, gran)
   local i, j, nf, ng, res, ffrun, ggrun, ibegin, jbegin, split;
   nf := Length(ff);    ng := Length(gg);
   res := []; i:=1; j:=1;

   while i <= nf or j <= ng do
       # handle domain=1 or range=1, which are always mergeable
       if (i<=nf and fdom(ff[i])=1) or (j<=ng and gran(gg[j])=1) then
	   if i<=nf and fdom(ff[i]) = 1 then
               Add(res, fid(ff[i]));
               i := i + 1;
	   fi;
	   if j<=ng and gran(gg[j]) = 1 then
               Add(res, gid(gg[j]));
               j := j + 1;
	   fi;

       # try to combine terms to get a match
       elif (i <= nf and j <= ng) then
	   if AnySyms(fdom(ff[i]), gran(gg[j])) then
               return false;
           elif fdom(ff[i]) = gran(gg[j]) then
               Add(res, compose(ff[i], gg[j]));
	   else
               ibegin := i;
               jbegin := j;
               if fdom(ff[i]) < gran(gg[j]) then
		   if CanTensorSplit(gg[j], fdom(ff[i])) then
		       split := TensorSplit(gg[j], fdom(ff[i]));
		       gg := Concat(gg{[1..j-1]}, split, gg{[j+1..ng]});
		       ng := ng + 1;
		   else
		       split := try_split(i, ff, fdom, gran(gg[j]));
		       if split=false then return false; fi;
		       i := split[1];
		       ff := split[2]; nf := Length(ff);
		   fi;
               else
		   if CanTensorSplit(ff[i], gran(gg[j])) then
		       split := TensorSplit(ff[i], gran(gg[j]));
		       if split=false then return false; fi;
		       ff := Concat(ff{[1..i-1]}, split, ff{[i+1..nf]});
		       nf := nf + 1;
		   else
		       split := try_split(j, gg, gran, fdom(ff[i]));
		       if split=false then return false; fi;
		       j := split[1];
		       gg := split[2]; ng := Length(gg);
		   fi;
               fi;
               Add(res, compose(ftensor(ff{[ibegin..i]}), gtensor(gg{[jbegin..j]})));
	   fi;
	   i := i+1;
	   j := j+1;
       else
	   return false;
       fi;
   od;
   target.val := res;
   return res;
end;

# *******************************************************************
# make sure ff[i] is not a diag function, where merging chains doesn't make sense (does it?)
fully_compat_tensor_chains := (ff,gg,fdom,gran) ->
    Length(ff) = Length(gg) and
    ForAll([1..Length(ff)], i -> ff[i].range()<>false and
                             fdom(ff[i])=gran(gg[i]));

merge_fc_tensor_chains := (ff,gg,combine) ->
   List([1..Length(ff)], i -> combine(ff[i], gg[i]));

# this assumes compat_domain_range compatibility
compat_tensor_chains := (f,g,fdom,gran) -> let(
    ff := Filtered(f, c->let(d:=fdom(c), IsSymbolic(d) or d > 1)),
    gg := Filtered(g, c->let(r:=gran(c), IsSymbolic(r) or r > 1)),
    fully_compat_tensor_chains(ff, gg, fdom, gran));

# this assumes compat_domain_range compatibility
merge_tensor_chains := function(ff, gg, compose, fidentity, gidentity, fdom, gran)
   local i, j, nf, ng, res;
   nf := Length(ff);    ng := Length(gg);
   res := []; i:=1; j:=1;

   while i <= nf or j <= ng do
       if (i <= nf and j <= ng) and fdom(ff[i]) = gran(gg[j]) then
       Add(res, compose(ff[i], gg[j]));
       i := i+1;
       j := j+1;
       else
       if i<=nf and fdom(ff[i]) = 1 then
           Add(res, fidentity(ff[i]));
           i := i + 1;
       elif j<=ng and gran(gg[j]) = 1 then
           Add(res, gidentity(gg[j]));
           j := j + 1;
       else Error("can't merge incompatible tensor chains");
       fi;
       fi;
   od;
   return res;
end;


text:
# This file has been generated by the GAP build system,
# do not edit manually!
GAParch=FAKE-GAP-ARCH
GAP_ABI=64
GAP_HPCGAP=no

GAP_VERSION="4.12dev"
GAP_BUILD_VERSION="4.12dev"
GAP_LIBTOOL_CURRENT=8
GAP_LIBTOOL_AGE=0
GAP_KERNEL_MAJOR_VERSION=8
GAP_KERNEL_MINOR_VERSION=0

GAP_BIN_DIR="/workspace/destdir/share/gap"
GAP_LIB_DIR="/workspace/destdir/share/gap"

GAP="/workspace/destdir/bin/gap"
GAC="/workspace/destdir/share/gap/gac"

GAP_CC="cc "
GAP_CXX="c++ -std=gnu++11 "
GAP_CFLAGS="-g -O2"
GAP_CXXFLAGS="-g -O2"
GAP_CPPFLAGS="-I/workspace/destdir/include/gap -DHAVE_CONFIG_H -fPIC "
GAP_LDFLAGS="-L/workspace/destdir/lib "
GAP_LIBS="-lgmp"

GAP_OBJS="build/obj/src/ariths.c.lo build/obj/src/bags.c.lo build/obj/src/blister.c.lo build/obj/src/bool.c.lo build/obj/src/calls.c.lo build/obj/src/code.c.lo build/obj/src/collectors.cc.lo build/obj/src/compiler.c.lo build/obj/src/costab.c.lo build/obj/src/cyclotom.c.lo build/obj/src/debug.c.lo build/obj/src/dt.c.lo build/obj/src/dteval.c.lo build/obj/src/error.c.lo build/obj/src/exprs.c.lo build/obj/src/ffdata.c.lo build/obj/src/finfield.c.lo build/obj/src/funcs.c.lo build/obj/src/gap.c.lo build/obj/src/gaptime.c.lo build/obj/build/gap_version.c.lo build/obj/src/gvars.c.lo build/obj/src/hookintrprtr.c.lo build/obj/src/info.c.lo build/obj/src/integer.c.lo build/obj/src/intfuncs.c.lo build/obj/src/intrprtr.c.lo build/obj/src/io.c.lo build/obj/src/iostream.c.lo build/obj/src/libgap-api.c.lo build/obj/src/listfunc.c.lo build/obj/src/listoper.c.lo build/obj/src/lists.c.lo build/obj/src/macfloat.c.lo build/obj/src/modules_builtin.c.lo build/obj/src/modules.c.lo build/obj/src/objcftl.c.lo build/obj/src/objects.c.lo build/obj/src/objfgelm.cc.lo build/obj/src/objpcgel.cc.lo build/obj/src/objset.c.lo build/obj/src/opers.cc.lo build/obj/src/permutat.cc.lo build/obj/src/plist.c.lo build/obj/src/pperm.cc.lo build/obj/src/precord.c.lo build/obj/src/profile.c.lo build/obj/src/range.c.lo build/obj/src/rational.c.lo build/obj/src/read.c.lo build/obj/src/records.c.lo build/obj/src/saveload.c.lo build/obj/src/scanner.c.lo build/obj/src/sctable.c.lo build/obj/src/set.c.lo build/obj/src/stats.c.lo build/obj/src/streams.c.lo build/obj/src/stringobj.c.lo build/obj/src/syntaxtree.c.lo build/obj/src/sysfiles.c.lo build/obj/src/sysroots.c.lo build/obj/src/sysstr.c.lo build/obj/src/system.c.lo build/obj/src/tietze.c.lo build/obj/src/tracing.c.lo build/obj/src/trans.cc.lo build/obj/src/trycatch.c.lo build/obj/src/vars.c.lo build/obj/src/vec8bit.c.lo build/obj/src/vecffe.c.lo build/obj/src/vecgf2.c.lo build/obj/src/vector.c.lo build/obj/src/weakptr.c.lo build/obj/src/julia_gc.c.lo build/obj/src/c_oper1.c.lo build/obj/src/c_type1.c.lo build/obj/src/compstat.c.lo"

JULIA=""
JULIA_CPPFLAGS="-I/workspace/destdir/include/julia -fPIC"
JULIA_LDFLAGS="-L/workspace/destdir/lib -L/workspace/destdir/lib/julia"
JULIA_LIBS="-Wl,-rpath,/workspace/destdir/lib -Wl,-rpath,/workspace/destdir/lib/julia -ljulia"

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


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

LocalConfig.bench := rec();

if IsBound(platforms.sse) then         LocalConfig.bench.SSE     := platforms.sse.benchSSE; fi;
if IsBound(platforms.intel) then       LocalConfig.bench.SMP_SSE := platforms.intel.benchCore2; fi;
if IsBound(platforms.benchScalar) then LocalConfig.bench.scalar  := platforms.benchScalar; fi;
if IsBound(platforms.avx) then         LocalConfig.bench.AVX     := platforms.avx.benchAVX; fi;

_dumpBench := function(obj, objPath, masterList, level)
	local newObj, newObjPath, fieldNames, name, size;

	if IsClass(obj) then
		if obj.__name__ = "DPBench" then
			if (level > 0) and IsBound(obj.sizes) and EndsWith(objPath, ")") then
				newObjPath := DropLast(Copy(objPath), 1);
				for size in obj.sizes do
					Append(masterList, [newObjPath::"["::String(size)::"])"]);
				od;
			else
				Append(masterList, [Copy(objPath)]);
			fi;
		fi;
		return;
	elif IsRec(obj) then
		fieldNames := UserRecFields(obj);
		for name in fieldNames do
			newObj := obj.(name);
			newObjPath := Copy(objPath)::"."::Copy(name);
			_dumpBench(newObj, newObjPath, masterList, level);
		od;
	elif IsFunc(obj) then
		newObj := obj();
		newObjPath := Copy(objPath)::"()";
		_dumpBench(newObj, newObjPath, masterList, level);
	fi;
end;

#F
#F DumpBenches(platform, level)
#F     platform : string, name of platform, eg., "SSE"
#F     level    : integer, level of detail
#F                0 : build benches that test a list of sizes
#F                1 : build an individual bench for each size
#F
#F Returns a list of all bench constructors for the specified platform
#F

DumpBenches := function(platform, level)
	local benchlist;

	benchlist := [];

	if IsBound(LocalConfig.bench.(platform)) and IsFunc(LocalConfig.bench.(platform)) then
		_dumpBench(LocalConfig.bench.(platform), "LocalConfig.bench."::platform, benchlist, level);
	fi;

	return Copy(benchlist);
end;

#F
#F DumpAllBenches(level = 0)
#F     level    : integer, optional (see DumpBenches)
#F
#F Returns a list of bench constructors for all supported platforms.
#F

DumpAllBenches := function(arg)
	local level, platform, platforms, benchlist;

	level := Cond(IsBound(arg[1]) and IsInt(arg[1]), arg[1], 0);

	##platforms := UserRecFields(LocalConfig.bench);
	platforms := ["SSE", "AVX"];

	benchlist := [];

	for platform in platforms do
		Append(benchlist, DumpBenches(platform, level));
	od;

	return Copy(benchlist);
end;


text:
# Copyright 2018-2019, Carnegie Mellon University
# See LICENSE for details

Declare(ScatHUnion, eUnion);

Class(ISumUnion, ISum, rec(
    toOperator := self >> (vec -> Checked(Length(vec)=self.dims()[2],
        Flat(List(TransposedMat(List([0..self.var.range-1],
        j->RulesStrengthReduce(SubstVars(Copy(self._children[1]), rec((self.var.id):=V(j)))).toOperator()(vec))), i->Filtered(i, k->k<>()))
    )))
));

Class(ISumReduction, ISum, rec(
    abbrevs := [ (v, rng, op, idval, isSaturated, expr) -> [v, rng, op, idval, isSaturated, expr] ],
    new := meth(self, var, domain, op, idval, isSaturated, expr)
        local res;
        Constraint(IsSPL(expr));
        # if domain is an integer (not symbolic) it must be positive
        Constraint(not IsInt(domain) or domain >= 0);
        var.isLoopIndex := true;
        res := SPL(WithBases(self, rec(_children := [expr], var := var, domain := domain, op := op, idval := idval, isSaturated := isSaturated)));
        res.dimensions := res.dims();
        return res;
    end,
    print := (self, i, is) >> Print(
        self.name, "(", self.var, ", ", self.domain, ", ", self.op, ", ", self.idval, ", ", self.isSaturated, ",\n",
        Blanks(i+is), self._children[1].print(i+is, is), "\n",
        Blanks(i), ")", self.printA(),
        When(IsBound(self._setDims), Print(".overrideDims(", self._setDims, ")"), Print(""))
    ),
    toOperator := self >> (vec -> Checked(Length(vec)=self.dims()[2],
        List(List(TransposedMat(List([0..self.var.range-1],
        j->RulesStrengthReduce(SubstVars(Copy(self._children[1]), rec((self.var.id):=V(j)))).toOperator()(vec))), i->Filtered(i, k->k<>())),
        e->_unwrap(FoldL(e, self.op, self.idval)))))
));

Class(GathH, SumsBase, BaseMat, rec(
    rChildren := self >> [self.N, self.n, self.base, self.stride],
    rSetChild := rSetChildFields("N", "n", "base", "stride"),
    new := (self, N, n, base, stride) >>
	SPL(WithBases(self, rec(
      	N:= N, n := n, base := base, stride := stride, dimensions := [n, N]))),
    dims := self >> self.dimensions,
    sums := self >> self,
    area := self >> self.n,
    isReal := self >> true,
    transpose := self >> ScatHUnion(self.N, self.n, self.base, self.stride),
    conjTranspose := self >> self.transpose(),
    inverse := self >> self.transpose(),
    toAMat := self >> let(
	   n := self.n,
       N := self.N,
       base := self.base,
       stride := self.stride,
       AMatMat(List([0..n-1], row -> let(
         idx := base+row*stride,
		 BasisVec(N, idx))))
    ),
    isIdentity := self >> (self.n=self.N) and self.base=0 and self.stride = 1,
    toOperator := self >> (vec-> Checked(Checked(Length(vec) = self.N, List([0..self.n-1], i->vec[_unwrap(self.base+i*self.stride+1)]))))
));

Class(ScatHUnion, SumsBase, BaseMat, rec(
    rChildren := self >> [self.N, self.n, self.base, self.stride],
    rSetChild := rSetChildFields("N", "n", "base", "stride"),
    new := (self, N, n, base, stride) >>

    SPL(WithBases(self, rec(
        N:= N, n := n, base := base, stride := stride, dimensions := [N, n]))),
    dims := self >> self.dimensions,
    sums := self >> self,
    area := self >> self.n,
    isReal := self >> true,
    transpose := self >> GathH(self.N, self.n, self.base, self.stride),
    conjTranspose := self >> self.transpose(),
    inverse := self >> self.transpose(),
    toAMat := self >> let(
	   n := self.n,
       N := self.N,
       base := self.base,
       stride := self.stride,
       TransposedAMat(AMatMat(List([0..n-1], row -> let(
         idx := base+row*stride,
		 BasisVec(N, idx))))
    )),
    isIdentity := self >> (self.n=self.N) and self.base=0 and self.stride = 1,
    toOperator := self >>(vec -> Checked(Length(vec) = self.n, List([0..self.N-1], j->let(i := (j-_unwrap(self.base))/_unwrap(self.stride), When(IsPosInt0(i) and i < _unwrap(self.n), vec[i+1], ())))))
));

Class(ScatHAcc, ScatHUnion);

Class(eT, SumsBase, BaseMat, rec(
    rChildren := self >> [self.N, self.base],
    rSetChild := rSetChildFields("N", "base"),
    new := (self, N, base) >> SPL(WithBases(self, rec(
      	N:= N, base := base, dimensions := [1, N]))),
    dims := self >> self.dimensions,
    sums := self >> self,
    area := self >> self.N,
    isReal := self >> true,
    transpose := self >> eUnion(self.N, self.base),
    conjTranspose := self >> self.transpose(),
    inverse := self >> self.transpose(),
    toAMat := self >> let(
	   n := self.n,
       N := self.N,
       base := self.base,
       stride := self.stride,
       AMatMat(List([0..n-1], row -> let(
         idx := base+row*stride,
		 BasisVec(N, idx))))
    ),
    isIdentity := self >> false,#(self.N=1),
    toOperator := self >>(vec -> Checked(Length(vec) = self.N, [vec[_unwrap(self.base)+1]]))
));

Class(eUnion, SumsBase, BaseMat, rec(
    rChildren := self >> [self.N, self.base],
    rSetChild := rSetChildFields("N", "base"),
    new := (self, N, base) >> SPL(WithBases(self, rec(
      	N:= N, base := base, dimensions := [N, 1]))),
    dims := self >> self.dimensions,
    sums := self >> self,
    area := self >> 1,
    isReal := self >> true,
    transpose := self >> eT(self.N, self.base),
    conjTranspose := self >> self.transpose(),
    inverse := self >> self.transpose(),
    toAMat := self >> let(
	   n := 1,
       N := self.N,
       base := self.base,
       stride := 1,
       TransposedAMat(AMatMat(List([0..n-1], row -> let(
         idx := base+row*stride,
		 BasisVec(N, idx))))
    )),
    isIdentity := False,
    toOperator := self >>(vec -> Checked(Length(vec) = 1, List([0..self.N-1], i->When(i = _unwrap(self.base), vec[1], ()))))
));

Class(Mul1, Blk1, rec(
    toOperator := self >> (vec -> Checked(Length(vec) = 1, [_unwrap(RulesStrengthReduce(vec[1]*_unwrap(self.element)))]))
));

Class(SUMUnion, SUM, rec(
    toOperator := self >> (vec -> Checked(Length(vec)=self.dims()[2],
        Flat(List(TransposedMat(List(self.children(), c->RulesStrengthReduce(c.toOperator()(vec)))), i->Filtered(i, k->k<>())))))
));

Class(Inductor,  BaseMat, rec(
    abbrevs := [ (N, i, op, initval) -> [N, i, op, initval] ],
    new := (self, N, i, op, initval) >> SPL(WithBases(self, rec(
        N := N, i := i, op := op, initval := initval, dimensions := [1, 1]))),
    print := (self, i, is) >> Print("Inductor(", self.N, ", ", self.i, ", ", self.op, ", ", self.initval, ")"),
    dims  := self >> self.dimensions,
    rChildren := self >> [self.N, self.i, self.op, self.initval],
    rSetChild := rSetChildFields("N", "i", "op", "initval"),
    isReal := True,
    recurrence := self >> ((j,v)->When(j=0, self.initval, self.op.at(self.recurrence()(j-1,v),v))),
    toOperator := self >> (vec -> Checked(Length(vec)=self.dims()[2], [List([0..self.N-1], i->_unwrap(RulesStrengthReduce(self.recurrence()(i,vec[1]))))[_unwrap(self.i)+1]]))
));