text:
# Arrays are better called lists in GAP. Lists may have elements of mixed types, e$
v := [ 10, 7, "bob", true, [ "inner", 5 ] ];
# [ 10, 7, "bob", true, [ "inner", 5 ] ]

# List index runs from 1 to Size(v)
v[1];
# 10

v[0];
# error

v[5];
# [ "inner", 5 ]

v[6];
# error

# One can assign a value to an undefined element
v[6] := 100;

# Even if it's not after the last: a list may have undefined elements
v[10] := 1000;
v;
# [ 10, 7, "bob", true, [ "inner", 5 ], 100,,,, 1000 ]

# And one can check for defined values
IsBound(v[10]);
# true

IsBound(v[9]);
# false

# Size of the list
Size(v);
# 10

# Appending a list to the end of another
Append(v, [ 8, 9]);
v;
# [ 10, 7, "bob", true, [ "inner", 5 ], 100,,,, 1000, 8, 9 ]

# Adding an element at the end
Add(v, "added");
v;
# [ 10, 7, "bob", true, [ "inner", 5 ], 100,,,, 1000, 8, 9, "added" ]


text:
# Built-in

Filtered([1 .. 100], IsPrime);
# [ 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97 ]

Filtered([1 .. 10], IsEvenInt);
# [ 2, 4, 6, 8, 10 ]

Filtered([1 .. 10], IsOddInt);
# [ 1, 3, 5, 7, 9 ]

ext: gap

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


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

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

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

  return GroupOfAutomFamily(AutomFamily(list, bind_vars));
end);

###############################################################################
##
#M  AutomatonGroup(<list>, <names>)
##
InstallMethod(AutomatonGroup, "for [IsList, IsList]", [IsList, IsList],
function(list, names)
  return AutomatonGroup(list, names, AG_Globals.bind_vars_autom_family);
end);

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

  return GroupOfAutomFamily(AutomFamily(list, names, bind_vars));
end);

###############################################################################
##
#M  AutomatonGroup(<string>)
#M  AutomatonGroup(<string>, <bind_vars>)
##
InstallMethod(AutomatonGroup, "for [IsString]", [IsString],
function(string)
  return AutomatonGroup(string, AG_Globals.bind_vars_autom_family);
end);

InstallMethod(AutomatonGroup, "for [IsString, IsBool]", [IsString, IsBool],
function(string, bind_vars)
  local s;
  s := AG_ParseAutomatonString(string);
  return AutomatonGroup(s[2], s[1], bind_vars);
end);

InstallMethod(AutomatonGroup, "for [IsMealyAutomaton]", [IsMealyAutomaton],
function(A)
  if not IsInvertible(A) then
    Error("Automaton <A> is not invertible");
  fi;
  return AutomatonGroup(AutomatonList(A), A!.states);
end);

InstallMethod(AutomatonGroup, "for [IsMealyAutomaton, IsBool]", [IsMealyAutomaton, IsBool],
function(A, bind_vars)
  if not IsInvertible(A) then
    Error("Automaton <A> is not invertible");
  fi;
  return AutomatonGroup(AutomatonList(A), A!.states, bind_vars);
end);

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

###############################################################################
##
#M  IsGroupOfAutomFamily(<G>)
##
InstallMethod(IsGroupOfAutomFamily, "for [IsAutomGroup]",
              [IsAutomGroup],
function(G)
  return G = GroupOfAutomFamily(G);
end);

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

###############################################################################
##
#M  __AG_SubgroupOnLevel(<G>, <gens>, <level>)
##
InstallMethod(__AG_SubgroupOnLevel, [IsAutomGroup,
                                    IsList and IsTreeAutomorphismCollection,
                                    IsPosInt],
function(G, gens, level)
  local overgroup;

  if IsEmpty(gens) or (Length(gens) = 1 and IsOne(gens[1])) then
    return TrivialSubgroup(G);
  fi;

  if HasIsGroupOfAutomFamily(G) and IsGroupOfAutomFamily(G) then
    overgroup := G;
  else
    overgroup := GroupOfAutomFamily(UnderlyingAutomFamily(G));
  fi;

  return SubgroupNC(overgroup, gens);
end);

InstallOtherMethod(__AG_SubgroupOnLevel, [IsAutomGroup, IsList and IsEmpty, IsPosInt],
function(G, gens, level)
  return TrivialSubgroup(G);
end);

InstallMethod(__AG_SubgroupOnLevel, [IsTreeAutomorphismGroup,
                                    IsList and IsInvertibleAutomCollection,
                                    IsPosInt],
function(G, gens, level)
  local overgroup;

  overgroup := GroupOfAutomFamily(FamilyObj(gens[1]));

  if Length(gens) = 1 and IsOne(gens[1]) then
    return TrivialSubgroup(overgroup);
  fi;

  return SubgroupNC(overgroup, gens);
end);

InstallMethod(__AG_SimplifyGroupGenerators, [IsList and IsInvertibleAutomCollection],
function(gens)
  local words, fam;

  if IsEmpty(gens) then
    return [];
  fi;

  fam := FamilyObj(gens[1]);
  words := FreeGeneratorsOfGroup(Group(List(gens, a -> a!.word)));

  if fam!.use_rws and not IsEmpty(words) then
    words := AG_ReducedForm(fam!.rws, words);
    if not IsEmpty(words) then
      words := FreeGeneratorsOfGroup(Group(words));
    fi;
  fi;

  if IsEmpty(words) then
    return [];
  fi;

  return List(words, w -> Autom(w, fam));
end);

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

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

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

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

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

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

  gens := GeneratorsOfGroup(G);

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

  return s;
end);

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

###############################################################################
##
#M  MihailovaSystem(G)
##
##  TODO XXX test it
##
InstallMethod(MihailovaSystem, "for [IsAutomatonGroup]", [IsAutomatonGroup],
function (G)
  local gens, mih, mih_gens, i;

  if not IsActingOnBinaryTree(G) then
    Error("MihailovaSystem(IsAutomGroup):\n  sorry, group is not acting on binary tree\n");
  fi;
  if not IsFractalByWords(G) then
    Info(InfoAutomGrp, 1, "given group is not IsFractalByWords");
    return fail;
  fi;

  gens := GeneratorsOfGroup(StabilizerOfFirstLevel(G));
  mih := AG_ComputeMihailovaSystemPairs(List(gens, a -> StatesWords(a)));

  if mih = fail then
    return fail;
  elif not mih[3] then
    return gens;
  fi;

  mih_gens := [];
  for i in [1..Length(gens)] do
    mih_gens[i] := AG_CalculateWord(mih[2][i], gens);
  od;
  return mih_gens;
end);

###############################################################################
##
#M  IsFractalByWords(G)
##
InstallMethod(IsFractalByWords, "for [IsAutomGroup]",
              [IsAutomGroup],
function (G)
  local freegens, stab, i, sym, f;

  sym := GroupWithGenerators(List(GeneratorsOfGroup(G), g -> Perm(g)));
  if not IsTransitive(sym, [1..DegreeOfTree(G)]) then
    Info(InfoAutomGrp, 1, "group is not transitive on first level");
    return false;
  fi;

  f := GroupWithGenerators(List(GeneratorsOfGroup(G), g -> Word(g)));
  stab := StabilizerOfFirstLevel(G);
  stab := List(GeneratorsOfGroup(stab), a -> StatesWords(a));

  for i in [1..DegreeOfTree(G)] do
    if f <> GroupWithGenerators(List(stab, s -> s[i])) then
      return false;
    fi;
  od;
  return true;
end);

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

  if CanEasilyTestSphericalTransitivity(G) and IsSphericallyTransitive(G) then
    Info(InfoAutomGrp, 3, "Size(G): infinity, G is spherically transitive");
    return infinity;
  fi;

  if IsFractalByWords(G) then
    Info(InfoAutomGrp, 3, "Size(G): infinity, G is fractal by words");
    return infinity;
  fi;

  if HasIsFractal(G) and IsFractal(G) then
    Info(InfoAutomGrp, 3, "Size(G): infinity, G is fractal");
    return infinity;
  fi;

  if IsAutomatonGroup(G) then
    A:=MealyAutomaton(AutomatonList(G));
    if IsMDTrivial(A) then
      Info(InfoAutomGrp, 3, "Size(G): automaton generating G is MD-trivial");
      lev:=LevelOfFaithfulAction(G,infinity);
      return Size(G);
    elif (DegreeOfTree(G)=2 or NumberOfStates(A)=2) and IsIRAutomaton(A) then
      Info(InfoAutomGrp, 3, "Size(G): infinity, G is generated by 2-letter or 2-state not MD-trivial IR-automaton ");
      return infinity;
    fi;
  fi;

  if IsAutomatonGroup(G) and LevelOfFaithfulAction(G,8)<>fail then
    return Size(G);
  fi;

  f := FindElementOfInfiniteOrder(G,10,10);

  if HasSize(G) or f <> fail then
    return Size(G);
  fi;

  Info(InfoAutomGrp, 1, "You can try to use IsomorphismPermGroup(<G>) or\n",
                        "   FindElementOfInfiniteOrder(<G>,<length>,<depth>) with bigger bounds");
  TryNextMethod();
end);

InstallOtherMethod(LevelOfFaithfulAction, "for [IsAutomGroup and IsSelfSimilar]",
              [IsAutomGroup and IsSelfSimilar,IsCyclotomic],
function(G,max_lev)
  local s,s_next,lev;
  if HasIsFinite(G) and not IsFinite(G) then return fail; fi;
  if HasLevelOfFaithfulAction(G) then return LevelOfFaithfulAction(G); fi;
  lev := 0; s := 1; s_next := Size(PermGroupOnLevel(G,1));
  while s<s_next and lev<max_lev do
    lev := lev+1;
    s := s_next;
    s_next := Size(PermGroupOnLevel(G,lev+1));
  od;
  if s = s_next then
    SetSize(G,s);
    SetLevelOfFaithfulAction(G,lev);
    return lev;
  else
    return fail;
  fi;
end);

InstallMethod(LevelOfFaithfulAction, "for [IsAutomGroup and IsSelfSimilar]",
              [IsAutomGroup and IsSelfSimilar],
function(G)
  return LevelOfFaithfulAction(G,infinity);
end);

InstallOtherMethod(IsomorphismPermGroup, "for [IsAutomGroup and IsSelfSimilar,IsCyclotomic]",
                   [IsAutomGroup and IsSelfSimilar, IsCyclotomic],
function (G, n)
  local H, lev;
  lev := LevelOfFaithfulAction(G, n);
  if lev <> fail then
    H := PermGroupOnLevel(G,LevelOfFaithfulAction(G));
    return AG_GroupHomomorphismByImagesNC(G, H, GeneratorsOfGroup(G), GeneratorsOfGroup(H));
  fi;
  return fail;
end);

InstallMethod(IsomorphismPermGroup, "for [IsAutomGroup]",
              [IsAutomGroup], SUM_FLAGS,
function (G)
  local H;
  H := AG_FiniteGroupId(G);
  return AG_GroupHomomorphismByImagesNC(G, H, GeneratorsOfGroup(G), GeneratorsOfGroup(H));
end);

BindGlobal("TestSelfSimilarity",
function(G)
  if CanEasilyTestSelfSimilarity(G) then
    IsSelfSimilar(G);
    return true;
  fi;

  if IsTrivial(G) then
    SetIsSelfSimilar(G, true);
    return true;
  fi;

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

  if Set(GeneratorsOfGroup(G)) = Set(GeneratorsOfGroup(GroupOfAutomFamily(UnderlyingAutomFamily(G)))) then
    SetIsSelfSimilar(G, true);
    return true;
  fi;

  return false;
end);

