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

# Transform symbol: Toeplitz
#
# Parameters:       <list>
#
# Definition:       Toeplitz(L) represents the (n X n) toeplitz matrix.
#                   L must contain the coefficients of the first row (from
#                   right to left) followed by all but first coefficents
# 		    on the first column specified. Toeplitz matrices have
#                   constant elements along all diagonals.
#
# 		    n = (len(L)+1) / 2.
#
# Example:          Transform( "Toeplitz", [1, 2, 3, 4, 5, 6, 7]) is the matrix
#                     [ [4, 3, 2, 1],
#                       [5, 4, 3, 2],
#                       [6, 5, 4, 3],
#                       [7, 6, 5, 4] ]
Class(Toeplitz, NonTerminal, DataNonTerminalMixin, rec(
    # above DataNonTerminalMixin must come after NonTerminal to pickup
    # NonTerminal.operations instead of DataNonTerminalMixin.operations
    # which is actually ClassBase.operations
    abbrevs := [
	L -> let(f:=toFunc(L),
	         [ Checked(IsOddInt(f.domain()), f), 0, f.domain()-1 ]),
	(L, l, r) -> let(f:=toFunc(L),
	         [ Checked(IsOddInt(f.domain()), f), l, r ] )
    ],

    dims := self >> Replicate(2, (self.params[1].domain()+1)/2),

    terminate := self >> let(
	coeffs := self.params[1],
	#lcoeffs := coeffs.lambda(),
	coeffs_ev := coeffs.lambda().tolist(),
	len := coeffs.domain(),
	l := self.params[2],
	r := self.params[3],
	#v := Ind(len),
        #toeplitz( Lambda(v, cond(leq(l,v,r), lcoeffs.at(v), 0)).tolist() )),
        toeplitz(coeffs_ev)),

    transpose := self >> let(nums := self.params[1],
        Toeplitz(fCompose(nums, J(nums.domain())))),

    isReal := self >> true,

    HashId := self >> let(n := self.params[1].domain(), l := self.params[2], r := self.params[3],
#	[n, When(IsInt(l) and IsInt(r), QuantizeQuarters(((r-l) mod n)/(n-1)), 1)]),
	n),

    hprint := self >> let(n := self.params[1].domain(), h := self.HashId,
	Print(self.name, "(", n, ")")), #, ", 0, ", ", Int((n-1)*h[2]), ")")),

    SmallRandom  := () -> [1 .. 1+2*Random([1..8])],

    setNonZero := meth(self, l, r)
        self.params[2] := l;
	self.params[3] := r;
    end
));

##
## Rules
##
RulesFor(Toeplitz, rec(
    ###################################################################
    ## Time domain methods
    ###################################################################

    #F Toeplitz_Base: (base case)
    #F
    #F NonTerm_Filter -> SPLMat
    #F Computes convolution by definition
    #F
    Toeplitz_Base := rec(
	info             := "Toeplitz -> Mat",
	forTransposition := false,
	limit            := 8,
	isApplicable     := (self, P) >> let(n := (P[1].domain()+1)/2, n <= self.limit),
	allChildren      := P -> [[ ]],
	rule := (P, C) -> ApplyFunc(Toeplitz, P).terminate()
    ),

    Toeplitz_SmartBase := rec(
	info             := "Toeplitz -> Lower-Triangular | Upper-Triangular | Full matrix",
	forTransposition := false,
	limit            := 5,
	isApplicable     := (self, P) >> let(n := (P[1].domain()+1)/2, n <= self.limit),
	allChildren      := P -> [[ ]],

	rule := (P, C) -> let(
	    coeffs := P[1],
	    coeffs_ev := coeffs.lambda().tolist(),
	    len := coeffs.domain(), # 7 (r 6,5,4,<3>,2,1,0 l)
	    halflen := (len-1) / 2, # 3
	    l := P[2],
	    r := P[3],
	    closed := leq(l, r),

	    A := leq(halflen, l),
	    B := leq(r, halflen),

	    lout := leq(len, l),
	    rout := leq(r, -1),

	    COND(logic_and(A, logic_or(closed, rout)), # lower triangular
		 toeplitz(Concat(Replicate(halflen,0), Drop(coeffs_ev, halflen))),
		 COND(logic_and(B, logic_or(closed, lout)), # upper triangular
		      toeplitz(Concat(Take(coeffs_ev, halflen+1), Replicate(halflen,0))),
		      toeplitz(coeffs_ev)))) # full
    ),

    Toeplitz_PreciseBase := rec(
	info             := "Toeplitz -> toeplitz of non-zero elements only",
	forTransposition := false,
	switch := false,
	isApplicable     := P -> let(len := P[1].domain(), len = 3),
	allChildren      := P -> [[ ]],
	rule := (P, C) -> let(
	    lcoeffs := P[1].lambda(),
	    len := P[1].domain(),
	    l := P[2],
	    r := P[3],
	    v := Ind(len),
	    toeplitz( Lambda(v, cond(leq(l,v,r), lcoeffs.at(v), 0)).tolist() ))
    ),

    #F Toeplitz_Blocking
    #F
    #F Toeplitz -> Blocks (Toeplitzes)
    #F
    #F Toeplitz is divided into square blocks of sizes given by all
    #F proper divisors of the input size.
    #F
    Toeplitz_Blocking := rec(
	info             := "Toeplitz_nk -> Toeplitz_k",
	forTransposition := false,

	isApplicable     := P -> let(N := (P[1].domain() + 1) / 2,
	    N > 2 and not IsPrime(N)),

	allChildren := P -> let(
	    len := P[1].domain(), N := (len+1)/2,
	    l := P[2], r:= P[3], nzlen := r-l+1, ratio := When(IsInt(nzlen), nzlen/len, 1),
	    List(DivisorPairs(N), d -> let(newlen := 2*d[1]-1,
		[ Toeplitz(FUnk(newlen), newlen - CeilingRat(newlen*ratio), newlen-1) ]))),

	rule := (P, C, Nonterms) -> let(
	    len := P[1].domain(),
	    N := (len + 1) / 2,
	    bksize := Rows(C[1]),
	    bklen  := 2*bksize - 1,
	    bkdim := N / bksize,
	    numblocks := 2 * bkdim - 1,
	    i := Ind(bkdim),
	    k := Ind(bkdim),
	    ofs := DataInd(TInt, N*2-bklen).setAttr("live_out"),
	    gfunc := fAdd(len, bklen, ofs),

	    l := P[2],
	    r := P[3],
	    open   := leq(r+1, l),
	    A := leq(l-bklen+1, ofs),
	    B := leq(ofs, r),
	    isNonZero := logic_or(logic_and(A, B), logic_and(open, logic_or(A, B))),

            nt := Nonterms[1].setData(fCompose(P[1], gfunc))
	                     .setNonZero(l-ofs, r-ofs),

	    ISum(i, i.range,
		Scat(fTensor(fBase(i), fId(bksize))) *
		ISumAcc(k, k.range,
		    Data(ofs, bksize*((bkdim-1)+i-k),
		       COND(isNonZero, C[1],
			    O(bksize))) *
		    Gath(fTensor(fBase(k), fId(bksize)))))
	)
    ),

    Toeplitz_BlockingDense := rec(
	info             := "Toeplitz_nk -> Toeplitz_k",
	forTransposition := false,

	isApplicable     := P -> let(N := (P[1].domain() + 1) / 2,
	    N > 2 and not IsPrime(N)),

	allChildren := P -> let(
	    len := P[1].domain(), N := (len+1)/2,
	    l := P[2], r:= P[3], nzlen := r-l+1, ratio := When(IsInt(nzlen), nzlen/len, 1),
	    List(DivisorPairs(N), d -> let(newlen := 2*d[1]-1,
		[ Toeplitz(FUnk(newlen), newlen - CeilingRat(newlen*ratio), newlen-1) ]))),

	rule := (P, C, Nonterms) -> let(
	    len := P[1].domain(),
	    N := (len + 1) / 2,
	    bksize := Rows(C[1]),
	    bklen  := 2*bksize - 1,
	    bkdim := N / bksize,
	    numblocks := 2 * bkdim - 1,
	    i := Ind(bkdim),
	    k := Ind(bkdim),
	    ofs := DataInd(TInt, N*2-bklen).setAttr("live_out"),
	    gfunc := fAdd(len, bklen, ofs),
	    isNonZero := When(P[3]-P[2]+1=len, V(1), leq(P[2]-bklen+1, ofs, P[3])),
            nt := Nonterms[1].setData(fCompose(P[1], gfunc))
	                     .setNonZero(P[2]-ofs, P[3]-ofs),

	    ISum(i, i.range,
		Scat(fTensor(fBase(i), fId(bksize))) *
		ISumAcc(k, k.range,
		    Data(ofs, bksize*((bkdim-1)+i-k),
		       C[1]) *
		    Gath(fTensor(fBase(k), fId(bksize)))))
	)
    ),
    ###################################################################
    ## Frequency domain methods
    ###################################################################

    #F Toeplitz_ExpandedConvolution
    #F
    #F Toeplitz -> Toeplitz_n -> submatrix Circulant(k)  where k>=2*n-1
    #F
    #F Computes Toeplitz matrix through convolution of expanded size.
    #F
    #F Van Loan, C. "Computational Frameworks for the Fast Fourier Transform",
    #F SIAM Philadelphia 1992, p208.
    #F
    Toeplitz_ExpandedConvolution := rec(
	info             := "Toeplitz_n -> submatrix Circulant(k), k>=2*n-1 ",
	forTransposition := false,
	switch           := false,
	isApplicable     := P -> P[1].domain() >= 5,
	allChildren      := function( P )
	    local conv_size,l,L,toep_size;
	    l := P[1].domain();
	    L := P[1].lambda().tolist();
	    conv_size := 2^LogInt(l, 2);
	    while conv_size < l do conv_size := conv_size*2; od;
	    return [[ Circulant(
		Concatenation(
		    List([1..(l+1)/2], x->L[(l+1)/2+x-1]),
		    List([1..conv_size-l], x->0),
		    List([1..(l-1)/2], x->L[x]))) ]];
	end,

	rule := (P, C) -> let(n := (P[1].domain() + 1)/2,
	    RI(n, Rows(C[1])) *
	    C[1] *
	    RI(Rows(C[1]), n)
	)
    )
));


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

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

###############################################################################
##
#M  SelfSimilarGroup(<list>, <bind_vars>)
##
InstallMethod(SelfSimilarGroup, "for [IsList, IsBool]", [IsList, IsBool],
function(list, bind_vars)
  if not AG_IsCorrectRecurList(list, true) then
    Error("in SelfSimilarGroup(IsList, IsBool):\n",
          "  given list is not a correct list representing self-similar group\n");
  fi;

  return GroupOfSelfSimFamily(SelfSimFamily(list, bind_vars));
end);

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

