text:
#############################################################################
##
#W  small.gi                 GAP group library             Hans Ulrich Besche
##                                               Bettina Eick, Eamonn O'Brien
##
##  This file contains the basic installations for the library of small
##  groups and the group identification routines.
##

#############################################################################
##
#F  SMALL_AVAILABLE_FUNCS
##
##  On every level of the small groups library one function is written into
##  this list. It will detect those sizes, which are contained in this
##  library level.
SMALL_AVAILABLE_FUNCS := [ ];

#############################################################################
##
#F  SMALL_AVAILABLE( size )
##
##  returns fail if the library of groups of <size> is not installed.
##  Otherwise a record with some information about the construction of the
##  groups of <size> is returned.
InstallGlobalFunction( SMALL_AVAILABLE, function( size )
    local l, r;

    if not IsPosInt( size ) then
      Error( "<size> must be a positive integer");
    fi;

    for l in [ 1 .. Length( SMALL_AVAILABLE_FUNCS ) ] do
        if IsBound( SMALL_AVAILABLE_FUNCS[ l ] ) then
            r := SMALL_AVAILABLE_FUNCS[ l ]( size );
            if r <> fail then
                return r;
            fi;
        fi;
    od;
    return fail;
end );

InstallGlobalFunction( SmallGroupsAvailable,
function(order)
    return SMALL_AVAILABLE(order) <> fail;
end);

InstallGlobalFunction( NumberSmallGroupsAvailable,
function(order)
    if order = 1024 then
        return true;
    else
        return SmallGroupsAvailable(order);
    fi;
end);

#############################################################################
##
#F  SMALL_GROUP_FUNCS
##
##  This list will contain all functions to construct/read the groups from
##  the library.
SMALL_GROUP_FUNCS := [ ];

#############################################################################
##
#F  CODE_SMALL_GROUP_FUNCS
##
##  This list will contain those functions used to read the code of the
##  groups of some sizes from the data files.
CODE_SMALL_GROUP_FUNCS := [ ];

#############################################################################
##
#F  NUMBER_SMALL_GROUPS_FUNCS
##
NUMBER_SMALL_GROUPS_FUNCS := [ ];

#############################################################################
##
#F  SELECT_SMALL_GROUPS_FUNCS
##
SELECT_SMALL_GROUPS_FUNCS := [ ];

#############################################################################
##
#V  SMALL_GROUP_LIB
##
##  This list will contain all data for the group construction read from the
##  small group library.
SMALL_GROUP_LIB := [ ];

#############################################################################
##
#V  PROPERTIES_SMALL_GROUPS
##
##  This list will contain all data for the group selection read from the
##  small group library.
PROPERTIES_SMALL_GROUPS := [ ];

#############################################################################
##
#F  SmallGroup(<size>,<i>)
##
##  returns the <i>th  group of  order <size> in the catalogue. It will return
##  an PcGroup, if the group is soluble and a permutation group otherwise.
##  If the groups of this size are not installed, it will return an error.
InstallGlobalFunction( SmallGroup, function( arg )
    local inforec, g, size, i, nid;

    if Length( arg ) = 1 then
        if not IsList( arg[1] ) or Length( arg[1] ) <> 2 then
            Error( "usage: SmallGroup( order, number )" );
        fi;
        size := arg[ 1 ][ 1 ];
        i    := arg[ 1 ][ 2 ];
    elif Length( arg ) = 2 then
        size := arg[ 1 ];
        i    := arg[ 2 ];
    else
        Error( "usage: SmallGroup( order, number )" );
    fi;
    if not IsPosInt( size ) or not IsPosInt( i ) then
        Error( "usage: SmallGroup( order, number )" );
    fi;
    inforec := SMALL_AVAILABLE( size );
    if inforec = fail then
        Error( "the library of groups of size ", size, " is not available" );
    fi;
    nid := i;
    if not SMALL_GROUPS_OLD_ORDER then
        if size = 3^7 then
            nid := SMALLGP_PERM3(i);
        elif size = 5^7 then
            nid := SMALLGP_PERM5(i);
        elif size = 7^7 then
            nid := SMALLGP_PERM7(i);
        elif size = 11^7 then
            nid := SMALLGP_PERM11(i);
        fi;
    fi;
    g := SMALL_GROUP_FUNCS[ inforec.func ]( size, nid, inforec );
    SetIdGroup( g, [ size, i ] );
    IsPGroup( g );
    return g;
end );

#############################################################################
##
#F  NumberSmallGroups(<size>)
##
##  returns the  number of groups of the order <size>.
InstallGlobalFunction( NumberSmallGroups, function( size )
    local inforec;

    if not IsPosInt( size ) then
        Error( "usage: NumberSmallGroups( order )" );
    fi;
    if size = 1024 then
        return 49487365422;
    fi;

    inforec := SMALL_AVAILABLE( size );
    if inforec = fail then
        Error( "the library of groups of size ", size, " is not available" );
    fi;

    if IsBound( inforec.number ) then
        return inforec.number;
    fi;
    return NUMBER_SMALL_GROUPS_FUNCS[ inforec.func ]( size, inforec ).number;
end );

#############################################################################
##
#F  SelectSmallGroups( argl, all, id )
##
InstallGlobalFunction( SelectSmallGroups, function( argl, all, id )
    local sizes, size, i, funcs, vals, gs, inforec, result, hasSizes, pos,
          idList;

    sizes := [ ];
    hasSizes := false;
    idList := fail;

    for i in [ 1 .. Length( argl ) ] do
        if i = 1 and argl[ i ] = Size then
            ;
        elif ( not hasSizes ) and IsList( argl[ i ] )
                              and Length( sizes ) = 1 then
            idList := argl[ i ];
        elif ( not hasSizes ) and IsList( argl[ i ] ) then
            Append( sizes, argl[ i ] );
        elif ( not hasSizes ) and IsPosInt( argl[ i ] ) then
            Add( sizes, argl[ i ] );
        elif ( not hasSizes ) and sizes <> [] and IsFunction( argl[i] ) then
            hasSizes     := true;
            funcs        := [ argl[ i ] ];
            vals         := [ [ ] ];
            pos          := 1;
        elif not hasSizes then
            Error( "usage: AllSmallGroups / OneGroup(\n",
                   "             Size, [ sizes ],\n",
                   "             function1, [ values1 ],\n",
                   "             function2, [ values2 ], ... )" );
        elif vals[ pos ] <> [ ] and IsFunction( argl[ i ] ) then
            pos          := pos + 1;
            funcs[ pos ] := argl[ i ];
            vals[ pos ]  := [ ];
        elif IsFunction( argl[ i ] ) then
            vals[ pos ]  := [ true ];
            pos          := pos + 1;
            funcs[ pos ] := argl[ i ];
            vals[ pos ]  := [ ];
        elif IsList( argl[ i ] ) and vals[ pos ] = [ ] then
            vals[ pos ]  := argl[ i ];
        elif IsList( argl[ i ] ) and IsInt( argl[ i ][ 1 ] ) and
             IsList( vals[ pos ][ 1 ] ) and IsInt( vals[ pos ][1][ 1 ] ) then
            Add( vals[ pos ], argl[ i ] );
        elif IsList( argl[ i ] ) and IsInt( argl[ i ][ 1 ] ) and
             IsInt( vals[ pos ][ 1 ] ) then
            vals[ pos ]  := [ vals[ pos ], argl[ i ] ];
        else
            Add( vals[ pos ], argl[ i ] );
        fi;
    od;

    if sizes <> [ ] and ( not IsBound( vals ) ) then
        funcs := [ ];
        vals  := [ ];
    elif vals[ pos ] = [ ] then
        vals[ pos ] := [ true ];
    fi;

    result := [ ];
    for size in sizes do
        inforec := SMALL_AVAILABLE( size );
        if inforec = fail then
            Error( "AllSmallGroups / OneGroup: groups of order ", size,
                   " not available" );
        fi;
        gs := SELECT_SMALL_GROUPS_FUNCS[ inforec.func ]
                             ( size, funcs, vals, inforec, all, id, idList );
        if all then
            Append( result, gs );
        elif gs <> fail then
            return gs;
        fi;
    od;

    if all then
        return result;
    else
        return fail;
    fi;
end );

#############################################################################
##
#F ID_AVAILABLE_FUNCS
##
ID_AVAILABLE_FUNCS := [ ];

#############################################################################
##
#F ID_AVAILABLE
##
InstallGlobalFunction( ID_AVAILABLE, function( size )
    local l, r;

    if not IsInt( size ) then return fail; fi;

    for l in [ 1 .. Length( ID_AVAILABLE_FUNCS ) ] do
        if IsBound( ID_AVAILABLE_FUNCS[ l ] ) then
            r := ID_AVAILABLE_FUNCS[ l ]( size );
            if r <> fail then
                return r;
            fi;
        fi;
    od;
    return fail;
end );

InstallGlobalFunction( IdGroupsAvailable,
function(order)
    return ID_AVAILABLE(order) <> fail;
end);

#############################################################################
##
#F  ID_GROUP_FUNCS
##
ID_GROUP_FUNCS := [ ];

#############################################################################
##
#M  IdGroup( G )
##
InstallMethod( IdGroup,
               "generic method for groups",
               true,
               [ IsGroup ],
               0,
function( G )
    local inforec, size, id;

    size := Size( G );
    if size = 1 then return [ 1, 1 ]; fi;

    inforec := ID_AVAILABLE( size );
    if inforec = fail then
        Error( "the group identification for groups of size ", size,
               " is not available" );
    fi;

    if not ( IsPcGroup( G ) or IsPermGroup( G ) ) then
        if IsSolvableGroup( G ) then
            G := Image( IsomorphismPcGroup( G ) );
        else
            G := Image( IsomorphismPermGroup( G ) );
        fi;
    fi;

    if Size( G ) > 1000 and IsPermGroup( G )
          and LargestMovedPoint( G ) > 100 and IsSolvableGroup( G ) then
        G := Image( IsomorphismPcGroup( G ) );
    fi;

    if IsPcGroup( G ) and HasParent( G ) and Size( Parent( G ) ) > 10000
       and Size( Parent( G ) ) / Size( G ) > 10 then
        G := PcGroupCode( CodePcGroup( G ), Size( G ) );
    fi;

    id := ID_GROUP_FUNCS[ inforec.func ]( G, inforec );

    if not SMALL_GROUPS_OLD_ORDER then
        if size = 3^7 then
          id:=First([1..NrSmallGroups(3^7)],x->SMALLGP_PERM3(x)=id);
        elif size = 5^7 then
          # note that the permutation is not an involution!
          id:=First([1..NrSmallGroups(5^7)],x->SMALLGP_PERM5(x)=id);
        elif size = 7^7 then
          id:=First([1..NrSmallGroups(7^7)],x->SMALLGP_PERM7(x)=id);
        elif size = 11^7 then
          id:=First([1..NrSmallGroups(11^7)],x->SMALLGP_PERM11(x)=id);
        fi;
    fi;
    return [ size, id ];
end );

#############################################################################
##
#V  ID_GROUP_TREE
##
##  Variable containing information for group identification
ID_GROUP_TREE := rec( fp := [ 1 .. 50000 ], next := [ ] );

#############################################################################
##
#F  ReadSmallLib( str, i, size, list )
##
##  universal reading function for data files
ReadSmallLib := function( str, i, size, list )
    local l, str2, str3, j;

    l := "abcdefghijklmnopqrstuvwxyz";
    str2 := Concatenation( str, String( size ) );
    str3 :=  [ ];
    for j in list do
        if j > 702 then
            Add( str3, l[ QuoInt( j - 27, 676 ) ] );
            j := 27 + ( j - 27 ) mod 676;
        fi;
        if j > 26 then
            Add( str3, l[ QuoInt( j - 1, 26 ) ] );
        fi;
        Add( str3, l[ ( j - 1 ) mod 26 + 1 ] );
    od;

    if Length( str2 ) > 8 then
        str3 := Concatenation( str2{[ 9..Length( str2 ) ]} , str3 );
        str2 := str2{[ 1..8 ]};
    elif Length( str3 ) = 0 then
        str3 := "z";
    elif Length( str3 ) > 3 then
        str2 := Concatenation( str2, str3{[ 1 .. Length( str3 ) - 3 ]} );
        str3 := str3{[ Length( str3 ) - 2 .. Length( str3 ) ]};
    fi;

    if str in [ "sml", "col", "prop", "nor" ] then
        READ_SMALL_FUNCS[ i ]( Concatenation( str2, ".", str3 ) );
    else
        READ_IDLIB_FUNCS[ i ]( Concatenation( str2, ".", str3 ) );
    fi;
end;

#############################################################################
##
#V  GAP3_CATALOGUE_ID_GROUP
##
##  List with the gap3-ids. Will be loaded before use.
GAP3_CATALOGUE_ID_GROUP := fail;

#############################################################################
##
#M  Gap3CatalogueIdGroup(<G>)
##
InstallMethod( Gap3CatalogueIdGroup,
               "for permgroups or pcgroups",
               true,
               [ IsGroup ],
               0,
function( G )
    if Size( G ) > 100 then
        Error( "Gap3CatalogueIdGroup: the group catalogue of gap-3.0 was\n",
               "limited to size 100" );
    fi;

    if GAP3_CATALOGUE_ID_GROUP = fail then
        ReadPackage( "smallgrp", "gap/gap3cat.g" );
    fi;

    if not IsBound( GAP3_CATALOGUE_ID_GROUP[ Size( G ) ] ) then
        return IdGroup( G );
    fi;

    return [ Size( G ),
             GAP3_CATALOGUE_ID_GROUP[ Size( G ) ][ IdGroup( G )[ 2 ] ] ];
end );