###############################################################################
##
#M  IsSphericallyTransitive(G)
##
InstallMethod(IsSphericallyTransitive, "for [IsAutomGroup]",
              [IsAutomGroup],
function (G)
  local x, rat_gens, abel_hom, lev;

  if DegreeOfTree(G)=1 then
    Info(InfoAutomGrp, 3, "IsSphericallyTransitive(G): true");
    Info(InfoAutomGrp, 3, "  G acts on 1-ary tree");
    return true;
  fi;

  if IsFractalByWords(G) then
    Info(InfoAutomGrp, 3, "IsSphericallyTransitive(G): true");
    Info(InfoAutomGrp, 3, "  G is fractal");
    return true;
  fi;

  if IsTrivial(G) then
    Info(InfoAutomGrp, 3, "IsSphericallyTransitive(G): false");
    Info(InfoAutomGrp, 3, "  G is trivial: G = ", G);
    return false;
  fi;

  if HasIsFinite(G) and IsFinite(G) then
    Info(InfoAutomGrp, 3, "IsSphericallyTransitive(G): false");
    Info(InfoAutomGrp, 3, "  IsFinite(G): G = ", G);
    return false;
  fi;

  if DegreeOfTree(G) = 2 and TestSelfSimilarity(G) and IsSelfSimilar(G) then
    if HasIsFinite(G) and IsFinite(G) = false then
      Info(InfoAutomGrp, 3, "IsSphericallyTransitive(G): true");
      Info(InfoAutomGrp, 3, "  <G> is infinite self-similar acting on binary tree");
      return true;
    fi;
    if PermGroupOnLevel(G,2) = Group((1,4,2,3)) then
      Info(InfoAutomGrp, 3, "IsSphericallyTransitive(G): true");
      Info(InfoAutomGrp, 3, "  any element which acts transitively on the first level acts spherically transitively");
      return true;
    fi;
  fi;

  for lev in [1..8] do
    if not IsTransitiveOnLevel(G,lev) then
      Info(InfoAutomGrp, 3, "IsSphericallyTransitive(G): false");
      Info(InfoAutomGrp, 3, "  the group does not act transitively on level ", lev);
      return false;
    fi;
  od;

  TryNextMethod();
end);

###############################################################################
##
#M  DiagonalPower(<G>, <n>)
##
InstallOtherMethod( DiagonalPower,
                    "for [IsAutomGroup and IsGroupOfAutomFamily, IsPosInt]",
                    [IsAutomGroup and IsGroupOfAutomFamily, IsPosInt],
function(G, n)
  return DiagonalPower(UnderlyingAutomFamily(G), n);
end);

###############################################################################
##
#M  MultAutomAlphabet(<G>, <n>)
##
InstallOtherMethod( MultAutomAlphabet,
                    "for [IsAutomGroup and IsGroupOfAutomFamily, IsPosInt]",
                    [IsAutomGroup and IsGroupOfAutomFamily, IsPosInt],
function(G, n)
  return MultAutomAlphabet(UnderlyingAutomFamily(G), n);
end);

###############################################################################
##
#M  \= (<G>, <H>)
##
InstallMethod(\=, "for [IsAutomGroup, IsAutomGroup]",
              IsIdenticalObj, [IsAutomGroup, IsAutomGroup],
function(G, H)
  local fgens1, fgens2, fam;

  if HasIsGroupOfAutomFamily(G) and HasIsGroupOfAutomFamily(H) then
    if IsGroupOfAutomFamily(G) <> IsGroupOfAutomFamily(H) then
      Info(InfoAutomGrp, 3, "G = H: false, exactly one is GroupOfAutomFamily");
      return false;
    fi;
    if IsGroupOfAutomFamily(G) then
      Info(InfoAutomGrp, 3, "G = H: true, both are GroupOfAutomFamily");
      return true;
    fi;
  fi;

  fgens1 := List(GeneratorsOfGroup(G), g -> Word(g));
  fgens2 := List(GeneratorsOfGroup(H), g -> Word(g));
  fam := UnderlyingAutomFamily(G);

  if fam!.rws <> fail then
    fgens1 := AG_ReducedForm(fam!.rws, fgens1);
    fgens2 := AG_ReducedForm(fam!.rws, fgens2);
  fi;

  if IsEmpty(fgens1) then
    return ForAll(fgens2, IsOne);
  elif IsEmpty(fgens2) then
    return ForAll(fgens1, IsOne);
  fi;

  if GroupWithGenerators(fgens1) = GroupWithGenerators(fgens2) then
    Info(InfoAutomGrp, 3, "G = H: true, by subgroups of free group");
    return true;
  fi;

  TryNextMethod();
end);

###############################################################################
##
#M  IsSubset (<G>, <H>)
##
InstallMethod(IsSubset, "for [IsAutomGroup, IsAutomGroup]",
              IsIdenticalObj, [IsAutomGroup, IsAutomGroup],
function(G, H)
  local h, fam, fgens1, fgens2;

  if HasIsGroupOfAutomFamily(G) and IsGroupOfAutomFamily(G) then
    Info(InfoAutomGrp, 3, "IsSubgroup(G, H): true");
    Info(InfoAutomGrp, 3, "  G is GroupOfAutomFamily");
    return true;
  fi;

  fgens1 := List(GeneratorsOfGroup(G), g -> Word(g));
  fgens2 := List(GeneratorsOfGroup(H), g -> Word(g));
  fam := UnderlyingAutomFamily(G);

  if fam!.rws <> fail then
    fgens1 := AG_ReducedForm(fam!.rws, fgens1);
    fgens2 := AG_ReducedForm(fam!.rws, fgens2);
  fi;

  if IsSubgroup(GroupWithGenerators(fgens1), GroupWithGenerators(fgens2)) then
    Info(InfoAutomGrp, 3, "IsSubgroup(G, H): true");
    Info(InfoAutomGrp, 3, "  by subgroups of free group");
    return true;
  fi;

  TryNextMethod();
end);

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

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

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

  fam := UnderlyingAutomFamily(G);

  if fam!.rws <> fail then
    fgens := AG_ReducedForm(fam!.rws, fgens);
    if IsEmpty(fgens) then
      return IsOne(g);
    fi;
    w := AG_ReducedForm(fam!.rws, w);
  fi;

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

  TryNextMethod();
end);

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

  if IsTrivial(G) then
    return One(G);
  elif IsAutomatonGroup(G) then
    return Autom(Random(rs, UnderlyingFreeGroup(G)), UnderlyingAutomFamily(G));
  else
    gens := GeneratorsOfGroup(G);
    F := FreeGroup(Length(gens));
    pi := GroupHomomorphismByImagesNC(F, G,  GeneratorsOfGroup(F), gens);
    return Random(rs, F)^pi;
  fi;
end);

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

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

###############################################################################
##
#M  TrivialSubmagmaWithOne(<G>)
##
InstallMethod(TrivialSubmagmaWithOne, "for [IsAutomGroup]",
              [IsAutomGroup],
function(G)
  return Subgroup(G, [One(G)]);
end);

###############################################################################
##
#M  IsAutomatonGroup(<G>)
##
InstallImmediateMethod(IsAutomatonGroup, IsAutomGroup, 0,
function(G)
  local fam;
  fam := UnderlyingAutomFamily(G);
  return fam!.numstates = 0 or
         GeneratorsOfGroup(G) = fam!.automgens{[1..fam!.numstates]};
end);

###############################################################################
##
#M  AutomatonList(<G>)
##
InstallMethod(AutomatonList, "for [IsAutomGroup]",
              [IsAutomGroup],
function(G)
  if IsAutomatonGroup(G) then
    return AutomatonList(GroupOfAutomFamily(UnderlyingAutomFamily(G)));
  else
    Error("Group <G> is not necessarily generated by automaton,");
  fi;
end);

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

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

###############################################################################
##
#M  UnderlyingAutomaton(<G>)
##
InstallMethod(UnderlyingAutomaton, "for [IsAutomGroup]",
              [IsAutomGroup],
function(G)
  local fam;
  fam := UnderlyingAutomFamily(G);
  return MealyAutomaton(AG_AddInversesList(fam!.automatonlist){[1..fam!.numstates+1]});
end);

#E

ext: gi

max_stars_repo_name: gap-packages/automgrp

max_issues_repo_name: gap-packages/automgrp

max_forks_repo_name: gap-packages/automgrp


text:
# Copyright (c) 2018-2020, Carnegie Mellon University
# See LICENSE for details
#
# n-qubit T and Tdag transforms
#

#F qTT( n, i ) - T Gate non-terminal
#F Definition: (2^n x 2^n)-matrix  that applies an n-point T Transform to the target qubits
#F Note:       qTT(1, 1) denotes the matrix[[1, 0], [0, e^(-ipi/4)]], and qTT(1, -1) denotes the matrix[[1, 0], [0, e^(ipi/4)]]
#F             qTT(_) is symmetric.
#F
#F qTT(n , i) -> an n-qubit T transform, or Tdag transform if i = -1
Class(qTT, TaggedNonTerminal, rec(
  abbrevs   := [ (n, i) -> Checked(IsPosInt(n), [n, When(i <= -1, -1, 1)]) ],
  dims      := self >> let(size := 2^self.params[1], [size, size]),
  terminate := self >> When(self.params[2] = -1, Tensor(Replicate(self.params[1], qTdag())), Tensor(Replicate(self.params[1], qT()))),
  isReal    := self >> false,
  rChildren := self >> self.params,
  from_rChildren := (self, rch) >> self.__bases__[1](rch[1], rch[2]),
  groups    := self >> [Replicate(self.params[1], 1)],
  recursive_def := (self, arch) >> self.__bases__[1](self.params[1], self.params[2]),
  SmallRandom := () -> Random([2..5]),
  LargeRandom := () -> Random([6..15]),
  normalizedArithCost := self >> Error("ArithCost not implemented"),
  TType := T_Complex(64),
));

NewRulesFor(qTT, rec(

    # qTT_BinSplit rule
    # qTT_BinSplit qTT_(k) -> (qTT_(k1) tensor I) (I tensor qTT_(k2))
    # k1 + k2 = k
    qTT_BinSplit := rec (
        forTransposition := false,
        minSize          := 2,
        applicable       := (self, nt) >> nt.params[1] > 1,
        children         := nt -> List( [1..nt.params[1] - 1],  i -> [ Tensor(qTT(i, nt.params[2]), qTT(nt.params[1]-i, nt.params[2])).withTags(nt.getTags()) ]  ),
        apply            := (nt, c, cnt) -> c[1],
        switch := true,
    ),

    #F qTT_Base: qTT(1, _) = qT() or qT().inverse() (i.e. qTdag) SPL object
    #F Directly represent as an implementable gate
    qTT_Base := rec(
        info             := "qTT_(1, _) -> qT()",
        forTransposition := false,
        applicable       := (self, nt) >> nt.params[1]=1,
        apply            := (nt, c, cnt) -> When(nt.params[2] = -1, qTdag(), qT()),
    )

));

max_stars_repo_name: spiral-software/spiral-package-quantum

max_issues_repo_name: spiral-software/spiral-package-quantum

max_forks_repo_name: spiral-software/spiral-package-quantum


text:
# Return the time passed in last function
time;


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

Class(ABuf, AGenericTag, rec(
    __call__ := (self, bs) >> WithBases(self, rec(bs:=bs)),
    print := (self) >> When(IsBound(self.bs), Print(self.name, "(", self.bs, ")"), Print(self.name)),
#D    isBuffer := true,
    operations := Inherit(PrintOps, rec(\= := (self, other) >> ObjId(other) = ObjId(self) and (Same(other, self) or IsBound(self.bs) and IsBound(other.bs) and self.bs = other.bs)))
));

Class(AVecMem, AVec, rec(
    __call__ := (self, v) >> WithBases(self, rec(v:=v)),
    print := (self) >> When(IsBound(self.v), Print(self.name, "(", self.v, ")"), Print(self.name)),
#D    isMem := true,
    operations := Inherit(PrintOps, rec(\= := (self, other) >> ObjId(other) = ObjId(self) and (Same(other, self) or IsBound(self.v) and IsBound(other.v) and self.v = other.v)))
));

Class(AVecMemL, AVec, rec(
    __call__ := (self, v) >> WithBases(self, rec(v:=v)),
    print := (self) >> When(IsBound(self.v), Print(self.name, "(", self.v, ")"), Print(self.name)),
#D    isMemL := true,
    operations := Inherit(PrintOps, rec(\= := (self, other) >> ObjId(other) = ObjId(self) and (Same(other, self) or IsBound(self.v) and IsBound(other.v) and self.v = other.v)))
));

Class(AVecMemR, AVec, rec(
    __call__ := (self, v) >> WithBases(self, rec(v:=v)),
    print := (self) >> When(IsBound(self.v), Print(self.name, "(", self.v, ")"), Print(self.name)),
#D    isMemR := true,
    operations := Inherit(PrintOps, rec(\= := (self, other) >> ObjId(other) = ObjId(self) and (Same(other, self) or IsBound(self.v) and IsBound(other.v) and self.v = other.v)))
));


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

