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

# We are interested to see if
# IMDPRDFT(l)*MDPRDFT(l) = Product(l)*MatSPL(Product(l));
# because we want to see if we can recover the original REAL data.

Declare(MDPRDFT);
Declare(IMDPRDFT);

# Same as (floor(n_t/2)+1)*2.
# RClength := (n) -> (n + 1) mod 2 + n + 1;
RClength := (n) -> PRDFT1(n).dims()[1];

IJmatrix := (n) -> DirectSum(I(1), J(n-1));

tensorIJmatrix := (l) -> When(Length(l) > 0,
                         Tensor(IJmatrix(l[1]), tensorIJmatrix(Drop(l, 1))),
                         Diag([1, -1]));

#F MDPRDFT(<dims>, [<exp>]) - multi-dimensional PRDFT non-terminal
#F   dims = [ <n_1>,.., <n_t> ] list of (positive) dimensions
#F   exp = root of unity exponent scaling (see DFT for exact definition)
#F
#F Definition : multidimensional matrix of size M x N, where
#F M = n_1*...*n_{t-1}*(floor(n_t/2)+1)*2
#F N = n_1*..*n_t
#F This matrix has real components, for an operator with
#F N real inputs,
#F M real outputs for interleaved real and imaginary components.
#F
#F Example (direct)  : MDPRDFT([2,4,4])
#F Example (inverse) : MDPRDFT([2,4,4], -1)
#F
#F In last dimension, t: real-to-complex DFT;
#F then dimensions t-1, ..., 1: complex-to-complex DFT.
#F

# To verify that components of m are real:
# Im(MatSPL(m)) = MatSPL(ApplyFunc(O, m.dims()));

Class(MDPRDFT, TaggedNonTerminal, rec(
    a_lengths := self >> self.params[1],
    a_exp := self >> self.params[2],

    # Components of self.a_lengths() that are equal to 1 are removed.
    abbrevs := [
        P     -> Checked(IsList(P),
                         ForAll(P, IsPosInt),
                         Product(P) > 1,
                         [ RemoveOnes(P), 1 ]),
        (P,k) -> Checked(IsList(P),
                         ForAll(P, IsPosInt),
                         IsInt(k),
                         Product(P) > 1,
                         Gcd(Product(P), k) = 1,
                         [ RemoveOnes(P), k mod Product(P) ])
        ],

    # dims() has 2 components, counting outputs and inputs:
    # dims()[1] is product of all but last dimension, and RClength on last.
    # dims()[2] is product of all dimensions given.
    dims := self >> let(a_lengths := self.a_lengths(),
                        [Product(DropLast(a_lengths, 1)) *
                         RClength(Last(a_lengths)),
                         Product(a_lengths)
                        ]),

terminate := self >>  let(
        a_lengths := self.a_lengths(),
        Ninputs := self.dims()[2],              # number of scalar inputs
        Nlast := Last(a_lengths),        # length of last dimension
        Nrest := DropLast(a_lengths, 1), # all dimension lengths but last
        rcNlast := RClength(Nlast),
        # interleave has a 2x1 matrix [1; 0] along diagonal for each input;
        # it is a block-diagonal matrix of size (2*Ninputs) x Ninputs
        # that takes real inputs and puts a 0 between each.
        interleave := Scat(fTensor(fId(Ninputs), fBase(2,0))),
        # Here we are defining mddft as multi-dimensional DFT that takes
        # input as interleaved real and imaginary data,
        # output as interleaved real and imaginary data.
        mddft := RC(ApplyFunc(MDDFT, self.params)),
        # gath is block-diagonal matrix consisting of Product(Nrest) blocks
        # of size rcNlast x (2*Nlast), so it has
        # input as Product(Nrest) blocks of length 2*Nlast,
        # output as Product(Nrest) blocks of length rcNlast.
        f := fAdd(2*Nlast, rcNlast, 0), # pad length rcNlast to 2*Nlast
        # could also set f := HStack(I(rcNlast), O(rcNlast, 2*Nlast-rcNlast))
        ff := fTensor(List(Nrest, N -> fId(N)), f),
        # could also set ff := Tensor(List(Nrest, N->I(N)), f), no Gath needed.
        gath := Gath(ff),
        # 1. interleave: take real input, put a 0 after each element.
        # 2. mddft: multi-dimensional DFT, interleaved complex input & output.
        # 3. gath: in last dimension, keep only rcNlast of the 2*Nlast components.
        mdprdft := gath * mddft * interleave,
        mdprdft),

    transpose := self >> IMDPRDFT(self.a_lengths(), -self.a_exp()),

    isReal := True,

    normalizedArithCost :=  (self) >> let(n := Product(self.a_lengths()),
                                        IntDouble(2 * n * d_log(n) / d_log(2)) )
));

# Conjugate in interleaved form:  (a, b) maps to (a, -b).
# That's multiplication by matrix [ [1, 0], [0, -1] ].
# That's fBase(2, 0) above -fBase(2, 1).
# Try adapting fId, replacing let(i := Ind(self.params[1]), Lambda(i,i)).
# Or dLin or dOmega.

#F IMDPRDFT(<dims>, [<exp>]) - multi-dimensional inverse PRDFT non-terminal
#F   dims = [ <n_1>,.., <n_t> ] list of (positive) dimensions
#F   exp = root of unity exponent scaling (see DFT for exact definition)
#F
#F Definition : multidimensional matrix of size N x M, where
#F N = n_1*..*n_t
#F M = n_1*...*n_{t-1}*(floor(n_t/2)+1)*2
#F This matrix has real components, for an operator with
#F M real inputs for interleaved real and imaginary components,
#F N real outputs.
#F
#F Example (direct)  : IMDPRDFT([2,4,4])
#F Example (inverse) : IMDPRDFT([2,4,4], -1)
#F
#F In dimensions 1, .., t-1: complex-to-complex DFT;
#F then in last dimension, t: complex-to-real DFT.
#F
Class(IMDPRDFT, TaggedNonTerminal, rec(
    a_lengths := self >> self.params[1],
    a_exp := self >> self.params[2],

    # Components of self.a_lengths() that are equal to 1 are removed.
    abbrevs := [
        P     -> Checked(IsList(P),
                         ForAll(P, IsPosInt),
                         Product(P) > 1,
                         [ RemoveOnes(P), 1 ]),
        (P,k) -> Checked(IsList(P),
                         ForAll(P, IsPosInt),
                         IsInt(k),
                         Product(P) > 1,
                         Gcd(Product(P), k) = 1,
                         [ RemoveOnes(P), k mod Product(P) ])
        ],

    # dims() has 2 components, counting outputs and inputs:
    # dims()[1] is product of all dimensions given.
    # dims()[2] is product of all but last dimension, and RClength on last.
    dims := self >> let(a_lengths := self.a_lengths(),
                        [Product(a_lengths),
                         Product(DropLast(a_lengths, 1)) *
                         RClength(Last(a_lengths))
                        ]),

    terminate := self >>  let(
        a_lengths := self.a_lengths(),
        Nlast := Last(a_lengths),        # length of last dimension
        Nrest := DropLast(a_lengths, 1), # all dimension lengths but last
        rcNlast := RClength(Nlast),
        newcomplex := Nlast - rcNlast/2,
        Nallbutlast := Product(Nrest),
        shuffle := Tensor( L(Nallbutlast*rcNlast/2, rcNlast/2), I(2)),
        unshuffle := Tensor( L(Nallbutlast*Nlast, Nallbutlast), I(2)),
        # This is the basic conjugation matrix.
        flipper := tensorIJmatrix(Nrest),
        newrows := 2*newcomplex*Nallbutlast,
        nzeroleft := When(Nlast mod 2 = 0,
                          Nallbutlast*(rcNlast - 2*newcomplex - 2),
                          Nallbutlast*(rcNlast - 2*newcomplex)),
        hpart := When(newrows > 0,
           HStack(
                  O(newrows, nzeroleft),
                  Tensor(J(newcomplex), flipper),
                  When(Nlast mod 2 = 0, O(newrows, 2*Nallbutlast), [ ])
                 ),
            [ ]),
        stacker := VStack(I(rcNlast*Nallbutlast), hpart),
        # Here we are defining mddft as multi-dimensional DFT that takes
        # input as interleaved real and imaginary data,
        # output as interleaved real and imaginary data.
        imddft := RC(ApplyFunc(MDDFT, self.params)),
        # extract real has 1x2 matrix [1, 0] along diagonal for each input;
        # it is a block-diagonal matrix of size (2*Noutputs) x Noutputs
        # that extracts every other component of the input.
        extractreal := Gath(fTensor(List(a_lengths, N->fId(N)), fBase(2, 0))),
        extend := unshuffle * stacker * shuffle,
        # 1. extend: in last dimension, extend the rcNlast interleaved complex input to 2*Nlast.
        # 2. imddft: inverse multi-dimensional DFT, interleaved complex input & output.
        # 3. extractreal: return the real parts of interleaved complex data.
        imdprdft := extractreal * imddft * extend,
        imdprdft),

    transpose := self >> MDPRDFT(self.a_lengths(), -self.a_exp()),

    isReal := True,

    normalizedArithCost :=  (self) >> let(n := Product(self.a_lengths()),
                                        IntDouble(2 * n * d_log(n) / d_log(2)) )
));

max_stars_repo_name: franzfranchetti/spiral-package-fftx

max_issues_repo_name: franzfranchetti/spiral-package-fftx

max_forks_repo_name: franzfranchetti/spiral-package-fftx


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

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

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

  # XXX
  return SemigroupOfSelfSimFamily(SelfSimFamily(list, bind_vars));
end);

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

  # XXX
  return SemigroupOfSelfSimFamily(SelfSimFamily(list, names));
end);

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

  #XXX
  return SemigroupOfSelfSimFamily(SelfSimFamily(list, names, bind_vars));
end);

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

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

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

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

###############################################################################
##
#M  SemigroupOfSelfSimFamily(<G>)
##
InstallMethod(SemigroupOfSelfSimFamily, "for [IsSelfSimSemigroup]",
                   [IsSelfSimSemigroup],
function(G)
  return SemigroupOfSelfSimFamily(UnderlyingSelfSimFamily(G));
end);

###############################################################################
##
#M  IsSemigroupOfSelfSimFamily(<G>)
##
InstallMethod(IsSemigroupOfSelfSimFamily, "for [IsSelfSimSemigroup]",
              [IsSelfSimSemigroup],
function(G)
  return G = SemigroupOfSelfSimFamily(G);
end);

###############################################################################
##
#M  IsSelfSimilar(<G>)
##
InstallMethod(IsSelfSimilar, "for [IsSelfSimSemigroup]",
              [IsSelfSimSemigroup],
function(G)
  local g, i, res;
  res := true;
  for g in GeneratorsOfSemigroup(G) do
    for i in [1..UnderlyingSelfSimFamily(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  UnderlyingSelfSimFamily(<G>)
##
InstallMethod(UnderlyingSelfSimFamily, "for [IsSelfSimSemigroup]",
              [IsSelfSimSemigroup],
function(G)
  return FamilyObj(GeneratorsOfSemigroup(G)[1]);
end);

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

InstallMethod(TopDegreeOfTree, "for [IsSelfSimSemigroup]",
              [IsSelfSimSemigroup],
function(G)
  return DegreeOfTree(UnderlyingSelfSimFamily(G));
end);

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

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

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

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

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

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

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

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

  return s;
end);

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

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

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

  TryNextMethod();
end);

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

  if HasIsGroupOfSelfSimFamily(G) and IsGroupOfSelfSimFamily(G) then
    return true;
  fi;

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

  fam := UnderlyingSelfSimFamily(G);

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

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

  TryNextMethod();
end);