#############################################################################
##
#F  Gap3CatalogueGroup(<size>,<i>)
##
InstallGlobalFunction( Gap3CatalogueGroup,
function( size, i )
    local p;

    if size > 100 then
        Error( "Gap3CatalogueIdGroup: the group catalogue of gap-3.0 was\n",
               "limited to size 100" );
    fi;

    if GAP3_CATALOGUE_ID_GROUP = fail then
        ReadPackage("smallgrp", "gap/gap3cat.g");
    fi;

    if not IsBound( GAP3_CATALOGUE_ID_GROUP[ size ] ) then
        return SmallGroup( size, i );
    fi;

    p := Position( GAP3_CATALOGUE_ID_GROUP[ size ], i );
    if p = fail then
        Error( "Gap3CatalogueGroup: there are just ",
               Length( GAP3_CATALOGUE_ID_GROUP[ size ] ),
               " groups of size ", size );
    fi;

    return SmallGroup( size, p );
end );

#############################################################################
##
#F  UnloadSmallGroupsData( )
##
##  will remove the data from the small groups library from memory.
InstallGlobalFunction( UnloadSmallGroupsData, function( )
    SMALL_GROUP_LIB := [ ];
    PROPERTIES_SMALL_GROUPS := [ ];
    GAP3_CATALOGUE_ID_GROUP := fail;
    ID_GROUP_TREE := rec( fp := [ 1 .. 50000 ], next := [ ] );
end );

#############################################################################
##
#M  FrattinifactorSize(<G>)
##
InstallMethod( FrattinifactorSize,
               "generic method for groups",
               true,
               [ IsGroup ],
               0,
function( G )
    return Size( G ) / Size( FrattiniSubgroup( G ) );
end );

#############################################################################
##
#M  FrattinifactorId(<G>)
##
InstallMethod( FrattinifactorId,
               "generic method for groups",
               true,
               [ IsGroup ],
               0,
function( G )
    local ff;

    ff := G / FrattiniSubgroup( G );
    if ID_AVAILABLE( Size( ff ) ) = fail then
        Error( "FrattinifactorId: IdGroup for groups of size ", Size( ff ),
               " not available" );
    fi;
    return IdGroup( ff );
end );

max_stars_repo_name: hulpke/smallgrp

max_issues_repo_name: hulpke/smallgrp

max_forks_repo_name: hulpke/smallgrp


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

Class(VDirectSum, DirectSum, rec(
    abbrevs := [(L, v) -> [L, v]],
    new := meth(self, L, v)
        if ForAll(L, IsIdentitySPL)
            then return I(Sum(L, Rows));
        else
            return SPL(WithBases( self,
            rec( _children := L,
            v := v,
            dimensions := [ Sum(L, t -> t.dimensions[1]),
                            Sum(L, t -> t.dimensions[2]) ] )));
        fi;
    end,
#    sums := meth(self)
#        nblocks, c, cols, rows, bkcols, bkrows, gcd;
#        nblocks  :=  self.numChildren();
#        c := self.children();
#        cols  :=  EvalScalar(Cols(self));
#        rows := EvalScalar(Rows(self));
#        bkcols := List(c, i->EvalScalar(Cols(i)));
#        bkrows := List(c, i->EvalScalar(Rows(i)));
#        gcd := self.v;
#        bkrows_by_gcd := List(bkrows, i->i/gcd);
#        psr_by_gcd := List(DropLast(_partSum(bkrows), 1), i ->idiv(i,gcd));
#    fi;
#    psc := DropLast(Concat([0], List(Drop(_partSum(bkcols), 1), i->i-overlap)), 1);
#    return SUM( Map([1..nblocks],
#        i -> Compose(
#                    Scat(fTensor(fAdd(rows/gcd, bkrows_by_gcd[i], Cond(i=1, 0, i=2, bkrows_by_gcd[i-1], psr_by_gcd[i])), fId(gcd))),
#                    self.child(i).sums(),
#                    Gath(fTensor(fAdd(cols/gcd, bkcols[i]/gcd, psc[i]/gcd), fId(gcd))))
#            ));
));

ext: gi


text:
for i in [1 .. 11] do
    if RemInt(i, 5) = 0 then
        Print(i, "\n");
        continue;
    fi;
    Print(i, ", ");
od;

# 1, 2, 3, 4, 5
# 6, 7, 8, 9, 10

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

########################################################################
#   (A x B) rules
NewRulesFor(TTensor, rec(
#   (A x B) -> L (B x I)(L(A x I))
    splitL_BxI__L_AxI := rec(
        info := "(A x B) -> L (B x I)(L(A x I))",
        forTransposition := false,
        applicable := nt -> let(P := nt.params,
            Rows(P[1]) = Cols(P[1])
            and Rows(P[2]) = Cols(P[2])
        ),
        children := nt -> let(P := nt.params,
            [[TCompose([TTensorI(P[2], P[1].dims()[1], AVec, AVec), TTensorI(P[1], P[2].dims()[2], APar, AVec)]).withTags(nt.getTags())]]
        ),
        apply := (t, C, nt) -> let(
            p:=nt[1].params[1][1],
            mn := Rows(nt[1]),
            n := p.params[2],
            DelayedPrm(L(mn, n)) * C[1]
        )

#D        isApplicable := P -> Rows(P[1]) = Cols(P[1]) and Rows(P[2]) = Cols(P[2]),
#D        allChildren := P -> [[TCompose([TTensorI(P[2], P[1].dims()[1], AVec, AVec), TTensorI(P[1], P[2].dims()[2], APar, AVec)], P[3])]],
#D        rule := (P, C, nt) -> let(p:=nt[1].params[1][1], mn := Rows(nt[1]), n := p.params[2], DelayedPrm(L(mn, n)) * C[1])
    ),
#   (A x B) -> ((A x I)L)(B x I) L
    AxI_L__BxI_splitL := rec(
        info := "(A x B) -> ((A x I)L)(B x I) L",
        forTransposition := false,

        applicable := nt -> let(P := nt.params,
            Rows(P[1]) = Cols(P[1])
            and Rows(P[2]) = Cols(P[2])
        ),

        children := nt -> let(P := nt.params,
            [[TCompose([TTensorI(P[1], P[2].dims()[1], AVec, APar), TTensorI(P[2], P[1].dims()[2], AVec, AVec)]).withTags(nt.getTags())]]
        ),

        apply := (t, C, nt) -> let(
            p := nt[1].params[1][1],
            mn := Rows(nt[1]),
            n := Rows(p.params[1]),
            C[1] * DelayedPrm(L(mn, n))
        )
#D        isApplicable := P -> Rows(P[1]) = Cols(P[1]) and Rows(P[2]) = Cols(P[2]),
#D        allChildren := P -> [[TCompose([TTensorI(P[1], P[2].dims()[1], AVec, APar), TTensorI(P[2], P[1].dims()[2], AVec, AVec)], P[3])]],
#D        rule := (P, C, nt) -> let(p:=nt[1].params[1][1], mn := Rows(nt[1]), n := Rows(p.params[1]), C[1] * DelayedPrm(L(mn, n)))
    ),
));


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

########################################################################
#   (A + B) rules
NewRulesFor(TDirectSum, rec(
#   (A + B) terminate
    A_dirsum_B_delayed := rec(
        forTransposition := false,
        children := t -> let(tags := t.getTags(),
            [[ t.params[1].withTags(tags), t.params[2].withTags(tags) ]]
        ),
#D        children := (self, t) >> let (tags:=GetTags(t),
#D            [[ AddTag(t.params[1], tags), AddTag(t.params[2], tags) ]]),
        apply := (t, C, Nonterms) -> DelayedDirectSum(C)
    )
));


text:
# First position is built-in
haystack := Eratosthenes(10000);;
needle := 8999;;
Position(haystack, needle);
# 1117

LastPosition := function(L, x)
  local old, new;
  old := 0;
  new := 0;
  while new <> fail do
    new := Position(L, x, old);
    if new <> fail then
      old := new;
    fi;
  od;
  return old;
end;

a := Shuffle(List([1 .. 100], x -> x mod 10));
# [ 0, 2, 4, 5, 3, 1, 0, 4, 8, 8, 2, 7, 6, 3, 3, 6, 4, 4, 3, 0, 7, 1, 8, 7, 2, 4, 7, 9, 4, 9, 4, 5, 9, 9, 6, 7, 8, 2, 3,
#   5, 1, 5, 4, 2, 0, 9, 6, 1, 1, 2, 2, 0, 5, 7, 6, 8, 8, 3, 1, 9, 5, 1, 9, 6, 8, 9, 2, 0, 6, 2, 1, 6, 1, 1, 2, 5, 3, 3,
#   0, 3, 5, 7, 5, 4, 6, 8, 0, 9, 8, 3, 7, 8, 0, 4, 9, 7, 0, 6, 5, 7 ]
Position(a, 0);
# 1
LastPosition(a, 0);
# 97


text:
#############################################################################
##
#M  GraphicSubgroupLattice(<G>) . . . . displays subgroup lattice graphically
#M  GraphicSubgroupLattice(<G>,<def>)  . . . . . . . . . . same with defaults
##
##  Displays a graphic poset which shows (parts of) the subgroup lattice of
##  the group <group>. Normally only the whole group and the trivial group are
##  shown (behaviour of "InteractiveLattice" in xgap3). Returns a
##  IsGraphicSubgroupLattice object. Calls DecideSubgroupLatticeType. See
##  there for details.
##
InstallMethod(GraphicSubgroupLattice,
    "for a group, and a record",
    true,
    [IsGroup, IsRecord],
    0,
function(G, def)
  local latticeType, graph, poset, cls, vmath, v2, v1, menu, m, knownArgs, cb;

  latticeType := DecideSubgroupLatticeType(G);
  # we do some heuristics to avoid the trivial group:
  # if we know all levels, we probably can calc. Size, if we shall generate
  # a vertex for the trivial subgroup, we should also know Size!
  if latticeType.knowsLevels or latticeType.trivial then
    # no trivial case:
    if Size(G) = 1 then
      return Error( "<G> must be non-trivial" );
    fi;
  fi;

  # create a new canvas
  poset := Canvas("GraphicSubgroupLattice");
  # and make it a GraphicSubgroupLattice
  SetFilterObj(poset, IsGraphicSubgroupLattice);

  poset!.latticeType := latticeType;

  if HasName(G) then
    SetTitle(poset, Concatenation(Title(poset)," of ", Name(G)));
  elif HasIdGroup(G) then
    SetTitle(poset, Concatenation(Title(poset)," of ", String(IdGroup(G))));
  fi;

  if IsBound(def.width) then SetWidth(poset, def.width); fi;
  if IsBound(def.height) then SetHeight(poset, def.height); fi;
  if IsBound(def.title) then SetTitle(poset, def.title); fi;

  # now the other filters, depending on type:
  if latticeType.knowsLevels then SetFilterObj(poset, KnowsAllLevels); fi;
  if latticeType.hasse then SetFilterObj(poset, HasseProperty); fi;
  if latticeType.canCompare then SetFilterObj(poset, CanCompareSubgroups); fi;

  if KnowsAllLevels(poset) then
    # we know how many levels so we can show the user
    Add(poset, FrancyMessage(Concatenation("There are ",
      String(Length(DivisorsInt(Size(G)))), " levels in this Group.")));
  else
    # we don't know how many levels so we try to compute them
    if latticeType.trivial then
      if CanComputeSize(G) and Size(G) <> infinity then
        Add(poset, FrancyMessage(Concatenation("There are ", String(Size(G)), " levels in this Group.")));
      else
        Add(poset, FrancyMessage("We don't know how many levels there are on this group."));
      fi;
    else
      # we really don't know
      Add(poset, FrancyMessage("We don't know how many levels there are on this group."));
    fi;
  fi;

  # create an undirected graph to represent the poset
  # by adding layers later, the undirected graph turns
  # into an Hasse diagram
  graph := Graph(GraphType.DIRECTED);
  Add(poset, graph);

  cls:=ConjugacyClassesSubgroups(G);

  # create one or two initial vertices (G itself and trivial subgroup):
  v2 := Shape(ShapeType.DIAMOND, "G");
  SetLayer(v2, -4);
  for m in latticeType.contextMenus do
    knownArgs := [poset, Representative(cls[Length(cls)])];
    if m.group = true then
      Add(knownArgs, G);
    fi;
    cb := Callback(m.func, knownArgs);
    Add(v2, Menu(m.name, cb));
  od;
  Add(graph, v2);

  if latticeType.trivial then
    v1 := Shape(ShapeType.DIAMOND, "1");
    SetLayer(v1, 0);
    for m in latticeType.contextMenus do
      knownArgs := [poset, Representative(cls[1])];
      if m.group = true then
        Add(knownArgs, G);
      fi;
      cb := Callback(m.func, knownArgs);
      Add(v1, Menu(m.name, cb));
    od;
    Add(graph, v1);

    # connect the two vertices
    Add(graph, Link(v1, v2));
  fi;

  # <G> is selected at first
  v2!.selected := true;

  # create menus:
  menu := Menu("Subgroup Lattice");
  for m in latticeType.menus do
    knownArgs := [poset];
    if m.group = true then
      Add(knownArgs, G);
    fi;
    cb := Callback(m.func, knownArgs);
    Add(menu, Menu(m.name, cb));
  od;
  Add(poset, menu);

  return Draw(poset);
end);

##
## without defaults record:
##
InstallOtherMethod(GraphicSubgroupLattice,
    "for a group",
    true,
    [IsGroup],
    0,
function(G)
  return GraphicSubgroupLattice(G, rec());
end);

#############################################################################
##
##  Decision function for subgroup lattice type:
##
#############################################################################