Class(VWrap, VWrapBase, rec(
    __call__ := (self,isa) >> Checked(IsSIMD_ISA(isa),
    WithBases(self, rec(operations:=PrintOps, isa:=isa))),

    wrap := (self,r,t,opts) >> let(
        isa := self.isa,
        v := isa.v,
        tt := When(t.isReal(),
            #paradigms.vector.breakdown.AxI_vec(TTensorI(t, v, AVec, AVec).withTags([AVecReg(isa)]), r),
            @_Base(paradigms.vector.sigmaspl.VTensor(r.node, v), r),
            paradigms.vector.breakdown.AxI_vec(
                TTensorI(TRC(t).withTags([AVecReg(isa)]), v, AVec, AVec).withTags([AVecReg(isa)]),
                paradigms.vector.breakdown.TRC_VRCLR(TRC(t).withTags([AVecReg(isa)]), r)
            )
        ),

#        Chain(Print("---> VWrap.wrap: ", r, "\n* ", t, "\n* ", tt, "\n"), DbgBreak("a"), tt)
        tt
    ),

    twrap := (self, t,opts) >> let(
        isa := self.isa,
        v := isa.v,
        tt := When(t.isReal(),
            paradigms.vector.sigmaspl.VTensor(t, v),
        #TTensorI(t, v, AVec, AVec).withTags([AVecReg(isa)]),
            TTensorI(TRC(t).withTags([AVecReg(isa)]), v, AVec, AVec).withTags([AVecReg(isa)])
        ),

#        Chain(Print("---> VWrap.twrap: ", t, "\n* ", tt, "\n"), DbgBreak("a"), tt)
        tt
    ),

    print := self >> Print(self.name, "(", self.isa, ")"),

));

Class( VWrapCpx_Mixin, rec(
    # as VWrap may change data type it must update data types
    complexToRealT := (self, t, wt) >> Cond(t.hasA("t_in") and t.hasA("t_out"),
	wt.setA(t_in => List(t.getA("t_in"), t -> t.realType()),
		t_out => List(t.getA("t_out"), t -> t.realType())),
	wt),

    opts := (self, t, opts) >> CopyFields(opts, rec(
        XType := StripList(List(Flat([opts.XType]), t -> t.realType())),
        YType := StripList(List(Flat([opts.YType]), t -> t.realType())),
    )),

));

Class(VWrapTRC, VWrapCpx_Mixin, VWrapBase, rec(
    __call__ := (self, isa) >> Checked(IsSIMD_ISA(isa),
    WithBases(self, rec(isa:=isa, operations:=PrintOps))),
    print := self >> Print(self.name, "(", self.isa, ")"),
    wrap := (self, r, t, opts) >> let(
        tt := paradigms.vector.breakdown.TRC_vect(self.complexToRealT(t, TRC(t).withTags([AVecReg(self.isa)])), r),
#        Chain(Print("---> VWrapTRC.wrap: ", r, "\n* ", t, "\n* ", tt, "\n"), tt)
        tt
    ),
    twrap := (self, t, opts) >> let(
        tt := self.complexToRealT(t, TRC(t).withTags([AVecReg(self.isa)])),
#        Chain(Print("---> VWrapTRC.twrap: ", t, "\n* ", tt, "\n"), tt)
        tt
    ),
));

Class(VWrapTRCcplx, VWrapCpx_Mixin, VWrapBase, rec(
    __call__ := (self, isa) >> Checked(IsSIMD_ISA(isa),
    WithBases(self, rec(isa:=isa, operations:=PrintOps))),
    print := self >> Print(self.name, "(", self.isa, ")"),
    wrap := (self, r, t, opts) >> paradigms.vector.breakdown.TRC_cplx(self.complexToRealT(t, TRC(t).withTags([AVecRegCx(self.isa)])), r),
    twrap := (self, t, opts) >> self.complexToRealT(t, TGrp(TRC(t).withTags([AVecRegCx(self.isa)]))),

#    wrap := (self, r, t) >> @_Base(self.twrap(t), r),
#    twrap := (self, t) >> paradigms.vector.sigmaspl.vRC(t.withTags([AVecRegCx(self.isa)]) )
));


text:
a := [2, 9, 4, 7, 5, 3];
Random(a);

max_stars_repo_name: mullikine/RosettaCodeData


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

############################################################################
##
#I  TestFlag(<n>, <i>)
##
##  tests if the i-th bit is set in binary representation the nonnegative
##  integer n
##
InstallGlobalFunction(TestFlag,
    function(n, i)
        return QuoInt(n, 2^i) mod 2 = 1;
    end);

############################################################################
##
#F  CodeCharacteristicPolynomial(<mat>, <q>)
##
##  given a matrix mat over GF(q), this computes a number which uniquely
##  identifies the characteristic polynomial of mat.
##
InstallGlobalFunction(CodeCharacteristicPolynomial,
    function(mat, q)

        local cf, z, sum, c;

        z := Z(q);
        if Length(mat) = 2 then
            cf := [ mat[1][1]*mat[2][2] - mat[1][2]*mat[2][1],
                    - mat[1][1] - mat[2][2],
                    z^0];
        elif Length(mat) = 3 then
            cf := [ - mat[1][1]*mat[2][2]*mat[3][3] - mat[1][2]*mat[2][3]*mat[3][1] - mat[1][3]*mat[2][1]*mat[3][2]
                    + mat[1][1]*mat[2][3]*mat[3][2] + mat[2][2]*mat[1][3]*mat[3][1] + mat[3][3]*mat[2][1]*mat[1][2],
                    mat[1][1]*mat[2][2] + mat[2][2]*mat[3][3] + mat[1][1]*mat[3][3]
                        - mat[1][2]*mat[2][1] - mat[2][3]*mat[3][2] - mat[1][3]*mat[3][1],
                    - mat[1][1] - mat[2][2] - mat[3][3],
                    z^0];
        else
            cf := CoefficientsOfUnivariatePolynomial(CharacteristicPolynomial(mat, 1));
        fi;

        sum := 0;
        for c in cf do
            if c = 0*z then
                sum := sum * q;
            else
                sum := sum * q + LogFFE(c, z) + 1;
            fi;
        od;

        return sum;
    end);

############################################################################
##
#F  FFMatrixByNumber(n, d, q)
##
##  computes a d x d matrix over GF(q) represented by the integer n
##
InstallGlobalFunction(FFMatrixByNumber,
    function(n, d, q)

        local z, m, i, j, k;

        z := Z(q);
        m := NullMat(d,d, GF(q));

        for i in [d, d-1..1] do
            for j in [d, d-1..1] do
                k := RemInt(n, q);
                n := QuoInt(n, q);
                if k > 0 then
                    m[i][j] := z^(k-1);
                fi;
            od;
        od;
        ConvertToMatrixRep(m, q);
        return m;
    end);

############################################################################
##
#F  CanonicalPcgsByNumber(<pcgs>, <n>)
##
##  computes the canonical pcgs wrt. pcgs represented by the integer n
##
InstallGlobalFunction(CanonicalPcgsByNumber,
    function(pcgs, n)

        local gens, cpcgs;

        gens := List(ExponentsCanonicalPcgsByNumber(RelativeOrders(pcgs), n),
            exp -> PcElementByExponents(pcgs, exp));
        cpcgs := InducedPcgsByPcSequenceNC(pcgs, gens);
        SetIsCanonicalPcgs(cpcgs, true);
        return cpcgs;
    end);

############################################################################
##
#F  OrderGroupByCanonicalPcgsByNumber(<pcgs>, <n>)
##
##  computes Order(Group(CanonicalPcgsByNumber(<pcgs>, <n>))) without
##  actually constructing the canonical pcgs or the group
##
InstallGlobalFunction(OrderGroupByCanonicalPcgsByNumber,
    function(pcgs, n)

        local ros, order, j;

        order := 1;
        ros := RelativeOrders(pcgs);
        n := RemInt(n, 2^Length(ros));
        for j in [1..Length(ros)] do
            if RemInt(n, 2) > 0 then
                order := order * ros[j];
            fi;
            n := QuoInt(n, 2);
        od;
        return order;
    end);

############################################################################
##
#F  ExponentsCanonicalPcgsByNumber(<ros>, <n>)
##
##  computes the list of exponent vectors of the
##  elements of CanonicalPcgsByNumber(<pcgs>, <n>)) without actually
##  constructing the canonical pcgs itself - it is enough to know the
##  relative orders of the pc elements
##
InstallGlobalFunction(ExponentsCanonicalPcgsByNumber,
    function(ros, n)

        local depths, len, d, exps, exp, j, cpcgs;
        depths := [];
        len := Length(ros);
        for j in [1..len] do
            d := RemInt(n, 2);
            n := QuoInt(n, 2);
            if d > 0 then
                Add(depths, j);
            fi;
        od;

        exps := [];
        for d in depths do
            exp := ListWithIdenticalEntries(len, 0);
            exp[d] := 1;
            for j in [d+1..len] do
                if not j in depths then
                    exp[j] := RemInt(n, ros[j]);
                    n := QuoInt(n, ros[j]);
                fi;
            od;
            Add(exps, exp);
        od;

        return exps;
    end);

############################################################################
##
#F  IsMatGroupOverFieldFam(famG, famF)
##
##  tests whether famG is the collections family of matrices over the field
##  whose family is famF
##
InstallGlobalFunction(IsMatGroupOverFieldFam, function(famG, famF)
    return CollectionsFamily(CollectionsFamily(famF)) = famG;
end);

############################################################################
##
#V  IRREDSOL_DATA.PRIME_POWERS
##
##  cache of proper prime powers, preset to all prime powers <= 65535
##
IRREDSOL_DATA.PRIME_POWERS :=
[ 4, 8, 9, 16, 25, 27, 32, 49, 64, 81, 121, 125, 128, 169, 243, 256, 289,
  343, 361, 512, 529, 625, 729, 841, 961, 1024, 1331, 1369, 1681, 1849, 2048,
  2187, 2197, 2209, 2401, 2809, 3125, 3481, 3721, 4096, 4489, 4913, 5041,
  5329, 6241, 6561, 6859, 6889, 7921, 8192, 9409, 10201, 10609, 11449, 11881,
  12167, 12769, 14641, 15625, 16129, 16384, 16807, 17161, 18769, 19321,
  19683, 22201, 22801, 24389, 24649, 26569, 27889, 28561, 29791, 29929,
  32041, 32761, 32768, 36481, 37249, 38809, 39601, 44521, 49729, 50653,
  51529, 52441, 54289, 57121, 58081, 59049, 63001 ];

############################################################################
##
#F  IsPPowerInt(q)
##
##  tests whether q is a positive prime power, caching new prime powers
##
InstallGlobalFunction(IsPPowerInt,
    function(q)
        if not IsPosInt(q) then
            return false;
        elif IsPrimeInt(q) then
            return true;
        elif q in IRREDSOL_DATA.PRIME_POWERS then
            return true;
        elif IsPrimePowerInt(q) then
            AddSet(IRREDSOL_DATA.PRIME_POWERS, q);
            return true;
        else
            return false;
        fi;
    end);

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

max_stars_repo_name: fingolfin/irredsol

max_issues_repo_name: fingolfin/irredsol

max_forks_repo_name: fingolfin/irredsol


text:
a := [ 8, 2, 5, 9, 1, 3, 6, 7, 4 ];
# Make a copy (with "b := a;", b and a would point to the same list)
b := ShallowCopy(a);

# Sort in place
Sort(a);
a;
# [ 1, 2, 3, 4, 5, 6, 7, 8, 9 ]

# Sort without changing the argument
SortedList(b);
# [ 1, 2, 3, 4, 5, 6, 7, 8, 9 ]
b;
# [ 8, 2, 5, 9, 1, 3, 6, 7, 4 ]


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