###############################################################################
##
#M  SelfSimilarGroup(<list>, <names>, <bind_vars>)
##
InstallMethod(SelfSimilarGroup,
              "for [IsList, IsList, IsBool]", [IsList, IsList, IsBool],
function(list, names, bind_vars)
  if not AG_IsCorrectRecurList(list, true) then
    Error("error in SelfSimilarGroup(IsList, IsList, IsBool):\n",
          "  given list is not a correct list representing self-similar group\n");
  fi;

  return GroupOfSelfSimFamily(SelfSimFamily(list, names, bind_vars));
end);

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

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

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

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

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

###############################################################################
##
#M  IsGroupOfSelfSimFamily(<G>)
##
InstallMethod(IsGroupOfSelfSimFamily, "for [IsSelfSimGroup]",
              [IsSelfSimGroup],
function(G)
  return G = GroupOfSelfSimFamily(G);
end);

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

###############################################################################
##
#M  __AG_SubgroupOnLevel(<G>, <gens>, <level>)
##
InstallMethod(__AG_SubgroupOnLevel, [IsSelfSimGroup,
                                    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 HasIsGroupOfSelfSimFamily(G) and IsGroupOfSelfSimFamily(G) then
    overgroup := G;
  else
    overgroup := GroupOfSelfSimFamily(UnderlyingSelfSimFamily(G));
  fi;

  return SubgroupNC(overgroup, gens);
end);

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

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

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

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

  return SubgroupNC(overgroup, gens);
end);