#############################################################################
##
#M  DecideSubgroupLatticeType(<grp>)  . . decides about the type of a lattice
##
##  This operation is called while creation of a new graphic subgroup lattice.
##  It has to decide about the type of the lattice. That means it has to
##  decide 5 questions:
##   1) Are all levels known right from the beginning?
##   2) Has the lattice the HasseProperty?
##   3) Can we test two subgroups for equality reasonably cheaply?
##   4) Shall we create a vertex for the trivial subgroup at the beginning?
##   5) What menu operations are possible?
##   6) What information is displayed on RightClick?
##  Returns a record. The first four entries are boolean values for  questions
##  1-4. Note that if the answer to 2 is true, then the answer to 3 must also
##  be true. The fifth and sixth entry are configuration lists as explained
##  in the configuration section of "ilatgrp.gi" for menu operations and
##  info displays respectively.
##
##  The following is the default "fallback" method suitable for reasonably
##  small finite groups.
##
InstallMethod(DecideSubgroupLatticeType,
    "for a group",
    true,
    [IsGroup],
    0,
function(G)
  local knowslevels, menus;
  menus:=[];
  Append(menus, MenuOpsForFiniteGroups);
  Append(menus, MenuOpsCommon);
  if Size(G) > 10^17 then # that is just heuristic!
    knowslevels := false;
  else
    knowslevels := Length(DivisorsInt(Size(G))) < 50;
  fi;
  return rec(
    knowsLevels  := knowslevels,
    hasse        := true, # we assume HasseProperty
    canCompare   := true, # we assume we can compare groups
    trivial      := true, # we want the trivial subgroup
    menus        := menus,
    contextMenus := ContextMenusForFiniteGroups,
  );
end);

## for finitely presented groups:
InstallMethod(DecideSubgroupLatticeType,
    "for a group",
    true,
    [IsGroup and IsFpGroup],
    0,

function(G)
  local menus;
  menus:=[];
  Append(menus, MenuOpsForFpGroups);
  Append(menus, MenuOpsCommon);
  return rec(
    knowsLevels  := false, # we create levels dynamically
    hasse        := false, # we assume HasseProperty
    canCompare   := false, # we assume we can compare groups
    trivial      := false, # we want the trivial subgroup
    menus        := menus,
    contextMenus := ContextMenusForFpGroups,
  );
end);

max_stars_repo_name: mcmartins/subgroup-lattice

max_issues_repo_name: mcmartins/subgroup-lattice

max_forks_repo_name: mcmartins/subgroup-lattice


text:
#
# Read("~/Workspace/groupsSB/epi/A2/case.2.gi");
#

type:="A";
rank:=2;
nr_pos_roots:=3;

Read("~/Workspace/groupsSB/epi/group.gi");

highest_weight:=[1,1];# corresponding highest weight module is the Lie algebra

Read("~/Workspace/groupsSB/epi/ext.gi");

#
# p>2
#

getRep:=function(k,t)
    local u,ua,ub,uab;
    ua:=ext_root_group(k,1,t);
    #ub:=ext_root_group(k,2,t);
    uab:=ext_root_group(k,3,t);
    u:=ua*uab;

    return u;
end;

#
# k:=3; # power of exterior algebra
#
handle:=function(k)
    local id_mat,u,e,h,h1,h2,ese,evh1,esh1,evh2,esh2,esh01,index_ev01,esh02,index_ev02,si,si1,si2,inter;
    u:=getRep(k,1);
    id_mat:=u^0;
    e:=u-id_mat;
    ese:=Eigenspaces(Rationals,e)[1]; #eigenvalues 0, centralizer

    # torus
    # centralizer of T1
    h1:=ext_e(k,cb0[7]);
    esh1:=Eigenspaces(Rationals,h1);
    evh1:=Eigenvalues(Rationals,h1);
    esh01:=[];
    index_ev01:=Position(evh1,0);
    esh01:=esh1[index_ev01];

    # centralizer of T2
    h2:=ext_e(k,cb0[8]);
    esh2:=Eigenspaces(Rationals,h2);
    evh2:=Eigenvalues(Rationals,h2);
    esh02:=[];
    index_ev02:=Position(evh2,0);
    esh02:=esh2[index_ev02];

    #
    # sum and intersection
    #
    #return [ese,esh01,esh02];
    si1:=SumIntersectionMat(Basis(ese),Basis(esh01));
    si2:=SumIntersectionMat(Basis(ese),Basis(esh02));
    si:=SumIntersectionMat(si1[2],si2[2]);
    inter:=si[2];
    return inter;
    #return [inter,esh01,esh02];
end;

testFixedPointsTrous:=function(k)
    local result1,result2,h1,esh1,evh1,esh01,index_ev01,h2,esh2,evh2,esh02,index_ev02;
    # torus
    # centralizer of T1
    h1:=ext_e(k,cb0[7]);
    esh1:=Eigenspaces(Rationals,h1);
    evh1:=Eigenvalues(Rationals,h1);
    esh01:=[];
    index_ev01:=Position(evh1,0);
    esh01:=esh1[index_ev01];

    result1:=List(Basis(esh01),b->ext_lin_comb(k,b));
    result1:=List(result1,i->cb0[7]^i);

    # centralizer of T2
    h2:=ext_e(k,cb0[8]);
    esh2:=Eigenspaces(Rationals,h2);
    evh2:=Eigenvalues(Rationals,h2);
    esh02:=[];
    index_ev02:=Position(evh2,0);
    esh02:=esh2[index_ev02];

    result2:=List(Basis(esh02),b->ext_lin_comb(k,b));
    result2:=List(result2,i->cb0[8]^i);

    return [result1,result2];
end;

max_stars_repo_name: iuliansimion/groupsSB

max_issues_repo_name: iuliansimion/groupsSB

max_forks_repo_name: iuliansimion/groupsSB


text:
#!/usr/bin/gap

# https://cp4space.wordpress.com/2020/05/10/minimalistic-quantum-computation/

Print("running universal.gap\n");;

r2 := Sqrt(2);;
ir2 := 1/r2;;

i := [[1, 0], [0, 1]];;
w := [[E(4), 0], [0, E(4)]];;
x := [[0, 1], [1, 0]];;
z := [[1, 0], [0, -1]];;
s := [[1, 0], [0, E(4)]];;
h := [[ir2, ir2], [ir2, -ir2]];;

Cliff1 := Group(w, s, h);; # Order 192
Pauli1 := Group(w, x, z);; # Order 32

xi := KroneckerProduct(x, i);;
ix := KroneckerProduct(i, x);;
zi := KroneckerProduct(z, i);;
iz := KroneckerProduct(i, z);;
si := KroneckerProduct(s, i);;
is := KroneckerProduct(i, s);;
hi := KroneckerProduct(h, i);;
ih := KroneckerProduct(i, h);;
wi := KroneckerProduct(w, i);;

cz := [
    [1, 0, 0, 0],
    [0, 1, 0, 0],
    [0, 0, 1, 0],
    [0, 0, 0, -1]];;

Cliff2 := Group(si, is, hi, ih, wi, cz);; # Order 92160
#for g in Cliff2 do Print(g, "\n"); od;
Pauli2 := Group(wi, xi, ix, zi, iz);;

xii := KroneckerProduct(xi, i);;
ixi := KroneckerProduct(i, xi);;
iix := KroneckerProduct(i, ix);;

zii := KroneckerProduct(zi, i);;
izi := KroneckerProduct(i, zi);;
iiz := KroneckerProduct(i, iz);;

sii := KroneckerProduct(si, i);;
isi := KroneckerProduct(i, si);;
iis := KroneckerProduct(i, is);;

hii := KroneckerProduct(hi, i);;
ihi := KroneckerProduct(i, hi);;
iih := KroneckerProduct(i, ih);;

wii := KroneckerProduct(wi, i);;

icz := KroneckerProduct(i, cz);;
czi := KroneckerProduct(cz, i);;

Tofolli := [
    [1, 0, 0, 0, 0, 0, 0, 0],
    [0, 1, 0, 0, 0, 0, 0, 0],
    [0, 0, 1, 0, 0, 0, 0, 0],
    [0, 0, 0, 1, 0, 0, 0, 0],
    [0, 0, 0, 0, 1, 0, 0, 0],
    [0, 0, 0, 0, 0, 1, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 1],
    [0, 0, 0, 0, 0, 0, 1, 0]];;

zhh := zii*ihi*iih;
hzh := izi*hii*iih;
hhz := iiz*hii*ihi;

xhh := xii*ihi*iih;
hxh := ixi*hii*iih;
hhx := iix*hii*ihi;

Gx := Group(Tofolli, xhh, hxh, hhx);
Gz := Group(Tofolli, zhh, hzh, hhz);

Print("Order(Gz) = ", Order(Gz), "\n");
Print("Order(Gx) = ", Order(Gx), "\n");

L7:= SimpleLieAlgebra("E", 7, Rationals);;
R7:= RootSystem(L7);;
W7:= WeylGroup(R7);
Print("Order(W7) = ", Order(W7), "\n");

f := IsomorphismGroups(Gz, W7);
Print("IsomorphismGroups(Gz, E7) = ", f, "\n");

L8:= SimpleLieAlgebra("E", 8, Rationals);;
R8:= RootSystem(L8);;
W8:= WeylGroup(R8);
Print("Order(W8) = ", Order(W8), "\n");

f := IsomorphismGroups(Gx, W8);
Print("IsomorphismGroups(Gx, E8) = ", f, "\n");

max_stars_repo_name: punkdit/bruhat

max_issues_repo_name: punkdit/bruhat

max_forks_repo_name: punkdit/bruhat


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

canReorder := e -> not IsBound(e.doNotReorder);
Class(RulesSums, RuleSet);

RewriteRules(RulesSums, rec(
# Drop_DPWrapper := Rule([DPWrapper, @(1), ...], e->@(1).val),
 Drop_GrpISum := Rule([Grp, @(1, ISum)], e -> ISum(@(1).val.var, @(1).val.domain, @(1).val.child(1))),
 Drop_GrpInplaceISum := Rule([Grp, [Inplace, @(1, ISum)]], e -> Inplace(ISum(@(1).val.var, @(1).val.domain, @(1).val.child(1)))),

 # ================================================================
 # Associativity of Compose and SUM
 # ================================================================
 ComposeAssoc := ARule( Compose, [ @(1,Compose) ],  e -> @(1).val.children() ),
 SUMAssoc     := ARule( SUM,     [ @(1,SUM)     ],  e -> @(1).val.children() ),
 SUMAccAssoc  := ARule( SUMAcc,  [ @(1,SUMAcc)  ],  e -> @(1).val.children() ),
 TensorAssoc  := ARule( Tensor,  [ @(1,Tensor)  ],  e -> @(1).val.children() ),
 RedundandSUMAcc := Rule([SUMAcc, @(1)], e -> @(1).val),

 # ================================================================
 # Identity Matrix (I)
 # ================================================================
 Compose_IR := ARule(Compose, [@(1), I], e->[@(1).val]),
 Compose_IL := ARule(Compose, [I, @(1)], e->[@(1).val]),
 # rule Prm(fId()) -> I may break vectorization
 fId_toI := Rule([@(1, [Gath, Scat]), [fId, @(2)]], e -> I(@(2).val)),

 # ================================================================
 # Zero Matrix (O)
 # ================================================================
 ZeroSUMAcc1 := ARule(SUMAcc, [O, @(1)], e -> [@(1).val]),
 ZeroSUMAcc2 := ARule(SUMAcc, [@(1), O], e -> [@(1).val]),
 ZeroSUM1    := ARule(SUM, [O, @(1)], e -> [@(1).val]),
 ZeroSUM2    := ARule(SUM, [@(1), O], e -> [@(1).val]),
 ZeroTensor := Rule([Tensor, ..., @(1,O), ...], e -> O(Rows(e), Cols(e))),

 # ================================================================
 # Combining Gath, Scat, and Prm
 # ================================================================
 ComposeGathGath := ARule(Compose, [ @(1, Gath), @(2, [Gath, Prm]) ], # o 1-> 2->
     e -> [ Gath(fCompose(@(2).val.func, @(1).val.func)) ]),

 ComposePrmPrm := ARule(Compose, [ @(1, Prm), @(2, Prm) ], # o 1-> 2->
     e -> [ Prm(fCompose(@(2).val.func, @(1).val.func)) ]),

 ComposeScatScat := ARule(Compose, [ @(1, [Scat, ScatAcc]), @(2, [Scat, ScatAcc]) ], # <-1 <-2 o
     e -> [ Cond(ObjId(@(1).val)=ScatAcc or ObjId(@(2).val)=ScatAcc,
                 ScatAcc(fCompose(@(1).val.func, @(2).val.func)),
                 Scat   (fCompose(@(1).val.func, @(2).val.func))) ]),

 ComposePrmScat  := ARule(Compose, [ @(1, Prm), @(2, [Scat, ScatAcc]) ], # 1-> <-2 o
     e -> [ ObjId(@(2).val)(fCompose(@(1).val.func.transpose(), @(2).val.func)) ]),

 ComposeDPrmDPrm   := ARule(Compose, [ @(1, DelayedPrm), @(2, DelayedPrm) ], # direct:  o 1-> 2->
     e -> [ DelayedPrm(fCompose(@(2).val.func, @(1).val.func)) ]),

 DropGathScat := ARule( Compose, [ @(1, Gath), @(2, Scat, x->x.transpose()=@(1).val) ],
     e -> []),

 # ================================================================
 # ScatGath
 # NOTE: YSV ScatGath is ugly and should be removed
 # ================================================================
 ComposeScat_ScatGath := ARule(Compose, [ @(1, Scat), @(2, ScatGath) ], # <-1 <-2 o
     e -> [ ScatGath(fCompose(@(1).val.func, @(2).val.sfunc), @(2).val.gfunc) ]),

 ComposeScatGath_Gath := ARule(Compose, [ @(1, ScatGath), @(2, Gath) ], # o 1-> 2->
     e -> [ ScatGath(@(1).val.sfunc, fCompose(@(2).val.func, @(1).val.gfunc)) ]),

 # ================================================================
 # Using distributivity to pull Prm, Gath, Scat into ISum, SUM, etc.
 # ================================================================
 PullInRight := ARule( Compose,
       [ @(1, [Prm, Scat, ScatAcc, TCast, PushR, PushLR, Conj, ConjL, ConjR, ConjLR, FormatPrm]),
         @(2, [RecursStep, Grp, BB, SUM, Buf, ISum, Data, COND, NoDiagPullin, NoDiagPullinLeft, NoDiagPullinRight, NeedInterleavedComplex]) ],
  e -> [ CopyFields(@(2).val, rec(
             _children :=  List(@(2).val._children, c -> @(1).val * c),
             dimensions := [Rows(@(1).val), Cols(@(2).val)] )) ]),

 PullInLeft := ARule( Compose,
       [ @(1, [RecursStep, Grp, BB, SUM, SUMAcc, Buf, ISum, ISumAcc, Data, COND, NoDiagPullin, NoDiagPullinLeft, NoDiagPullinRight, NeedInterleavedComplex]),
         @(2, [Prm, Gath, TCast, PushL, PushLR, Conj, ConjL, ConjR, ConjLR, FormatPrm]) ],
     e -> [ CopyFields(@(1).val, rec(
                _children := List(@(1).val._children, c -> c * @(2).val),
                dimensions := [Rows(@(1).val), Cols(@(2).val)] )) ]),

 ComposePushRPrmScat := ARule(Compose, [ @(1, PushR, e -> ObjId(e.child(1)) = Prm), @(2, Scat) ],
    e -> [ Scat(fCompose(@(1).val.child(1).func.transpose(), @(2).val.func)) ]),

 # Match Gath * PushL(Prm) structure. The normal Gath * PushL rule occasionally
 # does not apply properly -- not clear why! (MRT)
 ComposeGathPushLPrm := ARule(Compose, [ @(1, Gath), @(2, PushL, e -> ObjId(e.child(1)) = Prm) ],
    e -> [ Gath(fCompose(@(2).val.child(1).func, @(1).val.func)) ]),

 # ================================================================
 # Low Priority
 # ================================================================
 ComposeScatPrm  := ARule(Compose, [ @(1, [Scat, ScatAcc]), @(2, Prm) ],  # <-1 2-> o
     e -> [ ObjId(@(1).val)(fCompose(@(1).val.func, @(2).val.func.transpose())) ]),

 ComposePrmGath  := ARule(Compose, [ @(1, Prm), @(2, Gath) ], # o 1-> 2->
     e -> [ Gath(fCompose(@(2).val.func, @(1).val.func)) ]),

 # ================================================================
 # TCast
 # ================================================================
 ComposeConvScat  := ARule(Compose, [ @(1, TCast), @(2, [Scat, ScatAcc]) ],  # <-1 2-> o
     e -> [ @(2).val, TCast(Cols(@(2).val), @(1).val.params[2], @(1).val.params[3]) ]),

 ComposeGathConv  := ARule(Compose, [ @(1, Gath), @(2, TCast) ], # o 1-> 2->
     e -> [ TCast(Rows(@(1).val), @(2).val.params[2], @(2).val.params[3]), @(1).val ]),

 # ================================================================
 # ConjTranspose
 # ================================================================
 InertConjTranspose_@ := Rule(
     [InertConjTranspose, @.cond(e->not (IsBound(e.isInertConjTranspose) and e.isInertConjTranspose()))],
     e -> e.child(1).conjTranspose()),
));