Class(VecRaderMid, BaseMat, SumsBase, rec(
    sums := self >> self,
    isReal := False,
    dims := self >> self.dimensions,
    #-----------------------------------------------------------------------
    new := (self, p, k, r, v) >> SPL(WithBases(self,
        rec(p := p, k := k, r := r, v := v, dimensions := [v+_roundup(p-1,v), v+_roundup(p-1,v)]))),
    #-----------------------------------------------------------------------
    area := self >> self.dimensions[1] + 2,
    transpose := self >> self,
    #-----------------------------------------------------------------------
    print := (self,i,is) >> Print(self.name, "(", self.p, ", ", self.k, ", ", self.r, ", ", self.v, ")"),
    #-----------------------------------------------------------------------
    toAMat := self >> self.term().toAMat(),
    rChildren := self >> [],
    from_rChildren := (self, rch) >> self,
    #-----------------------------------------------------------------------
    term := meth(self)
        local dl, d, e0, ep, blk, p, k, root, v;
        p := self.p; k:=self.k; root :=self.r; v:= self.v;
        dl := Concat(TRaderMid.raderDiag(p, k, root), List([p-1..v], i->V(0.0)));
        e0 := Concat([V(1.0)], List([2..v], i->V(0.0)));
        ep := Concat([V(-1/(p-1))], dl{[1..v-1]});
        blk := VBlk([[V(e0),V(e0)],[V(e0),V(ep)]], v).setSymmetric();

        if p > 2*v then
            d := VDiag(fStretch(FList(TComplex, dl{[v..p-2]}), _roundup(p-1-v, v), p-1-v), v);
            return DelayedDirectSum(blk, d);
        else
            return blk;
        fi;
    end,
    #-----------------------------------------------------------------------
    stretchCT := meth(self, N)
        local dl, d, e0, ep, blk, p, k, root, v;
        p := self.p;
        k:=self.k;
        root:=self.r;
        v:= self.v;
        dl := fStretch(FList(TComplex, Concat([-1/(p-1)], TRaderMid.raderDiag(p, k, root))), _roundup(N, v), N).tolist();
        d := VDiag(FList(TComplex, Drop(dl, v)), v);
        e0 := Concat([V(1.0)], List([2..v], i->V(0.0)));
        ep := dl{[1..v]};
        blk := VBlk([[e0,e0],[e0,ep]], v).setSymmetric();
        return DelayedDirectSum(blk, d);
    end,
    #-----------------------------------------------------------------------
    stretchPFA := meth(self, part, N, func)
        local dl, d, e0, ep, blk, p, k, root, v;
        p := self.p;
        k:=self.k;
        root:=self.r;
        v:= self.v;
#        Error();
        dl := fStretch(fCompose(FList(TComplex, Concat([-1/(p-1)], TRaderMid.raderDiag(p, k, root))), func), _roundup(N/part, v), N/part).tolist();
        d := VDiag(FList(TComplex, Drop(dl, v)), v);
        e0 := Concat([V(1.0)], List([2..v], i->V(0.0)));
        ep := dl{[1..v]};
        blk := VBlk([[e0,e0],[e0,ep]], v).setSymmetric();
        return DelayedDirectSum(blk, d);
    end,
));


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

RulesFor(DCT2, rec(
    #F DCT2_DCT2and4: 1977
    #F
    #F   DCT2_n = perm * (DCT2_n/2 dirsum DCT4_n/2^perm) * (1 tensor DFT_2)^perm
    #F
    #F   Chen/Fralick/Smith:
    #F     A Fast Computational Algorithm for the
    #F     Discrete Cosine Transform, IEEE Trans. on Comm., 1977, pp. 1004--1009.
    #F   Wang:
    #F     Reconsideration of --above--
    #F     Circ., Systems, and Signal Proc., 1983, 121--123.
    #F   Rao/Yip:
    #F     Discrete Cosine Transform, Academic Press, 1990, pp. 53
    #F
    DCT2_DCT2and4 := rec(
	info         := "DCT2'_n --> DCT2'_n/2, DCT4_n/2",
	isApplicable := P -> P[1] > 2 and P[1] mod 2 = 0,
	allChildren  := P -> [[ DCT2(P[1]/2), DCT4(P[1]/2) ]],
	rule         := (P, C) ->
	    LIJ(P[1]) *
	    DirectSum(C[1], C[2] ^ J(P[1]/2)) *
	    Tensor(I(Int((P[1] / 2))), F(2)) ^ LIJ(P[1]) # this line was blocks1(P[1])
    ),

    #F DCT2_toRDFT:
    #F
    #F   DCT2_n = blocks * RDFT_n * perm
    #F
    #F for n even
    #F
    DCT2_toRDFT := rec (
	info         := "DCT2_n -> RDFT_n",
	isApplicable := P -> IsEvenInt(P[1]),
	allChildren  := P -> [[ SRDFT(P[1]) ]],

	rule := (P, C) -> let(i := Ind(P[1]/2-1),
            # first an X shaped matrix where opposite diagonal is lowered by 1
            LIJ(P[1]).transpose() *
	    DirectSum(I(1),
		      When(i.range > 0, IterDirectSum(i, i.range, Rot(fdiv(1+i, 2*P[1]))*J(2)), []),
		      Diag(Sqrt(1/2))) *
	    C[1] * K(P[1], 2)
	)
    ),
#rdft = d os dct p7
#dct = os^-1 d^-1 rdft * p7^-1
    DCT2_toRDFT_odd := rec(
	info             := "DCT2_n -> RDFT_n",
	isApplicable     := P -> IsOddInt(P[1]),
        allChildren      := P -> [[ PRDFT1(P[1]) ]],
	forTransposition := true,
	rule := (P, C) -> let(N := P[1],
	    OS(N, 2).transpose() *
	    Diag(List([0..N-1], i -> (-1)^i)) *
	    DirectSum(Mat([[1,0]]), LIJ(N-1).transpose()) *
	    C[1] *
	    perm7(N).transpose()
	)
    ),

    #F DCT2_PrimePowerInduction:
    #F
    #F   DCT2'_n = perm * sparse * (1 tensor DCT2'_n/3) *
    #F             ( dirsum of 3x3 blocks ) ^ perm * perm
    #F
    #F   Pueschel/Moura: Discrete Cosine and Sine Transforms
    #F
    DCT2_PrimePowerInduction := rec (
	info             := "DCT2_3n -> DCT2_n",
	isApplicable := P -> P[1] mod 3 = 0 and P[1] <> 3,
	allChildren  := P -> [[ DCT2(P[1]/3) ]],

        # the sparse matrix occuring in the rule; it has the form
        # [ [ I_n, Z_n ], [ Z_n, I_n ] ], Z_n has only 1's on the upper
        # diagonal; n = L/3
	sparsemat := function ( n )
	    local L, i;
            # diagonal 1's
	    L := List([1..2*n], i -> [i, i, 1]);
            # the other 1's
	    for i in [1..n-1] do  Add(L, [i, n+1+i, 1]);   od;
	    for i in [1..n-1] do  Add(L, [n+i, i+1, 1]);   od;
	    return Sparse(L);
	end,

        # the 3x3 blocks occuring; so to say the twiddle factors
        # note: the block for i = (n-1)/2 is a DCT of size 3 and can be done
        # in 4 adds, 2 mults; but here only in 5 adds and 2 mults
	block3 := function ( i, n )
	    local M, H, M1, ii;
	    ii := i + 1/2;
	    M := TransposedMat(
		[ [ 1,                1,                       1                      ],
		  [CosPi(  ii/(3*n)), CosPi(  (2*n-ii)/(3*n)), CosPi(  (2*n+ii)/(3*n))],
		  [CosPi(2*ii/(3*n)), CosPi(2*(2*n-ii)/(3*n)), CosPi(2*(2*n+ii)/(3*n))] ]
	    );
	    M := 3 * DiagonalMat([1, 1/2, 1/2]) * M^-1;
            # now M has 1's in the first row
	    H := [[1, 1, 1], [1, -1, 0], [1, 0, -1]];
	    M1 := M*H^-1;
            # M*H has the structure 1 dirsum 2x2 block
            # get the 2x2 block
	    M1 := Sublist(List(M1, r -> Sublist(r, [2,3])), [2,3]);
	    return DirectSum(I(1), Mat(M1)) * Mat(H);
	end,

	rule := (self, P, C) >> let(n := P[1]/3,
	    L(3*n, n) *
	    DirectSum(I(n), self.sparsemat(n)) *
	    Tensor(I(3), C[1]) *
	    (DirectSum(List([1..n], i -> self.block3(i, n))) ^ L(3*n, n)) *
	    IJ(3*n, n))
    ),

#F DCT2_PrimeFactor: 1985
#F
#F   DCT2_nm = P1 * (1_(n+m-1) dirsum (1_(nm-n-m+1) tensor F_2)) *
#F             P2 * (DCT2_n tensor DCT2_m) * P3,    gcd(n, m) = 1
#F
#F Yang/Narasimha: Prime Factor Decomposition of the Discrete Cosine Transform,
#F   Proc. ICASSP, pp. 772--775, 1985
#F
#F see also:
#F Feig/Linzer: Scaled DCT's on Input Sizes that Are Composite, IEEE Transactions on
#F   Signal Processing 43(1), pp. 43--50, 1995
#F
DCT2_PrimeFactor := rec (
  info             := "DCT2_nm -> DCT2_n tensor DCT2_m",
  isApplicable     := P -> not IsPrimePowerInt(P[1]),
  allChildren      := P -> List(DivisorPairsRP(P[1]), p -> [ DCT2(p[1]), DCT2(p[2]) ]),

  rule := function ( P, C )
    local n, n1, n2, dctperm1, dctperm2;

    n1 := C[1].dimensions[1];
    n2 := C[2].dimensions[1];
    n  := P[1];

    # first permutation
    dctperm1 := function ( n1, n2 )
      local n, g1, g2, L, i, j;

      n  := n1*n2;
      g1 := QuotientMod(1, n2, n1);
      g2 := n2 - QuotientMod(1, n1, n2);

      # indices
      L := [ ];
      for i in [0..n1 - 1] do
	for j in [0..n2 - 1] do
	  if (i + j) mod 2 = 0 then
	    Add(L, ((2*i + 1)*g1*n2 - (2*j + 1)*g2*n1) mod (4*n));
	  else
	    Add(L, ((2*i + 1)*g1*n2 + (2*j + 1)*g2*n1) mod (4*n));
	  fi;
	od;
      od;

      for i in [1..Length(L)] do
	if L[i] < 2*n then
	  L[i] := (L[i] - 1)/2;
	else
	  L[i] := (4*n - L[i] - 1)/2;
	fi;
      od;

      return PermList(L + 1);
    end;

    # second and third permutation
    dctperm2 := function ( n1, n2 )
      local n, K, L, i, j, q2;

      n  := n1*n2;

      L := [ ];
      K := [ ];

      # start with the border
      # i = 0
      for j in [0..n2 - 1] do
	Add(L, j);
	Add(K, j*n1);
      od;

      # j = 0
      for i in [1..n1 - 1] do
	Add(L, i*n2);
	Add(K, i*n2);
      od;

      # now the interior; two cases
      if n1 mod 2 = 0 then
	for i in [1..n1/2 - 1] do
	  for j in [1..n2 - 1] do
	    q2 := i*n2 + j*n1;
	    if q2 > n then
	      q2 := q2 - 2*n;
	      Add(L, (n1 - i)*n2 + (n2 - j));
	      Add(L, i*n2 + j);
	    else
	      Add(L, i*n2 + j);
	      Add(L, (n1 - i)*n2 + (n2 - j));
	    fi;
	    Add(K, AbsInt(i*n2 - j*n1));
	    Add(K, AbsInt(q2));
	  od;
	od;
	i := n1/2;
	for j in [1..(n2 - 1)/2] do
	  q2 := i*n2 + j*n1;
	  if q2 > n then
	    q2 := q2 - 2*n;
	    Add(L, (n1 - i)*n2 + (n2 - j));
	    Add(L, i*n2 + j);
	  else
	    Add(L, i*n2 + j);
	    Add(L, (n1 - i)*n2 + (n2 - j));
	  fi;
	  Add(K, AbsInt(i*n2 - j*n1));
	  Add(K, AbsInt(q2));
	od;
      else # n1 mod 2 <> 0
	for i in [1..(n1 - 1)/2] do
	  for j in [1..n2 - 1] do
	    q2 := i*n2 + j*n1;
	    if q2 > n then
	      q2 := q2 - 2*n;
	      Add(L, (n1 - i)*n2 + (n2 - j));
	      Add(L, i*n2 + j);
	    else
	      Add(L, i*n2 + j);
	      Add(L, (n1 - i)*n2 + (n2 - j));
	    fi;
	    Add(K, AbsInt(i*n2 - j*n1));
	    Add(K, AbsInt(q2));
	  od;
	od;
      fi;

      return [PermList(L + 1), PermList(K + 1)];
    end;

    # actual rule
    return
      Perm(dctperm2(n1, n2)[2]^-1, n) *
      DirectSum(
	I(n1 + n2 - 1),
	Tensor(I((n - n1 - n2 + 1)/2), F(2))
      ) *
      Perm(dctperm2(n1, n2)[1], n) *
      Tensor(C) *
      Perm(dctperm1(n1, n2), n);
  end
)
));