###############################################################################
##
#M  Random(<G>)
##
InstallMethodWithRandomSource(Random, "for a random source and [IsSelfSimSemigroup]",
              [IsRandomSource, IsSelfSimSemigroup],
function(rs, G)
  local w, monoid, F, gens, pi;
  if IsSelfSimilarSemigroup(G) then
    monoid := UnderlyingFreeMonoid(G);

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

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

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

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

InstallMethod(SphericalIndex, "for [IsSelfSimSemigroup]",
              [IsSelfSimSemigroup],
function(G)
  return SphericalIndex(GeneratorsOfSemigroup(G)[1]);
end);

InstallMethod(DegreeOfTree, "for [IsSelfSimSemigroup]",
              [IsSelfSimSemigroup],
function(G)
  return UnderlyingSelfSimFamily(G)!.deg;
end);

InstallMethod(TopDegreeOfTree, "for [IsSelfSimSemigroup]",
              [IsSelfSimSemigroup],
function(G)
  return UnderlyingSelfSimFamily(G)!.deg;
end);

###############################################################################
##
#M  IsSelfSimilarSemigroup(<G>)
##
##  Returns `true' if generators of <G> coincide with generators of the family
InstallMethod(IsSelfSimilarSemigroup, [IsSelfSimSemigroup],
function(G)
  local fam;
  fam := UnderlyingSelfSimFamily(G);
  return fam!.numstates = 0 or
         GeneratorsOfSemigroup(G) = fam!.recurgens{[1..fam!.numstates]};
end);

###############################################################################
##
#M  IsFiniteState(<G>)
##
InstallMethod(IsFiniteState, "for [IsSelfSimSemigroup]",
              [IsSelfSimSemigroup],
function(G)
  local states, MealyAutomatonLocal, aut_list, gens, images, H, g, hom_function, \
        inv_hom_function, hom, free_groups_hom, inv_free_groups_hom, inv_hom, \
        gens_in_freegrp, images_in_freegrp, preimages_in_freegrp, F, pi, pi_bar, \
        preimage_in_freegrp;

  MealyAutomatonLocal := function(g)
    local cur_state;
    if g!.word in states then return Position(states, g!.word); fi;
    Add(states, g!.word);
    cur_state := Length(states);
    aut_list[cur_state] := List([1..g!.deg], x -> MealyAutomatonLocal(Section(g, x)));
    Add(aut_list[cur_state], g!.perm);
    return cur_state;
  end;

  states := [];
  aut_list := [];
  gens := GeneratorsOfSemigroup(G);
  images := [];

  for g in gens do
    Add(images, MealyAutomatonLocal(g));
  od;

  H := AutomatonSemigroup(aut_list);

  images := UnderlyingAutomFamily(H)!.oldstates{images};

  SetIsomorphicAutomSemigroup(G, SemigroupByGenerators(UnderlyingAutomFamily(H)!.automgens{images}));
  SetUnderlyingAutomatonSemigroup(G, H);

# preimages of generators of G in UnderlyingFreeGroup(G)
  gens_in_freegrp := List(GeneratorsOfSemigroup(G), Word);

# preimages of generators of a subgroup of H isomorphic to G in UnderlyingFreeGroup(H)
  images_in_freegrp := List(UnderlyingAutomFamily(H)!.automgens{images}, Word);

  preimage_in_freegrp := function(x)
    local w;
    if IsOne(x!.word) then
      return states[ Position( UnderlyingAutomFamily(H)!.oldstates, UnderlyingAutomFamily(H)!.trivstate)];
    fi;
    w := LetterRepAssocWord(x!.word)[1];
    if w>0 then
      return states[ Position( UnderlyingAutomFamily(H)!.oldstates, w)];
    else
      return states[ Position( UnderlyingAutomFamily(H)!.oldstates, -w+UnderlyingAutomFamily(H)!.numstates)];
    fi;
  end;

# preimages of generators of H in UnderlyingFreeGroup(G)
  preimages_in_freegrp := List(GeneratorsOfSemigroup(H), x -> preimage_in_freegrp(x));

  if IsSelfSimilarSemigroup(G) then
    free_groups_hom :=
       GroupHomomorphismByImagesNC( Group(gens_in_freegrp), UnderlyingFreeGroup(H),
                                    gens_in_freegrp, images_in_freegrp );

    inv_free_groups_hom :=
       GroupHomomorphismByImagesNC( UnderlyingFreeGroup(H), UnderlyingFreeGroup(G),
                                    UnderlyingFreeGenerators(H), preimages_in_freegrp );

    hom_function := function(a)
      return Autom(Image(free_groups_hom, a!.word), UnderlyingAutomFamily(H));
    end;

    inv_hom_function :=  function(b)
      return SelfSim(Image(inv_free_groups_hom, b!.word), UnderlyingSelfSimFamily(G));
    end;

    hom := MappingByFunction(G, SemigroupByGenerators(UnderlyingAutomFamily(H)!.automgens{images}), hom_function, inv_hom_function);

    SetMonomorphismToAutomatonSemigroup(G, hom);
  else
    F := FreeGroup(Length(GeneratorsOfSemigroup(G)));
    #SetCoveringFreeGroup(G, F);

#        pi
#    F ----- ->  G --- ->  UnderlyingFreeGroup(H)
#      ------------- ->
#            pi_bar

    pi := GroupHomomorphismByImages(F,                     Group(gens_in_freegrp),
                                    GeneratorsOfGroup(F),  gens_in_freegrp);

    pi_bar := GroupHomomorphismByImages(F,                     UnderlyingFreeGroup(H),
                                        GeneratorsOfGroup(F),  images_in_freegrp);

    hom_function := function(g)
      return Autom(Image(pi_bar, PreImagesRepresentative(pi, g!.word)), UnderlyingAutomFamily(H));
    end;

    inv_hom_function :=  function(b)
      return SelfSim(Image(pi, PreImagesRepresentative(pi_bar, b!.word)), UnderlyingSelfSimFamily(G));
    end;

    hom := MappingByFunction(G, SemigroupByGenerators(UnderlyingAutomFamily(H)!.automgens{images}), hom_function, inv_hom_function);

    SetMonomorphismToAutomatonSemigroup(G, hom);
  fi;

  return true;
end);

###############################################################################
##
#M  IsomorphicAutomSemigroup(<G>)
##
InstallMethod(IsomorphicAutomSemigroup, "for [IsSelfSimSemigroup]",
              [IsSelfSimSemigroup],
function(G)
  if IsFiniteState(G) then return IsomorphicAutomSemigroup(G); fi;
end);

###############################################################################
##
#M  UnderlyingAutomatonSemigroup(<G>)
##
InstallMethod(UnderlyingAutomatonSemigroup, "for [IsSelfSimSemigroup]",
              [IsSelfSimSemigroup],
function(G)
  if IsFiniteState(G) then return UnderlyingAutomatonSemigroup(G); fi;
end);

###############################################################################
##
#M  MonomorphismToAutomatonSemigroup(<G>)
##
InstallMethod(MonomorphismToAutomatonSemigroup, "for [IsSelfSimSemigroup]",
              [IsSelfSimSemigroup],
function(G)
  if IsFiniteState(G) then return MonomorphismToAutomatonSemigroup(G); fi;
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-2021, Carnegie Mellon University
# See LICENSE for details

#F WHT( <log2(size)> ) - Walsh Hadamard Transform non-terminal
#F Definition: (2^n x 2^n)-matrix  DFT_2 tensor ... tensor DFT_2 (n-fold)
#F Note:       DFT_2 denotes the matrix [[1, 1], [1, -1]],
#F             WHT is symmetric.
#F Example:    WHT(3)
#F
#F Tweakable parameters: WHT_GeneralSplit.maxSize,
#F                       WHT_BinSplit.minSize
#F                       WHT_DDL.minSize
#F
#F WHT_GeneralSplit is a generalization of BinSplit, so it only makes
#F sense to have WHT_BinSplit.minSize := WHT_GeneralSplit.maxSize + 1
#F (not to have 2 rules enabled at the same time)
#F
#F By default, WHT_GeneralSplit is disabled.
#
Class(WHT, TaggedNonTerminal, rec(
  abbrevs   := [ n -> Checked(IsPosInt(n), [n]) ],
  dims      := self >> let(size := 2^self.params[1], [size, size]),
  terminate := self >> When(self.params[1] = 1, F(2), Tensor(Replicate(self.params[1], F(2)))),
  transpose := self >> Copy(self),
  isReal    := self >> true,
  SmallRandom := () -> Random([2..5]),
  LargeRandom := () -> Random([6..15]),
  normalizedArithCost := self >> self.params[1]*2^self.params[1],
  TType := T_Real(64)
));

# NOTE: A lot of WHT rules were moved to paradigms/common/wht.g, UGLY!
#
NewRulesFor(WHT, rec(
    #F WHT_Base: WHT_(2^1) = F_2
    #F
    WHT_Base := rec(
        info             := "WHT_(2^1) -> F_2",
        forTransposition := false,
        applicable       := (self, nt) >> nt.params[1]=1 and not nt.hasTags(),
        apply            := (nt, c, cnt) -> F(2),
    ),

    #F WHT_Base2: WHT_(2^1) = F_2
    # Added by Kyle 7-11-07
    # Need a better way to represent base cases for WHT so it is not fully expanded.
    WHT_Base2 := rec(
        info             := "WHT_(2^2) -> F_2",
        forTransposition := false,
        switch           := false,
        applicable       := (self, nt) >> nt.params[1]=2 and not nt.hasTags(),
        apply            := (nt, c, cnt) -> WHT(2),

    ),

    #F WHT_GeneralSplit: follows from definition
    #F
    #F   WHT_(2^k) =
    #F                        WHT_(2^k1) tensor I_(2^(k-k1)) *
    #F     I_(2^k1)    tensor WHT_(2^k2) tensor I_(2^(k-k1-k2)) *
    #F       ...
    #F       ...
    #F     I_(2^(k-kr)) tensor WHT_(2^kr)
    #F
    #F This rule has been restricted to small sizes
    #F
    WHT_GeneralSplit := rec (
        info             := "WHT_(2^k) -> (WHT_(2^k1) tensor I) .. (I tensor WHT_(2^kr))",
        forTransposition := false,
        switch           := false,
        maxSize          := 7,
        applicable       := (self, nt) >> nt.params[1] <> 1 and nt.params[1] <= self.maxSize,

        children := nt -> let(C := DropLast(OrderedPartitions(nt.params[1]), 1),
            List(C, p -> List(p, WHT))),

        # MRT 10-2007: .randomChildren does appear to be called in a certain case during RandomRuleTree()
        # but is it necessary? It only appears here!
        randomChildren := function ( Pin )
            local C, sum, rand, P;
            P := Pin[1];
            C := [];
            sum := 0;
            while sum < P do
                if sum = 0 then  rand := RandomList( [1..(P-1)] );
                else             rand := RandomList( [1..(P-sum)] );
                fi;
                Add( C, rand );
                sum := sum + rand;
            od;
            When(not (Sum(C)=P and Length(C) > 1),
                Error( "WHT_GeneralSplit.randomChildren has a bug"));
            return List(C, WHT);
        end,

        rule := function ( P, C )
            local P, left, right, i;
            right := 2^P / C[1].dimensions[1];
            left  := 1;
            P := [ Tensor(C[1], I(right)) ];

            # list of factors
            for i in [2..Length(C)-1] do
                left  := left * C[i-1].dimensions[1];
                right := right/C[i].dimensions[1];
                Add(P, Tensor(I(left), C[i], I(right)));
            od;
            left := left * C[Length(C)-1].dimensions[1];
            Add(P, Tensor(I(left), C[Length(C)]));

            return Compose(P); # return product
        end
    ),

    WHT_BinSplit_binloops := rec (
        info             := "WHT_(2^k) -> (WHT_(2^k1) tensor I) (I tensor WHT_(2^k2))",
        forTransposition := false,
        minSize          := 2,
        applicable       := (self, nt) >> nt.params[1] >= self.minSize and not nt.hasTags(),

        children := nt -> List([1..nt.params[1]-1], i -> [ WHT(i), WHT(nt.params[1]-i) ] ),

        apply := function(nt, c, cnt)
            local r1, r2, a, i, b;

            r1 := Rows(c[1]);
            r2 := Rows(c[2]);

            a := c[1];
            for i in [1..Log2Int(r2)] do
                a := Grp(Tensor(I(2), a));
            od;
            a := Grp(L(r1 * r2, r1) * a * L(r1 * r2, r2));

            b := c[2];
            for i in [1..Log2Int(r1)] do
                b := Grp(Tensor(I(2), b));
            od;

            return a * b;
        end,

        switch := true,
    )
));


text:
#
# Read("~/Workspace/Chevalley.gap/init.gi"); Read(Filename(home_dir,"load.gi")); Read(Filename(test_dir,"unichv.test.gi"));
#
# Read("~/Workspace/Chevalley.gap/init.gi"); Read(Filename(test_dir,"unichv.test.init.gi")); Read(Filename(test_dir,"unichv.test.gi"));
#
# Read(Filename(test_dir,"unichv.test.gi"));
#

sys:=Unipotent(ChevalleyAdj("F",4,GF(2)),[[1,1],[2,1]]);
Print("Created Unipotent ",sys,"\n");
Print("\t[",IsUnipotent(sys)=true,"] IsUnipotent(sys)=true\n");


text:
BindGlobal("MitM_rnams", ["name", "attributes", "content"]);

# OM elements - known as "omel" in the specification
BindGlobal("MitM_OMel", [ "OMS", "OMV", "OMI", "OMB", "OMSTR", "OMF",
                          "OMA", "OMBIND", "OME", "OMATTR", "OMR"]);

#
# ValidXSD: a collection of functions to check some XML string types
#
BindGlobal("MitM_ValidXSD", rec());

MitM_ValidXSD.Empty := function(str)
    if not IsEmpty(str) then
        return "must be empty";
    fi;
    return true;
end;

MitM_ValidXSD.Text := function(str)
    if '<' in str then
        return "XML strings cannot contain '<'";
    elif Number(str, c -> c = '&') > Number(str, c -> c = ';') then
        return "there is a '&' character without a closing ';'";
    fi;
    # Perhaps we could check some other simple XML things?
    return true;
end;

MitM_ValidXSD.NCName := function(str)
    local result, char;
    result := MitM_ValidXSD.Text(str);
    if result <> true then
        return result;
    elif IsEmpty(str) then
        return "must not be the empty string";
    elif IsDigitChar(str[1]) or str[1] = '-' or str[1] = '.' then
        return "must start with a letter or underscore";
    fi;
    for char in str do
        if not (IsAlphaChar(char) or IsDigitChar(char)
                or char = '_' or char = '-' or char = '.') then
            return Concatenation("must not contain the character \"",
                                 String(char), "\"");
        fi;
    od;
    return true;
end;

# Checking ID uniqueness is difficult, but we may want to do it in the future
MitM_ValidXSD.ID := MitM_ValidXSD.Text;

# Validating URIs is difficult, but we may want to do it in the future
MitM_ValidXSD.AnyURI := MitM_ValidXSD.Text;

MitM_ValidXSD.Base64Binary := function(str)
    local len, nreq, i;
    str := ShallowCopy(str);
    RemoveCharacters(str, " \n\t\r"); # ignore whitespace
    len := Length(str);
    if len mod 4 <> 0 then
        return "must have length divisible by 4";
    fi;
    nreq := 0;
    if str[len] = '=' then
        nreq := 1;
        if str[len - 1] = '=' then
            nreq := 2;
            if not str[len - 2] in "AQgw" then
                return "one of [AQgw] must come before '=='";
            fi;
        elif not str[len - 1] in "AEIMQUYcgkosw048" then
            return "one of [AEIMQUYcgkosw048=] must come before '='";
        fi;
    fi;
    for i in [1 .. len - nreq] do
        if str[i] = '=' then
            return "only 1 or 2 '=' characters allowed, and only at the end";
        elif not (IsAlphaChar(str[i]) or IsDigitChar(str[i])
                  or str[i] = '+' or str[i] = '/') then
            return Concatenation("cannot contain character '", str{[i]}, "'");
        fi;
    od;
    return true;
end;

#
# ValidAttr: each object type's valid attributes,
#            along with functions to check their values
#
BindGlobal("MitM_ValidAttr",
rec(
     OMOBJ := rec(cdbase := MitM_ValidXSD.AnyURI,
                  version := MitM_ValidXSD.Text),
     OMS := rec(name := MitM_ValidXSD.NCName,
                cd := MitM_ValidXSD.NCName,
                cdbase := MitM_ValidXSD.AnyURI),
     OMV := rec(name := MitM_ValidXSD.NCName),
     OMI := rec(),
     OMB := rec(),
     OMSTR := rec(),
     OMF := rec(dec := function(str)
                   if str = "INF" or str = "-INF" or str = "NaN" then
                       return true;
                   elif Float(str) = fail then
                       return Concatenation(str, " is not a valid float");
                   fi;
                   return true;
                end,
                hex := function(str)
                    local valid_chars, char, err;
                    if Length(str) <> 16 then
                        return "must be 16 characters long";
                    fi;
                    valid_chars := "0123456789ABCDEF";
                    for char in str do
                        if not char in valid_chars then
                            err := "contains non-hex character '";
                            Add(err, char);
                            Add(err, ''');
                            if char in "abcdef" then
                                Append(err, " (not capital)");
                            fi;
                            return err;
                        fi;
                    od;
                    return true;
                end),
     OMA := rec(cdbase := MitM_ValidXSD.AnyURI),
     OMBIND := rec(cdbase := MitM_ValidXSD.AnyURI),
     OME := rec(),
     OMATTR := rec(cdbase := MitM_ValidXSD.AnyURI),
     OMR := rec(),
     OMBVAR := rec(),
     OMATP := rec(cdbase := MitM_ValidXSD.AnyURI),
     common := rec(id := MitM_ValidXSD.ID)
));

#
# RequiredAttr: a list of required attributes for each object type
#
BindGlobal("MitM_RequiredAttr",
rec(
     OMOBJ := [],
     OMS := ["cd", "name"],
     OMV := ["name"],
     OMI := [],
     OMB := [],
     OMSTR := [],
     OMF := [],
     OMA := [],
     OMBIND := [],
     OME := [],
     OMATTR := [],
     OMR := [],
     OMBVAR := [],
     OMATP := []
));

#
# ValidCont: a function to check content for each object type
#
BindGlobal("MitM_ValidCont",
rec(
     OMOBJ := function(content)
       if Length(content) <> 1 then
         return "must be precisely one object";
       elif not (MitM_OMRec(content[1]) and MitM_Tag(content[1]) in MitM_OMel) then
         return "must be an OM element";
       fi;
       return MitM_IsValidOMRec(content[1]);
     end,

     OMS := MitM_ValidXSD.Empty,

     OMV := MitM_ValidXSD.Empty,

     OMI := function(content)
         local str, pos, valid_chars;
         if not (Length(content) = 1 and IsString(content[1])) then
             return "must be only a string";
         fi;
         str := ShallowCopy(content[1]);
         RemoveCharacters(str, " \n\t\r"); # ignore whitespace
         pos := 1;
         if str[pos] = '-' then
             Remove(str, 1);
         fi;
         if str[pos] = 'x' then
             Remove(str, 1);
             valid_chars := "0123456789ABCDEF";
         else
             valid_chars := "0123456789";
         fi;
         if ForAny(str, char -> not char in valid_chars) or Length(str) = 0 then
             return Concatenation(content[1], " is not an integer");
         fi;
         return true;
     end,

     OMB := function(content)
         if IsEmpty(content) then
             return true;
         elif not (Length(content) = 1 and IsString(content[1])) then
             return "must be only a string";
         fi;
         return MitM_ValidXSD.Base64Binary(content[1]);
     end,

     OMSTR := function(content)
         if IsEmpty(content) then
             return true;
         elif not (Length(content) = 1 and IsString(content[1])) then
             return "must be only a string";
         fi;
         return MitM_ValidXSD.Text(content[1]);
     end,

     OMF := MitM_ValidXSD.Empty,

     OMA := function(content)
         local item, result;
         if Length(content) = 0 then
             return "must not be empty";
         fi;
         for item in content do
             if not MitM_OMRec(item) then
                 return "must only contain OM elements";
             elif not MitM_Tag(item) in MitM_OMel then
                 return Concatenation("cannot contain ", MitM_Tag(item), " objects");
             fi;
             result := MitM_IsValidOMRec(item);
             if result <> true then
                 return result;
             fi;
         od;
         return true;
     end,

     OMBIND := function(content)
         local item, result;
         if not (Length(content) = 3
                 and ForAll(content, MitM_OMRec)
                 and MitM_Tag(content[1]) in MitM_OMel
                 and MitM_Tag(content[2]) = "OMBVAR"
                 and MitM_Tag(content[3]) in MitM_OMel) then
             return "must be [OM elm, OMBVAR, OM elm] (in that order)";
         fi;
         for item in content do
             result := MitM_IsValidOMRec(item);
             if result <> true then
                 return result;
             fi;
         od;
         return true;
     end,

     OME := function(content) return "not implemented"; end,

     OMATTR := function(content)
         local i, result;
         if Length(content) <> 2 then
             return "must contain precisely two objects";
         elif not (MitM_OMRec(content[1]) and MitM_Tag(content[1]) = "OMATP") then
             return "first object must be OMATP";
         elif not (MitM_OMRec(content[2]) and MitM_Tag(content[2]) in MitM_OMel) then
             return "second object must be an OM element";
         fi;
         for i in [1 .. Length(content)] do
             result := MitM_IsValidOMRec(content[i]);
             if result <> true then
                 return result;
             fi;
         od;
         return true;
     end,

     OMR := function(content) return "not implemented"; end,

     OMBVAR := function(content)
         local item, result;
         if IsEmpty(content) then
             return "must not be empty";
         fi;
         for item in content do
             if not (MitM_OMRec(item) and MitM_Tag(item) = "OMV") then
                 return "must only contain OMV objects";
                 # ... or attvar objects in the full spec
             fi;
             result := MitM_IsValidOMRec(item);
             if result <> true then
                 return result;
             fi;
         od;
         return true;
     end,

     OMATP := function(content)
         local i, result;
         if Length(content) = 0 then
             return "must not be empty";
         elif Length(content) mod 2 <> 0 then
             return "must contain an even number of objects";
         fi;
         for i in [1, 3 .. Length(content) - 1] do
             if not (MitM_OMRec(content[i]) and MitM_Tag(content[i]) = "OMS") then
                 return StringFormatted("item {} must be an OMS object", i);
             elif not (MitM_OMRec(content[i + 1])
                       and MitM_Tag(content[i + 1]) in MitM_OMel) then
                 # TODO: allow OMFOREIGN
                 return StringFormatted("item {} must be an OM element", i + 1);
             fi;
             result := MitM_IsValidOMRec(content[i]);
             if result <> true then
                 return result;
             fi;
             result := MitM_IsValidOMRec(content[i + 1]);
             if result <> true then
                 return result;
             fi;
         od;
         return true;
     end
));

#
# IsValidOMRec: the function we call on an object to check its validity
#
InstallGlobalFunction(MitM_IsValidOMRec,
function(tree)
    local rnam, attr, result;
    # Check that this is a proper object
    if not MitM_OMRec(tree) then
        return "<tree> must be an OM record";
    elif not MitM_Tag(tree) in RecNames(MitM_ValidAttr) then
        return Concatenation(MitM_Tag(tree), " is not a valid OM object name");
    fi;
    for rnam in NamesOfComponents(tree) do
        if not rnam in MitM_rnams then
            return Concatenation("invalid XML: ", rnam, " should not exist");
        fi;
    od;

    # Validate the attributes
    for attr in RecNames(MitM_Attributes(tree)) do
        if attr in RecNames(MitM_ValidAttr.(MitM_Tag(tree))) then
            result := MitM_ValidAttr.(MitM_Tag(tree)).(attr)(MitM_Attributes(tree).(attr));
        elif attr in RecNames(MitM_ValidAttr.common) then
            result := MitM_ValidAttr.common.(attr)(MitM_Attributes(tree).(attr));
        else
            return Concatenation(attr, " is not a valid attribute of ",
                                 MitM_Tag(tree), " objects");
        fi;
        if result <> true then
            return StringFormatted("{} attribute of {} object: {}",
                                   attr, MitM_Tag(tree), result);
        fi;
    od;

    # Check required attributes
    for attr in MitM_RequiredAttr.(MitM_Tag(tree)) do
        if not attr in RecNames(MitM_Attributes(tree)) then
            return Concatenation(MitM_Tag(tree), " objects must have the ",
                                 attr, " attribute");
        fi;
    od;

    # Check antirequisite attributes - just one of these
    if MitM_Tag(tree) = "OMF" then
        if MitM_Attributes(tree) = fail
               or not (IsBound(MitM_Attributes(tree).dec)
                       or IsBound(MitM_Attributes(tree).hex)) then
            return "OMF objects must have either the dec or the hex attribute";
        elif IsBound(MitM_Attributes(tree).dec) and IsBound(MitM_Attributes(tree).hex) then
            return "OMF objects cannot have both the dec and the hex attribute";
        fi;
    fi;

    # Validate the content
    if MitM_Content(tree) <> fail then
        result := MitM_ValidCont.(MitM_Tag(tree))(MitM_Content(tree));
    else
        result := MitM_ValidCont.(MitM_Tag(tree))([]);
    fi;
    if result <> true then
        return Concatenation(MitM_Tag(tree), " contents: ", result);
    fi;
    return true;
end);

max_stars_repo_name: markuspf/MathsInTheMiddle

max_issues_repo_name: markuspf/MathsInTheMiddle

max_forks_repo_name: markuspf/MathsInTheMiddle


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

Class(HCOLSumsGen, DefaultSumsGen, rec(
    OLCompose := (self, o, opts) >> Cond(o.isPermutation(), Prm(ApplyFunc(fCompose, List(Reversed(o.children()), c -> self(c, opts).func))),
	    OLCompose(Map(o.children(), child -> self(child, opts)))),
    Reduction := (self, o, opts) >> o,
    Induction := (self, o, opts) >> o,
    ScalarProd := (self, o, opts) >> let(i := Ind(o.N), x := var.fresh_t("r", TDouble),
        OLCompose(
            Reduction(o.N, (a,b)->add(a,b), V(0.0), False),
            PointWise(o.N, Lambda([x,i], mul(x,nth(o.scalars, i)))))),
    PointWise := (self, o, opts) >> o,
	COND := (self, o, opts) >> let(  #this is a recursive call.....must include check
		test := o.cond,
		ch := o.children(),
		COND(self(test, opts), self(ch[1], opts), self(ch[2], opts))
		),
    IterVStack := (self, o, opts) >> let(
    	ch := o.children()[1],
    	ISumUnion(o.var, o.domain,
	  OLCompose(
	    ScatHUnion(Rows(o), Rows(ch), o.var*Rows(ch), 1),
	    self(ch, opts)))),
    VStack := (self, o, opts) >> let(
      ch := o.children(),
      ApplyFunc(SUMUnion, List([1..Length(ch)], i->
	     OLCompose(
	       ScatHUnion(Rows(o), Rows(ch[i]),  Sum(List(ch{[1..i-1]}, child->child.dims()[1])), 1),
	       self(ch[i], opts)
	     )
      ))
    ),
    I := (self, o, opts) >> let (i := Ind(o.N()),
        ISumUnion(i, i.range, OLCompose(
          ScatHUnion(o.N(), 1, i, 1),
          GathH(o.N(), 1, i, o.N())
	))),
    IterDirectSum := (self, o, opts) >> let(
        ch := o.children(),
	  ISumUnion(o.var, o.domain, #o.var.range,
	    ApplyFunc(SUMUnion, List([1..Length(ch)], i->OLCompose(
              ScatHUnion(Rows(o), Rows(ch[i]), ch[1].dims()[1]*o.var, 1),
              self(ch[i], opts),
              GathH(Cols(o), Cols(ch[i]), ch[1].dims()[2]*o.var, 1)) )
	    )
	  )
        ),
    DirectSum := (self, o, opts) >> let(
        ch := o.children(),
        ApplyFunc(SUMUnion, List([1..Length(ch)], i->OLCompose(
            ScatHUnion(Rows(o), Rows(ch[i]), Sum(List(ch{[1..i-1]}, child->child.dims()[1])), 1),
            self(ch[i], opts),
            GathH(Cols(o), Cols(ch[i]), Sum(List(ch{[1..i-1]}, child->child.dims()[2])), 1))))),
    BinOp := (self, o, opts) >> When(o.N=1, o, let(i := Ind(o.N),
        ISumUnion(i, i.range, OLCompose(
          ScatHUnion(o.N, 1, i, 1),
          BinOp(1, o.op),
          GathH(2*o.N, 2, i, o.N)
        )))),
    Constant := (self, o, opts) >> o #Constant(o.value)
));

max_stars_repo_name: spiral-software/spiral-package-hcol

max_issues_repo_name: spiral-software/spiral-package-hcol

max_forks_repo_name: spiral-software/spiral-package-hcol


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

Declare(DWTper);

#F DWTper(<n>, <j>, [<h(z)>,<g(z)>])
#F DWTper(<n>, <j>, <L>, <V>)
#F   returns a 2-channel <j>-stage discrete wavelet transform with periodic extensions
#F   (circulant transforms) where the
#F   low pass and the high pass analysis filters given by polynomials
#F   [<h(z)>,<g(z)>]
#F
#F   <L> = [<Lh>,<Lg>] - coefficient lists for h and g
#F   <V> = [<vh>,<vg>] - valuation for h and g
#F
#F  <n> - size of the output of the reconstructed sequence
#F  <j> - number of filter bank stages
Class(DWTper, NonTerminal, rec(

  abbrevs := [
    function(n,j,M)
    local S,L,V, j;
    if IsList(M[1]) then
     S:=List([1,2], j->FillZeros([M[1][j],M[2][j]]));
     V:=List(S, j->j[2]);
     L:=List(S, j->j[1]);
     return [n,j,L,V];
    else
    Checked(IsPosInt(n), Checked(ForAll(M,i->IsPolynomial(i)),true));
      S:=List(M, i-> FillZeros(i));
      V:=List(S, i->i[2]);
      L:=List(S, i->i[1]);
    return [n,j,L,V];
    fi;
    end,
    function(n,j,L,V)
    local S,V,L,i,j;
     S:=List([1,2], j->FillZeros([L[j],V[j]]));
     V:=List(S, j->j[2]);
     L:=List(S, j->j[1]);
     return [n,j,L,V];
    end
  ],

  dims := self >> [self.params[1], self.params[1]],

  half1 := self >> let(
      n  := self.params[1],
      l1 := self.params[3][1],
      v1 := self.params[4][1],
      DownSample(n,2,0) * Circulant(n,l1,v1).terminate() ),

  half2 := self >> let(
      n  := self.params[1],
      l2 := self.params[3][2],
      v2 := self.params[4][2],
      DownSample(n,2,0) * Circulant(n,l2,v2).terminate() ),

  terminate := self >> let(
      n := self.params[1],
      j := self.params[2],
      L := self.params[3],
      V := self.params[4],
      # DWT as filter bank stage + downsampling
      res := Cond(j = 1, I(n),
	          DirectSum(DWTper(n/2,j-1,L,V).terminate(), I(n/2))) *
             VStack(self.half1(), self.half2()),
      When(self.transposed, res.transpose(), res)
  ),

  LiftingScheme := self >> self.lifting(),

  lifting := meth(self)
               local n,j,LS,h,g,he,ho,ge,go;
               n := self.params[1];
               j := self.params[2];
               h := DownsampleTwo([self.params[3][1],self.params[4][1]]);
               g := DownsampleTwo([self.params[3][2],self.params[4][2]]);
               # Lifting scheme does not converge for wavelets longer than 9
               # (need to be fixed)
               if (Maximum(List(h, i-> Length(i[1])))>9 or
                   Maximum(List(g, i-> Length(i[1])))>9) then return [[]];fi;
               he := Poly(h[1][1],h[1][2]);
               ho := Poly(h[2][1],h[2][2]);
               ge := Poly(g[1][1],g[1][2]);
               go := Poly(g[2][1],g[2][2]);
               LS := LiftingScheme([[he,ho],[ge,go]]);
               return LS;
             end,

  isReal := True,
));

#F RuleFilter_DWT: (base case) DWT -> Mat,
#F   Computes filter by definition
#F
RulesFor(DWTper, rec(
    #F RuleFilt_Mallat_2:
    #F
    #F DWTper(n,j,[h,g])->DWTper(n/2,j-1,[h,g]),DWTper(n,1,[h,g])
    #F
    #F Mallat rule recursive (single stage DWT)
    #F The I(n/2) is inefficient because it copies the output
    #F However, this rule allows application of Lifting, POlyphase, etc.
    #F on DWTper(n,1,[h,g])
    #F
    DWTper_Mallat_2 := rec(
	info             := "DWTper(n,j,[h,g])->DWTper(n,j-1,[h,g])",
	forTransposition := false,
	isApplicable     := P -> P[1] > 2 and P[1] mod 2 =0 and P[2] > 1,

	allChildren := P -> [[ DWTper(P[1]/2,P[2]-1,P[3],P[4]), DWTper(P[1],1,P[3],P[4]) ]],

	rule := (P, C) -> DirectSum(C[1],I(P[1]/2))*C[2]
    ),

    #F DWTper_Mallat:
    #F
    #F DWTper(n,j,[h,g])->DWTper(n/2,j-1,[h,g]), DSCirculant
    #F
    DWTper_Mallat := rec(
        info             := "DWTper(n,j,[h,g])->DWTper(n/2,j-1,[h,g])",
        forTransposition := false,
        isApplicable     := P -> P[1]>2 and P[1] mod 2 = 0,

        allChildren := P -> let(
	    n        := P[1],
	    dcirc1 := DSCirculant(n, P[3][1], P[4][1], 2, 0),
	    dcirc2 := DSCirculant(n, P[3][2], P[4][2], 2, 0),
	    When(P[2]=1, [[ dcirc1, dcirc2 ]],
		         [[ dcirc1, dcirc2, DWTper(P[1]/2, P[2]-1, P[3], P[4]) ]])),

         rule := (P, C) -> When(P[2]=1, VStack(C[1], C[2]),
                                        VStack(C[3]*C[1], C[2]))
    ),

    #F DWTper_Polyphase:
    #F
    #F DWTper(n,1,[h,g])->[ [Circulant(he), Circulant(ho)],
    #F                      [Circulant(ge), Circulant(go)] ]
    #F
    #F Single-stage periodic DWT into a matrix of circulants of downsampled filters
    #F
    DWTper_Polyphase := rec(
	info             := "DWTper(n,1,[h,g]) -> [[Circ(he), Circ(ho)], [Circ(ge), Circ(go)]]",
	forTransposition := false,
	isApplicable     := P -> P[1]>2 and P[1] mod 2 =0 and P[2]=1,
	allChildren := function(P)
	    local n,h,g,he,ho,ge,go;
	    n := P[1];
	    h := DownsampleTwo([P[3][1], P[4][1]]);
	    g := DownsampleTwo([P[3][2], P[4][2]]);
	    he := Circulant(n/2, h[1][1], h[1][2]);
	    ho := Circulant(n/2, h[2][1], h[2][2]);
	    ge := Circulant(n/2, g[1][1], g[1][2]);
	    go := Circulant(n/2, g[2][1], g[2][2]);
	    return [[ he, ho, ge, go ]];
	end,
	rule := (P, C) -> BlockMat( [[ C[1], C[2] ],
		                     [ C[3], C[4] ]] ) * L(P[1],2)
    ),

    #F RuleDWTper_Lifting:
    #F
    DWTper_Lifting := rec(
        info             := "DWTper(n,1,[h,g]) -> Lifting steps",
        forTransposition := false,
        isApplicable     := ( L ) -> L[1]>2 and L[1] mod 2 =0 and L[2]=1,

        allChildren := function ( P )
            local n,j,LS,Lc,scheme,step,pol;
            n := P[1];
            LS := Copy(HashLookupWav(HashTableWavelets, [P[3],P[4]]));
            Lc := List(LS, scheme->
                     List(scheme{[2..Length(scheme)]}, step->
                        let(pol := FillZeros(ListPoly(step)),
                            Circulant(n/2,pol[1],pol[2]))));

            for i in [1..Length(LS)] do
             # attach the indicator of the type of the first liftings step
              Lc[i][1].lift :=LS[i][1];

             # fuse in the constants/shifts in the last lifting step
             #if (LS[i][1]=0 or LS[i][1]=-2) then LS[i][4]:=LS[i][4]*LS[i][3];
             #else  LS[i][4]:=LS[i][4]*LS[i][2];
             #fi;

             pol := FillZeros(ListPoly(LS[i][4]));
             Lc[i][3] := Circulant(n/2,pol[1],pol[2]);
            od;
            return Lc;
        end,

        rule := function ( P, C, Nonterms )
            local n, i, ind, b, l, M, first, last, ind0, ind1, lstep, last_ls;
            n := P[1];
            b := Nonterms[1].lift;
            l:=Length(C);

	    ind0 := fTensor(fBase(2,0), fId(Rows(C[3])));
	    ind1 := fTensor(fBase(2,1), fId(Rows(C[3])));
	    lstep := (f,b) -> When(b=0,
		LStep(Scat(ind0) * f * Gath(ind1)),
		LStep(Scat(ind1) * f * Gath(ind0)));

            M := When(b < 0,
		SUM(Scat(ind0)*C[2]*Gath(ind1), Scat(ind1)*C[1]*Gath(ind0)),
		SUM(Scat(ind0)*C[2]*Gath(ind0), Scat(ind1)*C[1]*Gath(ind1)));

	    b := (b+2) mod 2; # make b positive
	    for i in [1 .. l-2] do
	        M := M * lstep(C[i+2], b);
		b := (b+1) mod 2;
            od;

	    first := M.child(1);
	    M := Inplace(Compose(Drop(M.children(),1)));
	    return first * M * L(n,2);
       end
    )
));


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

#P SPL Parametrized Symbols
#P ========================
#P
#P This package contains parametrized symbols - shortcut notation for common
#P idiomatic constructs.
#P

#F -----------------------------------------------------------------------------
#F I(N) : NxN identity matrix
#F -----------------------------------------------------------------------------
Class(I, Sym, rec(
    abbrevs := [ (x,y) -> Checked(x=y, IsPosInt0Sym(x), x),
                 (n)   -> Checked(IsPosInt0Sym(n), n) ],
    def := n -> Perm((), n),
    transpose := self >> self,
    conjTranspose := self >> self,
    inverse := self >> self,
    isIdentity := True,
    isSymmetric := True,
    doNotMeasure := true,
    printlatex := self >> Print(" \\one_{",self.params[1],"} "),
    area := self >> self.params[1],
    normalizedArithCost := self >> 0,
));

#F -----------------------------------------------------------------------------
#F O(M,N) : MxN matrix of zeros
#F -----------------------------------------------------------------------------
Class(O, BaseMat, rec(
    abbrevs := [ n -> Checked(IsPosInt(n), [n,n]) ],

    new := (self, r, c) >> SPL(WithBases(self, rec(
        params := [r, c], dimensions := [r, c]))),

    isReal := True,
    isPermutation := False,
    arithmeticCost := (self, costMul, costAddMul) >> costMul(0) - costMul(0),

    transpose     := self >> ObjId(self)(self.params[2], self.params[1]),
    conjTranspose := self >> ObjId(self)(self.params[2], self.params[1]),

    area := self >> Maximum(self.dimensions),

    toAMat := self >> AMatMat(NullMat(
	EvalScalar(self.params[1]), EvalScalar(self.params[2]))),

    print := Sym.print,
    dims  := self >> [ self.params[1], self.params[2] ],
    rChildren := Sym.rChildren,
    rSetChild := Sym.rSetChild
));

#F -----------------------------------------------------------------------------
#F RI(M,N) : MxN rectangular identity
#F           MIN(M,N) x MIN(M,N) identity padded with 0s
#F -----------------------------------------------------------------------------
Class(RI, Sym, rec(
    def := (r, c) -> Checked(IsPosInt(r), IsPosInt(c),
    Cond(r = c, I(r),
         r < c, Gath(fAdd(c,r,0)),
                Scat(fAdd(r,c,0)))),

    transpose := self >> ApplyFunc(self.__bases__[1], Reversed(self.params)),
    conjTranspose := self >> ApplyFunc(self.__bases__[1], Reversed(self.params)),
));

#F -----------------------------------------------------------------------------
#F F(N) - NxN Discrete Fourier Transform matrix
#F -----------------------------------------------------------------------------
Class(F, Sym, rec(
    def := size -> Checked(IsPosInt(size),
    Cond(size = 1, Mat([[1]]),
         size = 2, Mat([[1,1], [1,-1]]),
         Mat(Global.DFT(size)))),

    isReal    := self >> self.params[1] <= 2,
    isPermutation := False,
    transpose := self >> self,
    conjTranspose := self >> self,
    inverse := self >> let(n:=self.params[1], Cond(n=1, self, n=2, 1/2*F(2), Error("Inverse not supported"))),
    toAMat    := self >> DFTAMat(self.params[1]),
    printlatex := (self) >> Print(" F_{", self.params[1], "} ")
));

#F -----------------------------------------------------------------------------
#F R(a)     : 2x2 rotation matrix with angle a*pi, namely
#F
#F    cos(a) sin(a)
#F   -sin(a) cos(a)
#F
#F alternatively:
#F
#F R(c,s)   : specified by a pair such as c^2 + s^2 = 1
#F     <c>  <s>
#F    -<s>  <c>
#F
#F -----------------------------------------------------------------------------
Class(Rot, Sym, rec(
    abbrevs := [ (a) -> [cospi(a), sinpi(a)] ],

    def := (a,b) ->
    Mat( [[a, b], [-b, a]] ),

    transpose := self >>
        self.__bases__[1](self.params[1], -self.params[2]),

    conjTranspose := self >> self.transpose(),
    inverse := self >> self.transpose(),

    toAMat := self >> let(a := self.params[1], b := self.params[2],
        When(IsScalar(a) and IsScalar(b) and ScalarIsCos(a) and ScalarIsSin(b),
         RotationAMat(EvalScalar(ScalarCosArg(a))),
         AMatMat([[EvalScalar(a),  EvalScalar(b)],
              [EvalScalar(-b), EvalScalar(a)]]))),

    #-----------------------------------------------------------------------
    # Expansions
    expandLifting := self >>
        let(c := self.params[1],
            s := self.params[2],
            ev := (not IsScalar(c) or CanBeFullyEvaluated(c.val)) and
                  (not IsScalar(s) or CanBeFullyEvaluated(s.val)),
            Cond(ev and c = s,  c * F(2) * Perm((1,2), 2),
                 ev and c = -s, c * Perm((1,2), 2) * F(2),
                 Mat([[1,1,0], [0,-1,1]]) *
                 Diag([c-s, s, c+s]) *
                 Mat([[1,0], [1,1], [0,1]]))),

    expandWinograd:= self >>
        let(c := self.params[1],
            s := self.params[2],
            ev := (not IsScalar(c) or CanBeFullyEvaluated(c.val)) and
                  (not IsScalar(s) or CanBeFullyEvaluated(s.val)),
            Cond(ev and c = s,  c * F(2) * Perm((1,2), 2),
                 ev and c = -s, c * Perm((1,2), 2) * F(2),
                 Mat([[1, (1-c)/s], [0, 1]]) *
                 Mat([[1, 0], [-s, 1]]) *
                 Mat([[1, (1-c)/s], [0, 1]]))),

    expandDef := self >> self.obj,
));

Declare(toeplitz);

#F -----------------------------------------------------------------------------
#F toeplitz(elements) : toeplitz matrix given by elements of the first row and
#F                      column, starting from the upper right corner
#F -----------------------------------------------------------------------------
#F NOTE: use generating function not list (see 'Toeplitz')
Class(toeplitz, Sym, rec(
    abbrevs := [ arg -> When(Length(arg) > 1 or not IsList(arg[1]), [arg], arg) ],

    def := function(elements)
        local n;
        Constraint(IsList(elements) and Length(elements) > 0);
        Constraint(IsOddInt(Length(elements)));
        DoForAll(elements, x->Constraint(IsScalarOrNum(x)));
        n := (Length(elements) + 1) / 2;
        return Mat(ToeplitzMat(elements));
    end,

    transpose := self >> toeplitz(Reversed(self.params[1])),
    conjTranspose := self >> toeplitz(List(Reversed(self.params[1]), x->Global.Conjugate(x))),
    isReal := True,
    isPermutation := False,
    area := self >> Product(self.dimensions)
));


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

#############################################################################
##
##  AG_InverseLessThanForLetters(<w1>, <w2>)
##
##  Compares w1 and w2 according to lexicografic ordering given by
##  x1 < x1^-1 < x2 < x2^-1 < ...
##
BindGlobal("AG_InverseLessThanForLetters",
function(w1, w2)
  local i, er1, er2;

  if Length(w1) <> Length(w2) then
    return Length(w1) < Length(w2);
  fi;

  er1 := LetterRepAssocWord(w1);
  er2 := LetterRepAssocWord(w2);
  for i in [1..Length(er1)] do
    if AbsInt(er1[i]) <> AbsInt(er2[i]) then
      return AbsInt(er1[i]) < AbsInt(er2[i]);
    fi;
    if er1[i] <> er2[i] then
      return er1[i] > er2[i];
    fi;
  od;

  return false;
end);

#############################################################################
##
##  AG_ReducedListOfWordsByNielsen(<words_list>)
##  AG_ReducedListOfWordsByNielsenBack(<words_list>)
##  AG_ReducedListOfWordsByNielsen(<words_list>, <string>)
##  AG_ReducedListOfWordsByNielsenBack(<words_list>, <string>)
##
InstallMethod(AG_ReducedListOfWordsByNielsen, [IsAssocWordCollection],
function(words)
  return AG_ReducedListOfWordsByNielsen(words, \<);
end);

InstallMethod(AG_ReducedListOfWordsByNielsenBack, [IsAssocWordCollection],
function(words)
  return AG_ReducedListOfWordsByNielsen(words, \<);
end);

InstallMethod(AG_ReducedListOfWordsByNielsen, [IsAssocWordCollection, IsString],
function(words, string)
  return AG_ReducedListOfWordsByNielsen(words, AG_InverseLessThanForLetters);
end);

InstallMethod(AG_ReducedListOfWordsByNielsenBack, [IsAssocWordCollection, IsString],
function(words, string)
  return AG_ReducedListOfWordsByNielsen(words, AG_InverseLessThanForLetters);
end);

#############################################################################
##
#M  AG_ReducedListOfWordsByNielsen(<words_list>, <lt>)
##
InstallMethod( AG_ReducedListOfWordsByNielsen,
                    [IsAssocWordCollection, IsFunction],
function(words_list, lt)
  local result, transform, did_something, n, i, j, try_again, tmp;

  n := Length(words_list);
  result := ShallowCopy(words_list);
  transform := ShallowCopy(FreeGeneratorsOfFpGroup(FreeGroup(n)));
  did_something := false;
  try_again := true;

  for i in [1..n] do
    if not IsAssocWord(result[i]) then
      Print("error in AG_ReducedListOfWordsByNielsen(IsAssocWordCollection, IsFunction):\n");
      Print("  ", i, "-th element of list is not an associative word\n");
      return fail;
    fi;
  od;

  while try_again do
    try_again := false;

    for i in [1..n] do
    for j in [1..n] do

      if i = j then
        if lt(result[i]^-1, result[i]) then
          result[i] := result[i]^-1;
          transform[i] := transform[i]^-1;
          did_something := true;
          try_again := true;
        fi;
        continue;
      fi;

      if i > j and lt(result[i], result[j]) then
        tmp := result[i];
        result[i] := result[j];
        result[j] := tmp;
        tmp := transform[i];
        transform[i] := transform[j];
        transform[j] := tmp;
        did_something := true;
        try_again := true;
      fi;

      if lt(result[i]*result[j], result[i]) then
        result[i] := result[i]*result[j];
        transform[i] := transform[i]*transform[j];
        did_something := true;
        try_again := true;
      fi;
      if lt(result[i]*result[j], result[j]) then
        result[j] := result[i]*result[j];
        transform[j] := transform[i]*transform[j];
        did_something := true;
        try_again := true;
      fi;

      if lt(result[i]^-1*result[j], result[i]) then
        result[i] := result[i]^-1*result[j];
        transform[i] := transform[i]^-1*transform[j];
        did_something := true;
        try_again := true;
      fi;
      if lt(result[i]^-1*result[j], result[j]) then
        result[j] := result[i]^-1*result[j];
        transform[j] := transform[i]^-1*transform[j];
        did_something := true;
        try_again := true;
      fi;

      if lt(result[i]*result[j]^-1, result[i]) then
        result[i] := result[i]*result[j]^-1;
        transform[i] := transform[i]*transform[j]^-1;
        did_something := true;
        try_again := true;
      fi;
      if lt(result[i]*result[j]^-1, result[j]) then
        result[j] := result[i]*result[j]^-1;
        transform[j] := transform[i]*transform[j]^-1;
        did_something := true;
        try_again := true;
      fi;

      if lt(result[i]^-1*result[j]^-1, result[i]) then
        result[i] := result[i]^-1*result[j]^-1;
        transform[i] := transform[i]^-1*transform[j]^-1;
        did_something := true;
        try_again := true;
      fi;
      if lt(result[i]^-1*result[j]^-1, result[j]) then
        result[j] := result[i]^-1*result[j]^-1;
        transform[j] := transform[i]^-1*transform[j]^-1;
        did_something := true;
        try_again := true;
      fi;
    od;
    od;
  od;

  return [result, transform, did_something];
end);

#############################################################################
##
#M  AG_ReducedListOfWordsByNielsenBack(<words_list>, <lt>)
##
InstallMethod(AG_ReducedListOfWordsByNielsenBack,
                   [IsAssocWordCollection, IsFunction],
function(words_list, lt)
  local result, transform, did_something, n, i, j, try_again, tmp;

  n := Length(words_list);
  result := ShallowCopy(words_list);
  transform := ShallowCopy(FreeGeneratorsOfFpGroup(FreeGroup(n)));
  did_something := false;
  try_again := true;

  for i in [1..n] do
    if not IsAssocWord(result[i]) then
      Print("error in AG_ReducedListOfWordsByNielsenBack(IsAssocWordCollection, IsFunction):\n");
      Print("  ", i, "-th element of list is not an associative word\n");
      return fail;
    fi;
  od;

  while try_again do
    try_again := false;

    for i in [1..n] do
    for j in [1..n] do

      if i = j then
        if lt(result[i]^-1, result[i]) then
          result[i] := result[i]^-1;
          transform[i] := transform[i]^-1;
          did_something := true;
          try_again := true;
        fi;
        continue;
      fi;

      if i > j and lt(result[i], result[j]) then
        tmp := result[i];
        result[i] := result[j];
        result[j] := tmp;
        tmp := transform[i];
        transform[i] := transform[j];
        transform[j] := tmp;
        did_something := true;
        try_again := true;
      fi;

      if lt(result[i]^-1*result[j]^-1, result[j]) then
        result[j] := result[i]^-1*result[j]^-1;
        transform[j] := transform[i]^-1*transform[j]^-1;
        did_something := true;
        try_again := true;
      fi;
      if lt(result[i]^-1*result[j]^-1, result[i]) then
        result[i] := result[i]^-1*result[j]^-1;
        transform[i] := transform[i]^-1*transform[j]^-1;
        did_something := true;
        try_again := true;
      fi;

      if lt(result[i]*result[j]^-1, result[j]) then
        result[j] := result[i]*result[j]^-1;
        transform[j] := transform[i]*transform[j]^-1;
        did_something := true;
        try_again := true;
      fi;
      if lt(result[i]*result[j]^-1, result[i]) then
        result[i] := result[i]*result[j]^-1;
        transform[i] := transform[i]*transform[j]^-1;
        did_something := true;
        try_again := true;
      fi;

      if lt(result[i]^-1*result[j], result[i]) then
        result[i] := result[i]^-1*result[j];
        transform[i] := transform[i]^-1*transform[j];
        did_something := true;
        try_again := true;
      fi;
      if lt(result[i]^-1*result[j], result[j]) then
        result[j] := result[i]^-1*result[j];
        transform[j] := transform[i]^-1*transform[j];
        did_something := true;
        try_again := true;
      fi;

      if lt(result[i]*result[j], result[j]) then
        result[j] := result[i]*result[j];
        transform[j] := transform[i]*transform[j];
        did_something := true;
        try_again := true;
      fi;
      if lt(result[i]*result[j], result[i]) then
        result[i] := result[i]*result[j];
        transform[i] := transform[i]*transform[j];
        did_something := true;
        try_again := true;
      fi;

    od;
    od;

  od;

  return [result, transform, did_something];
end);

#############################################################################
##
##  AG_ComputeMihailovaSystemPairs(<pairs_list>)
##
InstallGlobalFunction(AG_ComputeMihailovaSystemPairs,
function(pairs)
  local result, i, nie, m, n, w, tmp,
        did_smth, npairs, transform,
        generate_full_group, nielsen_mihaylov, nielsen_low, rank,
        number_of_letters;

  if not IsDenseList(pairs) then
    Print("error in AG_ComputeMihailovaSystemPairs:  \n");
    Print("  argument is not an IsDenseList\n");
    return fail;
  fi;
  if not IsList(pairs[1]) then
    Print("error in AG_ComputeMihailovaSystemPairs:  \n");
    Print("  first element of list is not an IsList\n");
    return fail;
  fi;
  if Length(pairs[1]) <> 2 then
    Print("error in AG_ComputeMihailovaSystemPairs:  \n");
    Print("  can work only with pairs\n");
    return fail;
  fi;
  if not IsAssocWord(pairs[1][1]) then
    Print("error in AG_ComputeMihailovaSystemPairs:  \n");
    Print("  <arg>[1][1] is not IsAssocWord\n");
    return fail;
  fi;

  #############################################################################
  ##
  ##  generate_full_group
  ##
  generate_full_group := function(list, rank)
    local nie, i;
    nie := AG_ReducedListOfWordsByNielsen(list)[1];
    if Length(Difference(nie, [One(nie[1])])) <> rank then
      return false;
    fi;
    for i in [1..Length(nie)] do
      if Length(nie[i]) > 1 then return false; fi;
    od;
    return true;
  end;

  #############################################################################
  ##
  ##  rank
  ##
  rank := function(words)
    return Length(Difference(AG_ReducedListOfWordsByNielsen(words)[1],
                             [One(words[1])]));
  end;

  #############################################################################
  ##
  ##  number_of_letters
  ##
  number_of_letters := function(list)
    local letters, i, j;
    letters := [];
    for i in [1..Length(list)] do
      for j in [1..NumberSyllables(list[i])] do
        AddSet(letters, GeneratorSyllable(list[i], j));
      od;
    od;
    return Length(letters);
  end;

  #############################################################################
  ##
  ##  nielsen_mihaylov
  ##
  nielsen_mihaylov := function(words_list, m, n)
    local result, transform, did_something, try_again, nie, i, j, tf, pair,
          good_tf, good_pair, tmp;

    result := StructuralCopy(words_list);
    transform := ShallowCopy(FreeGeneratorsOfFpGroup(FreeGroup(m+n)));
    did_something := false;
    try_again := true;

    while try_again do
      try_again := false;

      nie := nielsen_low(result, m, n, \<);
      if nie[3] then
        did_something := true;
        try_again := true;
        result := nie[1];
        transform := AG_CalculateWords(nie[2], transform);
      fi;

      nie := AG_ReducedListOfWordsByNielsen(result{[m+1..m+n]});
      if nie[3] then
        did_something := true;
        try_again := true;
        result := Concatenation(result{[1..m]}, nie[1]);
        transform := Concatenation( transform{[1..m]},
                                    AG_CalculateWords(nie[2],
                                    transform{[m+1..m+n]}));
      fi;

      if rank(result{[m+1..m+n]}) = n then
        if List(result{[m+1..m+n]}, w -> Length(w)) = List([1..n], i -> 1) then
          ## ok
          try_again := false;
        else
          ##  try to minimize sum of lengths
          good_pair := false;
          for pair in ListX([m+1..m+n], [1..m], function(i,j) return [i,j]; end) do
            good_tf := false;
            for tf in [  [1,1,2,1],[2,1,1,1],[1,-1,2,1],[2,-1,1,1],
                        [1,1,2,-1],[2,1,1,-1],[1,-1,2,-1],[2,-1,1,-1]  ] do
              tmp := StructuralCopy(result);
              tmp[pair[1]] := tmp[pair[tf[1]]]^tf[2] * tmp[pair[tf[3]]]^tf[4];
              if  rank(tmp{[m+1..m+n]}) = n and
                  number_of_letters(tmp{[m+1..m+n]}) =
                    number_of_letters(result{[m+1..m+n]}) and
                  Sum(List(tmp{[m+1..m+n]}, w -> Length(w))) <
                    Sum(List(result{[m+1..m+n]}, w -> Length(w)))
              then
                good_tf := true;
                break;
              fi;
            od;
            if good_tf then
              good_pair := true;
              break;
            fi;
          od;
          if not good_pair then
            ##  give up
            return [result, transform, did_something];
          else
            result[pair[1]] :=  result[pair[tf[1]]]^tf[2] *
                                result[pair[tf[3]]]^tf[4];
            transform[pair[1]] := transform[pair[tf[1]]]^tf[2] *
                                  transform[pair[tf[3]]]^tf[4];
            try_again := true;
            did_something := true;
          fi;
        fi;
      else
        ##  try to make rank bigger
        for i in [1..m] do
          good_tf := false;
          pair := [m+1, i];
          for tf in [  [1,1,2,1],[2,1,1,1],[1,-1,2,1],[2,-1,1,1],
                      [1,1,2,-1],[2,1,1,-1],[1,-1,2,-1],[2,-1,1,-1]  ] do
            tmp := StructuralCopy(result);
            tmp[pair[1]] := tmp[pair[tf[1]]]^tf[2] * tmp[pair[tf[3]]]^tf[4];
            if   rank(tmp{[m+1..m+n]}) > rank(result{[m+1..m+n]}) and
                number_of_letters(tmp{[m+1..m+n]}) >=
                  number_of_letters(result{[m+1..m+n]})
            then
              good_tf := true;
              break;
            fi;
          od;
          if good_tf then
            good_pair := true;
            break;
          fi;
        od;
        if not good_pair then
          ##  give up
          return [result, transform, did_something];
        else
          result[pair[1]] :=  result[pair[tf[1]]]^tf[2] *
                              result[pair[tf[3]]]^tf[4];
          transform[pair[1]] := transform[pair[tf[1]]]^tf[2] *
                                transform[pair[tf[3]]]^tf[4];
          try_again := true;
          did_something := true;
        fi;
      fi;

    od;

    return [result, transform, did_something];
  end;

  #############################################################################
  ##
  ##  nielsen_low
  ##
  nielsen_low := function(words_list, m, n, lt)
    local result, transform, did_something, i, j, try_again, tmp, nie;

    result := ShallowCopy(words_list);
    transform := ShallowCopy(FreeGeneratorsOfFpGroup(FreeGroup(m+n)));
    did_something := false;
    try_again := true;

    while try_again do
      try_again := false;

      nie := AG_ReducedListOfWordsByNielsen(result{[1..m]});
      if nie[3] then
        result := Concatenation(AG_CalculateWords(nie[2], result{[1..m]}),
                                result{[m+1..m+n]});
        transform := Concatenation( AG_CalculateWords(nie[2], transform{[1..m]}),
                                    transform{[m+1..m+n]} );
        did_something := true;
        try_again := true;
      fi;

      for i in [1..m] do
        for j in [m+1..m+n] do
          if lt(result[i]^result[j], result[i]) then
            result[i] := result[i]^result[j];
            transform[i] := transform[i]^transform[j];
            did_something := true;
            try_again := true;
          fi;
          if lt(result[i]^(result[j]^-1), result[i]) then
            result[i] := result[i]^(result[j]^-1);
            transform[i] := transform[i]^(transform[j]^-1);
            did_something := true;
            try_again := true;
          fi;
          if lt((result[i]^-1)^result[j], result[i]) then
            result[i] := (result[i]^-1)^result[j];
            transform[i] := (transform[i]^-1)^transform[j];
            did_something := true;
            try_again := true;
          fi;
          if lt((result[i]^-1)^(result[j]^-1), result[i]) then
            result[i] := (result[i]^-1)^(result[j]^-1);
            transform[i] := (transform[i]^-1)^(transform[j]^-1);
            did_something := true;
            try_again := true;
          fi;
        od;
      od;
    od;

    return [result, transform, did_something];
  end;

  #############################################################################
  ##
  ##  MihailovaSystem body
  ##
  n := Length(FreeGeneratorsOfWholeGroup(Group(pairs[1][1])));
  m := Length(pairs) - n;
  npairs := StructuralCopy(pairs);
  transform := StructuralCopy(FreeGeneratorsOfFpGroup(FreeGroup(n+m)));
  did_smth := false;
  if not generate_full_group(List(pairs, p -> p[1]), n)
      or not generate_full_group(List(pairs, p -> p[2]), n)
  then
    Print("error in AG_ComputeMihailovaSystemPairs:  \n");
    Print("  projections do not generate full free group\n");
    return fail;
  fi;

  ##  if rank equals number of pairs then just make one coordinate nicer
  if m = 0 then
    nie := AG_ReducedListOfWordsByNielsen(List(npairs, p -> p[1]), "r");
    if nie[3] then
      tmp := List(npairs, p -> []);
      for i in [1..m+n] do
        tmp[i][1] := AG_CalculateWord(nie[2][i], List(npairs, p -> p[1]));
        tmp[i][2] := AG_CalculateWord(nie[2][i], List(npairs, p -> p[2]));
      od;
      npairs := StructuralCopy(tmp);
      transform := StructuralCopy(nie[2]);
      did_smth := true;
    fi;
    return [npairs, transform, did_smth];
  fi;

  ##  else try to do as much as possible

  ##  1. Apply Nielsen to first coordinate
  nie := AG_ReducedListOfWordsByNielsen(List(npairs, p -> p[1]), "r");
  if nie[3] then
    tmp := StructuralCopy(npairs);
    for i in [1..m+n] do
      tmp[i][1] := AG_CalculateWord(nie[2][i], List(npairs, p -> p[1]));
      tmp[i][2] := AG_CalculateWord(nie[2][i], List(npairs, p -> p[2]));
    od;
    npairs := StructuralCopy(tmp);
    transform := StructuralCopy(nie[2]);
    did_smth := true;
  fi;

  ##  2. Now apply nielsen_mihaylov to the second coordinate
  nie := nielsen_mihaylov(List(npairs, p -> p[2]), m, n);
  if nie[3] then
    tmp := StructuralCopy(npairs);
    for i in [1..m+n] do
      tmp[i][1] := AG_CalculateWord(nie[2][i], List(npairs, p -> p[1]));
      tmp[i][2] := AG_CalculateWord(nie[2][i], List(npairs, p -> p[2]));
    od;
    npairs := StructuralCopy(tmp);
    tmp := StructuralCopy(transform);
    for i in [1..m+n] do
      tmp[i] := AG_CalculateWord(nie[2][i], transform);
    od;
    transform := StructuralCopy(tmp);
    did_smth := true;
  fi;

  ##  3. Try to get nice generators on first coordinate
  nie := AG_ReducedListOfWordsByNielsenBack(List(npairs{[m+1..m+n]}, p -> p[1]), "r");
  if nie[3] then
    tmp := StructuralCopy(npairs);
    for i in [1..n] do
      tmp[m+i][1] := AG_CalculateWord(nie[2][i], List(npairs{[m+1..m+n]}, p -> p[1]));
      tmp[m+i][2] := AG_CalculateWord(nie[2][i], List(npairs{[m+1..m+n]}, p -> p[2]));
    od;
    npairs := StructuralCopy(tmp);
    tmp := StructuralCopy(transform);
    for i in [1..n] do
      tmp[m+i] := AG_CalculateWord(nie[2][i], transform{[m+1..m+n]});
    od;
    transform := StructuralCopy(tmp);
    did_smth := true;
  fi;

  return [npairs, transform, did_smth];
end);

#############################################################################
##
#M  AG_ReducedByNielsen(<words_list>)
##
InstallMethod(AG_ReducedByNielsen,
              "for [IsList and IsAssocWordCollection]",
              [IsList and IsAssocWordCollection],
function(words)
  if AG_Globals.use_inv_order_in_apply_nielsen then
    return AG_ReducedListOfWordsByNielsen(words, "back")[1];
  else
    return AG_ReducedListOfWordsByNielsen(words)[1];
  fi;
end);

#############################################################################
##
#M  AG_ReducedByNielsen(<autom_list>)
##
InstallMethod(AG_ReducedByNielsen,
              "for [IsList and IsAutomCollection]",
              [IsList and IsAutomCollection],
function(automs)
  local words;
  if IsEmpty(automs) then
    return [];
  fi;
  words := AG_ReducedByNielsen(List(automs, a -> a!.word));
  return List(words, w -> Autom(w, automs[1]));
end);

#E


text:
# Here we use generators : the given formula doesn't need one, but the alternate
# non-squares function is better done with a generator.

# The formula is implemented with exact floor(sqrt(n)), so we use
# a trick: multiply by 100 to get the first decimal digit of the
# square root of n, then add 5 (that's 1/2 multiplied by 10).
# Then just divide by 10 to get floor(1/2 + sqrt(n)) exactly.
# It looks weird, but unlike floating point, it will do the job
# for any n.
NonSquaresGen := function()
	local ns, n;
	n := 0;
	ns := function()
		n := n + 1;
		return n + QuoInt(5 + RootInt(100*n), 10);
	end;
	return ns;
end;

NonSquaresAlt := function()
	local ns, n, q, k;
	n := 1;
	q := 4;
	k := 3;
	ns := function()
		n := n + 1;
		if n = q then
			n := n + 1;
			k := k + 2;
			q := q + k;
		fi;
		return n;
	end;
	return ns;
end;

gen := NonSquaresGen();
List([1 .. 22] i -> gen());
# [ 2, 3, 5, 6, 7, 8, 10, 11, 12, 13, 14, 15, 17, 18, 19, 20, 21, 22, 23, 24, 26, 27 ]

a := NonSquaresGen();
b := NonSquaresAlt();

ForAll([1 .. 1000000], i -> a() = b());
# true

ext: gap

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


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

#F Random SPLs
#F -----------
#F

#F RandomPermSPL ( <n> )
#F   returns a random spl of type "perm" or "symbol" representing
#F   a <n> x <n> permutation matrix.
#F

RandomPermSPL := function ( n )
  local type;

  type := Random(["perm", "J", "L"]);
  if type = "perm" then
    return Perm(Random(SymmetricGroup(n)), n);
  elif type = "J" then
    return SPLSymbol("J", n);
  elif type = "L" then
    return SPLSymbol("L", [n, Random(DivisorsInt(n))]);
  else
    Error("should not happen");
  fi;
end;

#F RandomSymbolSPL ( <n> )
#F   returns a random spl of type "symbol" representing
#F   a <n> x <n> matrix.
#F

RandomSymbolSPL := function ( n )
  local type, d;

  if n = 1 then
    type := Random(["I", "J"]);
  elif n = 2 then
    type := Random(["I", "J", "F", "Rot"]);
  elif IsPrime(n) then
    type := Random(["I", "J"]);
  else
    type := Random(["I", "J", "L", "T"]);
  fi;
  if type in ["I", "J", "F"] then
    return SPLSymbol(type, n);
  elif type = "Rot" then
    return SPLSymbol("Rot", Random(Rationals));
  elif type in ["L", "T"] then
    d := Random([2..Tau(n)-1]);
    return SPLSymbol(type, [n, DivisorsInt(n)[d]]);
  else
    Error("should not happen");
  fi;
end;

#F RandomTerminalSPL ( <n> )
#F   returns a random terminal spl representing a <n> x <n> matrix.
#F

RandomTerminalSPL := function ( n )
  local type, randscalar, L, i, P, p;

  # function to return a somewhat random scalar <> 0
  randscalar := function ( )
    local L, den;

    L :=
      Concatenation(
        [-5..-1], [1..5], List([1/3,2/5,3/7], CosPi), List([2..5], E)
      );
    return Random(L)/Random([1..5]);
  end;

  # choose type, mats only for n <= 3, give better odds to
  # composed types
  if n <= 3 then
    type :=
      Random( [
	"mat", "sparse", "diag", "perm", "symbol",
        "compose", "directSum", "tensor",
        "compose", "directSum", "tensor",
        "scalarMultiple", "conjugate"
      ] );
  else
    type :=
      Random( [
	"sparse", "diag", "perm", "symbol",
        "compose", "directSum", "tensor",
        "scalarMultiple", "conjugate"
      ] );
  fi;

  # dispatch on types
  if type = "mat" then

    return Mat(List([1..n], i -> List([1..n], j -> randscalar())));

  elif type = "sparse" then

    # choose n entries
    # make sure there are no doubles
    L := [ ];
    P := [ ];
    for i in [1..n] do
      repeat
        p := [Random([1..n]), Random([1..n])];
      until not p in P;
      AddSet(P, p);
      Add(L, [p[1], p[2], randscalar()]);
    od;
    if ForAll(L, t -> t[1] <> n or t[2] <> n) then

      # add entry in lower right corner
      Add(L, [n, n, randscalar()]);
    fi;
    return Sparse(L);

  elif type = "diag" then

    return Diag(List([1..n], i -> randscalar()));

  elif type = "perm" then

    return RandomPermSPL(n);

  elif type = "symbol" then

    return RandomSymbolSPL(n);

  elif type = "compose" then

    # choose 2 factors
    i := Random([2]);
    return Compose(List([1..i], i -> RandomTerminalSPL(n)));

  elif type = "directSum" then

    if n <= 3 then
      return RandomTerminalSPL(n);
    fi;

    # choose partition into 2 summands
    i := Random([2]);
    return
      DirectSum(List(Random(OrderedPartitions(n, i)), RandomTerminalSPL));

  elif type = "tensor" then

    if IsPrime(n) or n = 1 then
      return RandomTerminalSPL(n);
    else
      i := Random(DivisorsInt(n));
      return Tensor(RandomTerminalSPL(i), RandomTerminalSPL(n/i));
    fi;

  elif type = "scalarMultiple" then

    return Scale(randscalar(), RandomTerminalSPL(n));

  elif type = "conjugate" then

    return Conjugate(RandomTerminalSPL(n), RandomPermSPL(n));

  else
    Error("should not happen");
  fi;
end;


text:
################################################################################
##
#W TY_data.gi              GroupTheoretical Package
##
#W Paul Bruillard, Cesar Galindo, Siu-Hung Ng, Julia Plavnik, Eric Rowell,
#W Zhenghan Wang
##
## Installation file for near_group_data functions of the GroupTheoretical Package
##
#Y Copyright (C) 2016, Battelle Memorial Institute
##
################################################################################

################################################################################
##
#F near_group_data(<group>,<int>). . . . . . compute the data
## . . . . . . . . . . . . . . . . . . . . . for a neargroup category k=0 is TY.
##
InstallGlobalFunction(near_group_data, function(G,k)
  local Simples,m,N,rank,x,FPdimC,FPdim,S,T,g_idx,h_idx,g,h;

  if k < 0 or not(IsInt(k)) then
    Error("k must be a non-negative integer\n");
  fi;

  if k > 0 then

    Error("k > 0 is not currently supported. We need to unravel Siehler's paper 'Near-group categories\n");

    if Order(G) > k+1 then
      Error("near group categories only admit monoidal structure for |G| <= k+1");
    fi;

    if Order(G) = k+1 and not(IsCyclic(G) and Size(Set(FactorsInt(Order(G)+1))) = 1) then
      Error("for |G| = k+1, the near group category only admits monoidal structure if G  is the multiplicative group of a finite field.");
    fi;

    if Order(G) = 1 and ( (k mod 4) = 2 or (k mod 4) = 3) then
      Error("For |G| = 1, the near group category does not admit  monoidal structure if k = 2 or 3 mod 4.");
    fi;
  fi;

  Simples := ShallowCopy(AsSet(G));
  m:="m";
  Add(Simples,m);
  rank:=Size(Simples);
  N:=List([1..rank],x->NullMat(rank,rank));
  # compute the fusion rules
  for g_idx in [1..rank] do
    g:=Simples[g_idx];
    for h_idx in [1..rank] do
      h:=Simples[h_idx];
      if g_idx = rank then
        if h_idx = rank then
          # m\otimes m = Sum_{k\in G}k+km
          N[g_idx][h_idx]:=List([1..rank],x->1);
          N[g_idx][h_idx][rank] := k;
        else
          N[g_idx][h_idx][rank] := 1; # m\otimes h = m
        fi;
      else
        if h_idx = rank then
          N[g_idx][h_idx][rank] := 1; # g\otimes m = m
        else
          # g\otimes h = gh
          N[g_idx][h_idx][Position(Simples,g*h)] := 1;
        fi;
      fi;
    od;
  od;

  FPdimC := 2*Order(G);
  FPdim := List([1..rank],x->1);
  FPdim[rank] := Sqrt(Order(G));
  return [Simples, N, FPdim, FPdimC];
end);

################################################################################
##
#F TY_data(<group>,<bicharacter>,<num>,<list>,<sign>). . . . compute the data for TY(G,chi,tau,delta,Tsign)
##
InstallGlobalFunction(TY_data, function(G,chi,tau,delta,Tsign)
  local Simples,m,N,rank,x,FPdimC,FPdim,S,T,sigma1,sigma3_1,sigma3,sigma3_squared,sigma3_squared_trace,g_idx,h_idx,g,h;

  if not(IsInt(Tsign)) or Tsign^2 <> 1 then
    Error("Tsign must be +/-1");
  fi;

  if Size(delta)<>Order(G)+1 or false in List(delta,x->(IsInt(x) and x^2=1)) then
    Error("delta must be a list of length |G|+1 and consist entirely of +/- 1");
  fi;

  if (1/tau)^2 <> Order(G) then
    Error("tau must be a choice of squareroot of 1/|G|");
  fi;

  Simples := ShallowCopy(AsSet(G));
  m:="m";
  Add(Simples,m);
  rank:=Size(Simples);
  N:=List([1..rank],x->NullMat(rank,rank));
  # compute the fusion rules
  for g_idx in [1..rank] do
    g:=Simples[g_idx];
    for h_idx in [1..rank] do
      h:=Simples[h_idx];
      if g_idx = rank then
        if h_idx = rank then
          # m\otimes m = Sum_{k\in G}k+km
          N[g_idx][h_idx]:=List([1..rank],x->1);
          N[g_idx][h_idx][rank] := 0;
        else
          N[g_idx][h_idx][rank] := 1; # m\otimes h = m
        fi;
      else
        if h_idx = rank then
          N[g_idx][h_idx][rank] := 1; # g\otimes m = m
        else
          # g\otimes h = gh
          N[g_idx][h_idx][Position(Simples,g*h)] := 1;
        fi;
      fi;
    od;
  od;

  FPdimC := 2*Order(G);
  FPdim := List([1..rank],x->1);
  FPdim[rank] := Sqrt(Order(G));

  if not(IsAbelian(G) and Exponent(G)=2) then
    # not braided. Return the fusion rules and simples (arxiv:math/0011037v1)
    Display("G is not an elementary abelian 2-group and so the associated category does not admit a braiding.");
    return [Simples, N, FPdim, FPdimC];
  else
    # the category admits a braiding.
    #
    # need to compute the S-matrix and T-matrix
    #
    # First we compute T
    sigma1:=List([1..rank-1],x->delta[x]*Sqrt(chi[x][x]));
    sigma3_1:=delta[1]*Sqrt(tau*Sum(sigma1));
    T:=List([1..rank],x->0);
    for x in [1..rank-1] do
      T[x] := chi[x][x];
    od;
    T[rank]:=Tsign/sigma3_1;
    S:=NullMat(rank,rank);
    # sigma1 = sigma2 (Seihler section 2.3)
    # since delta[x] = pm1 we have sigma1[g]*sigma2[g] = chi[g][g]
    sigma3:=List([1..rank-1],x->sigma3_1*sigma1[x]*chi[x][x]);
    sigma3_squared:=List(sigma3,x->x^2);
    sigma3_squared_trace:=Sum(sigma3);
    for g in [1..rank-1] do
      for h in [1..rank-1] do
        # group elements are invertible
        #S[g][h] := chi[g][h]*chi[h][g]*FPdim[g]*FPdim[h];
        S[g][h] := chi[g][h]*chi[h][g];
      od;
      #S[g][m] := sigma1[g]*sigma2[g]*FPdim[g]*FPdim[rank];
      # group elements are invertible
      S[g][rank] := chi[g][g]*FPdim[rank];
      S[rank][g] := S[g][rank];
    od;
    S[rank][rank] := sigma3_squared_trace*FPdim[rank]*FPdim[rank];
    Display("VERIFY S and T\n");
    return [Simples,S,T,N,FPdim,FPdimC];
  fi;
end);
#E near_group_data.gi  . . . . . . . . . . . . . . . . . . . . . . . . ends here

max_stars_repo_name: pnnl/GroupTheoretical

max_issues_repo_name: pnnl/GroupTheoretical

max_forks_repo_name: pnnl/GroupTheoretical


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

Import(paradigms.vector, platforms.sse);

#	Tests SIMD functionality
TestSIMD := function()
	if SSE_4x32f.active then
		doSimdDft(5, SSE_4x32f);
		doSimdDft(5, SSE_4x32f, rec(oddSizes:=false, svct:=false));
		doSimdDft([2..20], SSE_4x32f);
	fi;
end;