# ================================================================
# COND Rules
# ================================================================

# NOTE: this replication should not be necessary
# Rule([@(0,SUMAcc), ..., @(1,COND), ...],
#     e -> let( ch := e.children(),
#         left := Left(...),
#         right := Right(...),
#         leftch := ch{[1..left]},
#         rightch := ch{[right..Length(ch)]},
#         OP := ObjId(@(0).val),
#         COND(@(1).val.cond,
#         List(@(1).val.children(),
#              c -> OP(Concatenation(Copy(leftch), [c], Copy(rightch)))))));

# ================================================================
# OL Rules
# ================================================================

Class(OLRules, RuleSet);
RewriteRules(OLRules, rec(
  Gath_fCross := Rule([Gath, @(1, fCross)], e -> Cross(List(@(1).val.rChildren(),x->Gath(x)))),

  fCross_fCross := ARule(fCompose, [@(1,fCross), @(2,fCross)],
      e-> [fCross(List([1..Length(@(1).val.rChildren())],i->fCompose(@(1).val.rChildren()[i], @(2).val.rChildren()[i])))]),

  f2DTensor_f2DTensor := ARule(fCompose, [@(1,f2DTensor), @(2,f2DTensor)],
      e-> [f2DTensor(List([1..Length(@(1).val.rChildren())],i->fCompose(@(1).val.rChildren()[i], @(2).val.rChildren()[i])))]),

  f4DTensor_f4DTensor := ARule(fCompose, [@(1,f4DTensor), @(2,f4DTensor)],
      e-> [f4DTensor(List([1..Length(@(1).val.rChildren())],i->fCompose(@(1).val.rChildren()[i], @(2).val.rChildren()[i])))]),

  f2DTr_f2DTensor := ARule(fCompose, [@(1,f2DTr), @(2,f2DTensor)],
      e-> let(f:=f2DTensor(Reversed(@(2).val.rChildren())),
              [f, ApplyFunc(ObjId(@(1).val), Reversed(f.advdomain()[1]))])),

  f2DTr_f2DTr := ARule(fCompose, [@(1,f2DTr), @(2,f2DTr)],
      e-> []),

  f2DTensor_f2DTrExplicit := ARule(fCompose, [@(1,f2DTensor), @(2, [f2DTrExplicitL, f2DTrExplicitR])],
      e-> let(f:=f2DTensor(Reversed(@(1).val.rChildren())),
              [ApplyFunc(ObjId(@(2).val), Reversed(f.advrange()[1])), f])),

  f2DTr_f2DTrExplicit := ARule(fCompose, [@(1,f2DTr), @(2,[f2DTrExplicitL, f2DTrExplicitR])],
      e-> []),

  f2DTrExplicitRL := ARule(fCompose, [@(1,f2DTrExplicitR), @(2,f2DTrExplicitL)],
      e-> []),

  ParSeq_ParSeq := Rule([@(1, ParSeq), ..., @(2, ParSeq, x -> x.fb_cnt=@(1).val.fb_cnt), ...],
      e-> ParSeq(e.fb_cnt, Flat(List(e.children(), x -> Cond( x _is ParSeq and x.fb_cnt=e.fb_cnt, x.children(), x))))),
));


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

#F CondPat(<obj>, <pat1>, <result1>, ..., <patN>, <resultN>, [<resultElse>])
#F
#F CondPat is a conditional control-flow construct similar to Cond, but instead
#F of normal boolean conditions it uses patterns that are matched against its
#F first argument <obj>. CondPat returns the result that corresponds to the
#F first matched pattern. If nothing matches <resultElse> is returned.
#F
#F Note: resultElse is optional
#F
#F CondPat is implemented as a function with delayed evaluation, and will
#F only evaluate the result that corresponds to the matched patterns, all
#F other results will not be evaluated, for example:
#F
#F CondPat(mul(X,2), [mul, @, 2], Print("yes"), [add, @, 2], Print("no"));
#F
#F will only execute Print("yes") -- as expected.
#F

CondPat := UnevalArgs(function(arg)
    local usage, i, pat, obj;

	usage := "CondPat(<obj>, <pat1>, <result1>, ..., <patN>, <resultN>, <resultElse>)";

	if Length(arg) < 3 then
		Error(usage);
	fi;

    obj := Eval(arg[1]);
    i := 2;
    while i <= Length(arg)-1 do
        pat := Eval(arg[i]);
		if PatternMatch(obj, pat, empty_cx()) then
			return Eval(arg[i+1]);
		fi;
		i := i+2;
    od;
    if i=Length(arg) then
		return Eval(Last(arg));
    else
		Error("CondPat: No 'else' result, and no patterns match");
    fi;
end);


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

Import(simt);

_isSIMTTag := tag -> IsBound(tag.isSIMTTag) and tag.isSIMTTag;
_toSIMTDim := (tags, n) -> Cond(
    ObjId(tags[1]) = ASIMTLoopDim, tags[1],
    IsBound(tags[1].params[1]) and not IsInt(tags[1].params[1]), ApplyFunc(ObjId(tags[1]), [tags[1].params[1](n)]),
    tags[1](n));

NewRulesFor(TIterHStack, rec(
    TIterHStack_SIMT := rec(
        applicable := nt->nt.hasTags() and _isSIMTTag(nt.firstTag()),
        forTransposition := false,
        children := nt -> [[ nt.params[1].withTags(Drop(nt.getTags(), 1)), InfoNt(nt.params[2]) ]],
        apply := (nt, c, cnt) -> SIMTIterHStack1(_toSIMTDim(nt.getTags(), cnt[2].params[1].range), cnt[2].params[1], cnt[2].params[1].range, c[1])
    )
));

NewRulesFor(TIterVStack, rec(
    TIterVStack_SIMT := rec(
        applicable := nt->nt.hasTags() and _isSIMTTag(nt.firstTag()),
        forTransposition := false,
        children := nt -> [[ nt.params[1].withTags(Drop(nt.getTags(), 1)), InfoNt(nt.params[2]) ]],
        apply := (nt, c, cnt) -> SIMTIterVStack(_toSIMTDim(nt.getTags(), cnt[2].params[1].range), cnt[2].params[1], cnt[2].params[1].range, c[1])
    )
));

NewRulesFor(MDDFT, rec(
    MDDFT_RowCol_3D_SIMT := rec(
        info :="MDDFT(n_1,n_2,...,n_t) = (RC(MDDFT(n_1,n_2,...,n_{t-1}) x I(.))) (I(n_1*n_2*...*n_{t-1}) x PRDFT1(n_t))",
        applicable := nt->nt.hasTags() and Length(nt.getTags()) = 2 and ForAll(nt.getTags(), _isSIMTTag) and Length(nt.params[1]) = 3,
        children  := nt -> let(a_lengths := nt.params[1],
                               a_exp := nt.params[2],
                               tags := nt.getTags(),
                               taglen := Length(tags),
                               rtags1 := tags{[2]},
                               [[ MDDFT(DropLast(a_lengths, 1), a_exp).withTags(rtags1),
                                  DFT(Last(a_lengths), a_exp) ]]),

        apply := (nt, C, cnt) -> let(
            a := nt.params[1],
            n1 := a[1],
            n2 := a[2],
            n3 := a[3],
            SIMTTensor(_toSIMTDim(nt.getTags(){[1]}, n3), C[1], I(n3)) *
                                     SIMTTensor(_toSIMTDim(nt.getTags(){[1]}, n1), I(n1),
                                         SIMTTensor(_toSIMTDim(nt.getTags(){[2]}, n2), I(n2),
                                             C[2])))
    ),
    MDDFT_tSPL_RowCol_SIMT := rec(
        info :="MDDFT(n_1,n_2,...,n_t) = (RC(MDDFT(n_1,n_2,...,n_{t-1}) x I(.))) (I(n_1*n_2*...*n_{t-1}) x PRDFT1(n_t))",
        applicable := nt->nt.hasTags() and ForAll(nt.getTags(), _isSIMTTag) and Length(nt.params[1]) > 1,
        children  := nt -> let(a_lengths := nt.params[1],
                               a_exp := nt.params[2],
                               tags := nt.getTags(),
                               [[ TTensorI(MDDFT(DropLast(a_lengths, 1), a_exp), Last(a_lengths), AVec, AVec).withTags(tags),
                                  FoldR(DropLast(a_lengths, 1), (a,b)->TTensorI(a, b, APar, APar), DFT(Last(a_lengths), a_exp)).withTags(tags) ]]),
        apply := (nt, C, cnt) -> C[1] * C[2]
    ),
    MDDFT_tSPL_Pease_SIMT := rec(
        applicable := nt->nt.hasTags() and ForAll(nt.getTags(), _isSIMTTag) and Length(nt.params[1]) > 1,
        children  := nt -> let(a_lengths := nt.params[1],
                               a_exp := nt.params[2],
                               tags := nt.getTags(),
                               [ [TCompose(List(nt.params[1], i->TTensorI(DFT(i, a_exp), Product(nt.params[1])/i, AVec, APar))).withTags(tags) ]]),
        apply := (nt, C, cnt) -> C[1]
    )
));