text:
#############################################################################
####
##
#W  anupqopt.gi                ANUPQ package                    Werner Nickel
#W                                                                Greg Gamble
##
##  Install file for functions to do with option manipulation.
##
#Y  Copyright (C) 2001  Lehrstuhl D fuer Mathematik,  RWTH Aachen,  Germany
##

#############################################################################
##
#V  PQ_FUNCTION . . . . . . . . . internal functions called by user functions
##
##  A record whose fields are (function)  names  and  whose  values  are  the
##  internal functions called by the functions with those names.
##
InstallValue( PQ_FUNCTION,
              rec( Pq                   := PQ_EPI_OR_PCOVER,
                   PqDescendants        := PQ_DESCENDANTS,
                   StandardPresentation := PQ_EPIMORPHISM_STANDARD_PRESENTATION,
                   PqDescendantsTreeCoclassOne := PqDescendantsTreeCoclassOne
                  )
             );

#############################################################################
##
#V  ANUPQoptions  . . . . . . . . . . . . . . . . . . . .  admissible options
##
##  is a record of lists of names of admissible {\ANUPQ} options,  such  that
##  each field is either the name of a (``key'') {\ANUPQ}  function  and  the
##  corresponding value is the list of option names that are  admissible  for
##  the function.
##
InstallValue( ANUPQoptions,
              rec( # options for `Pq' and `PqEpimorphism'
                   Pq  := [ "Prime",
                            "ClassBound",
                            "Exponent",
                            "Metabelian",
                            "OutputLevel",
                            "Relators",
                            "Identities",
                            "GroupName",
                            "SetupFile",
                            "PqWorkspace",
                            "RedoPcp" ],

                   # options for `PqDescendants'
                   PqDescendants
                       := [ "ClassBound",
                            "OrderBound",
                            "Relators",
                            "GroupName",
                            "StepSize",
                            "PcgsAutomorphisms",
                            "RankInitialSegmentSubgroups",
                            "SpaceEfficient",
                            "CapableDescendants",
                            "AllDescendants",
                            "Exponent",
                            "Metabelian",
                            "SubList",
                            "BasicAlgorithm",
                            "CustomiseOutput",
                            "SetupFile",
                            "PqWorkspace" ],

                   # options for `[Epimorphism][Pq]StandardPresentation'
                   StandardPresentation
                       := [ "Prime",
                            "pQuotient",
                            "ClassBound",
                            "Relators",
                            "GroupName",
                            "PcgsAutomorphisms",
                            "Exponent",
                            "Metabelian",
                            "OutputLevel",
                            "StandardPresentationFile",
                            "SetupFile",
                            "PqWorkspace" ],

                   # options for `PqDescendantsTreeCoclassOne'
                   PqDescendantsTreeCoclassOne
                       := [ "ClassBound",
                            "OrderBound",
                            "Relators",
                            "GroupName",
                            "StepSize",
                            "PcgsAutomorphisms",
                            "RankInitialSegmentSubgroups",
                            "SpaceEfficient",
                            "CapableDescendants",
                            "AllDescendants",
                            "Exponent",
                            "Metabelian",
                            "SubList",
                            "BasicAlgorithm",
                            "CustomiseOutput",
                            "TreeDepth",
                            "SetupFile",
                            "PqWorkspace" ],

                   PqList
                       := [ "SubList" ],
                   PqPcPresentation
                       := [ "Prime",
                            "ClassBound",
                            "Exponent",
                            "Metabelian",
                            "OutputLevel",
                            "Relators",
                            "Identities",
                            "GroupName" ],
                   PqNextClass
                       := [ "QueueFactor" ],
                   PqEvaluateIdentities
                       := [ "Identities" ],
                   PqDoExponentChecks
                       := [ "Bounds" ],
                   PqDisplayStructure
                       := [ "Bounds" ],
                   PqDisplayAutomorphisms
                       := [ "Bounds" ],
                   PqSPComputePcpAndPCover
                       := [ "Prime",
                            "ClassBound",
                            "Exponent",
                            "Metabelian",
                            "OutputLevel",
                            "Relators",
                            "GroupName" ],
                   PqSPSavePresentation
                       := [ "ClassBound",
                            "PcgsAutomorphisms",
                            "StandardPresentationFile" ],
                   PqPGSetDescendantToPcp
                       := [ "Filename" ],
                   PqPGSupplyAutomorphisms
                       := [ "NumberOfSolubleAutomorphisms",
                            "RelativeOrders" ],
                   PqPGConstructDescendants
                       := [ "ClassBound",
                            "OrderBound",
                            "StepSize",
                            "PcgsAutomorphisms",
                            "RankInitialSegmentSubgroups",
                            "SpaceEfficient",
                            "CapableDescendants",
                            "AllDescendants",
                            "Exponent",
                            "Metabelian",
                            "BasicAlgorithm",
                            "CustomiseOutput" ],
                   PqAPGDegree
                       := [ "Exponent" ],
                   PqAPGPermutations
                       := [ "PcgsAutomorphisms",
                            "SpaceEfficient",
                            "PrintAutomorphisms",
                            "PrintPermutations" ],
                   PqAPGOrbits
                       := [ "PcgsAutomorphisms",
                            "SpaceEfficient",
                            "CustomiseOutput" ],
                   PqAPGOrbitRepresentatives
                       := [ "PcgsAutomorphisms",
                            "SpaceEfficient",
                            "CapableDescendants",
                            "AllDescendants",
                            "Exponent",
                            "Metabelian",
                            "CustomiseOutput",
                            "Filename" ],
                   PqAPGSingleStage
                       := [ "StepSize",
                            "PcgsAutomorphisms",
                            "RankInitialSegmentSubgroups",
                            "SpaceEfficient",
                            "CapableDescendants",
                            "AllDescendants",
                            "Exponent",
                            "Metabelian",
                            "BasicAlgorithm",
                            "CustomiseOutput" ]
                  )
             );

#############################################################################
##
#F  AllANUPQoptions() . . . . . . . .  lists all options of the ANUPQ package
##
##  lists all the {\GAP}  options  defined  for  functions  of  the  {\ANUPQ}
##  package.
##
InstallGlobalFunction( AllANUPQoptions, function()
  return Set( Concatenation(
                  List( RecNames(ANUPQoptions), fld -> ANUPQoptions.(fld) )
                  ) );
end );

#############################################################################
##
#V  ANUPQGlobalOptions . . . . .  options that can be set globally by PqStart
##
##  A list of the options that `PqStart' can set and thereby  make  available
##  to any function  interacting  with  the  {\ANUPQ}  process  initiated  by
##  `PqStart'.
##
InstallValue( ANUPQGlobalOptions, [ "Prime", "Exponent", "Relators" ] );

#############################################################################
##
#V  ANUPQoptionChecks . . . . . . . . . . . the checks for admissible options
##
##  A record whose fields are the names of admissible ANUPQ options and whose
##  values are one-argument functions that return `true' when given  a  value
##  that is a valid value for the option, and `false' otherwise.
##
InstallValue( ANUPQoptionChecks,
              rec( Prime := x -> IsInt(x) and IsPrimeInt(x),
                   pQuotient  := IsPcGroup and IsPGroup,
                   ClassBound := IsPosInt,
                   OrderBound := IsPosInt,
                   Exponent   := IsPosInt,
                   Metabelian := IsBool,
                   GroupName  := IsString,
                   Identities := x -> IsList(x) and ForAll(x, IsFunction),
                   OutputLevel := x -> x in [0..3],
                   Relators   := x -> IsList(x) and ForAll(x, IsString),
                   StandardPresentationFile := IsString,
                   SetupFile  := IsString,
                   PqWorkspace := IsPosInt,
                   StepSize := x -> IsPosInt(x) or
                                    (IsList(x) and ForAll(x, IsPosInt)),
                   PcgsAutomorphisms := IsBool,
                   BasicAlgorithm := IsBool,
                   RankInitialSegmentSubgroups := x -> x = 0 or IsPosInt(x),
                   SpaceEfficient := IsBool,
                   CapableDescendants := IsBool,
                   AllDescendants := IsBool,
                   SubList := x -> IsPosInt(x) or
                                   (IsSet(x) and ForAll(x, IsInt)
                                    and IsPosInt(x[1])),
                   CustomiseOutput := IsRecord,
                   Bounds := x -> IsSet(x) and 2 = Length(x) and
                                  ForAll(x, IsPosInt),
                   QueueFactor := IsPosInt,
                   RedoPcp := IsBool,
                   PrintAutomorphisms := IsBool,
                   PrintPermutations := IsBool,
                   NumberOfSolubleAutomorphisms := x -> x = 0 or IsPosInt(x),
                   RelativeOrders := x -> IsList(x) and ForAll(x, IsPosInt),
                   Filename := IsString,
                   TreeDepth := IsPosInt
                   )
             );

#############################################################################
##
#V  ANUPQoptionTypes . . . . . .  the types (in words) for admissible options
##
##  A record whose fields are the names of admissible ANUPQ options and whose
##  values are valid types of the options in plain words.
##
InstallValue( ANUPQoptionTypes,
              rec( Prime := "prime integer",
                   pQuotient  := "pc p-group",
                   ClassBound := "positive integer",
                   OrderBound := "positive integer",
                   Exponent   := "positive integer",
                   Metabelian := "boolean",
                   GroupName  := "string",
                   Identities := "list of functions",
                   OutputLevel := "integer in [0..3]",
                   Relators   := "list of strings",
                   StandardPresentationFile := "string",
                   SetupFile  := "string",
                   PqWorkspace := "positive integer",
                   StepSize := "positive integer or positive integer list",
                   PcgsAutomorphisms := "boolean",
                   BasicAlgorithm := "boolean",
                   RankInitialSegmentSubgroups := "nonnegative integer",
                   SpaceEfficient := "boolean",
                   CapableDescendants := "boolean",
                   AllDescendants := "boolean",
                   SubList
                       := "pos've integer or increasing pos've integer list",
                   CustomiseOutput := "record",
                   Bounds := "pair of increasing positive integers",
                   QueueFactor := "positive integer",
                   RedoPcp := "boolean",
                   PrintAutomorphisms := "boolean",
                   PrintPermutations := "boolean",
                   NumberOfSolubleAutomorphisms := "nonnegative integer",
                   RelativeOrders := "list of positive integers",
                   Filename := "string",
                   TreeDepth := "positive integer"
                   )
             );

#############################################################################
##
#F  PQ_OTHER_OPTS_CHK( <funcname>, <interactive> ) . check opts belong to f'n
##
##  checks the `OptionsStack'  only  has  recognised  options  for  (generic)
##  function <funcname> and if not and if  `ANUPQWarnOfOtherOptions  =  true'
##  (see~"ANUPQWarnOfOtherOptions") `Info's  the  non-<funcname>  options  at
##  `InfoANUPQ' level 1.
##
##  The argument <interactive> is only relevant for those functions that have
##  both an interactive and non-interactive form, namely those with fields in
##  `PQ_FUNCTION', for which some options need to be excluded.
##
InstallGlobalFunction(PQ_OTHER_OPTS_CHK, function(funcname, interactive)
local optnames, excopts, generic, interactivestr;
  if ANUPQWarnOfOtherOptions and ValueOption("recursive") = fail and
     not IsEmpty(OptionsStack) then
    excopts := [];
    if funcname in RecNames(PQ_FUNCTION) then
      if interactive then
        excopts := ["PqWorkspace", "SetupFile"];
        interactivestr := "interactive";
      else
        interactivestr := "non-interactive";
        if funcname = "Pq" then
          excopts  := ["RedoPcp"];
        fi;
      fi;
      generic := "generic ";
    else
      generic := "";
    fi;
    optnames := Difference( RecNames( OptionsStack[ Length(OptionsStack) ] ),
                            Difference( ANUPQoptions.(funcname), excopts ) );
    if funcname = "Pq" then
      optnames := Difference( optnames, ["PqEpiOrPCover"] );
    fi;
    if not IsEmpty(optnames) then
      Info( InfoANUPQ + InfoWarning, 1,
            "ANUPQ Warning: Options: ", optnames, " ignored" );
      if IsSubset(excopts, optnames) then
        Info( InfoANUPQ + InfoWarning, 1,
              "(invalid for ", interactivestr, " call of generic function: `",
              funcname, "')." );
      else
        Info( InfoANUPQ + InfoWarning, 1,
              "(invalid for ", generic, "function: `", funcname, "')." );
      fi;
    fi;
  fi;
end);