InstallMethod(__AG_SimplifyGroupGenerators, "for [IsList and IsInvertibleSelfSimCollection]",
                          [IsList and IsInvertibleSelfSimCollection],
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 IsEmpty(words) then
      return [];
    fi;
    words := FreeGeneratorsOfGroup(Group(words));
  fi;

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

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

###############################################################################
##
#M  Display(<G>)
##
InstallMethod(Display, "for [IsSelfSimGroup]",
              [IsSelfSimGroup],
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 [IsSelfSimGroup]", [IsSelfSimGroup],
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 [IsSelfSimGroup]",
              [IsSelfSimGroup],
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  IsFractalByWords(G)
##
InstallMethod(IsFractalByWords, "for [IsSelfSimGroup]",
              [IsSelfSimGroup],
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 [IsSelfSimGroup]", [IsSelfSimGroup],
function (G)
  local f;
  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 IsSelfSimilarGroup(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 [IsSelfSimGroup and IsSelfSimilar, IsCyclotomic]",
              [IsSelfSimGroup 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 [IsSelfSimGroup and IsSelfSimilar]",
              [IsSelfSimGroup and IsSelfSimilar],
function(G)
  return LevelOfFaithfulAction(G, infinity);
end);

################################################################################
##
#O  IsomorphismPermGroup (<G>)
#O  IsomorphismPermGroup (<G>, <max_lev>)
##
##  For a given finite group <G> generated by initial automata or by elements defined by
##  wreath recursion
##  computes an isomorphism from <G> into a finite permutational group.
##  If <G> is not known to be self-similar (see "IsSelfSimilar") the isomorphism is based on the
##  regular representation, which works generally much slower. If <G> is self-similar
##  there is a level of the tree (see "LevelOfFaithfulAction"), where <G> acts faithfully.
##  The corresponding representation is returned in this case. If <max_lev> is given
##  it finds only the first <max_lev> quotients by stabilizers and if all of them have
##  different size it returns `fail'.
##  If <G> is infinite and <max_lev> is not specified it will loop forever.
##
##  For example, consider a subgroup $\langle a, b\rangle$ of Grigorchuk group.
##  \beginexample
##  gap> Grigorchuk_Group := AutomatonGroup("a=(1,1)(1,2),b=(a,c),c=(a,d),d=(1,b)");
##  < a, b, c, d >
##  gap> f := IsomorphismPermGroup(Group(a, b));
##  MappingByFunction( < a, b >, Group(
##  [ (1,2)(3,5)(4,6)(7,9)(8,10)(11,13)(12,14)(15,17)(16,18)(19,21)(20,22)(23,
##      25)(24,26)(27,29)(28,30)(31,32), (1,3)(2,4)(5,7)(6,8)(9,11)(10,12)(13,
##      15)(14,16)(17,19)(18,20)(21,23)(22,24)(25,27)(26,28)(29,31)(30,32)
##   ]), function( g ) ... end, function( b ) ... end )
##  gap> Size(Image(f));
##  32
##  gap> H := SelfSimilarGroup("a=(a*b,1)(1,2), b=(1,b*a^-1)(1,2), c=(b, a*b)");
##  < a, b, c >
##  gap> f1 := IsomorphismPermGroup(H);
##  MappingByFunction( < a, b, c >, Group([ (1,3)(2,4), (1,3)(2,4), (1,2)
##   ]), function( g ) ... end, function( b ) ... end )
##  gap> Size(Image(f1));
##  8
##  gap> PreImagesRepresentative(f1, (1,3,2,4));
##  a*c
##  gap> (a*c)^f1;
##  (1,3,2,4)
##  \endexample
##
InstallOtherMethod(IsomorphismPermGroup, "for [IsSelfSimilarGroup, IsCyclotomic]",
                   [IsSelfSimGroup 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);

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

  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 [IsSelfSimGroup and IsGroupOfSelfSimFamily, IsPosInt]",
                    [IsSelfSimGroup and IsGroupOfSelfSimFamily, IsPosInt],
function(G, n)
  return DiagonalPower(UnderlyingSelfSimFamily(G), n);
end);

###############################################################################
##
#M  MultAutomAlphabet(<G>, <n>)
##
InstallOtherMethod( MultAutomAlphabet,
                    "for [IsSelfSimGroup and IsGroupOfSelfSimFamily, IsPosInt]",
                    [IsSelfSimGroup and IsGroupOfSelfSimFamily, IsPosInt],
function(G, n)
  return MultAutomAlphabet(UnderlyingSelfSimFamily(G), n);
end);

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

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

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

  if fam!.rws <> fail then
    fgens1 := AG_ReducedForm(fam!.rws, fgens1);
    fgens2 := AG_ReducedForm(fam!.rws, fgens2);
  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 [IsSelfSimGroup, IsSelfSimGroup]",
              IsIdenticalObj, [IsSelfSimGroup, IsSelfSimGroup],
function(G, H)
  local h, fam, fgens1, fgens2;

  if HasIsGroupOfSelfSimFamily(G) and IsGroupOfSelfSimFamily(G) then
    Info(InfoAutomGrp, 3, "IsSubgroup(G, H): true");
    Info(InfoAutomGrp, 3, "  G is GroupOfSelfSimFamily");
    return true;
  fi;

  fgens1 := List(GeneratorsOfGroup(G), g -> Word(g));
  fgens2 := List(GeneratorsOfGroup(H), g -> Word(g));
  fam := UnderlyingSelfSimFamily(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 [IsSelfSim, IsSelfSimGroup]",
              [IsSelfSim, IsSelfSimGroup],
function(g, G)
  local fam, fgens, w;

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

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

  fam := UnderlyingSelfSimFamily(G);

  if fam!.rws <> fail then
    fgens := AG_ReducedForm(fam!.rws, fgens);
    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);

###############################################################################
##
#O  Random(<G>)
##
##  Returns a random element of a group (semigroup) <G>. The operation is based
##  on the generator of random elements in free groups and semigroups.
##
##  \beginexample
##  gap> Basilica := AutomatonGroup( "u=(v,1)(1,2), v=(u,1)" );
##  < u, v >
##  gap> Random( Basilica );
##  v*u^-3
##  \endexample
##
InstallMethodWithRandomSource(Random, "for a random source and [IsSelfSimGroup]",
              [IsRandomSource, IsSelfSimGroup],
function(rs, G)
  local F, gens, pi;

  if IsTrivial(G) then
    return One(G);
  elif IsSelfSimilarGroup(G) then
    return SelfSim(Random(rs, UnderlyingFreeGroup(G)), UnderlyingSelfSimFamily(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 [IsSelfSimGroup]",
              [IsSelfSimGroup],
function(G)
  local f;
  if HasIsGroupOfSelfSimFamily(G) and IsGroupOfSelfSimFamily(G) then
    return UnderlyingFreeGroup(G);
  fi;
  f := Subgroup(UnderlyingFreeGroup(G), UnderlyingFreeGenerators(G));
  if f = UnderlyingFreeGroup(G) then
    SetIsGroupOfSelfSimFamily(G, true);
  fi;
  return f;
end);

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

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

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

###############################################################################
##
#M  RecurList(<G>)
##
InstallMethod(RecurList, "for [IsSelfSimGroup]",
              [IsSelfSimGroup],
function(G)
  if IsSelfSimilarGroup(G) then
    return RecurList(GroupOfSelfSimFamily(UnderlyingSelfSimFamily(G)));
  else
    Error("Group <G> is not necessarily self-similar");
  fi;
end);

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

###############################################################################
##
#M  IsFiniteState(<G>)
##
InstallMethod(IsFiniteState, "for [IsSelfSimGroup]",
              [IsSelfSimGroup],
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, MealyAutomatonLocalFinite;

# if we do not know much, we compare just words in free group
  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;

# if we do know that the groups is finite, we compare actual elements of the group
  MealyAutomatonLocalFinite := function(g)
    local cur_state;
    if g in states then return Position(states, g); fi;
    Add(states, g);
    cur_state := Length(states);
    aut_list[cur_state] := List([1..g!.deg], x -> MealyAutomatonLocalFinite(Section(g, x)));
    Add(aut_list[cur_state], g!.perm);
    return cur_state;
  end;

  if IsTrivial(G) then return true; fi;

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

  if HasIsFinite(G) and IsFinite(G) then
    for g in gens do
      Add(images, MealyAutomatonLocalFinite(g));
    od;
    states := List(states, Word);
  else
    for g in gens do
      Add(images, MealyAutomatonLocal(g));
    od;
  fi;

  H := AutomatonGroup(aut_list);

  if IsTrivial(H) then
    SetIsTrivial( G, true);
    return true;
  fi;

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

  SetIsomorphicAutomGroup(G, GroupWithGenerators(UnderlyingAutomFamily(H)!.automgens{images}));
  SetUnderlyingAutomatonGroup(G, H);

# preimages of generators of G in UnderlyingFreeGroup(G)
  gens_in_freegrp := List(GeneratorsOfGroup(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;
    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([1..Length(GeneratorsOfGroup(H))], x->states[Position(UnderlyingAutomFamily(H)!.oldstates, x)]);
  preimages_in_freegrp := List(GeneratorsOfGroup(H), x -> preimage_in_freegrp(x));

  if IsSelfSimilarGroup(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 := GroupHomomorphismByFunction(G, GroupWithGenerators(UnderlyingAutomFamily(H)!.automgens{images}), hom_function, inv_hom_function);

    SetMonomorphismToAutomatonGroup(G, hom);
  else
    F := FreeGroup(Length(GeneratorsOfGroup(G)));

#        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 := GroupHomomorphismByFunction(G, GroupWithGenerators(UnderlyingAutomFamily(H)!.automgens{images}), hom_function, inv_hom_function);

    SetMonomorphismToAutomatonGroup(G, hom);
  fi;

  return true;
end);

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

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

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

###############################################################################
##
#M  IsContracting( <G> )
##
InstallMethod(IsContracting, "for [IsSelfSimilarGroup]",
              [IsSelfSimilarGroup],
function(G)
  local res;
  if not IsFiniteState(G) then
#   every contracting self-similar group is finite-state
    return false;
  fi;

  res := IsContracting(GroupOfAutomFamily(UnderlyingAutomFamily(UnderlyingAutomatonGroup(G))));

  UnderlyingSelfSimFamily(G)!.use_contraction := true;
  UnderlyingAutomFamily(UnderlyingAutomatonGroup(G))!.use_contraction := true;

  return res;
end);

###############################################################################
##
#M  GroupNucleus( <G> )
##
InstallMethod(GroupNucleus, "for [IsSelfSimilarGroup]",
              [IsSelfSimilarGroup],
function(G)
  local H;
  if not IsFiniteState(G) then
#   every contracting self-similar group is finite-state
    Error("Group <G> is not finite-state");
  fi;

  if not IsContracting(G) then
    Error("Group <G> is not contracting");
  fi;

  H := GroupOfAutomFamily( UnderlyingAutomFamily( UnderlyingAutomatonGroup(G)));

  return List( GroupNucleus(H), x -> PreImagesRepresentative( MonomorphismToAutomatonGroup(G), x));
end);

###############################################################################
##
#M  UseContraction( <G> )
##
InstallMethod(UseContraction, "for [IsSelfSimGroup]", true,
              [IsSelfSimGroup],
function(G)
  if not IsSelfSimilarGroup(G) then
    Print("Error in UseContraction(<G>): The method is implemented only for IsSelfSimilarGroup\n");
    return fail;
  fi;

  if not HasIsContracting(G) then
    Print("Error in UseContraction(<G>): It is not known whether the group <G> is contracting\n");
    return fail;
  elif not IsContracting(G) then
    Print("Error in UseContraction(<G>): The group <G> is not contracting");
    return fail;
  fi;
  #  IsContracting returns either true or false or an error (it can not return fail)

  UnderlyingSelfSimFamily(G)!.use_contraction := true;
  UnderlyingAutomFamily( UnderlyingAutomatonGroup(G))!.use_contraction := true;

  return true;
end);

###############################################################################
##
#M  DoNotUseContraction( <G> )
##
InstallMethod(DoNotUseContraction, "for [IsSelfSimGroup]", true,
              [IsSelfSimGroup],
function(G)
  UnderlyingAutomFamily(G)!.use_contraction := false;

  if HasUnderlyingAutomatonGroup(G) then
    UnderlyingAutomFamily( UnderlyingAutomatonGroup(G))!.use_contraction := false;
  fi;
  return true;
end);

###############################################################################
##
#M  FindNucleus( <G> )
##
InstallMethod(FindNucleus, "for [IsSelfSimilarGroup, IsCyclotomic]",
              [IsSelfSimilarGroup, IsCyclotomic],
function(G, max_nucl)
  local H, nuclH, nuclG;
  if not IsFiniteState(G) then
#   every contracting self-similar group is finite-state
    Error("Group <G> is not finite-state");
  fi;

  if HasIsContracting(G) and not IsContracting(G) then
    Error("Group <G> is not contracting");
  fi;

  H := GroupOfAutomFamily( UnderlyingAutomFamily( UnderlyingAutomatonGroup( G )));

  if HasIsContracting(H) and not IsContracting(H) then
    Error("Group <G> is not contracting");
  fi;

  nuclH := FindNucleus(H, max_nucl);

  if nuclH=fail then return fail; fi;

  nuclG := [];
  Add(nuclG, List( GeneratingSetWithNucleus(H), x -> PreImagesRepresentative( MonomorphismToAutomatonGroup( G ), x )));
  Add(nuclG, List( GroupNucleus(H), x -> PreImagesRepresentative( MonomorphismToAutomatonGroup( G ), x )));
  Add(nuclG, GeneratingSetWithNucleusAutom(H));

  SetGroupNucleus(G, nuclG[1]);
  SetGeneratingSetWithNucleus(G, nuclG[2]);
  SetGeneratingSetWithNucleusAutom(G, nuclG[3]);
  SetContractingLevel(G, ContractingLevel(H));

  return nuclG;
end);

InstallMethod(FindNucleus, "for [IsSelfSimilarGroup]", true,
              [IsSelfSimilarGroup],
function(G)
  return FindNucleus(G, infinity);
end);

InstallMethod(GeneratingSetWithNucleus, "for [IsSelfSimilarGroup]", true,
              [IsSelfSimilarGroup],
function(G)
  if IsContracting(G) then return GeneratingSetWithNucleus(G); fi;
end);

InstallMethod(GeneratingSetWithNucleusAutom, "for [IsSelfSimilarGroup]", true,
              [IsSelfSimilarGroup],
function(G)
  if IsContracting(G) then return GeneratingSetWithNucleusAutom(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:
lpbp := (1,2,3);
rpbp := (1,3,4);

gpbp := Group([lpbp, rpbp]);

lfbp := (1,2,3)(4,5,6);
rfbp := (1,6,7)(4,3,8);

gfbp := Group([lfbp,rfbp]);

bs := Blocks(gfbp, [1..8]);

hom := ActionHomomorphism(gfbp, bs, OnSets);

im := Image(hom);

iso := IsomorphismGroups(im, gfbp);

k := Kernel(im);

ext: gap

max_stars_repo_name: dvberkel/ultimate-skweb

max_issues_repo_name: dvberkel/ultimate-skweb

max_forks_repo_name: dvberkel/ultimate-skweb


text: Flat([[1], 2, [[3,4], 5], [[[]]], [[[6]]], 7, 8, []]);

max_stars_repo_name: mullikine/RosettaCodeData

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text:
#
# usage:
# Read("~/Workspace/Chevalley.gap/init.gi"); Read(Filename(home_dir,"load.gi")); Read(Filename(f4_dir_char2,"_init_char2.gi")); Read(Filename(f4_dir_char2,"sA2.gi"));
#
# Read(Filename(f4_dir_char2,"sA2.gi"));
#

label:="\\tilde A_2";
orb_nr:=Position(labels,label);
#orb:=AllClasses(orbs)[14];
orb:=AllClasses(orbs)[orb_nr];
info:=infos[orb_nr];

#
# ------------------------------------------------------------------------
#

Print("Consider the class ",Label(orb)," in characteristic ",Characteristic(ring(sys)),":\n");
Print("\tBorel representative \n\t",coefficients(BorelRep(orb)),"\n");
Print("\tconnected C_U(u) \n\t",coefficients(info[1]),"\n");
Print("\tconnected double C_U(u) in Levi \n\t",coefficients(info[2]),"\n");
Print("\tconnected double C_U(u) \n\t",coefficients(info[3]),"\n");

Z0:=info[2];
# move Z0 to match the representative
Z0:=FromPositiveBorel(orb,Z0);
# adjust ordering of roots
Z0:=Unipotent(chevalleyAdj(Z0),coefficients(Z0),[1..Length(positiveRoots(sys))]);
Print("Z0 in position of Representative:\n\tZ0=",coefficients(Z0),"\n");

Print("Component group is C2.\n");
Print("On Z0:\n");
n:=[8,10];
aZ0:=ApplyRootsReflections(Z0,n);
u:=Unipotent(chevalleyAdj(Z0),[[3,-1]]);
aZ0:=Conj(aZ0,u);

Print("\ta^Z0=",coefficients(aZ0),"\n");

#
# ------------------------------------------------------------------------
#

tmp:=handle14char2();
Print("\nComponent group check:\n");
Print("\ta^u= \n",coefficients(tmp[1]),"\n");
Print("\tconnected center= \n",coefficients(tmp[2]),"\n");


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

Class(SVECodegen, DefaultCodegen, rec(
    GathH :=  ( self, o, y, x, opts ) >> assign(nth(y, 0), sve_gath(nth(x, o.base).toPtr(TVectSVE), o.n, opts.pg, o.stride)),
#    ScatH :=  ( self, o, y, x, opts ) >> sve_scat(nth(y, o.base).toPtr(TVectSVE), o.n, opts.pg, o.stride, nth(x, 0)),
    ScatH :=  ( self, o, y, x, opts ) >> sve_scat(nth(y, o.base).toPtr(TReal), o.n, opts.pg, o.stride, nth(x, 0)),

    VGath_L :=  ( self, o, y, x, opts ) >> let(_v := var.fresh_t("svex2_", TVectSVEx2),
        decl([_v], chain(
            assign(_v, sve_ld2(nth(x, o.base).toPtr(TReal), o.n, opts.pg)),
            assign(nth(y, 0), velem(_v, 0)),
            assign(nth(y, 1), velem(_v, 1))
        ))
    ),

    VScat_L :=  ( self, o, y, x, opts ) >> let(_v := var.fresh_t("svex2_", TVectSVEx2),
        decl([_v], chain(
            assign(velem(_v, 0), nth(x, 0)),
            assign(velem(_v, 1), nth(x, 1)),
            sve_st2(nth(y, o.base).toPtr(TReal), o.n, opts.pg, _v)
        ))
    )

));

max_stars_repo_name: spiral-software/spiral-package-ffte

max_issues_repo_name: spiral-software/spiral-package-ffte

max_forks_repo_name: spiral-software/spiral-package-ffte


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

@1 := @(1);
@2 := @(2);
zero   := @.cond(e -> e=0 or  (IsValue(e) and (e.v=0  or e.v=0.0)));
one    := @.cond(e -> e=1 or  (IsValue(e) and (e.v=1  or e.v=1.0)));
negone := @.cond(e -> e=-1 or (IsValue(e) and (e.v=-1 or e.v=-1.0)));

IS := objid -> When(IsList(objid),
    o -> ObjId(o) in objid,
    o -> Same(ObjId(o), objid)
);

IS_A := (o, objid) -> Same(ObjId(o), objid);

CHILDREN_ARE := (objid) -> (o -> ForAll(o.children(), IS(objid)));

NOT_ALL_CHILDREN_ARE := (objid) -> (o -> not ForAll(o.children(), IS(objid)));

NOT_ALL_ARGS_ARE := (objid) -> (o -> not ForAll(o.args, IS(objid)));

countRealAdds := c -> Sum(Collect(c, @(1, [add,sub], x->IsRealT(x.t))), op -> Length(op.args)-1);
countRealMuls := c -> Sum(Collect(c, @(1, mul, x->IsRealT(x.t))), op -> Length(op.args)-1);

countVectMuls := c -> Sum(Collect(c, @(1, mul, x->IsVecT(x.t))), op -> Length(op.args)-1);
countVectAdds := c -> Sum(Collect(c, @(1, [add,sub], x->IsVecT(x.t))), op -> Length(op.args)-1);


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

X := var("X", TPtr(TReal));
Y := var("Y", TPtr(TReal));
X1 := var("X1", TPtr(TReal));
Y1 := var("Y1", TPtr(TReal));
XY1 := var("XY1", TPtr(TReal));

# defined in cgstrat.gi
Declare(ApplyCodegenStrat);

RemoveAssignAcc := code -> SubstTopDownRulesNR(code, rec(
    assign_acc_toassign := Rule(assign_acc, e -> assign(Copy(e.loc), e.loc+e.exp)),
    assign_donothing := Rule(assign, e -> e),
));

#F CodeSums(<sums>, <opts>)
#F    Creates code for Sigma-SPL <sums>.
#F    Sigma-SPL rewrite rules are not applied.
#F    Compilation proceeds according to <opts>.compileStrategy
#F

CodeSums := function(sums, opts)
    local X, Y, code;

        if IsList(sums.dims()[2]) then
            X := List([1..Length(sums.dims()[2])],
                x -> var(Concatenation("X",String(x)), TPtr(TReal)));
        else
            X := Cond(IsBound(opts.X), opts.X, var("X", TPtr(TReal)));
        fi;

        if IsList(sums.dims()[1]) then
            Y := List([1..Length(sums.dims()[1])],
                x -> var(Concatenation("Y",String(x)), TPtr(TReal)));
        else
            Y := Cond(
                opts.inplace, X,
                IsBound(opts.Y), opts.Y,
                var("Y", TPtr(TReal)));
        fi;

        code := opts.codegen(Formula(sums), Y, X, opts);
        code.ruletree := Cond(IsBound(sums.ruletree), sums.ruletree, rec());
        return code;
end;

#F See CodeSums
#F
CodeSumsOpts := CodeSums;

#F CodeRuleTree(<rt>, <opts>)
#F
#F <opts> flags used:
#F   opts.formulaStrategies.sigmaSpl    Sigma-SPL rewriting strategy
#F   opts.formulaStrategies.rc          RC(.) rewriting strategy
#F   opts.generateComplexCode == bool   if set to false, then RC rewriting strategy is applied
#F
CodeRuleTree := function(rt, opts)
    local sums, code;
    sums := SumsRuleTree(rt, opts);
    code := CodeSums(sums, opts);
    return code;
end;

#F See CodeRuleTree
#F
CodeRuleTreeOpts := CodeRuleTree;

#F RealSums(<sums>)
#F    Convert complex Sigma-SPL formula to a real formula.
#F
RealSums := sums -> StandardSumsRules(RC(sums));

#F CodeSPL(<unroll>, <spl>, opts)
#F
CodeSPL := function(spl, opts)
    local sums;
    sums := SumsSPL(spl, opts);
    sums := ApplyStrategy(sums, opts.formulaStrategies.sigmaSpl, UntilDone, opts);
    sums := ApplyStrategy(sums, opts.formulaStrategies.preRC, UntilDone, opts);
    if not spl.isReal() and not opts.generateComplexCode then
        sums := ApplyStrategy(RC(sums), opts.formulaStrategies.rc, UntilDone, opts); fi;
    sums := ApplyStrategy(sums, opts.formulaStrategies.postProcess, UntilDone, opts);
    return CodeSums(sums, opts);
end;

#F PrintCode(<funcname>, <code>, <opts>)
#F    Prints the code using opts.unparser
#F
PrintCode := (funcname, code, opts) -> opts.unparser.gen(funcname, code, opts);

_ExportCodeRuleTree := function(ruletree, file, funcname, opts)
    local code;
    Constraint(IsRuleTree(ruletree));
    if opts.verbosity > 0 then Print("Generating code...\n"); fi;
    code := CodeRuleTree(ruletree, opts);
    PrintTo(file, opts.unparser.gen(funcname, code, opts));
end;

#F ExportCodeRuleTree(<ruletree>, <funcname>, <opts>)
#F    Generates and exports C code to a file for <ruletree>.
#F    The C transform function will have the name <funcname>.
#F    Code will be saved in <funcname>.c
#F
ExportCodeRuleTree := (ruletree, funcname, opts) ->
    _ExportCodeRuleTree(ruletree, Concat(funcname, ".c"), funcname, opts);

#F PrintCodeRuleTree(<ruletree>, <opts>)
#F    Generates C code and prints it out.
#F
PrintCodeRuleTree := (ruletree, opts) ->
    _ExportCodeRuleTree(ruletree, "*stdout*", "sub", opts);

#F ImplementRuleTree(<ruletree>, <file>, <opts>)
#F
ImplementRuleTree := (ruletree, file, opts) ->
    _ExportCodeRuleTree(ruletree, file,
    When(IsBound(opts.subName), opts.subName, ruletree.node.name), opts);

#F VerifyMatrixRuleTree(<ruletree>, <opts>)
#F
VerifyMatrixRuleTree := function(ruletree, opts)
    local code, mat;
    Constraint(IsRuleTree(ruletree));
    if opts.verbosity > 0 then Print("Generating code...\n"); fi;
    code := CodeRuleTree(ruletree, opts);
    if opts.verbosity > 0 then Print("Computing reference matrix...\n"); fi;
    mat := When(ruletree.node.isReal() or opts.dataType = "complex" or opts.generateComplexCode,
            MatSPL(ruletree.node),
        RCMatCyc(MatSPL(ruletree.node)));
    if opts.verbosity > 0 then Print("Running code and computing the norm...\n"); fi;
    return InfinityNormMat(CMatrix(code, opts) - mat);
end;

#F VerifyMatrixCode(<code>, <definition-matrix>, <opts>)
#F
VerifyMatrixCode := (code, def_matrix, opts) -> Checked(IsCommand(code),
    InfinityNormMat(CMatrix(code, opts) - def_matrix)
);

FailedTrees := [];
InaccurateTrees := [];

CMeasureRuleTree := function(rt, opts)
    local res, c, mfunc, tol;

    c := CodeRuleTree(rt, opts);

    mfunc := When(IsBound(opts.profile) and IsBound(opts.profile.meas), opts.profile.meas, CMeasure);

    if opts.faultTolerant then
		res := Try(mfunc(c, opts));
    else
		res := [true, mfunc(c, opts)];
    fi;

    if res[1]=false then
        Add(FailedTrees, rt);
        return 1e20;
    else
        if IsBound(opts.verifyDP) and opts.verifyDP then
            tol := VerifyMatrixRuleTree(rt, opts);
            if tol > opts.verifyTolerance then
                Add(InaccurateTrees, rec(rt:=rt, opts := opts, tol :=tol));
            fi;
        fi;
        return res[2];
    fi;
end;


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

########################################################################
#   vector constructs
########################################################################

IsSIMD_ISA := s -> IsRec(s) and IsBound(s.isSIMD_ISA) and s.isSIMD_ISA;
IsISA      := x -> IsRec(x) and IsBound(x.isISA) and x.isISA;

Class(AVecReg, AGenericTag, rec(
    isReg := true,
    isRegCx := false,
    isVec := true,
    updateParams := meth(self)
        Checked(IsSIMD_ISA(self.params[1]));
        Checked(Length(self.params)=1);
        self.v := self.params[1].v;
        self.isa := self.params[1];
    end,
    container := (self, spl) >> paradigms.vector.sigmaspl.VContainer(spl, self.isa)
));

Class(AVecRegCx, AVecReg, rec(
    updateParams := meth(self)
        Checked(IsSIMD_ISA(self.params[1]));
        Checked(Length(self.params)=1);
        self.v := self.params[1].v/2;
        self.isa := self.params[1];
    end,
    container := (self, spl) >> paradigms.vector.sigmaspl.VContainer(spl, self.isa.cplx()),
    isRegCx := true
));

# AMultiVec - list of ISAs, must be list of AVecReg tags in the future
#

Class(AMultiVec, AGenericTag, rec(
    isVec := true,
    updateParams := meth(self)
        Checked(ForAll(self.params, IsSIMD_ISA));
        Checked(Length(self.params)>=1);
    end,
));

Class(AISA, AGenericTag, rec(
    updateParams := meth(self)
        Checked(IsISA(self.params[1]));
        Checked(Length(self.params)=1);
        self.isa := self.params[1];
    end,
    # it's not a vectorized code, maybe defferent kind of containers?
    # containers do not go along well with OL though
    container := (self, spl) >> paradigms.vector.sigmaspl.VContainer(spl, self.isa)
));


text:
#
# cls 40: A_2
#
handle40char2:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[40];
    info:=infos[40];

    rep:=Representative(o);
    #pr[17]=[0,0,1,1,1,0,0]
    #pr[15]=[0,1,1,1,0,0,0]
    #pr[32]=[1,1,1,2,1,0,0]
    tmp1:=ApplyRootsReflections(rep,[17,15,32]);
    uu:=Unipotent(chevalleyAdj(rep),[[3,1]]);
    tmp1:=Conj(tmp1,uu);

    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[17,15,32]);
    uu:=Unipotent(chevalleyAdj(cent),[[3,1]]);
    tmp2:=Conj(tmp2,uu);

    return [tmp1,tmp2];
end;

#
# cls 38: A_2A_1
#
handle38char2:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[38];
    info:=infos[38];

    rep:=Representative(o);
    #pr[44]=[1,1,1,2,1,1,1]
    #pr[43]=[1,1,1,2,2,1,0]
    #pr[62]=[1,2,3,4,3,2,1]
    tmp1:=ApplyRootsReflections(rep,[44,43,62]);
    uu:=Unipotent(chevalleyAdj(rep),[[1,-1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);

    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[44,43,62]);
    uu:=Unipotent(chevalleyAdj(cent),[[1,-1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);

    return [tmp1,tmp2];
end;

#
# cls 30: D_4(a_1)
#
handle30char2:=function(orbs,infos)
    local tmp1,tmp2,tmp3,tmp4,o,info,rep,cent,uu;
    o:=Classes(orbs)[30];
    info:=infos[30];

    rep:=Representative(o);
    #pr[4] =[0,0,0,1,0,0,0]
    #pr[42]=[1,1,2,2,1,1,0]
    #pr[43]=[1,1,1,2,2,1,0]
    tmp1:=ApplyRootsReflections(rep,[4,42,43]);
    uu:=Unipotent(chevalleyAdj(rep),[[3,-1],[5,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);
    tmp1:=ConjugateByTori(tmp1,[4,1,2,6],[-1,-1,-1,-1]);

    #pr[23]=[0,1,0,1,1,1,0]
    #pr[24]=[0,0,1,1,1,1,0]
    tmp3:=ApplyRootsReflections(rep,[23,24]);
    uu:=Unipotent(chevalleyAdj(rep),[[3,1],[4,1],[5,-1],[10,1],[11,-1]],Ordering(rep));
    tmp3:=Conj(tmp3,uu);
    tmp3:=ConjugateByTori(tmp3,[2,4],[-1,-1,-1,-1]);

    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[4,42,43]);
    uu:=Unipotent(chevalleyAdj(cent),[[3,-1],[5,1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);
    tmp2:=ConjugateByTori(tmp2,[4,1,2,6],[-1,-1,-1,-1]);

    tmp4:=ApplyRootsReflections(cent,[23,24]);
    uu:=Unipotent(chevalleyAdj(cent),[[3,1],[4,1],[5,-1],[10,1],[11,-1]],Ordering(cent));
    tmp4:=Conj(tmp4,uu);
    tmp4:=ConjugateByTori(tmp4,[2,4],[-1,-1,-1,-1]);

    return [tmp1,tmp3,tmp2,tmp4];
end;

#
# cls 27: D_4(a_1)A_1 [2,5,9,10,7]
#
handle27char2:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[27];
    #o:=orbs;
    info:=infos[27];
    #info:=infos;

    rep:=Representative(o);

    #pr[20]=[1,1,1,1,0,0,0]
    #pr[21]=[1,0,1,1,1,0,0]
    tmp1:=ApplyRootsReflections(rep,[20,21]);
    uu:=Unipotent(chevalleyAdj(rep),[[3,-1],[4,-1],[5,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);
    tmp1:=ConjugateByTori(tmp1,[1,4],[-1,-1,-1,-1]);

    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[20,21]);
    uu:=Unipotent(chevalleyAdj(cent),[[3,-1],[4,-1],[5,1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);
    tmp2:=ConjugateByTori(tmp2,[1,4],[-1,-1,-1,-1]);

    return [tmp1,tmp2];
end;

#
# cls 26: A_3A_2^{(2)} from [1,3,4,6,7,37] we get [31,20,40,30,27,37] has no component group
#
handle26char2:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[27];
    #o:=orbs;
    info:=infos[27];
    #info:=infos;

    rep:=Representative(o);

    #pr[20]=[1,1,1,1,0,0,0]
    #pr[21]=[1,0,1,1,1,0,0]
    tmp1:=ApplyRootsReflections(rep,[20,21]);
    uu:=Unipotent(chevalleyAdj(rep),[[3,-1],[4,-1],[5,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);
    tmp1:=ConjugateByTori(tmp1,[1,4],[-1,-1,-1,-1]);

    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[20,21]);
    uu:=Unipotent(chevalleyAdj(cent),[[3,-1],[4,-1],[5,1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);
    tmp2:=ConjugateByTori(tmp2,[1,4],[-1,-1,-1,-1]);

    return [tmp1,tmp2];
end;

##
## cls 26: A_3A_2 [3,4,2,6,7] (not finished - check it)
##
#handle26char3:=function(orbs,infos)
#    local tmp1,tmp2,o,info,rep,cent,uu;
#    o:=Classes(orbs)[26];
#    info:=infos[26];
#
#    rep:=Representative(o);
#
#    #pr[23]=[0,1,0,1,1,1,0]
#    #pr[24]=[0,0,1,1,1,1,0]
#    #pr[45]=[0,1,1,2,2,1,1]
#    tmp1:=ApplyRootsReflections(rep,[23,24,45]);
#    uu:=Unipotent(chevalleyAdj(rep),[[2,-1],[3,1],[6,-1]],Ordering(rep));
#    tmp1:=Conj(tmp1,uu);
#    tmp1:=ConjugateByTori(tmp1,[4],[-1]);
#
#    cent:=FromPositiveBorel(o,info[2]);
#    tmp2:=ApplyRootsReflections(cent,[23,24,45]);
#    uu:=Unipotent(chevalleyAdj(cent),[[6,-1]],Ordering(cent));
#    tmp2:=Conj(tmp2,uu);
#    tmp2:=ConjugateByTori(tmp2,[4],[-1]);
#
#    return [tmp1,tmp2];
#end;

#
# cls 25: A_4 [1,2,3,4]
#
handle25char3:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[25];
    info:=infos[25];

    rep:=Representative(o);

    #pr[29]=[0,1,1,1,1,1,0]
    #pr[27]=[1,0,1,1,1,1,0]
    #pr[51]=[1,1,2,2,2,1,1]
    #pr[61]=[1,2,2,4,3,2,1]
    tmp1:=ApplyRootsReflections(rep,[29,27,51,61]);
    uu:=Unipotent(chevalleyAdj(rep),[[2,-1],[4,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);
    tmp1:=ConjugateByTori(tmp1,[1,2,3,4],[-1,-1,-1,-1]);

    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[29,27,51,61]);
    uu:=Unipotent(chevalleyAdj(cent),[[2,-1],[4,1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);
    tmp2:=ConjugateByTori(tmp2,[1,2,3,4],[-1,-1,-1,-1]);

    return [tmp1,tmp2];
end;

#
# cls 20: D_5(a_1) (not finished)? cred ca e gata
#
handle20char2:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[20];
    info:=infos[20];

    rep:=Representative(o);

    #pr[59]=[1,1,2,3,3,2,1]
    #pr[58]=[1,2,2,3,2,2,1]
    tmp1:=ApplyRootsReflections(rep,[59,58]);
    uu:=Unipotent(chevalleyAdj(rep),[[4,-1],[11,-1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);
    tmp1:=ConjugateByTori(tmp1,[1,2,4,5],[-1,-1,-1,-1]);

    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[59,58]);
    uu:=Unipotent(chevalleyAdj(cent),[[4,-1],[11,-1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);
    tmp2:=ConjugateByTori(tmp2,[1,2,4,5],[-1,-1,-1,-1]);

    return [tmp1,tmp2];
end;

#
# cls 21: A_4A_1
#
handle21char2:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[21];
    info:=infos[21];

    rep:=Representative(o);

    #pr[36]=[0,1,1,1,1,1,1]
    #pr[34]=[1,0,1,1,1,1,1]
    #pr[46]=[1,1,2,2,2,1,0]
    #pr[61]=[1,2,2,4,3,2,1]
    tmp1:=ApplyRootsReflections(rep,[36,34,46,61]);
    uu:=Unipotent(chevalleyAdj(rep),[[2,-1],[4,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);
    tmp1:=ConjugateByTori(tmp1,[1,2,3,4],[-1,-1,-1,-1]);

    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[36,34,46,61]);
    uu:=Unipotent(chevalleyAdj(cent),[[2,-1],[4,1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);
    tmp2:=ConjugateByTori(tmp2,[1,2,3,4],[-1,-1,-1,-1]);

    return [tmp1,tmp2];
end;

#
# cls 6: E_6(a_1) (not finished)
#
handle6char3:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[6];
    info:=infos[6];

    #pr[47]=[1,1,2,2,1,1,1]
    #pr[48]=[1,1,1,2,2,1,1]
    #pr[49]=[0,1,1,2,2,2,1]

    rep:=Representative(o);
    rep:=Unipotent(alg,coefficients(rep),[1..Length(pr)]);#ordered rep
    tmp1:=ApplyRootsReflections(rep,[47,48,49]);
    tmp1:=ConjNegUa(tmp1,4,1);
    uu:=Unipotent(chevalleyAdj(rep),[[1,1],[8,1],[6,-1]],Ordering(rep));# this works for orderd rep
    tmp1:=Conj(tmp1,uu);
    tmp1:=ConjugateByTori(tmp1,[3,7],[-1,-1]);
    tmp1:=Unipotent(chevalleyAdj(rep),coefficients(tmp1),Ordering(Representative(o)));

    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[47,48,49]);
    tmp2:=ConjNegUa(tmp2,4,1);
    #tmp2:=ConjugateByTorus(cent,4,-1);# h_1 --> h_4

    return [tmp1,tmp2];
end;

#
# cls 14: E_6(a_3) for ([1,8,9,11,14,19] in E6) [1,9,10,12,16,22] in E7
#
handle14char2:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[14];
    info:=infos[14];

    rep:=Representative(o);
    tmp1:=ConjNegUa(rep,3,1);
    tmp1:=ConjNegUa(tmp1,2,1);
    uu:=Unipotent(chevalleyAdj(rep),[[5,1],[6,1],[8,1],[9,1],[15,1],[16,1],[22,1],[23,1],[29,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);

    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ConjNegUa(cent,3,1);
    tmp2:=ConjNegUa(tmp2,2,1);
    uu:=Unipotent(chevalleyAdj(cent),[[5,1],[6,1],[8,1],[9,1],[15,1],[16,1],[22,1],[23,1],[29,1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);

    return [tmp1,tmp2];
end;

#
# cls 4: E_7(a_3) [1,7,9,10,12,16,22] in E7
#
handle4char2:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[4];
    info:=infos[4];

    rep:=Representative(o);
    tmp1:=ConjNegUa(rep,3,1);
    tmp1:=ConjNegUa(tmp1,2,1);
    uu:=Unipotent(chevalleyAdj(rep),[[5,1],[6,1],[8,1],[9,1],[15,1],[16,1],[22,1],[23,1],[29,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);

    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ConjNegUa(cent,3,1);
    tmp2:=ConjNegUa(tmp2,2,1);
    uu:=Unipotent(chevalleyAdj(cent),[[5,1],[6,1],[8,1],[9,1],[15,1],[16,1],[22,1],[23,1],[29,1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);

    return [tmp1,tmp2];
end;


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

# Top-Level Functions for SPIRAL
# ==============================
# BWS, MP, from 02/15/00

#F Creating Transforms
#F -------------------
#F

#F Transforms
#F   the list of symbols for valid transforms, e.g., "DFT".
#F
Transforms := NonTerminalListSPL;

#F Transform ( <symbol>, <parameters> )
#F   alias for SPLNonTerminal, returns the transform defined by
#F   <symbol> and <parameters>, e.g. Transform( "DFT", 8 ).
#F
#Transform := SPLNonTerminal;

#F Info ( <symbol> )
#F   prints information on how to create the transform <symbol>
#F
Info := Doc;

#F toSPL( <obj> )
#F   Converts ruletree or nonterminal (using random ruletree) into an SPL,
#F   If <obj> is an SPL it is returned as is.
toSPL := x ->
    Cond(IsNonTerminal(x), SPLRuleTree(RandomRuleTree(x)),
	 IsSPL(x), x,
	 IsRuleTree(x), SPLRuleTree(x),
	 Error("<x> must be an SPL or a RuleTree"));


text:
#
# MathInTheMiddle: Math-in-the-Middle functionality for GAP
#
# Type export to JSON for consumation by MMT import
#
# TODO: We need at some point advertise our abilities via
#       OpenMath directly

# Make GAP Type graph as a record
InstallGlobalFunction(MitM_TypesInfo, function()
    local res, lres, i, j, f, ff, a, meths, mpos, objs, m, mres, n, t, notcovered, v, op, loc,
          flags, flagslist;

    objs := NewDictionary(IsObject, true);
    res := [ rec( name := "IsObject", type := "Category", implied := [] ) ];

    Print("exporting operations...\c");
    for i in [1..Length(OPERATIONS)] do
        op := OPERATIONS[i];

        lres := rec();
        AddDictionary(objs, op, lres);

        lres.type := TypeOfOperation(op);
        lres.name := NameFunction(op);

        # Locations in which this operation is declared
        # There can be multiple locations with different filters
        # We convert the GAP representation as a list of pairs
        # into a list of records for the export
        lres.locations := List( GET_DECLARATION_LOCATIONS(op)
                              , x -> rec( file := x[1], line := x[2] ) );

        # This is way too complicated
        # if IsAttribute(op) then
        # elif IsProperty(op) then
        # else
        lres.filters := [];

        # This is a list if argument filters for every declaration
        # for op
        flagslist := GET_OPER_FLAGS(op);
        if flagslist = fail then
        # TODO, and filter? Synonym?
            Print("Op: ", op, " filters need some looking after\n");
        else
            # TODO: clean up
            for flags in flagslist do
                ff := Concatenation(List(flags, x -> List(TRUES_FLAGS(x), z -> FILTERS[z])));
                ff := List(ff, NameFunction);
            od;
        fi;
        lres.methods := rec( 0args := [], 1args := [], 2args := [],
                             3args := [], 4args := [], 5args := [],
                             6args := [] );

        for a in [1..6] do
            meths := MethodsOperation(op, a);

            for j in meths do
                mres := rec( filters := List(j.argFilt, x -> List(TRUES_FLAGS(x)
                                                                 ,y -> NameFunction(FILTERS[y])))
                           , rank := j.rank
                           , comment := j.info );
                # Methods are not bound to global variables directly (usually...)
                # but some methods might be global functions
                AddDictionary(objs, j.func, mres);
                Add(lres.methods.(Concatenation(String(a),"args")), mres);
            od;
        od;

        Add(res, lres);
    od;
    Print("   done\n");

    Print("exporting global functions...\c");
    for f in [1..Length(GLOBAL_FUNCTION_NAMES)] do
        lres := rec();
        lres.type := "Function";
        lres.name := GLOBAL_FUNCTION_NAMES[f];
        lres.location := rec();
        AddDictionary(objs, ValueGlobal(GLOBAL_FUNCTION_NAMES[f]), lres);
        Add(res, lres);
    od;
    Print("   done\n");

    Print("collecting global variable references...\c");
    notcovered := [];
    for n in NamesGVars() do
        if IsBoundGlobal(n) then
            v := ValueGlobal(n);
            t := LookupDictionary(objs, v);
            if t <> fail then
                if IsBound(t.aka) then
                    Add(t.aka, n);
                else
                    t.aka := [n];
                fi;
            else
                v := ValueGlobal(n);
                if IsFilter(v) then
                    lres := rec();
                    lres.type := "GAP_AndFilter";
                    ff := FLAGS_FILTER(v);
                    if ff <> false then
                        ff := TRUES_FLAGS(FLAGS_FILTER(v));
                        ff := List(ff, function(f)
                                      if IsBound(FILTERS[f]) then
                                          return NameFunction(FILTERS[f]);
                                      else
                                          return "<<unknown>>";
                                      fi;
                                  end);
                        lres.conjunction_of := ff;
                        lres.name := (NameFunction(v));
                        lres.aka := [n];
                        AddDictionary(objs, v, lres);
                        Add(res, lres);
                    else
                #        Print("strange: ", n, " ", v);
                    fi;
                elif IsFunction(v) and not '_' in n and ForAny("abcdefghijklmnopqrstuvwxyz", c -> c in n) then
                    # Defined using BindGlobal, not InstallGlobalFunction
                    # If it contains '_' or is all-caps it is probably internal
                    lres := rec();
                    lres.type := "Function";
                    lres.name := n;
                    lres.location := rec();
                    AddDictionary(objs, v, lres);
                    Add(res, lres);
                else
                    Add(notcovered, n);
                fi;
            fi;
        fi;
    od;
    Print("   done\n");
    return [res, notcovered];
end);

# Write the graph of type info to JSon file
InstallGlobalFunction(MitM_TypesToJson,
function(file)
    local fd, n, typeinfo;

    fd := IO_File(file, "w");
    if fd = fail then
        Error("Opening file ", file, " failed");
    fi;
    typeinfo := MitM_TypesInfo();
    n := IO_Write(fd, GapToJsonString(typeinfo[1]));
    IO_Close(fd);

    return n;
end);

max_stars_repo_name: markuspf/MathsInTheMiddle

max_issues_repo_name: markuspf/MathsInTheMiddle

max_forks_repo_name: markuspf/MathsInTheMiddle


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

#
# Functions
#
IsASPFreqBasis := x -> IsRec(x) and IsBound(x.isASPFreqBasis) and x.isASPFreqBasis;
IsASPTimeBasis := x -> IsRec(x) and IsBound(x.isASPTimeBasis) and x.isASPTimeBasis;

Class(BasisOps, PrintOps, rec(
    \= := (b1, b2) -> Cond(
        not IsASPFreqBasis(b1) or not IsASPFreqBasis(b2), false,
        ObjId(b1)=ObjId(b2) and b1.rChildren() = b2.rChildren())
));

#
# Base classes
#
Class(ASPTimeBasis, rec(
    isAtomic := true, # NOTE: these probably should not be atomic, maybe use TimeComp(Freq_E) for composed bases
    isASPTimeBasis := true,
    doNormalize := false,
    operations := PrintOps,
    print := self >> Print(self.__name__, When(self.doNormalize, ".norm()")),
    norm := self >> CopyFields(self, rec(doNormalize:=true)),
    rChildren := self >> [],
    from_rChildren := (self, rch) >> self,
    unparse := "ConstClass",
));

Class(ASPComposedTimeBasis, ASPTimeBasis, rec(
    composedFrom := self >> Error("must be defined in the subclass"),
    toX := (self, alg) >> Checked(ObjId(alg) = XN_skew, IsEvenInt(alg.n),
        let(C := When(self.doNormalize, self.composedFrom().norm(),
                                        self.composedFrom()),
            Tensor(C.to1X(XN_skew(2, alg.a, alg.rot)), I(alg.n/2))))
));

Class(ASPFreqBasis, rec(
    doNormalize := false,
    isASPFreqBasis := true,
    dim := self >> Error("must be defined in the subclass"),
    from_rChildren := (self, rch) >> ApplyFunc(ObjId(self), rch),
    unparse := "RewritableObjectExp",
));

Class(ASPFreqBasis1d,  ASPFreqBasis, rec(
    __call__ := (self, scale1d) >> WithBases(self, rec(
            scale1d := scale1d, operations := BasisOps
    )),
    rChildren := self >> [self.scale1d],
    rSetChild := rSetChildFields("scale1d"),
    print := self >> Print(self.name, "(", self.scale1d, ")", When(self.doNormalize, ".norm()")),
    dim := self >> 1,
    hashAs := self >> ObjId(self)(1)
));

Class(ASPFreqBasis1d2d,  ASPFreqBasis, rec(
    __call__ := (self, scale1d, scale2d) >> WithBases(self, rec(
            scale1d := scale1d, scale2d := scale2d, operations := BasisOps
    )),
    rChildren := self >> [self.scale1d, self.scale2d],
    rSetChild := rSetChildFields("scale1d", "scale2d"),
    print := self >> Print(self.name, "(", self.scale1d, ", ", self.scale2d, ")",
        When(self.doNormalize, ".norm()")),

    norm := self >> CopyFields(self, rec(doNormalize:=true)),

    dim := self >> 2,
    hashAs := self >> ObjId(self)(1, 1),

    to1X_2d0 := (self, rot) >> I(2),
    to1X_2d := (self, a, rot) >> Error("Must be implemented in a subclass"),
    to1X := (self, alg) >> let(t := ObjId(alg), s1 := self.scale1d, s2 := self.scale2d, Cond(
        alg.n=1, 1/s1 * I(1),

        alg.n=2 and (t in [XN_min_1, XN_min_1U]),
            1/s2 * self.to1X_2d0(alg.rot),

        alg.n=2 and (t in [XN_skew, XN_plus_1]), let(
            a := alg.a, #EvalScalar(alg.a),
            aa := Cond(self.doNormalize, cond(gt(a, 0.5), 1-a, a), a), #Cond(self.doNormalize and a > 1/2, 1-a, a),
            1/s2 * self.to1X_2d(aa, alg.rot)),

        Error("Algebra <alg> must have dimension of 1 or 2"))),

    from1X_2d0 := (self, rot) >> I(2),
    from1X_2d := (self, a, rot) >> Mat(MatSPL(self.to1X_2d(a, rot))^-1),
    from1X := (self, alg) >> let(t := ObjId(alg), s1 := self.scale1d, s2 := self.scale2d, Cond(
        alg.n=1, s1 * I(1),

        alg.n=2 and (t in [XN_min_1, XN_min_1U]),
            s2 * self.from1X_2d0(alg.rot),

        alg.n=2 and (t in [XN_skew, XN_plus_1]), let(
            a := alg.a, #EvalScalar(alg.a),
            aa := Cond(self.doNormalize, cond(gt(a, 0.5), 1-a, a), a),
            s2 * self.from1X_2d(aa, alg.rot)),

        Error("Algebra <alg> must have dimension of 1 or 2"))),
));

#
# TIME BASIS Classes
#
Declare(Freq_1, Freq_1H, Freq_T, Freq_S, Freq_TH, Freq_THU,  Freq_E, Freq_EH,  Freq_H, Freq_HH);

Declare(Time_TX, Time_SX, Time_THX, Time_THUX,  Time_EX, Time_EHX,  Time_HX, Time_HHX);

# ASP time basis: standard "time" basis (1, x, ..., x^{n-1})
Class(Time_TX,  ASPTimeBasis, rec(
    toX := (self, alg) >> I(alg.n)
));

# ASP time basis: symmetric "time" basis (-x^{n/2}, ..., 1, x, ..., x^{n/2-1})
Class(Time_SX,  ASPTimeBasis, rec(
    toX := (self, alg) >> let(t := ObjId(alg), hf := Int(alg.n/2),
        Cond(t in [XN_min_1, XN_min_1U],
                Z(alg.n, hf),
             t = XN_plus_1 or (t = XN_skew and alg.a = 1/4),
                 DirectSum(I(Int((alg.n+1)/2)), -I(hf)) * Z(alg.n, hf),
             t = XN_skew,
                 Tensor(Freq_S(1,1).to1X(XN_skew(2, alg.a, alg.rot)), I(alg.n/2)),
             Error("Unrecognized algebra <alg>")))
));

Class(Time_THX,  ASPComposedTimeBasis, rec(composedFrom := self >> Freq_TH(1,1)));
Class(Time_THUX, ASPComposedTimeBasis, rec(composedFrom := self >> Freq_THU(1,1)));
Class(Time_EX,   ASPComposedTimeBasis, rec(composedFrom := self >> Freq_E(1,1))); # e_a(x^n) (*) t_n
Class(Time_EHX,  ASPComposedTimeBasis, rec(composedFrom := self >> Freq_EH(1,1)));
Class(Time_HX,   ASPComposedTimeBasis, rec(composedFrom := self >> Freq_H(1,1))); # h_a(x^n) (*) t_n
Class(Time_HHX,  ASPComposedTimeBasis, rec(composedFrom := self >> Freq_HH(1,1)));

#
# FREQUENCY BASIS Classes
#

# ASP frequency basis: (1)
Class(Freq_1,  ASPFreqBasis1d, rec(
    timeBasis := self >> Time_TX
));

Class(Freq_1H,  ASPFreqBasis1d, rec(
    timeBasis := self >> Error("not implemented"),
    to1X := (self, alg) >> Error("not implemented"),
));

# ASP frequency basis: (1, x), also implies 2d spectral components, ie real DFT
Class(Freq_T, ASPFreqBasis1d2d, rec(
    timeBasis := self >> Time_TX,
    to1X := (self, alg) >>  let(t := ObjId(alg), s1 := self.scale1d, s2 := self.scale2d, Cond(
        alg.n = 2, 1/s2 * I(2),
        alg.n = 1, 1/s1 * I(1),
        Error("Algebra <alg> must have dimension of 1 or 2"))),
    from1X := (self, alg) >>  let(t := ObjId(alg), s1 := self.scale1d, s2 := self.scale2d, Cond(
        alg.n = 2, s2 * I(2),
        alg.n = 1, s1 * I(1),
        Error("Algebra <alg> must have dimension of 1 or 2"))),
));

# ASP frequency basis: (x^-1, 1), also implies 2d spectral components, ie real DFT
Class(Freq_S, ASPFreqBasis1d2d, rec(
    timeBasis := self >> Time_SX,
    to1X_2d0 := (self, rot) >> J(2),
    to1X_2d := (self, a, rot) >> Mat(_stripval([[ 2*cospi(2*a*rot), 1],
                                                [             -1  , 0]])),
));

# ASP frequency basis: (x^1/2, x^3/2) for type 2--4 transforms, also implies 2d spectral components, ie real DFT
Class(Freq_TH, ASPFreqBasis1d2d, rec(
    timeBasis := self >> Time_THX,
    to1X_2d := (self, a, rot) >> let(c := cospi(2*a*rot), q := 1/(2*cospi(a*rot)),
        Mat(_stripval(q * [[1+2*c, 1],
                           [-1,    1]])))
));

Class(Freq_THU, ASPFreqBasis1d2d, rec(
    timeBasis := self >> Time_THUX,
    to1X_2d := (self, a, rot) >> let(c := cospi(2*a*rot), q := 1/(2*cospi(a * rot)),
        Mat(_stripval(q^2 * [[1+2*c, 1],
                             [-1,    1]])))
));

# ASP frequency basis: (1, (x-cos(a)) / sin(a)), also implies 2d spectral components, ie real DFT
Class(Freq_E,  ASPFreqBasis1d2d, rec(
    timeBasis := self >> Time_EX,
    to1X_2d := (self, a, rot) >> let(c := cospi(2*rot*a), s := sinpi(2*rot*a),
        Mat(_stripval([[1, -c/s],
                       [0,  1/s]]))),
    from1X_2d := (self, a, rot) >> let(c := cospi(2*rot*a), s := sinpi(2*rot*a),
        Mat(_stripval([[1,  c],
                       [0,  s]])))
));

Class(Freq_EH,  ASPFreqBasis1d2d, rec(
    timeBasis := self >> Time_EX,
    to1X_2d := (self, a, rot) >> let(c := r->cospi(r*2*rot*a), s := r->sinpi(r*2*rot*a),
        Mat(_stripval(1/s(1) * [[ s(3/2), -c(3/2)],
                                [-s(1/2),  c(1/2)]])))
));

# ASP frequency basis: (1 + (x-cos(a)) / sin(a), 1 - (x-cos(a))/sin(a)), ie Hartley transform
Class(Freq_H,  ASPFreqBasis1d2d, rec(
    timeBasis := self >> Time_HX,
    to1X_2d := (self, a, rot) >> let(c := cospi(2*rot*a), s := sinpi(2*rot*a),
        Mat(_stripval(1/2 * [[1-c/s, 1+c/s],
                             [  1/s,  -1/s]])))
));

Class(Freq_HH,  ASPFreqBasis1d2d, rec(
    timeBasis := self >> Time_EX,
    to1X_2d := (self, a, rot) >> let(s := sinpi(2*rot*a),
        cas := r -> cospi(r*2*rot*a)+sinpi(r*2*rot*a), cms := r -> cospi(r*2*rot*a)-sinpi(r*2*rot*a),
        Mat(_stripval(1/(2*s) * [[-cms(3/2), -cas(3/2)],
                                 [ cms(1/2),  cas(1/2)]])))
));


text:
n := 0;
repeat
    n := n + 1;
    Print(n, "\n");
until RemInt(n, 6) = 0;


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

NewRulesFor(TFCall, rec(
    TFCall_tag := rec(
        forTransposition := false,
        children := nt -> [[ nt.params[1].withTags(nt.getTags()) ]],
        apply := (nt, c, cnt) -> c[1]
    )
));

NewRulesFor(TDeviceCall, rec(
    TFCall_tag := rec(
        forTransposition := false,
        children := nt -> [[ nt.params[1].withTags(nt.getTags()) ]],
        apply := (nt, c, cnt) -> c[1]
    )
));

NewRulesFor(TSparseMat, rec(
    TSparseMat_base := rec(
        forTransposition := false,
        children := nt -> [[ ]],
        applicable := (self, nt) >> When(IsBound(self.max_rows), nt.params[1][1] <= self.max_rows, true),
        apply := function(nt, c, cnt)
            local P, mat, i, j, rdata;
            P := nt.params;
            mat := MatSPL(ApplyFunc(O, P[1]));

            for i in [0..P[1][1]-1] do
                rdata := Drop(Filtered(P[2], ii->ii[1] = i)[1], 1);
                for j in rdata do
                    mat[i+1][j[1]+1] := j[2];
                od;
            od;

            return Mat(mat);
        end
    ),
    TSparseMat_VStack := rec(
        forTransposition := false,
        applicable := nt -> nt.params[1][1] > 1,
        children := nt -> [ List([0..nt.params[1][1]-1], i->TSparseMat([1, nt.params[1][2]], [let(row := Filtered(nt.params[2], j->j[1] = i)[1], [0] :: Drop(row, 1))])) ],
        apply := (nt, c, cnt) -> BB(VStack(c))
    ),
));

NewRulesFor(TIterHStack, rec(
    TIterHStack_tag := rec(
        forTransposition := false,
        children := nt -> [[ nt.params[1].withTags(nt.getTags()), InfoNt(nt.params[2]) ]],
        apply := (nt, c, cnt) -> IterHStack(cnt[2].params[1], cnt[2].params[1].range, c[1])
    )
));

NewRulesFor(TIterVStack, rec(
    TIterVStack_tag := rec(
        forTransposition := false,
        children := nt -> [[ nt.params[1].withTags(nt.getTags()), InfoNt(nt.params[2]) ]],
        apply := (nt, c, cnt) -> IterVStack(cnt[2].params[1], cnt[2].params[1].range, c[1])
    )
));

NewRulesFor(TNoDiagPullinLeft, rec(
    TNoDiagPullinLeft_tag := rec(
        forTransposition := false,
        children := nt -> [[ nt.params[1].withTags(nt.getTags()) ]],
        apply := (nt, c, cnt) -> NoDiagPullinLeft(c[1])
    )
));

NewRulesFor(TNoDiagPullinRight, rec(
    TNoDiagPullinRight_tag := rec(
        forTransposition := false,
        children := nt -> [[ nt.params[1].withTags(nt.getTags()) ]],
        apply := (nt, c, cnt) -> NoDiagPullinRight(c[1])
    )
));

NewRulesFor(TNoPullLeft, rec(
    TNoPullLeft_tag := rec(
        forTransposition := false,
        children := nt -> [[ nt.params[1].withTags(nt.getTags()) ]],
        apply := (nt, c, cnt) -> NoPullLeft(c[1])
    )
));

NewRulesFor(TNoPullRight, rec(
    TNoPullRight_tag := rec(
        forTransposition := false,
        children := nt -> [[ nt.params[1].withTags(nt.getTags()) ]],
        apply := (nt, c, cnt) -> NoPullRight(c[1])
    )
));

NewRulesFor(TColMajor, rec(
    TColMajor_MDDFT_2D := rec(
        forTransposition := false,
        switch := false,
        applicable := (self, nt) >> ObjId(nt.params[1]) = MDDFT and Length(nt.params[1].params[1]) = 2,

        children := nt -> [[ nt.params[1].withTags(nt.getTags()) ]],
        apply := (nt, c, cnt) -> let(n1 := nt.params[1].params[1][1], n2 := nt.params[2][1].params[1][2],
            prm := L(n1*n2, n2),
            iprm := prm.transpose(),
            iprm * c[1] * prm)
    ),
    TColMajor_MDDFT_3D := rec(
        forTransposition := false,
        switch := false,
        applicable := (self, nt) >> ObjId(nt.params[1]) = MDDFT and Length(nt.params[1].params[1]) = 3,

        children := nt -> [[ nt.params[1].withTags(nt.getTags()) ]],
        apply := (nt, c, cnt) -> let(n1 := nt.params[1].params[1][1], n2 := nt.params[1].params[1][2], n3 := nt.params[1].params[1][3],
            # prm := Tensor(I(n1), L(n2*n3, n3)) * L(n1*n2*n3, n2*n3),
            prm := Prm(fCompose(L(n1*n2*n3, n2*n3), fTensor(fId(n1), L(n2*n3, n3)))),
            # iprm := L(n1*n2*n3, n1) * Tensor(I(n1), L(n2*n3, n2)),
            iprm := prm.transpose(),
            iprm * c[1] * prm)
    ),
    TColMajor_MDDFT_tensorFlip := rec(
        forTransposition := false,
        applicable := (self, nt) >> ObjId(nt.params[1]) = MDDFT,

        children := nt -> [[ ApplyFunc(MDDFT, [Reversed(nt.params[1].params[1])]::Drop(nt.params[1].params, 1)).withTags(nt.getTags()) ]],
        apply := (nt, c, cnt) -> c[1]
    ),
));

max_stars_repo_name: franzfranchetti/spiral-package-fftx

max_issues_repo_name: franzfranchetti/spiral-package-fftx

max_forks_repo_name: franzfranchetti/spiral-package-fftx


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

ImportAll(spiral.formgen);
ImportAll(spiral.transforms);
ImportAll(spiral.compiler);

doScalarDft := function(arg)
	local szs, userOpts, isSingle, opts, direction, func, file, t, rt, ir, N, cycles, flops;

	if not (Length(arg) in [1..2]) then
		Error("usage: doScalarDft(sizes, [userOptions])");
	fi;

	szs      := When(IsList(arg[1]), arg[1], [ arg[1] ]);
	userOpts := When(IsBound(arg[2]), arg[2], rec());

	isSingle := IsBound(userOpts.precision) and userOpts.precision = "single";

	opts := CopyFields(SpiralDefaults, rec(
				precision  := When(isSingle, "single", "double"),
				TRealCtype := When(isSingle, "float", "double"),
			));

	opts.globalUnrolling := 512;

	direction := When(IsBound(userOpts.transInverse) and userOpts.transInverse, 1, -1);

	for N in szs do
		if IsBound(userOpts.functionNameRoot) then
			func := ReplaceAll(userOpts.functionNameRoot, "%sz1", StringInt(N));
		else
			func := When(direction=1,"i","")::"DFT"::When(isSingle, "_SC", "_DC")::StringInt(N);
		fi;

		file := func::".c";

		t    := DFT(N, direction);
		rt   := RandomRuleTree(t, opts);
		ir   := CodeRuleTree(rt, opts);

		PrintTo(file, PrintCode(func, ir, opts));

		cycles := CMeasure(ir, opts);
		flops:= t.normalizedArithCost();
		PrintLine(t, "  ", cycles, " [cyc]  ", _compute_gflops(flops, cycles), " [Gf/s]");
	od;

	return true;
end;

doScalar2DDft := function(arg)
	local szs, userOpts, isSingle, opts, direction, func, file, t, rt, ir, dims, cycles, flops;

	if not ((Length(arg) in [1..2]) and IsList(arg[1]) and Length(arg[1]) > 0) then
		Error("usage: doScalar2DDft(sizes, [userOptions])");
	fi;

	szs      := arg[1];
	if Length(szs) = 2 and IsInt(szs[1]) and IsInt(szs[2]) then
		szs := [ szs ];
	elif not ForAll(szs, dims -> Length(dims) = 2 and IsInt(dims[1]) and IsInt(dims[2])) then
		Error("sizes must be one or more pairs of integers");
	fi;

	userOpts := When(IsBound(arg[2]), arg[2], rec());

	isSingle := IsBound(userOpts.precision) and userOpts.precision = "single";

	opts := CopyFields(SpiralDefaults, rec(
				precision  := When(isSingle, "single", "double"),
				TRealCtype := When(isSingle, "float", "double"),
			));

	opts.globalUnrolling := 128;

	direction := When(IsBound(userOpts.transInverse) and userOpts.transInverse, 1, -1);

	for dims in szs do
		if IsBound(userOpts.functionNameRoot) then
			func := ReplaceAll(userOpts.functionNameRoot, "%sz1", StringInt(dims[1]));
			func := ReplaceAll(func, "%sz2", StringInt(dims[2]));
		else
			func := When(direction=1,"i","")::"DFT"::When(isSingle, "_SC", "_DC")::StringInt(dims[1])::"x"::StringInt(dims[2]);
		fi;

		file := func::".c";

		t    := MDDFT(dims, direction);
		rt   := RandomRuleTree(t, opts);
		ir   := CodeRuleTree(rt, opts);

		PrintTo(file, PrintCode(func, ir, opts));

		cycles := CMeasure(ir, opts);
		flops:= t.normalizedArithCost();
		PrintLine(t, "  ", cycles, " [cyc]  ", _compute_gflops(flops, cycles), " [Gf/s]");
	od;

	return true;
end;

Class(ScriptGenScalar, ScriptGenBase, rec(

	_arch := () -> "Scalar",

	_init := meth(self)
		self._setTransform(SGKEY_FFT);
		self._setType(SGKEY_SPCX);
		self._setSize(16);
		self._setFilename("");
	end,

	_validTransforms := meth(arg)
		return [SGKEY_FFT, SGKEY_IFFT, SGKEY_FFT_2D, SGKEY_IFFT_2D];
	end,

	_validSizes := meth(arg)
		local self, szs, xform, type;
		self  := arg[1];
		szs   := [];
		xform := self.getSettingsValue(SGKEY_TRANSFORM);
		type  := self.getSettingsValue(SGKEY_DATATYPE);

		if xform in [SGKEY_FFT, SGKEY_IFFT] then
			szs := Filtered([2..512], i->ForAll(Factors(i), j->j<=19));
			Append(szs, Filtered([513..1024], i->ForAll(Factors(i), j->j<=19) and IsInt(i/16)));
			Append(szs, List([11..16], i->2^i));
		elif xform in [SGKEY_FFT_2D, SGKEY_IFFT_2D] then
			szs := List(Filtered([2..360], i->ForAll(Factors(i), j->j<=19)), n -> [n, n]);
		elif xform = SGKEY_WHT then
			if type = SGKEY_DPCX then
				szs := List([2..10], i->2^i);
			else
				szs := List([4..10], i->2^i);
			fi;
		fi;

		return szs;
	end,

	_validTypes := meth(arg)
		return [SGKEY_SPCX, SGKEY_DPCX];
	end,

	getScriptChoices := (self) >> [SGSTR_RUNRANDOMALL],

	_genScript := meth(self, runType)
		local scrstr, xform, type, szs, tempopts, optrec, optstr, funcstr;

		xform	 := self.getSettingsValue(SGKEY_TRANSFORM);
		type	 := self.getSettingsValue(SGKEY_DATATYPE);
		szs      := self.getSettingsValue(SGKEY_SIZE);

		tempopts := "scriptOpts";

		optrec	 := rec(
						precision  := When(type=SGKEY_SPCX, "single", "double"),
						searchType := runType,
					);
		if IsBound(self._settings.(SGKEY_FUNCNAME)) then
			optrec.functionNameRoot := self._settings.(SGKEY_FUNCNAME);
		fi;
		if xform in [SGKEY_IFFT, SGKEY_IFFT_2D] then
			optrec.transInverse := true;
		fi;

		if xform in [SGKEY_FFT, SGKEY_IFFT] then
			funcstr := "doScalarDft("::String(szs)::", "::tempopts::")";
		elif xform in [SGKEY_FFT_2D, SGKEY_IFFT_2D] then
			funcstr := "doScalar2DDft("::String(szs)::", "::tempopts::")";
		else
			return "";
		fi;

		optstr := StringPrint(optrec);

		scrstr := "Import(scriptgen);;\n";
		Append(scrstr, tempopts::" := "::optstr::";;\n");
		Append(scrstr, funcstr::";;\n");

		return scrstr;
	end,
)); # Class ScriptGenScalar

SetScriptGenConstructor(ScriptGenScalar);


text:
method vfClipMan vfClipMan.mLoad <alias=vfClipMan_mLoad>( )
{
//	this->vForm.mCreateWin()
	ustr empus
	str  emps
	ustr ustmp
	uint comp
	comp as this
		comp.AutoLang=1
		comp.Border=$fbrdSizeable
		comp.Bottom=0
		comp.Caption=empus
		comp.Enabled=1
		comp.FormStyle=$fsChild
		comp.Height=521
		comp.HelpTopic=empus
		comp.Hint=empus
		comp.HorzAlign=$alhLeft
		comp.IconName=empus
		comp.Left=0
		comp.Name="fClipMan"
		comp.Right=0
		comp.StartPos=$spDesigned
		comp.Style=emps
		comp.TabOrder=0
		comp.Tag=0
		comp.Top=0
		comp.TopMost=0
		comp.VertAlign=$alvTop
		comp.Visible=1
		comp.Width=610
		comp.WindowState=$wsNormal
		comp.OnPosChanged.Set( this, fClipMan_PosChanged )
		comp.OnCreate.Set( this, fClipMan_Create )
		comp.OnCloseQuery.Set( this, fClipMan_Close )
		comp.OnLanguage.Set( this, fClipMan_Lang )
		comp as this.Tray0
		comp.Owner = this
			comp.AutoLang=1
			comp.Caption=ustmp.fromutf8("progname")
			comp.Image=ustmp.fromutf8("main\\mclip")
			comp.Name="Tray0"
			comp.RBtnPopupMenu=.pmTray
			comp.Tag=0
			comp.Visible=1
			comp.OnMouse.Set( this, fClipMan_TrayMouse )
		comp as this.pmTray
		comp.Owner = this
			comp.AutoLang=1
			comp.Name="pmTray"
			comp.Tag=0
			comp as this.miTrayExit
			comp.Owner = this.pmTray
				comp.AutoCheck=0
				comp.AutoLang=1
				comp.Caption=ustmp.fromutf8("exit")
				comp.Checked=0
				comp.Ellipsis=0
				comp.Enabled=1
				comp.Image=empus
				comp.Name="miTrayExit"
				comp.RadioCheck=0
				comp.Separator=0
				comp.ShortKey=empus
				comp.Tag=0
				comp.Visible=1
				comp.OnClick.Set( this, fClipMan_Exit )
			comp as this.mPref
			comp.Owner = this.pmTray
				comp.AutoCheck=0
				comp.AutoLang=1
				comp.Caption=ustmp.fromutf8("preferences")
				comp.Checked=0
				comp.Ellipsis=1
				comp.Enabled=1
				comp.Image=empus
				comp.Name="mPref"
				comp.RadioCheck=0
				comp.Separator=0
				comp.ShortKey=empus
				comp.Tag=0
				comp.Visible=1
				comp.OnClick.Set( this, fClipMan_Settings )
			comp as this.mAbout
			comp.Owner = this.pmTray
				comp.AutoCheck=0
				comp.AutoLang=1
				comp.Caption=ustmp.fromutf8("about")
				comp.Checked=0
				comp.Ellipsis=1
				comp.Enabled=1
				comp.Image=empus
				comp.Name="mAbout"
				comp.RadioCheck=0
				comp.Separator=0
				comp.ShortKey=empus
				comp.Tag=0
				comp.Visible=1
				comp.OnClick.Set( this, fClipMan_About )
		comp as this.Menu0
		comp.Owner = this
			comp.AutoLang=1
			comp.Name="Menu0"
			comp.Tag=0
			comp as this.MenuItem0
			comp.Owner = this.Menu0
				comp.AutoCheck=0
				comp.AutoLang=1
				comp.Caption=ustmp.fromutf8("tools")
				comp.Checked=0
				comp.Ellipsis=0
				comp.Enabled=1
				comp.Image=empus
				comp.Name="MenuItem0"
				comp.RadioCheck=0
				comp.Separator=0
				comp.ShortKey=empus
				comp.Tag=0
				comp.Visible=1
				comp as this.MenuItem1
				comp.Owner = this.MenuItem0
					comp.AutoCheck=0
					comp.AutoLang=1
					comp.Caption=ustmp.fromutf8("settings")
					comp.Checked=0
					comp.Ellipsis=0
					comp.Enabled=1
					comp.Image=empus
					comp.Name="MenuItem1"
					comp.RadioCheck=0
					comp.Separator=0
					comp.ShortKey=empus
					comp.Tag=0
					comp.Visible=1
					comp.OnClick.Set( this, fClipMan_Settings )
	comp as this
		comp.ClientHeight=476
		comp.ClientWidth=592
		comp.Menu=.Menu0

	return this
}

method vfClipMan vfClipMan.init( )
{
   this.pTypeId = vfClipMan
   return this
}
func init_vfClipMan <entry>()
{
   regcomp( vfClipMan, "vfClipMan", vForm, $vForm_last,
      %{ %{$mLoad,     vfClipMan_mLoad}},
      0->collection )

}

max_stars_repo_name: novostrim/macroclip

max_issues_repo_name: novostrim/macroclip

max_forks_repo_name: novostrim/macroclip