NewRulesFor(MDPRDFT, rec(
    MDPRDFT_3D_SIMT := rec(
        info :="MDDFT(n_1,n_2,...,n_t) = (RC(MDDFT(n_1,n_2,...,n_{t-1}) x I(.))) (I(n_1*n_2*...*n_{t-1}) x PRDFT1(n_t))",
        applicable := nt->nt.hasTags() and Length(nt.getTags()) = 2 and ForAll(nt.getTags(), _isSIMTTag) and Length(nt.a_lengths()) = 3,
        # Generally, children is list of ordered lists of all possible children.
        # 2 children, both with same parameter a_exp():
        # [1] MDDFT on all but last dimension;
        # [2] PRDFT1 on last dimension, the fastest varying dimension.
        children  := nt -> let(a_lengths := nt.a_lengths(),
                               a_exp := nt.a_exp(),
                               tags := nt.getTags(),
                               taglen := Length(tags),
                               rtags1 := tags{[2]},
                               [[ MDDFT(DropLast(a_lengths, 1), a_exp).withTags(rtags1),
                                  PRDFT1(Last(a_lengths), a_exp) ]]),

        # nonterminal, children, children non-terminals
        # NOT YET CORRECT
        apply := (nt, C, cnt) -> let(a_lengths := nt.a_lengths(),
                                     Clast := C[2],
                                     Crest := C[1],
                                     Nlastcomplex := Clast.dims()[1]/2,
                                     Nrest := DropLast(a_lengths, 1),
                                     Nallbutlast := Product(Nrest),
                                     # Apply Clast on last dimension,
                                     # then Crest on all dimensions but last.
                                     RC(SIMTTensor(_toSIMTDim(nt.getTags(){[1]}, Nlastcomplex), Crest, I(Nlastcomplex))) *
                                     SIMTTensor(_toSIMTDim(nt.getTags(){[1]}, Nrest[1]), I(Nrest[1]),
                                         SIMTTensor(_toSIMTDim(nt.getTags(){[2]}, Nrest[2]), I(Nrest[2]),
                                             Clast)))
    ),
    MDPRDFT_tSPL_RowCol_SIMT := rec(
        info :="MDDFT(n_1,n_2,...,n_t) = (RC(MDDFT(n_1,n_2,...,n_{t-1}) x I(.))) (I(n_1*n_2*...*n_{t-1}) x PRDFT1(n_t))",
        applicable := nt->nt.hasTags() and ForAll(nt.getTags(), _isSIMTTag) and Length(nt.params[1]) > 1,
        children  := nt -> let(a_lengths := nt.params[1],
                               a_exp := nt.params[2],
                               tags := nt.getTags(),
                               mddft := MDDFT(DropLast(a_lengths, 1), a_exp),
                               prdft := PRDFT1(Last(a_lengths), a_exp),
                               [[ TTensorI(mddft, prdft.dims()[1]/2, AVec, AVec).withTags(tags),
                                  FoldR(DropLast(a_lengths, 1), (a,b)->TTensorI(a, b, APar, APar), prdft).withTags(tags) ]]),
        apply := (nt, C, cnt) -> RC(C[1]) * C[2]
    ),
    MDPRDFT_tSPL_Pease_SIMT := rec(
        applicable := nt->nt.hasTags() and ForAll(nt.getTags(), _isSIMTTag) and Length(nt.params[1]) > 1,
        children  := nt -> let(a_lengths := nt.params[1],
                               a_exp := nt.params[2],
                               tags := nt.getTags(),
                               prdft := PRDFT1(Last(a_lengths), a_exp),
                               rcdim := Rows(prdft),
                               [ [ TCompose(List(DropLast(nt.params[1], 1), i->TRC(TTensorI(DFT(i, a_exp), rcdim * Product(DropLast(nt.params[1], 1))/(2*i), AVec, APar)))::
                                           [ TGrp(TCompose([TL(rcdim * Product(DropLast(nt.params[1], 1)) / 2, rcdim / 2, 1, 2),
                                             TTensorI(PRDFT1(Last(a_lengths), a_exp), Product(DropLast(nt.params[1], 1)), APar, APar)])) ]).withTags(tags) ]] ),
        apply := (nt, C, cnt) -> C[1]
    )

));

NewRulesFor(IMDPRDFT, rec(
    IMDPRDFT_3D_SIMT := rec(
        info :="IMDDFT(n_1,n_2,...,n_t) = (RC(IMDDFT(n_1,n_2,...,n_{t-1}) x I(.))) (I(n_1*n_2*...*n_{t-1}) x PRDFT1(-n_t))",
        applicable := nt->nt.hasTags() and Length(nt.getTags()) = 2 and ForAll(nt.getTags(), _isSIMTTag) and Length(nt.a_lengths()) = 3,
        # Generally, children is list of ordered lists of all possible children.
        # 2 children, both with same parameter a_exp():
        # [1] IMDDFT on all but last dimension;
        # [2] IPRDFT1 on last dimension, the fastest varying dimension.
        children  := nt -> let(a_lengths := nt.a_lengths(),
                               a_exp := nt.a_exp(),
                               tags := nt.getTags(),
                               taglen := Length(tags),
                               rtags1 := tags{[2]},
                               [[ MDDFT(DropLast(a_lengths, 1), a_exp).withTags(rtags1),
                                  IPRDFT1(Last(a_lengths), a_exp) ]]),

        # nonterminal, children, children non-terminals
        # NOT YET CORRECT
        apply := (nt, C, cnt) -> let(a_lengths := nt.a_lengths(),
                                     Clast := C[2],
                                     Crest := C[1],
                                     Nlastcomplex := Clast.dims()[2]/2,
                                     Nrest := DropLast(a_lengths, 1),
                                     Nallbutlast := Product(Nrest),
                                     # Apply Crest on all dimensions but last,
                                     # then Clast on last dimension.
                                     SIMTTensor(_toSIMTDim(nt.getTags(){[1]}, Nrest[1]), I(Nrest[1]),
                                         SIMTTensor(_toSIMTDim(nt.getTags(){[2]}, Nrest[2]), I(Nrest[2]),
                                             Clast)) *
                                     RC(SIMTTensor(_toSIMTDim(nt.getTags(){[1]}, Nlastcomplex), Crest, I(Nlastcomplex))))
    ),
    IMDPRDFT_tSPL_RowCol_SIMT := rec(
        info :="IMDDFT(n_1,n_2,...,n_t) = (RC(IMDDFT(n_1,n_2,...,n_{t-1}) x I(.))) (I(n_1*n_2*...*n_{t-1}) x PRDFT1(-n_t))",
        applicable := nt->nt.hasTags() and ForAll(nt.getTags(), _isSIMTTag) and Length(nt.params[1]) > 1,
        children  := nt -> let(a_lengths := nt.params[1],
                               a_exp := nt.params[2],
                               tags := nt.getTags(),
                               mddft := MDDFT(DropLast(a_lengths, 1), a_exp),
                               iprdft := IPRDFT1(Last(a_lengths), a_exp),
                               [[ FoldR(DropLast(a_lengths, 1), (a,b)->TTensorI(a, b, APar, APar), iprdft).withTags(tags),
                                   TTensorI(mddft, iprdft.dims()[2]/2, AVec, AVec).withTags(tags) ]]),
        apply := (nt, C, cnt) -> C[1] * RC(C[2])
    ),
    IMDPRDFT_tSPL_Pease_SIMT := rec(
        applicable := nt->nt.hasTags() and ForAll(nt.getTags(), _isSIMTTag) and Length(nt.params[1]) > 1,
        children  := nt -> let(a_lengths := nt.params[1],
                               a_exp := nt.params[2],
                               tags := nt.getTags(),
                               iprdft := IPRDFT1(Last(a_lengths), a_exp),
                               rdim := Rows(iprdft),
                               cdim := Cols(iprdft),
                               [ [ TCompose([ TGrp(TCompose([
                                             TTensorI(IPRDFT1(Last(a_lengths), a_exp), Product(DropLast(a_lengths, 1)), APar, APar),
                                             TL(cdim * Product(DropLast(a_lengths, 1)) / 2, Product(DropLast(a_lengths, 1)), 1, 2),
                                       ])) ] ::
                                       Reversed(List(DropLast(a_lengths, 1), i->TRC(TTensorI(DFT(i, a_exp), cdim * Product(DropLast(a_lengths, 1))/(2*i), APar, AVec))))
                                    ).withTags(tags) ]] ),
        apply := (nt, C, cnt) -> C[1]
    )

));

NewRulesFor(TRC, rec(
    TRC_SIMT := rec(
        forTransposition := false,
        applicable := nt->nt.hasTags() and _isSIMTTag(nt.firstTag()),
        children := nt -> [[ nt.params[1].withTags(nt.getTags()) ]],
        apply := (nt, c, cnt) -> RC(c[1])
    )
));

NewRulesFor(TL, rec(
#   TL(N,n,l,r,[]) -> I_l x L(N,n) x I_r
    L_SIMT := rec(
        forTransposition := false,
        applicable := nt->nt.hasTags() and _isSIMTTag(nt.firstTag()),
        apply := (nt, c, cnt) -> let(
            c1 := When(nt.params[3]=1, [], [I(nt.params[3])]),
            c2 := When(nt.params[4]=1, [], [I(nt.params[4])]),
            Tensor(Concat(c1, [ L(nt.params[1], nt.params[2]) ], c2))
        )
    )
));

NewRulesFor(TTensorInd, rec(
#   base cases
#   I x A
    TTensorInd_SIMT := rec(
        info := "IxA base",
        forTransposition := false,
        applicable := nt -> nt.hasTags() and _isSIMTTag(nt.firstTag()) and IsParPar(nt.params),
        children := nt -> [[ nt.params[1].withTags(Drop(nt.getTags(), 1)), InfoNt(nt.params[2]) ]],
        apply := (nt, c, cnt) -> SIMTIDirSum(_toSIMTDim(nt.getTags(), cnt[2].params[1].range), cnt[2].params[1], cnt[2].params[1].range, c[1])
    )
));

NewRulesFor(TTensorI, rec(
#   base cases
#   I x A
    IxA_SIMT := rec(
        info := "IxA base",
        forTransposition := false,
        applicable := nt -> nt.hasTags() and _isSIMTTag(nt.firstTag()) and  IsParPar(nt.params) and nt.params[2] > 1,
        children := nt -> [[ nt.params[1].withTags(Drop(nt.getTags(), 1)) ]],
        apply := (nt, c, cnt) -> #When(nt.params[2] > 1,
            SIMTTensor(_toSIMTDim(nt.getTags(), nt.params[2]), I(nt.params[2]), c[1])#,
            #c[1])
    ),
#   A x I
    AxI_SIMT := rec(
        info := "AxI base",
        forTransposition := false,
        applicable := nt -> nt.hasTags() and _isSIMTTag(nt.firstTag()) and IsVecVec(nt.params) and nt.params[2] > 1,
        children := nt -> [[ nt.params[1].withTags(Drop(nt.getTags(), 1)) ]],
        apply := (nt, c, cnt) -> #When(nt.params[2] > 1,
            SIMTTensor(_toSIMTDim(nt.getTags(), nt.params[2]), c[1], I(nt.params[2]))#,
            #c[1])
    ),
#   A x I
    trivialLoop := rec(
        info := "AxI base",
        forTransposition := false,
        applicable := nt -> nt.hasTags() and _isSIMTTag(nt.firstTag()) and nt.params[2] = 1,
        children := nt -> [[ nt.params[1].withTags(DropLast(nt.getTags(), 1)) ]],
        apply := (nt, c, cnt) -> c[1]
    ),
#   L (I x A)
    L_IxA_SIMT := rec(
        forTransposition := false,

        # these config parameters need to be moved into the opts...
        mem := 1024*96,
        mem_per_pt := 2*8*2,
        max_threads := 2048,
        max_kernel := 18 * 18,
        _peelof := (self,n,m) >> Maximum(Filtered(self.mem_per_pt * Filtered(n*DivisorsInt(m), e-> e<self.max_threads),
            f -> f < When(n >= self.max_kernel, self.mem/2, self.mem)))/(self.mem_per_pt*n),

        applicable := nt -> nt.hasTags() and _isSIMTTag(nt.firstTag()) and IsVecPar(nt.params) and nt.params[2] > 1,
        children := (self, nt) >> let(n := Rows(nt.params[1]), m:= nt.params[2], peelof := self._peelof(n,m), remainder := m/peelof,
            [[ TCompose([TL(Rows(nt)/peelof, n, 1, peelof),
                TTensorI(
                    TCompose([TL(Rows(nt.params[1]) * peelof, Rows(nt.params[1])), TTensorI(nt.params[1], peelof, APar, APar)]),
                    remainder, APar, APar)]).withTags(nt.getTags()) ]]),
        apply := (nt, c, cnt) -> c[1]
    ),
#   (I x A) L
    IxA_L_SIMT := rec(
        forTransposition := false,

        # these config parameters need to be moved into the opts...
        mem := 1024*96,
        mem_per_pt := 2*8*2,
        max_threads := 2048,
        max_kernel := 18 * 18,
        _peelof := (self,n,m) >> Maximum(Filtered(self.mem_per_pt * Filtered(n*DivisorsInt(m), e-> e<self.max_threads),
            f -> f < When(n >= self.max_kernel, self.mem/2, self.mem)))/(self.mem_per_pt*n),

        applicable := nt -> nt.hasTags() and _isSIMTTag(nt.firstTag()) and IsParVec(nt.params) and nt.params[2] > 1,
        children := (self, nt) >> let(n := Rows(nt.params[1]), m:= nt.params[2], peelof := self._peelof(n,m), remainder := m/peelof,
            [[ TCompose([
                TTensorI(
                    TCompose([TTensorI(nt.params[1], peelof, APar, APar), TL(Rows(nt.params[1]) * peelof, peelof)]), remainder, APar, APar),
                TL(Rows(nt)/peelof, (Rows(nt)/peelof)/n, 1, peelof),
                ]).withTags(nt.getTags()) ]]),
        apply := (nt, c, cnt) -> c[1]
    )

));

max_stars_repo_name: mikefranusich/spiral-package-fftx

max_issues_repo_name: mikefranusich/spiral-package-fftx

max_forks_repo_name: mikefranusich/spiral-package-fftx


text:
#############################################################################
##
#W  extend.gi                                                   Karel Dekimpe
#W                                                               Bettina Eick
##

#############################################################################
##
#F FiniteSubgroupsAreCyclic( G )
##
FiniteSubgroupsAreCyclic := function( G )
    local fin, rep, flg;
    fin := FiniteSubgroupClasses( G );
    rep := List( fin, Representative );
    flg := ForAll( rep, IsCyclic );
    return flg;
end;

#############################################################################
##
#F AllActionsHolonomy( G )
##
## Computes all actions of G on Z which can lead to a torsion-free
## extension centralizing Fit(G). Thus all finite subgroups of G and Fit(G)
## centralizes Z and the kernel of the action has index at most 2.
##
AllActionsHolonomy := function( G )
    local mats, acts, fins, reps, gens, U, hom, F, low, H, i;

    # first we add the trivial action
    mats := List( Igs(G), x -> IdentityMat(1));
    acts := [mats];

    # determine the subgroup known to centralize Z
    U := FittingSubgroup(G);

    # determine all subgroups of index 2 in G/U
    hom := NaturalHomomorphismByNormalSubgroup( G, U );
    F   := Image( hom );
    low := LowIndexSubgroupClasses( F, 2 );
    low := List( low, Representative );
    low := List( low, x -> PreImage( hom, x ) );

    # for each subgroup we get one action
    for H in low do
        mats := List( Igs(G), x -> IdentityMat(1));
        for i in [1..Length(Igs(G))] do
            if not Igs(G)[i] in H then
                mats[i] := - mats[i];
            fi;
        od;
        Add( acts, mats );
    od;
    return acts;
end;