#############################################################################
##
#F  VALUE_PQ_OPTION( <optname> ) . . . . . . . . . enhancement of ValueOption
#F  VALUE_PQ_OPTION( <optname>, <defaultval> )
#F  VALUE_PQ_OPTION( <optname>, <datarec> )
#F  VALUE_PQ_OPTION( <optname>, <defaultval>, <datarec> )
##
##  If the value <optval> of <optname> is not `fail' and it is  an  ok  value
##  for <optname> then <optval> is returned; if <optval> is not an  ok  value
##  an error is signalled. If <optval> is `fail' and <datarec> is  given  and
##  <datarec>.(<optname>) is already  bound  then  that  value  is  returned;
##  otherwise, if  <optval>  is  `fail'  and  a  default  value  <defaultval>
##  different  from  `fail'  is  supplied  then  <defaultval>  is   returned.
##  Supplying a <defaultval> of `fail' is special; it indicates  that  option
##  <optname> must have a value i.e. <optval> is not allowed to be `fail' and
##  if it is an error is signalled. If  a  <datarec>  argument  is  supplied,
##  which must be a record, then the return value, if not `fail' and a  legal
##  value, is also stored in `<datarec>.(<optname>)'.
##
##  *Note:* <defaultval> cannot be a record.
##
InstallGlobalFunction(VALUE_PQ_OPTION, function(arg)
local optname, optval, len;
  optname := arg[1];
  optval := ValueOption(optname);
  len := Length(arg);
  if optval = fail then
    if 1 = len then
      return optval;
    elif IsRecord( arg[len] ) and IsBound( arg[len].(optname) ) then
      # return the previously recorded value
      return arg[len].(optname);
    elif not IsRecord(arg[2]) then
      if arg[2] = fail then
        Error("you must supply a value for option: \"", optname, "\"\n");
      fi;
      optval := arg[2];
    fi;
  elif not ANUPQoptionChecks.(optname)(optval) then
    Error("\"", optname, "\" value must be a ",
          ANUPQoptionTypes.(optname), "\n");
  fi;
  if (optval <> fail) and (2 <= len) and IsRecord(arg[len]) then
    arg[len].(optname) := optval;
  fi;
  return optval;
end);

#############################################################################
##
#F  PQ_OPTION_CHECK(<basefn>,<datarec>) . check optns present/setable if nec.
##
##  If `<basefn> = "Pq"' (i.e. this check is  carried  out  if  the  function
##  called is `Pq', `PqEpimorphism' or  `PqPCover')  check  that  the  option
##  `Prime' has been  passed  or,  in  a  special  `PqPCover'  case,  can  be
##  determined from the `<datarec>.group'  which  must  be  present.  In  the
##  special `PqPCover' case, the options `Prime' and `ClassBound'  determined
##  are saved in <datarec>. If `Prime' is not supplied in the cases where  it
##  needs to be, an error is emitted.
##
InstallGlobalFunction(PQ_OPTION_CHECK, function(basefn, datarec)
local optname, out;
  if basefn = "Pq" then
    if datarec.calltype = "interactive" then
      if VALUE_PQ_OPTION("RedoPcp", false) then
        PQ_UNBIND(datarec, ["Prime",  "ClassBound", "Exponent", "Metabelian",
                            "pCover", "pQuotient",  "pQepi"] );
      fi;
    fi;
    if ValueOption("PqEpiOrPCover") = "pCover" and
       HasIsPGroup(datarec.group) and IsPGroup(datarec.group) then
      if VALUE_PQ_OPTION("Prime", datarec) = fail then
        if not HasPrimePGroup(datarec.group) then
          Error( "supplied group is not known to be a p-group or p unknown.\n",
                 "Option `Prime' must be supplied" );
        else
          datarec.Prime := PrimePGroup(datarec.group);
        fi;
      fi;
      if VALUE_PQ_OPTION("ClassBound", datarec) = fail and
         HasPClassPGroup(datarec.group) then
        datarec.ClassBound := PClassPGroup(datarec.group);
      fi;
    else
      VALUE_PQ_OPTION("Prime", fail, datarec);
    fi;
    VALUE_PQ_OPTION("ClassBound", 63, datarec);
  elif basefn = "StandardPresentation" then
    if VALUE_PQ_OPTION("Prime", datarec) = fail and
       VALUE_PQ_OPTION("pQuotient", datarec) = fail then
      if IsPcGroup(datarec.group) and IsPGroup(datarec.group) then
        datarec.Prime := PrimePGroup(datarec.group);
      else
        Error( "since group of process is not a pc p-group, a prime or\n",
               "p-quotient (pc group) of the group of the process ",
               "must be supplied\n" );
      fi;
    fi;
  fi;
end);

#############################################################################
##
#F  PQ_CUSTOMISE_OUTPUT(<datarec>, <subopt>, <suboptstring>, <suppstrings>)
##
##  writes the required output to the `pq' binary for the sub-option <subopt>
##  of  the  option  `CustomiseOutput',  the  value  of  that  option  having
##  previously been stored in `<datarec>.des.CustomiseOutput'; <suboptstring>
##  is part of the comment written to the `pq' binary for the sub-option  and
##  <suppstrings> is a list of such comments for the supplementary  questions
##  asked by the `pq' binary for the sub-option <subopt>.
##
InstallGlobalFunction( PQ_CUSTOMISE_OUTPUT,
function(datarec, subopt, suboptstring, suppstrings)
local optrec, isOptionSet, i;
  optrec := datarec.des.CustomiseOutput;
  if IsEmpty(suppstrings) then
    isOptionSet := IsBound( optrec.(subopt) ) and optrec.(subopt) in [1, true];
    ToPQ_BOOL(datarec, isOptionSet, suboptstring);
  elif IsBound( optrec.(subopt) ) and IsList( optrec.(subopt) ) then
    ToPQ(datarec, [ 0 ], [ "  #customise ", suboptstring ]);
    for i in [1 .. Length(suppstrings)] do
      isOptionSet := IsBound( optrec.(subopt)[i] ) and
                     optrec.(subopt)[i] in [1, true];
      ToPQ_BOOL(datarec, isOptionSet, suppstrings[i]);
    od;
  else
    ToPQ(datarec, [ 1 ], [ "  #default ", suboptstring ]);
  fi;
end);

#############################################################################
##
#F  PQ_APG_CUSTOM_OUTPUT(<datarec>, <subopt>, <suboptstring>, <suppstrings>)
##
##  writes the required output to the `pq' binary for the sub-option <subopt>
##  of the option `CustomiseOutput',  as  required  by  an  Advanced  p-Group
##  Generation Menu item, the value of that  option  having  previously  been
##  stored in `<datarec>.des.CustomiseOutput'; <suboptstring> is part of  the
##  comment written to the `pq' binary for the sub-option  and  <suppstrings>
##  is a list of such comments for the supplementary questions asked  by  the
##  `pq' binary for the sub-option <subopt>.
##
InstallGlobalFunction( PQ_APG_CUSTOM_OUTPUT,
function(datarec, subopt, suboptstring, suppstrings)
local optrec, optlist, isOptionSet, i;
  optrec := datarec.des.CustomiseOutput;
  if not( IsRecord(optrec) and IsBound( optrec.(subopt) ) and
          IsList( optrec.(subopt) ) ) then
    optlist := [];
    datarec.des.CustomiseOutput.(subopt) := optlist;
  else
    optlist := optrec.(subopt);
  fi;
  for i in [1 .. Length(suppstrings)] do
    isOptionSet := IsBound( optlist[i] ) and optlist[i] in [1, true];
    ToPQ_BOOL(datarec, isOptionSet, suppstrings[i]);
  od;
end);

#############################################################################
##
#F  SET_ANUPQ_OPTIONS( <funcname>, <fnname> )  . set options from OptionStack
##
##  When called by a function with name  <funcname>  sets  the  options  from
##  `OptionsStack'    checking    that    they    are     a     subset     of
##  `ANUPQoptions.<fnname>'. Both <funcname> and <fnname> should be strings.
##
InstallGlobalFunction( SET_ANUPQ_OPTIONS, function( funcname, fnname )
    local optrec, optnames, opt;

    # there should be options
    if IsEmpty( OptionsStack ) then
        # no options??
        optrec := rec();
        Info( InfoANUPQ, 1, funcname, " called with no options!" );
    else
        optrec := ShallowCopy( OptionsStack[ Length( OptionsStack ) ] );
        optnames := Set( REC_NAMES(optrec) );
        SubtractSet( optnames, Set( ANUPQoptions.(fnname) ) );
        Info( InfoANUPQ, 2, funcname, " called with options: ",
                            OptionsStack[ Length( OptionsStack ) ] );
        if 0 < Length(optnames) then
            # it's not an error to have unknown options,
            # function may have been called recursively and the
            # options may be intended for some other function
            Info( InfoWarning + InfoANUPQ, 2,
                  funcname, " called with unknown options: ", optnames);
        fi;
        for opt in optnames do
            Unbind( optrec.(opt) );
        od;
    fi;
    return optrec;
end );

#############################################################################
##
#F  ANUPQoptError( <funcname>, <illegal> )  . . . . . create an error message
##
##  creates an error message  for  the  function  with  name  <funcname>.  If
##  <illegal> is a string it is taken to be  the  first  line  of  the  error
##  message. Otherwise <illegal> should be alist of illegal options (strings)
##  found. The error message (string) returned also gives the list  of  valid
##  options together with the value types expected for function <funcname>.
##
InstallGlobalFunction( ANUPQoptError, function( funcname, illegal )
    local Optstring, Valstring, errmsg, optname;

    Optstring := optname -> Concatenation("\"", optname, "\"");
    Valstring := optval  -> Concatenation("<",  optval,  ">");

    if IsString(illegal) then
        errmsg := illegal;
    else # IsList(illegal)
        errmsg := Concatenation("Illegal ", funcname, " option");
        if Length(illegal) > 1 then
            Append(errmsg, "s");
        fi;
        Append(errmsg, ": ");
        Append(errmsg, JoinStringsWithSeparator( List(illegal, Optstring) ));
    fi;
    Append(errmsg, Concatenation(".\nValid ", funcname, " options:\n"));
    for optname in ANUPQoptions.(funcname) do
        Append(errmsg, "    ");
        Append(errmsg, Optstring(optname));
        if ANUPQoptionChecks.(optname) <> IsBool then
            Append(errmsg, ", ");
            Append(errmsg, Valstring( ANUPQoptionTypes.(optname) ));
        fi;
        Append(errmsg, "\n");
    od;
    return errmsg;
end );

#############################################################################
##
#F  ANUPQextractOptions( <funcname>, <args> ) . . . . . . . . extract options
##
##  extracts options from  <args>  for  function  with  name  <funcname>  and
##  returns a record suitable for use with `PushOptions'.  Abbreviations  are
##  allowed for option names so long as  each  abbreviates  a  unique  option
##  name.
##
InstallGlobalFunction( ANUPQextractOptions, function(funcname, args)
    local   Match, error, optrec, i, optname;

    # allow to give only a prefix
    Match := function( argi )
        local matches;
        if IsString(argi) then
            matches := Filtered(ANUPQoptions.(funcname),
                                optname -> 0 < Length(argi) and
                                           Length(argi) <= Length(optname) and
                                           optname{[1..Length(argi)]} = argi);
            if 1 = Length(matches) then
                return matches[1];
            fi;
            error := Concatenation( "argument: \"", argi,
                                    "\" doesn't abbreviate a unique option");
        else
            error := Concatenation( "argument: ", String(argi),
                                    " is not a (non-null) string (option name)"
                                    );
        fi;
        return fail;
    end;

    # extract options from args
    optrec := rec();
    i := 1;
    while i <= Length(args)  do
        optname := Match( args[i] );
        if optname = fail then
            Error( ANUPQoptError( funcname, error ) );
        elif ANUPQoptionChecks.(optname) = IsBool then
            optrec.(optname) := true;
            i := i + 1;
        elif i = Length(args) then
            # all remaining options are non-boolean and expect a value to
            # follow
            Error( ANUPQoptError(
                       funcname,
                       Concatenation( "Expected value for option: ", args[i] )
                       ) );
        else
            # checking values are ok is done later
            optrec.(optname) := args[i + 1];
            i := i + 2;
        fi;
    od;
    return optrec;

end );