#############################################################################
##
#F AllActionsForTorsionFreeExtension( G )
##
## Computes all actions of G on Z which can lead to a torsion-free
## extension centralizing Fit(G). Thus all finite subgroups of G and Fit(G)
## centralizes Z and the kernel of the action has index at most 2.
##
AllActionsForTorsionFreeExtension := function( G )
    local mats, acts, fins, reps, gens, U, hom, F, low, H, i;

    # first we add the trivial action
    mats := List( Igs(G), x -> IdentityMat(1));
    acts := [mats];

    # determine the subgroup known to centralize Z
    fins := FiniteSubgroupClasses( G );
    reps := List( fins, Representative );
    gens := Flat( List( reps, Igs ) );
    Append( gens, GeneratorsOfGroup( FittingSubgroup(G) ) );
    U := NormalClosure( G, Subgroup( G, gens ) );
    if Index( G, U ) = 1 then return acts; fi;

    # determine all subgroups of index 2 in G/U
    hom := NaturalHomomorphismByNormalSubgroup( G, U );
    F   := Image( hom );
    low := LowIndexSubgroupClasses( F, 2 );
    low := List( low, Representative );
    low := List( low, x -> PreImage( hom, x ) );

    # for each subgroup we get one action
    for H in low do
        mats := List( Igs(G), x -> IdentityMat(1));
        for i in [1..Length(Igs(G))] do
            if not Igs(G)[i] in H then
                mats[i] := - mats[i];
            fi;
        od;
        Add( acts, mats );
    od;
    return acts;
end;

#############################################################################
##
#F ExpandedTail( tail, d )
##
ExpandedTail := function( t, d )
    local r, i;
    r := List( [1..d], x -> 0 );
    for i in [1..Length(t)] do
        if IsBound( t[i] ) then r[i] := t[i][1][1]; fi;
    od;
    return r;
end;

#############################################################################
##
#F HasTorsionFreeExtension( G, CR )
##
## Check if G has a torsion-free extension. We assume that all finite
## subgroups of G are cyclic (otherwise, no torsion-free extension exists).
## CR ist a record describing the second cohomology group of G by Z.
##
HasTorsionFreeExtension := function( G, CR )
    local fins, reps, repp, cc, d, subs, U, g, e, p, w, v, i, t, r, pr,
          f, V, spcs, sub;

    # compute the finite subgroups of prime order up to conjugacy
    fins := FiniteSubgroupClasses( G );
    reps := List( fins, Representative );
    repp := Filtered( reps, x -> IsPrime( Size(x) ) );
    if Length( repp ) = 0 then return true; fi;

    # extract second cohomology
    cc := CR.twocohom;

    # for each finite subgroup of prime order compute the corresponding
    # relation on cc.gcc
    subs := [];
    for U in repp do
        g := GeneratorsOfGroup( U )[1];
        e := Exponents( g );
        p := Order( g );
        w := rec( word := e, tail := [] );
        v := ShallowCopy( w );
        for i in [2..p] do
            v := CollectedTwoCR( CR, v, w );
        od;
        r := ExpandedTail( v.tail, Length( CR. enumrels ) );
        sub := NullspaceIntMod( cc.factor.prei, r, p );
        Add( subs, sub );
    od;

    e := Lcm( List( repp, Size ) );
    d := Length( cc.factor.prei );
    return SizeOfUnionMod( subs, e ) < e^d;
end;

#############################################################################
##
#F HasMinimalCentreExtension( G, CR )
##
## Check if there is an extension of G that has minimal centre; that is,
## Z(Fit(E)) = Z. This is only implemented for groups G with abelian
## Fitting subgroup.
##
HasMinimalCentreExtension := function( G, CR )
    local cc, d, t, F, f, g, l, c, i, j, v, w, a, b;

    # get cohomology and generic elements
    cc := CR.twocohom;
    d  := Length( CR.enumrels );
    t  := List( [1..Length(cc.gcc)], x -> Indeterminate( Integers, x ) );
    t  := t * cc.gcc;

    # get Fitting subgroup generators
    F := FittingSubgroup( G );

    # the abelian case is more effective - consider this first
    if IsAbelian( F ) then
        f := GeneratorsOfGroup( F );
        l := Length( f );

        # find tails of commutators
        c := NullMat( l, l );
        for i in [1..l] do
            for j in [i+1..l] do
                v := rec( word := Exponents( f[i]^-1 ), tail := [] );
                w := rec( word := Exponents( f[j]^-1 ), tail := [] );
                v := CollectedTwoCR( CR, v, w );
                w := rec( word := Exponents( f[i] ), tail := [] );
                v := CollectedTwoCR( CR, v, w );
                w := rec( word := Exponents( f[j] ), tail := [] );
                v := CollectedTwoCR( CR, v, w );
                b := ExpandedTail( v.tail, d );
                a := b * t;
                c[i][j] := a;
                c[j][i] := -a;
            od;
        od;

        # get determinant
        c := GenericDeterminantMat( c );
        return Length( ExtRepPolynomialRatFun(c) ) <> 0;
    fi;

    # the non-abelian case
    f := MinimalGeneratingSet( F );
    g := GeneratorsOfGroup( Centre( F ) );
    if Length(f) < Length( g ) then return false; fi;

    # find tails of commutators
    c := NullMat( Length(f), Length(g) );
    for i in [1..Length(f)] do
        for j in [1..Length(g)] do
            v := rec( word := Exponents( f[i]^-1 ), tail := [] );
            w := rec( word := Exponents( g[j]^-1 ), tail := [] );
            v := CollectedTwoCR( CR, v, w );
            w := rec( word := Exponents( f[i] ), tail := [] );
            v := CollectedTwoCR( CR, v, w );
            w := rec( word := Exponents( g[j] ), tail := [] );
            v := CollectedTwoCR( CR, v, w );
            b := ExpandedTail( v.tail, d );
            c[i][j] := b * t;
        od;
    od;

    # get determinant of rxr submats
    Error("not yet implemented");
    c := Determinant( c );
    c := ExtRepPolynomialRatFun( c );
    c := Filtered( c, x -> IsInt(x) );
    return ForAny( c, x -> x <> 0 );
end;

#############################################################################
##
#F HasExtensionOfType( G, torfree, mincent )
##
## Let G be a pcp group. This function checks if G has a torsion-free
## extension or an extension with Z(Fit(H)) = Z$ with the free abelian
## module Z.
##
InstallGlobalFunction( HasExtensionOfType, function( G, torfree, mincent )
    local mats, CR, found, i;

    # if both flags are false, then there is nothing to do
    if not torfree and not mincent then return true; fi;

    # first check if G is torsion-free or has non-cyclic finite subgrps
    Print("    check that all finite subgroups are cyclic\n");
    if not FiniteSubgroupsAreCyclic( G ) then return false; fi;

    # if G is torsion-free, then all extensions of G are torsion-free
    if Length( FiniteSubgroupClasses(G) ) = 1 then torfree := false; fi;
    if not torfree and not mincent then return true; fi;

    # now loop over actions - the trivial action first
    Print("    consider trivial action \n");
    mats := List( Igs(G), x -> IdentityMat( 1 ) );
    CR := CRRecordByMats( G, mats );
    CR.twocohom := TwoCohomologyCR( CR );
    if mincent then
        found := HasMinimalCentreExtension( G, CR );
        #if found <> HasMinimalCentreExtensionByCRRec2( G, CR ) then
        #    Error("something wrong with min cent");
        #fi;
    fi;
    if (not mincent) or (found and torfree) then
        found := HasTorsionFreeExtension( G, CR );
    fi;
    if found then return true; fi;

    # now consider the remaining actions
    mats := AllActionsForTorsionFreeExtension( G );
    for i in [2..Length(mats)] do
        Print("    consider non-trivial action number ",i-1,"\n");
        CR := CRRecordByMats( G, mats[i] );
        CR.twocohom := TwoCohomologyCR( CR );
        if mincent then
            found := HasMinimalCentreExtension( G, CR );
            #if found <> HasMinimalCentreExtensionByCRRec2( G, CR ) then
            #    Error("something wrong with min cent");
            #fi;
        fi;
        if (not mincent) or (found and torfree) then
            found := HasTorsionFreeExtension( G, CR );
        fi;
        if found then return true; fi;
    od;

    # we have not found a suitable extension
    return false;
end );

max_stars_repo_name: alex-konovalov/aclib

max_issues_repo_name: alex-konovalov/aclib

max_forks_repo_name: alex-konovalov/aclib


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

Class(InertTranspose, BaseContainer, SumsBase, rec(
    dims := self >> Reversed(self.child(1).dims()),
    toAMat := self >> AMatMat(TransposedMat(MatSPL(self.child(1))))
));

Class(InertConjTranspose, BaseContainer, SumsBase, rec(
    dims := self >> Reversed(self.child(1).dims()),
    toAMat := self >> AMatMat(Global.Conjugate(TransposedMat(MatSPL(self.child(1)))))
));

InertConjTranspose.conjTranspose := self >> self.child(1);


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

benchAVX := function()
    if LocalConfig.cpuinfo.SIMD().hasAVX() then
        return rec(
            4x64f := rec(
                wht := _defaultSizes(s->doSimdWht(s, AVX_4x64f, rec(verify := true, oddSizes := false, svct := true)), List([4..8], i->2^i)),
                1d := rec(
                    dft_ic := rec(
                        small := rec(
                            real := _defaultSizes(s->doSimdDft(s, AVX_4x64f, rec(verify := true, globalUnrolling := 10000, tsplBluestein:=true,
                              PRDFT:=true, URDFT:= true, CT := true, minCost := true, tsplPFA := true, tsplCT := true, oddSizes := true, PD := true, PFA := true,
                              RealRader := true, Rader := true, realVect := true, cplxVect := false, stdTTensor := false, interleavedComplex := true)), [2..64]),
                            cmplx := _defaultSizes(s->doSimdDft(s, AVX_4x64f, rec(verify := true, globalUnrolling := 10000, tsplBluestein:=true,
                              PRDFT:=false, URDFT:= true, CT := true, minCost := true, tsplPFA := true, tsplCT := true, oddSizes := true, PD := true, PFA := true,
                              RealRader := true, Rader := true, realVect := false, cplxVect := true, stdTTensor := false, interleavedComplex := true)), [2..64])
                        ),
                        medium := _defaultSizes(s->doSimdDft(s, AVX_4x64f, rec(globalUnrolling := 128, tsplRader:=false, tsplBluestein:=false, tsplPFA:=false, oddSizes:=false,
                                interleavedComplex := true, cplxVect := false, realVect := true)),
                              List([4..16], i->2^i)),
                        medium_cx := _defaultSizes(s->doSimdDft(s, AVX_4x64f, rec(
                                globalUnrolling := 128,
                                tsplRader       := false,
                                tsplBluestein   := false,
                                tsplPFA         := false,
                                URDFT           := true,
                                PRDFT           := false,
                                oddSizes        := false,
                                cplxVect        := true,
                                realVect        := false,
                                svct            := false,
                                splitL          := true,
                                interleavedComplex := true)),
                              List([2..16], i->2^i))
                    ),

                    dft_sc := rec(
                        small :=  _defaultSizes(s -> doSimdDft(s, AVX_4x64f, rec(
				       verify     := true,
				       stdTTensor := false,
				       svct       := true,
				       flipIxA    := true,
				       realVect   := true,
				       oddSizes   := true,
				       interleavedComplex := false,
				       globalUnrolling := 10000
				  )), List([1 .. 8], i -> 2^i)), # tsplRader:=false, tsplBluestein:=false, tsplPFA:=false,

                        medium := _defaultSizes(s->doSimdDft(s, AVX_4x64f, rec(
				       verify        := true,
				       tsplRader     := false,
				       tsplBluestein := false,
				       tsplPFA       := false,
				       oddSizes      := false,
				       globalUnrolling := 10000,
				       interleavedComplex := false
				  )), 16 * [1 .. 8]),
                    ),

                    trdft := _defaultSizes(s -> doSimdSymDFT(TRDFT, s, AVX_4x64f, rec(
				  verify     := true,
				  PRDFT      := false,
				  URDFT      := true,
				  cplxVect   := true,
                                  stdTTensor := false,
				  realVect   := true,
				  svct       := true, # <-- added from dft_sc
				  #flipIxA    := true, #<-- added from dft_sc, do we need this?
				  useDeref   := true,
                                  oddSizes   := true, # need this for TRDFT(16)
				  propagateNth := true,
                                  globalUnrolling :=10000,
                                  interleavedComplex := true,
			      )), 8 * [1 .. 32]),
                )
            ),
            8x32f := rec(
                wht := _defaultSizes(s->doSimdWht(s, AVX_8x32f, rec(verify := true, oddSizes := false, svct := true)), List([6..8], i->2^i)),
                1d := rec(
                    dft_ic := rec(
                        small := rec(
                            real := _defaultSizes(s -> doSimdDft(s, AVX_8x32f, rec(
					  verify     := true,
					  PRDFT      := true,
					  URDFT      := true,
					  CT         := true,
					  minCost    := true,
					  tsplPFA    := true,
					  tsplCT     := true,
					  oddSizes   := true,
					  PD         := true,
					  PFA        := true,
					  RealRader  := true,
					  Rader      := true,
					  realVect   := true,
					  cplxVect   := false,
					  stdTTensor := false,
					  tsplBluestein := true,
					  globalUnrolling := 10000,
					  interleavedComplex := true)), [2..64]),

                            cmplx := _defaultSizes(s->doSimdDft(s, AVX_8x32f, rec(
					  verify := true, globalUnrolling := 10000, tsplBluestein:=true,
					  PRDFT:=false, URDFT:= true, CT := true, minCost := true, tsplPFA := true,
					  tsplCT := true, oddSizes := true, PD := true, PFA := true,
					  RealRader := true, Rader := true, realVect := false, cplxVect := true,
					  stdTTensor := false, interleavedComplex := true)), [2..64])
                        ),

                        medium := _defaultSizes(s->doSimdDft(s, AVX_8x32f, rec(
				      tsplRader     := false,
				      tsplBluestein := false,
				      tsplPFA := false,
				      oddSizes := false,
				      cplxVect := false,
				      realVect := true,
				      globalUnrolling := 128,
                                      interleavedComplex := true,
				  )), List([6..16], i->2^i)),

                        medium_cx := _defaultSizes(s->doSimdDft(s, AVX_8x32f, rec(
                                    globalUnrolling := 128,
                                    tsplRader       := false,
                                    tsplBluestein   := false,
                                    tsplPFA         := false,
                                    PRDFT           := false,
                                    URDFT           := true,
                                    CT              := false,
                                    PD              := false,
                                    oddSizes        := false,
                                    cplxVect        := true,
                                    realVect        := false,
                                    svct            := true,
                                    flipIxA         := true,
                                    interleavedComplex := true
				)), List([4..16], i->2^i))
                    ),

                    dft_sc := rec(
                        small :=  _defaultSizes(s -> doSimdDft(s, AVX_8x32f, rec(
				       verify     := true,
				       stdTTensor := false,
				       svct       := true,
				       flipIxA    := true,
				       realVect   := true,
				       oddSizes   := true,
				       interleavedComplex := false,
				       globalUnrolling := 10000
				  )), List([1 .. 8], i -> 2^i)),

                        medium := _defaultSizes(s->doSimdDft(s, AVX_8x32f, rec(
				       verify        := true,
				       tsplRader     := false,
				       tsplBluestein := false,
				       tsplPFA       := false,
				       oddSizes      := false,
				       globalUnrolling := 10000,
				       interleavedComplex := false
				  )), 64*[1..4]),
                    ),

                    trdft := _defaultSizes(s -> doSimdSymDFT(TRDFT, s, AVX_8x32f, rec(
				  verify     := true,
				  PRDFT      := false,
				  URDFT      := true,
				  cplxVect   := true,
                                  stdTTensor := false,
				  realVect   := true,
				  svct       := true, # <-- added from dft_sc
				  #flipIxA    := true, #<-- added from dft_sc, do we need this?
				  useDeref   := true,
                                  oddSizes   := true, # need this for TRDFT(16)
				  propagateNth := true,
                                  globalUnrolling :=10000,
                                  interleavedComplex := true,
			      )), 16 * [1 .. 32]),
                )
            )
        );
    else
        return false;
    fi;
end;


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

Class(LegacySumsGen, HierarchicalVisitor, rec(

   __call__ := meth(arg)

        Constraint(Length(arg) >= 2);

        #BWD: I believe LegacySumsGen.visit only expects Formula objects
        if not PatternMatch(arg[2], Formula, empty_cx()) then
            arg[2] := Formula(arg[2]);
        fi;

        return ApplyFunc(arg[1].visit, arg{[2..Length(arg)]});

   end,

	Formula := meth(self, o, opts)
		local spl, sums;
		spl := _children(o)[1];

    	if IsBound(spl.memo_sums) then
    		sums := spl.memo_sums;
    	else
    	sums := spl.sums();
    		spl._sums := sums;
    	fi;
    	sums.root := spl;
    	#sums.sums := self >> self;
    	return sums;
	end,
));

Class(DefaultSumsGen, HierarchicalVisitor, rec(

    __call__ := meth(arg)
        local res;
        res := ApplyFunc(arg[1].visit, arg{[2..Length(arg)]});
        trace_log.addExpansion(ObjId(arg[2]), arg[2], res, var);
        return res;
   end,

    _directSumSums := (self, obj, roverlap, coverlap, opts) >>
        let(cspans  :=  BaseOverlap.spans(roverlap, List(obj.children(), Cols)),
            rspans  :=  BaseOverlap.spans(coverlap, List(obj.children(), Rows)),
            nblocks  :=  obj.numChildren(),
            C  :=  Cols(obj),
            R  :=  Rows(obj),
            OP  :=  When(coverlap > 0, SUMAcc, SUM),
            OP( Map([1..nblocks],
                 i -> Compose( #W(N, bksize_cols, fTensor(fBase(nblocks, ind), fId(bksize)))
                      WID(R, 1 + rspans[i][2] - rspans[i][1], rspans[i][1]-1),
                      self(obj.child(i), opts),
                      RID(C, 1 + cspans[i][2] - cspans[i][1], cspans[i][1]-1) )))),

    Formula     := (self, o, opts) >> self(_children(o)[1], opts),
    NonTerminal := (self, o, opts) >> o,
    PermClass   := (self, o, opts) >> Prm(o),

    Rot := (self, o, opts) >> self(ExpandRotationsSPL(o), opts),
    Sym := (self, o, opts) >> self(o.getObj(), opts),

    O  := (self, o, opts) >> o,
    L  := (self, o, opts) >> Prm(o),
    F  := (self, o, opts) >> self(o.getObj(), opts),
    I  := (self, o, opts) >> Prm(fId(o.params[1])),

    RC := (self, o, opts) >> CopyFields(o, rec(_children := [self(o.child(1), opts)])),

    Perm := (self, o, opts) >> let(
	N := o.size,
	plist_r := Permuted([1..N], o.element^-1) - 1, # we want 0-based perms
	#plist_w := Permuted([1..N], o.element) - 1,    # we want 0-based perms
	Gath(FList(N, plist_r).setRange(N))
    ),

    Diag := (self, o, opts) >> o,
    RCDiag := (self, o, opts) >> o,
    Mat := (self, o, opts) >> Blk(o.element),
    Sparse := (self, o, opts) >> Blk(MatSparseSPL(o)),

    B := (self, o, opts) >> Prm(fB(o.size, o.l)),
    B2 := (self, o, opts) >> Prm(o),
    Scale := (self, o, opts) >> Inherit(o, rec(_children := [self(o.child(1), opts)])),

    Conjugate := (self, o, opts) >> let(
	O := o.child(1), perm := o.child(2), Tperm := perm.transpose(),
	self(Compose(Tperm, O, perm), opts)
    ),

    Compose := (self, o, opts) >> Cond(o.isPermutation(),
	Prm(ApplyFunc(fCompose, List(Reversed(o.children()), c -> self(c, opts).func))),
	Compose(Map(o.children(), c -> self(c, opts)))
    ),

    Cross := (self, o, opts) >> Cond(o.isPermutation(),
	Prm(ApplyFunc(fCross, List(o.children(), c -> self(c, opts).func))),
	Cross(Map(o.children(), c -> self(c, opts)))
    ),

    ICompose := (self, o, opts) >> CopyFields(o, rec(
        _children := List(o.children(), c -> self(c, opts))
    )),

    DirectSum := (self, o, opts) >> Cond(o.isPermutation(),
	Prm(ApplyFunc(fDirsum, List(o.children(), c->self(c, opts).func))),
        self._directSumSums(o, 0, 0, opts)
    ),

    DelayedDirectSum := (self, o, opts) >> CopyFields(o, rec(
        _children := List(o.children(), c -> self(c, opts))
    )),

    RowDirectSum := (self, o, opts) >> self._directSumSums(o, o.overlap, 0, opts),
    ColDirectSum := (self, o, opts) >> self._directSumSums(o, 0, o.overlap, opts),

    RowTensor := (self, o, opts) >> let(
        A := o.child(1),
	i := Ind(o.isize),
	ISum(i, i.range,
	    Scat(fTensor(fBase(o.isize, i), fId(Rows(A)))) *
	    SumsSPL(A, opts) *
	    Gath(fAdd(Cols(o), Cols(A), i * (Cols(A) - o.overlap))))
    ),

    ColTensor := (self, o, opts) >> let(
        A := o.child(1),
	i := Ind(o.isize),
	ISumAcc(i, i.range,
	    Scat(fAdd(Rows(o), Rows(A), i * (Rows(A) - o.overlap))) *
	    SumsSPL(A, opts) *
	    Gath(fTensor(fBase(o.isize, i), fId(Cols(A)))))
    ),

    Tensor := meth(self, o, opts)
    	local ch, col_prods, row_prods, col_prods_rev, row_prods_rev, i, j1, j2, j1b, j2b, bkcols, bkrows, prod, term;

    	if o.isPermutation() then
	    return Prm(ApplyFunc(o.fTensor, List(o.children(), c->self(c, opts).func)));
	elif ForAll(o.children(), x->ObjId(x)=Diag) then
	    return Diag(ApplyFunc(diagTensor, List(o.children(), c->c.element)));
	fi;

	ch := o.children();
	col_prods := ScanL(ch, (x,y)->x*Cols(y), 1);
	col_prods_rev := Drop(ScanR(ch, (x,y)->x*Cols(y), 1), 1);

	row_prods := ScanL(ch, (x,y)->x*Rows(y), 1);
	row_prods_rev := Drop(ScanR(ch, (x,y)->x*Rows(y), 1), 1);

	prod := [];
	for i in [1..Length(ch)] do
	    if not IsIdentitySPL(ch[i]) then
		bkcols := Cols(ch[i]);
		bkrows := Rows(ch[i]);
		j1 := Ind(col_prods[i]);
		j2 := Ind(row_prods_rev[i]);
		j1b := When(j1.range = 1, 0, fBase(j1));
		j2b := When(j2.range = 1, 0, fBase(j2));

		term := Scat(o.scatTensor(i)(Filtered([j1b, fId(bkrows), j2b], x->x<>0))) *
		        self(ch[i], opts) *
			Gath(o.gathTensor(i)(Filtered([j1b, fId(bkcols), j2b], x->x<>0)));

		if j2.range <> 1 then
		    term := ISum(j2, j2.range, term);
		fi;

		if j1.range <> 1 then
		    term := ISum(j1, j1.range, term);
		fi;

		Add(prod, term);
	    fi;
	od;
	return Compose(prod);
    end,

    IDirSum := (self, o, opts) >> let(
	bkcols := Cols(o.child(1)),
	bkrows := Rows(o.child(1)),
	nblocks := o.domain,
	cols := Cols(o),
	rows := Rows(o),

	ISum(o.var, o.domain,
	    Compose(Scat(fTensor(fBase(nblocks, o.var), fId(bkrows))),
		    self(o.child(1), opts),
		    Gath(fTensor(fBase(nblocks, o.var), fId(bkcols)))))
    ),

    HStack := (self, o, opts) >> let(
	ch   := List(o.children(), e -> self(e, opts)),
	l    := Length(ch),
	it0  := Scat(fId(Rows(o))) * ch[1] * Gath(H(Cols(o), Cols(ch[1]), 0, 1)),
	ofs  := ScanL(ch, (p,x)->p+Cols(x), 0),
	its  := List([2..l], i ->
	           ScatAcc(fId(Rows(o))) * ch[i] * Gath(H(Cols(o), Cols(ch[i]), ofs[i], 1))),
	SUM( [it0] :: its )
    ),

    VStack := (self, o, opts) >> let(
	ch   := List(o.children(), e -> self(e, opts)),
	l    := Length(ch),
	ofs  := ScanL(ch, (p,x)->p+Rows(x), 0),
	SUM(List([1..l], i ->
	        Scat(H(Rows(o), Rows(ch[i]), ofs[i], 1)) * ch[i] * Gath(fId(Cols(o)))))
    ),

    HStack1 := (self, o, opts) >> let(
    	ch   := List(o.children(), e -> self(e, opts)),
    	l    := Length(ch),
    	it0  := Scat(fId(Rows(o))) * ch[1] * Gath(H(Cols(o), Cols(ch[1]), 0, 1)),
    	ofs  := ScanL(ch, (p,x)->p+Cols(x), 0),
    	its  := List([2..l], i ->
    	           Scat(fId(Rows(o))) * ch[i] * Gath(H(Cols(o), Cols(ch[i]), ofs[i], 1))),
    	SUM( [it0] :: its )
    ),

    VStack := (self, o, opts) >> let(
	ch   := List(o.children(), e -> self(e, opts)),
	l    := Length(ch),
	ofs  := ScanL(ch, (p,x)->p+Rows(x), 0),
	SUM(List([1..l], i ->
	        Scat(H(Rows(o), Rows(ch[i]), ofs[i], 1)) * ch[i] * Gath(fId(Cols(o)))))
    ),

    IterHStack := (self, o, opts) >> let(
	bkcols := Cols(o.child(1)),
	bkrows := Rows(o.child(1)),
	nblocks := o.domain,
	cols := Cols(o), rows := Rows(o),
	j := Ind(nblocks-1),
	ch := self(o.child(1), opts),

	SUM(
	    Scat(fId(bkrows)) *
	    SubstVars(Copy(ch), rec((o.var.id) := V(0))) *
	    Gath(fTensor(fBase(nblocks, 0), fId(bkcols))),

	    ISum(j, j.range,
		ScatAcc(fId(bkrows)) *
		SubstVars(Copy(ch), rec((o.var.id) := j+1)) *
		Gath(fTensor(fBase(nblocks, j+1), fId(bkcols)))
	    )
	)
    ),

    IterHStack1 := (self, o, opts) >> let(
    	bkcols := Cols(o.child(1)),
    	bkrows := Rows(o.child(1)),
    	nblocks := o.domain,
    	cols := Cols(o), rows := Rows(o),
    	j := Ind(nblocks),
    	ch := self(o.child(1), opts),

	    ISum(j, j.range,
		Scat(fId(bkrows)) *
		SubstVars(Copy(ch), rec((o.var.id) := j)) *
		Gath(fTensor(fBase(nblocks, j), fId(bkcols)))
    	)
    ),

    DPWrapper := (self, o, opts) >> self(o.child(1), opts),

    IterVStack := (self, o, opts) >> let(
	bkcols := Cols(o.child(1)),
	bkrows := Rows(o.child(1)),
	nblocks := o.domain,
	cols := Cols(o), rows := Rows(o),
	ISum(o.var, o.domain,
	    Scat(fTensor(fBase(nblocks, o.var), fId(bkrows))) *
	    self(o.child(1), opts) *
	    Gath(fId(bkcols)))
    ),

    ISum := (self, o, opts) >> CopyFields(o, rec(_children := [self(o._children[1], opts)])),

    RecursStep := ~.ISum,
    VPerm      := (self, o, opts) >> o,
    BlockVPerm := (self, o, opts) >> o,
    ExpDiag    := (self, o, opts) >> o,

    VTensor    := (self, o, opts) >> ObjId(o)(self(o.child(1), opts), o.vlen),

    SumsBase := meth(self, o, opts)
       local children, i, res;
       res := Copy(o);
       ## MRT, aug08
       # IsSPL -> IsSPL and IsFunction
       # joint SPL/Sigma objects like fId were improperly triggering the
       # SPL condition, resulting in illegal rewrites like:
       #   Scat(fId(2)) -> Scat(Prm(fId(2))

       # YSV: This is still buggy, because VTensor(L) --> VTensor(L)
       children := Map(res.rChildren(), c -> Cond(IsSPL(c) and not IsFunction(c), self(c, opts), c));
       for i in [1..Length(children)] do
           res.rSetChild(i, children[i]);
       od;
       return res;
    end,

     Split := (self, o, opts) >> o,
     Glue := (self, o, opts) >> o,
     VIxL := (self, o, opts) >> o,

));