#E  anupqopt.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:
# Copyright (c) 2018-2021, Carnegie Mellon University
# See LICENSE for details

#F ==========================================================================
#F ParCell(<num_spus>, <pkSize>) - Cell parallelization tag
#F    .params[1] is num_spus
#F    .params[2] is packetSize
#F
Class(ParCell, AGenericTag, rec(isCell := true));

#F ==========================================================================
#F ParCell_auto (<num_spus>) - Cell parallelization tag
#F    .params[1] is num_spus
#F Highest possible packet size is automatically chosen
Class(ParCell_auto, AGenericTag, rec(isCell := true));

#F ==========================================================================
#F ParDMPCell_old(<num_spus>) - Cell DMP parallelization tag
#F    .params[1] is num_spus
Class(ParCellDMP_old, AGenericTag, rec(isCell := true));

#F ==========================================================================
#F ParDMPCell(<num_spus>) - Cell DMP parallelization tag
#F    .params[1] is num_spus
#F    .params[2] is v
Class(ParCellDMP, AGenericTag, rec(isCell := true));

Class(StickyL, AGenericTag);


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

Declare(DPBench);

#F DPBench(rec(experiment1 := opts1, ...), <dpopts>)
#F
#F d := DPBench(rec(default := SpiralDefaults), rec())
#F
#F DPBench Interface:
#F
#F     .build(transforms, opts, bopts, name, dpopts),
#F
#F     .generateCfiles := [true, false],  default = true
#F     .matrixVerify   := [true, false],  default = false
#F     .fftwVerify     := [true, false],  default = false
#F     .quickVerify     := [true, false],  default = false
#F
#F     .fileTransform(exp, t, opts),    .c filenames:   override in a subclass
#F     .funcTransform(exp, t, opts),    function names: override in a subclass
#F     .txtFileName(exp, runMethod),    timing file name
#F
#F     .runAll()                        run all transforms from every experiment's opts.benchTransforms
#F     .resumeAll()                     same as .runAll() but try to reload hash from disk
#F
#F     .run(transforms)                 run given transforms in all experiments
#F     .resume(transforms)              same as .run() but try to relaod hash from disk
#F
#F      .runRandomAll(),  .runRandom()    run a random ruletree
#F      .runExhaustiveAll(), .runExhaustive()   run all ruletrees
#F
#F     .generateCode(transforms, exp)
#F     .generateProductionCode(transforms, exp)   same as .generateCode() but will use opts.production()
#F
#F     .entries(transforms, exp)
#F     .times(transforms, exp)
#F     .alltimes(transforms)
#F     .speedup(transforms, baselineExp)
#F
#F     .flopcyc(transforms)      # FLoating point Operations Per Cycle
#F     .mflops(transforms, mhz)
#F     .scaledTimes(transforms, scaleFunc)
#F
#F      getResult(i)                    returns "rec(t, opts, rt)" for the <i>th experiment run by runAll()
Class(DPBench, rec(
    ##
    ## Private methods
    ##
    _checkExp := exp ->
        Cond(not IsRec(exp),        Error("Experiments record <exp> must be a record"),
         NumRecFields(exp) < 1, Error("Experiments record is empty"),
         not ForAll(UserRecFields(exp), f -> IsRec(exp.(f))),
                 Error("Each entry must be a valid Spiral options record (ie. SpiralDefaults)"),
         exp),

    _startHashFile := (self, hfile, exp, d) >> PrintTo(hfile,
        "<# DPBench experiment '", exp, "'\n",
        " # Started ", d[2], " ", d[3], " ", d[1], "  ", d[4], ":", d[5], "\n",
        " # Transforms: ", Cond(IsBound(self.transforms), self.transforms, ""), "#> \n\n",
        "ImportAll(spiral); Import(paradigms.common, paradigms.smp, platforms.sse, platforms.avx, paradigms.vector, nontransforms.ol); \n",
        "ImportAll(platforms.scalar); \n",
        "ImportAll(paradigms.vector); \n\n",
        "hash := HashTableDP(); \n"
    ),

    _loadHash := meth(self, hfile, opts)
        local b, bkdowns, ns, result;
	# Create a package that contains mappings from breakdown names to the corresponding
	# objects. This is needed because global names map to global breakdown rule objects
	# which might have different settings from the ones in opts.breakdownRules
	# Note: package is essentially a namespace that is always on top of imports
	#       so imports within a hash file will be superceded by it
	bkdowns := ConcatList(UserRecFields(opts.breakdownRules), f->opts.breakdownRules.(f));
        ns := tab();
	for b in bkdowns do ns.(b.name) := b; od;

        result := READ(hfile, ns);
        if result = false or not IsBound(ns.hash) then return false;
        else return ns.hash;
        fi;
    end,

    _saveHash := meth(self, hfile, exp, date, hash)
        local bucket, e;
        var.print := var.printFull;
        self._startHashFile(hfile, exp, date);
        for bucket in hash.entries do
            for e in bucket do
                if e.data<>[] then
                    AppendTo(hfile, "HashAdd(hash, ", e.key, ", [", e.data[1], "]);\n");
                fi;
            od;
        od;
        var.print := var.printShort;
    end,

    # merging two hash files by taking fastest entries, returns resulting hash
    _mergeHashes := meth(self, src_file1, src_file2, opts)
        local h1, h2;
        h1 := self._loadHash(src_file1, opts);
        h2 := self._loadHash(src_file2, opts);
        Checked(h1<>false and h2<>false,
            HashWalk(h1, function(key, data)
                local d;
                d := HashLookup(h2, key);
                if d=false or (d[1].measured>data[1].measured and data[1].measured>0) then
                    HashAdd(h2, key, data);
                fi;
            end));
        return h2;
    end,

    _reloadAllHashes := meth(self)
       local hash, exp, e;
       for e in UserRecFields(self.exp) do
           exp := self.exp.(e);
           hash := self._loadHash(exp.hashFile, exp);
           if (self.verbosity>0) then
              PrintLine(When(hash=false, "Could not load ", "Loaded "), e, " (", exp.hashFile, ")");
           fi;
           if hash <> false then
               exp.hashTable := hash;
           fi;
       od;
    end,

    _generateCode := meth(self, transforms, exp, opts)
         local entries, e, c, r, t;
         for t in transforms do
             e := self.entries([HashAsSPL(t)], exp)[1];
             if e = false then Error("Transform ", t, " not found in hashTable for experiment '", exp, "'"); fi;
             r := ApplyRuleTreeSPL(e.ruletree, t, opts);
             c := CodeRuleTree(r, opts);
             PrintLine(t, " -> ", self.prodFileTransform(exp, t, opts));
             PrintTo(self.prodFileTransform(exp, t, opts), PrintCode(self.prodFuncTransform(exp, t, opts), c, opts));
         od;
    end,

    ##
    ## Public methods
    ##
    __call__ := meth(self, experiments, dpopts)
        local e, exp;
        self._checkExp(experiments);
        exp:=rec();
        for e in UserRecFields(experiments) do
           exp.(e) := CopyFields(experiments.(e));
           if not IsBound(exp.(e).hashTable) then
               exp.(e).hashTable := HashTableDP(); fi;
           if not IsBound(exp.(e).hashFile) then
               exp.(e).hashFile := Concat(e, ".hash"); fi;
        od;

        return WithBases(self,
            rec(dpopts:=dpopts, ran:=false, exp:=exp, transforms:=[], verbosity:=1, callbacks:=[]));
    end,

    resume := meth(self, transforms)
        if not ForAll(UserRecFields(self.exp), e -> ForAll(self.entries(transforms, e), e->e<>false))
            then self._reloadAllHashes();
        fi;
        self.run(transforms);
    end,

    generateCfiles := true,
    measureFinal := true,

	_fileRoot := meth(self, exp, t, opts)
		if IsBound(opts.vector) and IsBound(opts.vector.conf) and IsBound(opts.vector.conf.functionNameRoot) then
			return opts.vector.conf.functionNameRoot;
		else
			return Concat(exp, "_", Drop(CodeletName(CodeletShape(t)), 1));
		fi;
	end,

	_freq := meth(self, opts)
		if IsBound(opts.vector) and IsBound(opts.vector.conf) and IsBound(opts.vector.conf.target) and IsBound(opts.vector.conf.target.freq) then
			return opts.vector.conf.target.freq;
		elif IsBound(LocalConfig.cpuinfo.freq) then
			return LocalConfig.cpuinfo.freq;
		else
			return 1000;
		fi;
	end,

    fileTimer     := (self, exp, t, opts) >> Concat(self._fileRoot(exp, t, opts), ".timer"),
    fileVerifierf := (self, exp, t, opts) >> Concat(self._fileRoot(exp, t, opts), ".verifier"),
    fileStub      := (self, exp, t, opts) >> Concat(self._fileRoot(exp, t, opts), ".h"),
    fileTransform := (self, exp, t, opts) >> Concat(self._fileRoot(exp, t, opts), ".c"),
    funcTransform := (self, exp, t, opts) >> self._fileRoot(exp, t, opts),
    prodFileTransform := (self, exp, t, opts) >> self.fileTransform(exp, t, opts), # used in generateProductionCode
    prodFuncTransform := (self, exp, t, opts) >> self.funcTransform(exp, t, opts),

    txtFileName   := (exp, runMethod) -> Concat(exp, ".", SubString(runMethod, 5), ".txt"),

    verify := meth(self, opts, ruletree, code)
        local mat;
        mat := When(ruletree.node.isReal() or opts.dataType = "complex" or opts.generateComplexCode,
                    MatSPL(ruletree.node),
                    RCMatCyc(MatSPL(ruletree.node)));
        return VerifyMatrixCode(code, mat, opts);
    end,

    #NOTE: Slightly hacked in. Check for opts.profile being bound etc. Look at VerifyMatrixRuleTree.
    verifyfftw := (self, opts, code) >> opts.profile.verifyfftw(code, opts),
    verifyquick := (self, opts, code) >> opts.profile.verifyquick(code, opts),

    resumeAll := meth(self)
        local e;
        for e in UserRecFields(self.exp) do
            if (self.verbosity>0) then
               PrintLine("Resuming ", e);
            fi;
            self.resume(self.exp.(e).benchTransforms);
        od;
    end,

    _runAll := meth(self, runMethod)
        local e;
        for e in UserRecFields(self.exp) do
            PrintLine("Running ", e);
            self.(runMethod)(self.exp.(e).benchTransforms);
        od;
    end,

    allTrees := (self) >> let(exp := self.exp.(UserRecFields(self.exp)[1]),
    List(exp.benchTransforms, t ->
        ApplyRuleTreeSPL( HashLookup(exp.hashTable, HashAsSPL(t))[1].ruletree,
        t, exp))),

    runAll            := (self) >> self._runAll("run"),
    runExhaustiveAll  := (self) >> self._runAll("runExhaustive"),
    runRandomAll      := (self) >> self._runAll("runRandom"),
    runRandomSaveAll  := (self) >> self._runAll("runRandomSave"),
    pickRandomSaveAll := (self) >> self._runAll("pickRandomSave"),

    run            := (self, transforms) >> self._run(transforms, "_runDP", true),
    runExhaustive  := (self, transforms) >> self._run(transforms, "_runExhaustive", true),
    runRandom      := (self, transforms) >> self._run(transforms, "_runRandom", false),
    runRandomSave  := (self, transforms) >> self._run(transforms, "_runRandomSave", true),
    pickRandomSave := meth(self, transforms)
        local generateCfiles;

        #NOTE: once generateCfiles quits measuring things, this can be removed
        generateCfiles := self.generateCfiles;
        self.generateCfiles := false;
        self._run(transforms, "_pickRandomSave", true);
        self.generateCfiles := generateCfiles;
    end,

    # Find best using DP
    _runDP         := (self, e, t, opts) >> TimedAction(DP(t, self.dpopts, opts)),

    # Find best using an exhaustive search
    outputExhaustive := false,
    _runExhaustive := meth(self, e, t, opts)
       local r, searchTime, mincycles, mintree, rt, c, compiletime, cm, measuretime;

       r := AllRuleTrees(t, opts);
       searchTime := 0;
       mincycles := 10^100;

       for rt in r do
          [c,  compiletime] := TimedAction(CodeRuleTreeOpts(rt, opts));
          [cm, measuretime] := TimedAction(CMeasure(c, opts));
          if self.outputExhaustive then _seqPerfStatsGflops(Concat(e, ".Exhaustive-all.txt"), t, self._freq(opts), self.artcost(rt), cm, compiletime+measuretime); fi;
          if (cm < mincycles) then
              mincycles := cm; mintree := Copy(rt);
          fi;
          searchTime:=searchTime+compiletime+measuretime;
       od;
       HashDelete(opts.hashTable,t);
       HashAdd(opts.hashTable, t, [rec(ruletree:=mintree, measured:=mincycles)]);
       return([mintree, searchTime, c, mincycles]);
    end,

    # Run a Random ruletree. Useful for quick, dirty, non-comprehensive tests.
    _runRandom := meth(self, e, t, opts)
       local r, c, rrtime, codetime, runtime, cycles;

       [r, rrtime]      := TimedAction(RandomRuleTree(t, opts));
       [c, codetime]    := TimedAction(CodeRuleTreeOpts(r, opts));
       [cycles, runtime]:= TimedAction(CMeasure(c, opts));

       return([r, (rrtime+codetime+runtime), c, cycles]);
     end,

    # Run random search and save result in hash.
    _runRandomSave := meth(self, e, t, opts)
       local r, c, rrtime, codetime, runtime, cycles;

       [r, rrtime]      := TimedAction(RandomRuleTree(t, opts));
       [c, codetime]    := TimedAction(CodeRuleTreeOpts(r, opts));
       [cycles, runtime]:= TimedAction(CMeasure(c, opts));

       HashDelete(opts.hashTable, t);
       HashAdd(opts.hashTable, t, [rec(ruletree:=r, measured:=runtime)]);

       return([r, (rrtime+codetime+runtime), c, cycles]);
    end,

    # Pick random and save result in hash (NOT measured).
    _pickRandomSave := meth(self, e, t, opts)
       local r, searchtime;
       [r, searchtime] := TimedAction(RandomRuleTree(t, opts));

       HashDelete(opts.hashTable,t);
       HashAdd(opts.hashTable,t,[rec(ruletree:=r)]);

       return([r, searchtime, false, -1]);
    end,

    _run := meth(self, transforms, runMethod, useHash)
        local code, t, e, outf, res, opts, cycles, hentry,  date, i, searchTime, acc, optsForFile, f, ruletree, randomRes;

        Constraint(ForAll(transforms, IsSPL));

        for f in self.callbacks do
			f(self);
		od;

        for e in UserRecFields(self.exp) do
            opts := self.exp.(e);
            date := Date();

            for t in transforms do
				code := false;
				t := SumsUnification(t, opts);
				if useHash then
					# For run methods that use hash tables
					hentry := HashLookup(opts.hashTable, HashAsSPL(t));
					if hentry = false or hentry = [] then
						res := self.(runMethod)(e, HashAsSPL(t), opts);
						if res[1] = [] then
							Error("DP did not find any ruletrees for <t> (", t, ")");
						fi;
						self._saveHash(opts.hashFile, e, date, opts.hashTable);
						hentry := HashLookup(opts.hashTable, HashAsSPL(t))[1];
						hentry.searchTime := res[2];
						searchTime := res[2];
						if Length(res) >= 3 then
							# _runDP doesn't return code
							code := res[3];
						fi;
					else
						hentry := hentry[1];
						searchTime := -1;
					fi;
					hentry.spectree := ApplyRuleTreeSPL(hentry.ruletree, t, opts);
                    #NOTE: exhaustive search will not update cycles
					cycles := When(IsBound(hentry.measured), hentry.measured, 0);
					if not t in self.transforms then
						Add(self.transforms, t);
					fi;
					ruletree := hentry.spectree;
				else
					# For run methods that don't use hash tables
					randomRes  := self.(runMethod)(e, HashAsSPL(t), opts);
					ruletree   := ApplyRuleTreeSPL(randomRes[1], t, opts);
					searchTime := randomRes[2];
					code       := randomRes[3];
					cycles     := randomRes[4];
				fi;

				#HACK: It's a pain to do this in a cleaner way
				if runMethod = "_pickRandomSave" then
					return;
				fi;

				if self.generateCfiles then
					#NOTE: Should also output stub.h as <filename.h>
					compiler.CMEASURE_CURRENT_TREE := ruletree;
					compiler.CMEASURE_LAST_CODE := false;
					if (code=false or HashAsSPL(t)<>t) then
						code := CodeRuleTree(ruletree, opts);
					fi;
					compiler.CMEASURE_LAST_CODE := code;
					if (self.measureFinal and (runMethod <> "_runRandom"))then
						cycles := CMeasure(code, opts);
					fi;

					if useHash then
						hentry.measured := cycles;
					fi;

					opts.fileinfo := rec(
						cycles  := cycles,
						flops   := self.artcost(ruletree),
						file    := self.fileTransform(e,t,opts),
						algorithm := ruletree
					);
					PrintTo(self.fileTransform(e,t,opts), PrintCode(self.funcTransform(e,t,opts), code, opts));
					Unbind(opts.fileinfo);
				fi;

				if self.matrixVerify or self.fftwVerify or self.quickVerify then
					if code=false then
						code := CodeRuleTree(ruletree, opts);
					fi;
					if self.matrixVerify then
						acc := self.verify(opts, ruletree, code);
					elif self.fftwVerify then
                        acc := self.verifyfftw(opts, code);
                    else
                        acc := self.verifyquick(opts, code);
                    fi;
                    if IsBound(opts.outputVecStatistics) and opts.outputVecStatistics then
                        _seqPerfStatsGflopsAccCount(self.txtFileName(e, runMethod), t,
                                [self.artcost(ruletree), code.countedArithCost(opts.vector.isa.countrec)*opts.vector.vlen], cycles, searchTime, acc,
                                    code.countOps(opts.vector.isa.countrec), opts.vector.isa.countrec);
                    else
                        _seqPerfStatsGflopsAcc(self.txtFileName(e, runMethod), t,
                                [self.artcost(ruletree), self.countedArithCost(code, opts)], cycles, searchTime, acc);
                    fi;
                else
					if (opts.verbosity>-1) then
                       if IsBound(opts.outputVecStatistics) and opts.outputVecStatistics then
                            _seqPerfStatsGflopsCount(self.txtFileName(e, runMethod), t,
                                    [self.artcost(ruletree), code.countedArithCost(opts.vector.isa.countrec)*opts.vector.vlen], cycles, searchTime,
                                        code.countOps(opts.vector.isa.countrec), opts.vector.isa.countrec);
                       else
                            _seqPerfStatsGflops(self.txtFileName(e, runMethod), t, self._freq(opts),
                                    [self.artcost(ruletree), self.countedArithCost(code, opts)], cycles, searchTime);
                       fi;
                    fi;
                fi;
            od;
        od;
        self.ran := true;
    end,

    generateCode := (self, transforms, exp) >> self._generateCode(transforms, exp, self.exp.(exp)),

    generateProductionCode := (self, transforms, exp) >> self._generateCode(transforms, exp, self.exp.(exp).production()),

    generateAllCode := self >> DoForAll(UserRecFields(self.exp), exp ->
        self._generateCode(self.exp.(exp).benchTransforms, exp, self.exp.(exp))),

    generateAllProductionCode := self >> DoForAll(UserRecFields(self.exp), exp ->
        self._generateCode(self.exp.(exp).benchTransforms, exp, self.exp.(exp).production())),

    entries := (self, transforms, exp) >>
        When(not IsBound(self.exp.(exp)), Error("No such experiment '",exp, "'"),
         Map(transforms,
         x -> let(lookup := MultiHashLookup(Concatenation([self.exp.(exp).hashTable], self.exp.(exp).baseHashes), x),
                When(lookup = false, false,
                #Error("Transform '", x, "' not found in the '", exp, "' table"),
                lookup[1])))),

    times := (self, transforms, exp) >> Map(self.entries(transforms, exp), x->x.measured),

    alltimes := (self, transforms) >>
        let(tr := When(IsList(transforms), transforms, [transforms]),
            Map(UserRecFields(self.exp), e -> Concatenation([e], self.times(tr, e)))),

    speedup := (self, transforms, baselineExp) >>
        let(b := self.times(transforms, baselineExp),
        Map(self.alltimes(transforms),
            times -> Map([1..Length(times)],
                         i -> When(not IsInt(times[i]), times[i], times[i] / b[i-1])))),

    flopcyc := (self, transforms) >>
        self.scaledTimes(transforms, e -> self.acost(e) / e.measured),

    mflops := (self, transforms, mhz) >>
        self.scaledTimes(transforms, e -> self.acost(e) * mhz / e.measured),

    # FLoating point Operations Per Cycle
    scaledTimes := (self, transforms, scaleFunc) >>
        let(tr := When(IsList(transforms), transforms, [transforms]),
            Map(UserRecFields(self.exp),
                e -> Concatenation([e], Map(self.entries(tr, e), scaleFunc)))),

    # Use NonTerminal.normalizedArithCost() if it is there, otherwise return 0
    acost := entry -> let(
        t := entry.ruletree.node,
        When(IsBound(t.normalizedArithCost), t.normalizedArithCost(), 0)
    ),

    artcost := ruletree -> let(
        t := ruletree.node,
        When(IsBound(t.normalizedArithCost), t.normalizedArithCost(), 0)
    ),

    countedArithCost := (self, c, opts) >>
        When(IsRec(c) and
             IsBound(c.countedArithCost) and
             IsBound(opts.vector) and IsBound(opts.vector.isa) and IsBound(opts.vector.isa.countrec) and IsBound(opts.vector.vlen),
             c.countedArithCost(opts.vector.isa.countrec)*opts.vector.vlen, 0),

    matrixVerify := false,
    setMatrixVerify := self >> CopyFields(self, rec(matrixVerify := true)),

    fftwVerify   := false,
    quickVerify   := false,

    getOpts := self >> self.exp.(UserRecFields(self.exp)[1]),

    _getResult := meth(arg)
        local self, t, exp, lookup, rt;
        self := arg[1];
        exp := UserRecFields(self.exp)[1];
        t := When(Length(arg) >= 2, self.exp.(exp).benchTransforms[arg[2]], self.exp.(exp).benchTransforms[1]);
        lookup := MultiHashLookup(Concatenation([ self.exp.(exp).hashTable ], self.exp.(exp).baseHashes), HashAsSPL(t));
        rt := When(lookup = false, false, lookup[1].ruletree);
        return rec(opts := self.exp.(exp), t := t, rt:= rt);
    end,

    getResult := meth(arg)
        local self, l, rt, opts, c;
        self := arg[1];
        l := When(Length(arg)>1, self._getResult(arg[2]), self._getResult());
        rt := l.rt;
        opts := l.opts;
        c := CodeRuleTree(rt, opts);
        return CopyFields(l, rec(c := c, opcount := c.countOps(opts.vector.isa.countrec)));
    end,

    _runAll := meth(self, runMethod)
        local e;
        for e in UserRecFields(self.exp) do
            PrintLine("Running ", e);
            self.(runMethod)(self.exp.(e).benchTransforms);
        od;
    end,

    build := function(arg)
        local transforms, opts, bopts, name, dpr;

        transforms := When(IsList(arg[1]), arg[1], [arg[1]]);
        opts := When(Length(arg) >= 2, arg[2], SpiralDefaults);
        opts.benchTransforms := transforms;
        bopts := When(Length(arg) >= 3, arg[3], rec());
        name := When(Length(arg) >= 4, arg[4], "spiral");
        dpr := When(Length(arg) >= 5, arg[5], rec(verbosity := 0, timeBaseCases:=true));

        return CopyFields(DPBench(rec((name) := opts), dpr), bopts);

    end

));

# Class(NewDPBench, DPBench, rec(
#     measureNtimes := 10,
#     remeasure := (self, c, opts) >> List([1..self.measureNtimes],
#         i -> CMeasure(c, opts)),

#     runHooks := [
#         meth(self, t, c, hentry, opts)
#             hentry.remeasure := self.remeasure(c, opts);
#             PrintLine("remeasure : ", hentry.remeasure);
#         end
#     ]
# ));

#opts := SpiralDefaults; opts2 := CopyFields(opts, rec(declareConstants := true));
#d := NewDPBench(rec(default:=opts,pullconst:=opts2), rec(timeBaseCases := false, verbosity := 0));
sampleDPBench := DPBench(rec(default := SpiralDefaults), rec(timeBaseCases := false, verbosity := 0));
sampleDPBench.transforms := [DFT(2), DFT(3), DFT(4)];