text:
Logical := function(a, b)
    return [ a or b, a and b, not a ];
end;

Logical(true, true);
# [ true, true, false ]

Logical(true, false);
# [ true, false, false ]

Logical(false, true);
# [ true, false, true ]

Logical(false, false);
# [ false, false, true ]


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

# ==========================================================================
# RowDirectSum
# ==========================================================================
Class(RowDirectSum, BaseOverlap, rec(
    #-----------------------------------------------------------------------
    abbrevs := [
        function(arg)
          arg:=Flat(arg);
          return [arg[1], arg{[2..Length(arg)]}];
        end ],
    #-----------------------------------------------------------------------
    new := meth(self, overlap, spls)
       return self._new(0, overlap, spls);
    end,
    #-----------------------------------------------------------------------
    _dims := meth(self)
        local ovdim;
    ovdim := self.ovDim(self.overlap,
                        List(self._children, t -> t.dimensions[2]));
    return [ Sum(self._children, t -> t.dimensions[1]),
             ovdim[3] - ovdim[2] + 1 # max - min + 1
    ];
    end,
    dims := self >> When(IsBound(self._setDims), self._setDims, let(d := Try(self._dims()),
	    When(d[1], d[2], [errExp(TInt), errExp(TInt)]))),
    #-----------------------------------------------------------------------
    toAMat := meth(self)
        return
      DirectSumAMat(List(self._children, AMatSPL)) *
      AMatSPL(
          Sparse(
          self._rowoverlap(
              List(self._children, t -> t.dimensions[2]),
              self.overlap
          )
          )
      );
    end,
    #-----------------------------------------------------------------------
    arithmeticCost := meth(self, costMul, costAddMul)
        return Sum(List(self.children(), x -> x.arithmeticCost(costMul, costAddMul)));
    end
));


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

#Concatenation(Conf("spiral_dir"), Conf("path_sep"), "libgen", Conf("path_sep"),

#F InitLibgen(<opts>)
#F  opts must have a .libgen field.
#F  Example: opts := InitLibgen(LibgenDefaults);
#F
InitLibgen := function(opts)
    local base_opts, bench;
    opts := Copy(opts);

    base_opts := CopyFields(opts, rec(
            globalUnrolling := opts.libgen.basesUnrolling,
            breakdownRules  := opts.libgen.basesBreakdowns,
            hashFile        := opts.libgen.basesHashFile,
            benchTransforms := opts.libgen.bases));

    bench:= CopyFields(DPBench(rec(bases := base_opts), rec(verbosity:=0)),
                       rec(generateCfiles := false, generateSums := true, matrixVerify := false));
    bench.resumeAll();

    opts.libgen.baseBench := bench;
    opts.baseHashes := Concatenation(opts.baseHashes, [ CreateRecursBaseHash(bench.exp.bases.hashTable) ]);
    return opts;
end;

_bases := function(max)
   local t, facs, sizesTwoPower, sizesGeneral;
   sizesTwoPower := List([1..Log2Int(max)], x->2^x);
   facs := [2, 3, 5, 7, 11];
   sizesGeneral := Filtered([2 .. max], i -> ForAll(FactorsInt(i), f -> f in (facs)));
   t := [];
   Append(t, List(sizesTwoPower, PRDFT3));
   Append(t, List(sizesTwoPower, n -> PRDFT3(n, -1).transpose()));
   Append(t, List(sizesTwoPower, n -> DFT(n, -1 mod n)));
#  Append(t, List(sizesTwoPower, URDFT));
   Append(t, List(sizesGeneral, PRDFT));
   Append(t, List(sizesGeneral, DFT));
   return t;
end;

LibgenDefaultsMixin := rec(
    codegen := RecCodegen,
    verbosity := 1,
    hashTable := HashTableDP(),

    formulaStrategies := rec(
	sigmaSpl    := [ StandardSumsRules, HfuncSumsRules ],
	preRC       := [],
	rc          := [ StandardSumsRules, HfuncSumsRules ],
        postProcess := [
	    (s, opts) -> BlockSums(opts.globalUnrolling, s),
	    (s, opts) -> Process_fPrecompute(s, opts),
            RecursStepTerm
        ]
    ),

    libgen := rec(
	codeletTab := CreateCodeletHashTable(),
	terminateStrategy := [ HfuncSumsRules ], # used in SumsCodelet

	bases := _bases(34),
	basesUnrolling := 2^16,
        basesHashFile := let(p:=Conf("path_sep"),
            Concat(Conf("spiral_dir"), p, "spiral", p, "libgen", p, "realbases.hash")),
	basesBreakdowns := rec(
	    DFT    := [DFT_Base, DFT_CT, DFT_GoodThomas, DFT_PD, DFT_Rader],
	    PRDFT  := [PRDFT1_Base2, PRDFT1_CT, PRDFT1_PF, PRDFT_Rader, PRDFT_PD],
	    PRDFT3 := [PRDFT3_Base1, PRDFT3_Base2, PRDFT3_CT, PRDFT3_OddToPRDFT1])
    )
);

LibgenDefaults := CopyFields(SpiralDefaults, LibgenDefaultsMixin, rec(
    globalUnrolling := 66,
    breakdownRules := rec(
	DFT    := [DFT_Base, DFT_CT],
	PRDFT  := [PRDFT1_Base2, PRDFT1_CT],
	PRDFT3 := [PRDFT3_Base2, PRDFT3_CT],
	PDCT4  := [PDCT4_Base2, PDCT4_CT],
	PDST4  := [PDST4_Base2, PDST4_CT],
	InfoNt := [Info_Base]),
    libgen := Copy(LibgenDefaultsMixin.libgen)
));

CplxLibgenDefaults := CopyFields(CplxSpiralDefaults, LibgenDefaultsMixin, rec(
    globalUnrolling := 34,
    breakdownRules := rec(
	DFT    := [DFT_Base, DFT_CT],
	PRDFT  := [PRDFT1_Base2, PRDFT1_CT],
	PRDFT3 := [PRDFT3_Base2, PRDFT3_CT],
	InfoNt := [Info_Base]),
    libgen := CopyFields(Copy(LibgenDefaultsMixin.libgen), rec(
	bases := List([1..5], x->DFT(2^x)),
        basesHashFile := let(p:=Conf("path_sep"),
            Concat(Conf("spiral_dir"), p, "spiral", p, "libgen", p, "cplxbases.hash")),
	basesBreakdowns := rec(
	    DFT    := [DFT_Base, DFT_PRDFT],
	    PRDFT  := [PRDFT1_Base1, PRDFT1_Base2, PRDFT1_CT, PRDFT1_PF, PRDFT_PD, PRDFT_Rader],
	    PRDFT3 := [PRDFT3_Base1, PRDFT3_Base2, PRDFT3_CT, PRDFT3_OddToPRDFT1])
    ))
));


text:
IsEvenInt(n);
IsOddInt(n);


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

apply_nt := (NT, input) ->
    List(TransposedMat(MatSPL(NT) * TransposedMat([input]))[1], EvalScalar);

dft_a := [0, 0, 1/2, 1/2];
dft_b := [0, 1/2, 0, 1/2];

dft_algo := [ [DFT1, DFT1], # f1 -> f1, f1
              [DFT2, DFT1], # f2 -> f2, f1
	      [DFT1, DFT3], # f3 -> f1, f3
	      [DFT2, DFT3]  # f4 -> f2, f3
	    ];
twid := (type,N,m,i) -> List(
    Twid(N,m,1,dft_a[type], dft_b[type],i).lambda().tolist(), EvalScalar);

# F_N -> Tensor(F_m, I_n) * T(mn, n) * Tensor(I_m, F_n) * L(mn, m)
#
[n, inputs, dfts, inputs2, twids, tinputs2, outputs, out, tdfts] := [0,0,0,0,0,0,0,0,0];

dft_projection := function(type, N, m, input)
    n := N/m;
    inputs := List([0..m-1], x->part(input,m,x));
    dfts    := List(inputs, x->apply_nt(dft_algo[type][2](Length(x)), x));
    inputs2 := TransposedMat(dfts);
    twids := List([0..n-1], i -> twid(type,N,m,i));
    tinputs2 := List([1..n], i -> listmul(inputs2[i], twids[i]));
    outputs := List(tinputs2, x->apply_nt(dft_algo[type][1](Length(x)), x));
    return List([inputs, dfts, inputs2, tinputs2, outputs], x->List(x,csym));
end;

tdft_projection := function(type, N, m, input)
    n := N/m;
    inputs := List([0..n-1], x->part(input,n,x));
    dfts    := List(inputs, x->apply_nt(dft_algo[type][2](Length(x)), x));
    twids := List([0..n-1], i->twid(Cond(type=2,3,type=3,2,type),N,m,i));
    tdfts := List([1..n], i -> listmul(dfts[i], twids[i]));
    tinputs2 := TransposedMat(tdfts);
    outputs := List(tinputs2, x->apply_nt(dft_algo[type][1](Length(x)), x));
    return List([inputs, dfts, tdfts, tinputs2, outputs], x->List(x,csym));
end;

Declare(inp, p1, p2, imid, dft1, dft2);
dft1_radproj := function(N, input)
    local mid; #,imid, p1, p2, inp, out;
    mid := DFT_Rader.raderMid(N,1,PrimitiveRootMod(N));
    input := apply_nt(RR(N), input);
    p1 := [input[1]]; inp := Drop(input,1);
    dft1 := apply_nt(DFT(N-1,-1), inp);
    imid := apply_nt(mid, Concat(p1, dft1));
    p2 := [imid[1]]; imid := Drop(imid, 1);
    dft2 := apply_nt(DFT(N-1,-1), imid);
    out := apply_nt(RR(N).transpose(), Concat(p2, dft2));
    return List([inp, dft1, imid, dft2, out], csym);
end;
# left transform in Rader rule:
# jDFT1 -> IRDFT1 dirsum jIRDFT2
# jDFT1 -> (jDFT1   dirsum jDFT3  )^L (R^2 dirsum C^2 ... I^2)
# jDFT3 -> (jDFT1 D dirsum jDFT1 D)^L (R^2 dirsum C^2 ... R^2)
#
#
rcexp := l -> ConcatList(l,
    e -> let(ee := Complex(e), [ReComplex(ee), ImComplex(ee)]));

# DFT_mn = CRT(m,n,1,1)^T (OS(m,n) F_m x I_n) (I_m x OS(n,m) F_n) CRT(m,n,1,1)
pf_projection := function(N, m, input)
    local i, j;
    n := N/m;
    inputs := List([0..m-1], j -> List([0..n-1], i->input[1+ (j*n + i*m) mod N]));
    dfts    := List(inputs, x->apply_nt(OS(n,m)*DFT(n), x));
    inputs2 := TransposedMat(dfts);
    outputs := List(inputs2, x->apply_nt(OS(m,n)*DFT(m), x));
    out := [1..N];
    for j in [0..n-1] do
       for i in [0..m-1] do
           out[1+ (j*m + i*n) mod N] := outputs[1+j][1+i];
       od;
    od;
    return List([inputs, dfts, inputs2, outputs], x->List(x,csym));
end;


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

#F qHT( n ) - Hadamard Gate non-terminal
#F Definition: (2^n x 2^n)-matrix  that applies an n-point Hadamard Transform to the target qubits
#F Note:       qHT(1) denotes the matrix 1/sqrt(2)* [[1, 1], [1, -1]],
#F             qHT(_) is symmetric.
#F
#F qHT( n ) -> an n-qubit Hadamard transform
Class(qHT, TaggedNonTerminal, rec(
  abbrevs   := [ n -> Checked(IsPosInt(n), [n]) ],
  dims      := self >> let(size := 2^self.params[1], [size, size]),
  terminate := self >> Tensor(Replicate(self.params[1], qH())),
  isReal    := self >> true,
  groups    := self >> [Replicate(self.params[1], 1)],
  connected := self >> false,
  rChildren := self >> self.params,
  from_rChildren := (self, rch) >> self.__bases__[1](rch[1]),
  recursive_def := (self, arch) >> self.__bases__[1](self.params[1]),
  SmallRandom := () -> Random([2..5]),
  LargeRandom := () -> Random([6..15]),
  normalizedArithCost := self >> Error("ArithCost not implemented"),
  TType := T_Real(64)
));

##
#F qHT Breakdown Rules
##
NewRulesFor(qHT, rec(

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

    #F qHT_Base: qHT(1) = qH() SPL object
    #F Directly represent as an implementable gate
    qHT_Base := rec(
        info             := "qHT_(1) -> qH()",
        forTransposition := false,
        applicable       := (self, nt) >> nt.params[1]=1,
        apply            := (nt, c, cnt) -> qH(),
    )

));

max_stars_repo_name: spiral-software/spiral-package-quantum

max_issues_repo_name: spiral-software/spiral-package-quantum

max_forks_repo_name: spiral-software/spiral-package-quantum