text:
#############################################################################
##
#W  loewy.gi             GAP 4 package SingerAlg                Thomas Breuer
##
##  This file contains implementations of GAP functions for studying
##  the Loewy structure of Singer algebras A[q,z].
##

##  for those who worked with an earlier version of the package ...
LoewyLayersData:= function( arg )
    Print( "use LoewyStructureInfo not LoewyLayersData\n" );
    return CallFuncList( LoewyStructureInfoGAP, arg );
    end;

if IsPackageMarkedForLoading( "JuliaInterface", "" ) then

#############################################################################
##
#M  LoewyStructureInfoJulia( <A> )
##
##  Call a Julia function, convert its result to a record,
##  but keep the record components in Julia.
##
InstallMethod( LoewyStructureInfoJulia,
    [ "IsSingerAlgebra" ],
    A -> CallFuncList( LoewyStructureInfoJulia,
                       ParametersOfSingerAlgebra( A ) ) );

InstallMethod( LoewyStructureInfoJulia,
    [ "IsPosInt", "IsPosInt" ],
    { q, z } -> LoewyStructureInfoJulia( q, OrderMod( q, z ), z ) );

InstallMethod( LoewyStructureInfoJulia,
    [ "IsPosInt", "IsPosInt", "IsPosInt" ],
    { q, n, z } -> CallFuncList( Julia.SingerAlg.LoewyStructureInfo,
                       List( [ q, n, z ], GAPToJulia ) ) );

fi;

if IsPackageMarkedForLoading( "JuliaInterface", "" ) then

#############################################################################
##
#M  LoewyStructureInfo( <A> )
##
##  We prefer 'LoewyStructureInfoJulia' to 'LoewyStructureInfoGAP'.
##  For that, we convert the Julia dictionary into a GAP record.
##
InstallMethod( LoewyStructureInfo,
    [ "IsSingerAlgebra" ],
    A -> JuliaToGAP( IsRecord, LoewyStructureInfoJulia( A ), true ) );

InstallMethod( LoewyStructureInfo,
    [ "IsPosInt", "IsPosInt" ],
    { q, z } -> LoewyStructureInfo( q, OrderMod( q, z ), z ) );

InstallMethod( LoewyStructureInfo,
    [ "IsPosInt", "IsPosInt", "IsPosInt" ],
    { q, n, z } -> JuliaToGAP( IsRecord, LoewyStructureInfoJulia( q, n, z ),
                               true ) );

fi;

#############################################################################
##
#F  LoewyStructureInfoGAP( <A> )
#F  LoewyStructureInfoGAP( <q>[, <n>], <z> )
##
InstallMethod( LoewyStructureInfoGAP,
    [ "IsSingerAlgebra" ],
    A -> CallFuncList( LoewyStructureInfoGAP,
                       ParametersOfSingerAlgebra( A ) ) );

InstallMethod( LoewyStructureInfoGAP,
    [ "IsPosInt", "IsPosInt" ],
    { q, z } -> LoewyStructureInfoGAP( q, OrderMod( q, z ), z ) );

InstallMethod( LoewyStructureInfoGAP,
    [ "IsPosInt", "IsPosInt", "IsPosInt" ],
    function( q, n, z )
    local islessorequal, monomials, layers, degrees, predecessors, m, i,
          mon, lambda, pred, j, mm;

    if n mod OrderMod( q, z ) <> 0 then
      Error( "<z> must divide <q>^<n> - 1" );
    fi;

    islessorequal:= function( mon1, mon2 )
      local i;

      for i in [ 1 .. n ] do
        if mon2[i] < mon1[i] then
          return false;
        fi;
      od;
      return true;
    end;

    monomials:= [ 0 * [ 1 .. n ] ];
    layers:= [ 0 ];
    degrees:= [ 0 ];    # just for speedup: avoid comparisons

    predecessors:= [ 0 ];

    m:= n * (q-1);
    for i in [ 1 .. z ] do
      mon:= CoefficientsQadicReversed( i, z, q, n );
      lambda:= 0;
      pred:= 1;
      mm:= Sum( mon );
      for j in [ 2 .. i ] do
        if lambda < layers[j] and degrees[j] < mm
                              and islessorequal( monomials[j], mon ) then
          lambda:= layers[j];
          pred:= j;
        fi;
      od;
      monomials[i+1]:= mon;
      layers[i+1]:= lambda + 1;
      degrees[i+1]:= mm;
      predecessors[i+1]:= pred;

      if lambda = 0 then
        if mm < m then
          m:= mm;
        fi;
      fi;
    od;

    # Extract information about one longest chain.
    i:= Length( monomials );
    pred:= [];
    while i > 0 do
      Add( pred, i );
      i:= predecessors[i];
    od;

    return rec( monomials:= monomials,
                layers:= layers,
                chain:= pred,
                m:= m,
                LL:= lambda + 2,
                parameters:= [ q, n, z ] );
    end );

#############################################################################
##
#F  DimensionsLoewyFactorsGAP( <A> )
#F  DimensionsLoewyFactorsJulia( <A> )
##
##  Return the GAP object or the Julia object, respectively.
##
BindGlobal( "DimensionsLoewyFactorsGAP", function( A )
    local data, v, i;

    data:= LoewyStructureInfoGAP( A );
    v:= ListWithIdenticalEntries( data.LL, 0 );
    for i in data.layers do
      v[ i+1 ]:= v[ i+1 ] + 1;
    od;

    return v;
    end );

if IsPackageMarkedForLoading( "JuliaInterface", "" ) then

BindGlobal( "DimensionsLoewyFactorsJulia",
    A -> Julia.SingerAlg.LoewyVector( LoewyStructureInfoJulia( A ) ) );

fi;

#############################################################################
##
#M  DimensionsLoewyFactors( <A> )
##
##  If Julia is available then use the Julia functionality,
##  otherwise use the GAP code.
##
InstallMethod( DimensionsLoewyFactors,
    [ "IsSingerAlgebra" ],
    A -> DimensionsLoewyFactorsGAP( A ) );

if IsPackageMarkedForLoading( "JuliaInterface", "" ) then

InstallMethod( DimensionsLoewyFactors,
    [ "IsSingerAlgebra" ],
    A -> JuliaToGAP( IsList, DimensionsLoewyFactorsJulia( A ), true ) );

fi;

#############################################################################
##
#F  LoewyVectorAbbreviated( <v> )
#F  LoewyVectorExpanded( <v> )
##
InstallGlobalFunction( LoewyVectorAbbreviated, function( v )
    local result, x, m, i;

    if not IsDenseList( v ) then
      Error( "<v> must be a dense list" );
    elif Length( v ) = 0 then
      return [];
    fi;
    result:= [];
    x:= v[1];
    m:= 1;
    for i in [ 2 .. Length( v ) ] do
      if v[i] = x then
        m:= m + 1;
      elif m = 1 then
        Add( result, x );
        x:= v[i];
      else
        Add( result, [ x, m ] );
        x:= v[i];
        m:= 1;
      fi;
    od;

    if m = 1 then
      Add( result, x );
    else
      Add( result, [ x, m ] );
    fi;

    return result;
    end );

InstallGlobalFunction( LoewyVectorExpanded, function( v )
    local result, x, i;

    if not IsDenseList( v ) then
      Error( "<v> must be a dense list" );
    fi;
    result:= [];
    for x in v do
      if IsList( x ) then
        for i in [ 1 .. x[2] ] do
          Add( result, x[1] );
        od;
      else
        Add( result, x );
      fi;
    od;

    return result;
    end );

#############################################################################
##
#F  SingerAlgebra( <q>[, <n>], <z>[, <R>] )
#F  SingerAlgebra( <arec>[, <R>] )
##
##  The algebra returned by this function claims to be a structure constants
##  algebra (via the filter 'IsSCAlgebraObjCollection'),
##  but the structure constants table is not created until one asks for
##  something that needs it.
##  Currently a special method for 'CanonicalBasis' makes sure that the
##  relevant data get stored.
##
InstallGlobalFunction( SingerAlgebra, function( arg )
    local arec, paras, q, n, e, R, z, zero, filter, Fam, A;

    if Length( arg ) = 1 and IsRecord( arg[1] ) then
      # The record is assumed to belong to the database of Singer algebras.
      arec:= arg[1];
      q:= arec.q;
      n:= arec.n;
      if IsBound( arec.z ) then
        z:= arec.z;
      elif IsBound( arec.e ) then
        z:= ( q^n - 1 ) / arec.e;
      else
        Error( "<arec> must contain one of the components z or e" );
      fi;
      R:= Rationals;
    elif Length( arg ) = 2 and IsRecord( arg[1] ) and IsRing( arg[2] ) then
      # The record is assumed to belong to the database of Singer algebras.
      arec:= arg[1];
      q:= arec.q;
      n:= arec.n;
      if IsBound( arec.z ) then
        z:= arec.z;
      elif IsBound( arec.e ) then
        z:= ( q^n - 1 ) / arec.e;
      else
        Error( "<arec> must contain one of the components z or e" );
      fi;
      R:= arg[2];
    elif Length( arg ) = 1 and IsList( arg[1] ) then
      # The parameters must be q, n, e[, R].
      paras:= arg[1];
      if Length( paras ) = 3 then
        q:= paras[1];
        n:= paras[2];
        e:= paras[3];
        R:= Rationals;
      elif Length( paras ) = 4 then
        q:= paras[1];
        n:= paras[2];
        e:= paras[3];
        R:= paras[4];
      else
        Error( "usage: SingerAlgebra( <list> ) with list of length 3 or 4" );
      fi;
      z:= ( q^n - 1 ) / e;
    elif Length( arg ) = 2 then
      # The parameters must be q, z.
      q:= arg[1];
      z:= arg[2];
      n:= OrderMod( q, z );
      R:= Rationals;
    elif Length( arg ) = 3 and IsInt( arg[3] ) then
      # The parameters must be q, n, z.
      q:= arg[1];
      n:= arg[2];
      z:= arg[3];
      R:= Rationals;
    elif Length( arg ) = 3 and IsRing( arg[3] ) then
      # The parameters must be q, z, R.
      q:= arg[1];
      z:= arg[2];
      R:= arg[3];
      n:= OrderMod( q, z );
    elif Length( arg ) = 4 then
      # The parameters must be q, n, z, R.
      q:= arg[1];
      n:= arg[2];
      z:= arg[3];
      R:= arg[4];
    else
      Error( "usage: SingerAlgebra( <arec>[, <R>] ) or ",
             "SingerAlgebra( <q>[, <n>], <z>[, <R>] ) or ",
             "SingerAlgebra( <list> )" );
    fi;

    if not ( IsPosInt( q ) and IsPosInt( n ) and IsPosInt( z ) ) then
      Error( "<q>, <n>, <z> must be positive integers" );
    elif z <> 1 and PowerMod( q, n, z ) <> 1 then
      Error( "<z> must divide <q>^<n> - 1" );
    elif q = 1 then
      Error( "<q> must be an integer > 1" );
    fi;

    # Create the algebra as far as necessary.
    # (Do not set a 'Name' value, in order to use the different methods
    # installed for 'ViewObj' and 'PrintObj'.)
    zero := Zero( R );
    filter:= IsSCAlgebraObj;
    if IsAdditivelyCommutativeElementFamily( FamilyObj( zero ) ) then
      filter:= filter and IsAdditivelyCommutativeElement;
    fi;
    Fam:= NewFamily( "SCAlgebraObjFamily", filter );
    if Zero( ElementsFamily( FamilyObj( R ) ) ) <> fail then
      SetFilterObj( Fam, IsFamilyOverFullCoefficientsFamily );
    else
      Fam!.coefficientsDomain:= R;
    fi;
    Fam!.zerocoeff := zero;
    SetCharacteristic( Fam, Characteristic( R ) );
    SetCoefficientsFamily( Fam, ElementsFamily( FamilyObj( R ) ) );

    A:= Objectify( NewType( CollectionsFamily( Fam ),
                   IsSingerAlgebra and IsAttributeStoringRep ),
                   rec() );
    SetLeftActingDomain( A, R );
    SetParametersOfSingerAlgebra( A, [ q, n, z ] );
    SetGeneratingSubsetOfCanonicalBasisOfSingerAlgebra( A, [ 1 .. z+1 ] );
    SetDimension( A, z+1 );
    Fam!.fullSCAlgebra:= A;
    SetIsFullSCAlgebra( A, true );

    return A;
    end );

#############################################################################
##
#M  Display( <A> )
#M  ViewObj( <A> )
#M  PrintObj( <A> )
#M  DisplayString( <A> )
#M  ViewString( <A> )
#M  PrintString( <A> )
#M  String( <A> )
##
##  According to the Section "View and Print" in the GAP Reference Manual,
##  we need methods (only) for
##  'String' (which then covers 'PrintString' and 'PrintObj') and
##  'ViewString' (which then covers 'ViewObj').
##  (We do not need 'DisplayString', which would cover 'Display'.)
##
##  We want to *view* Singer algebras as 'A[q,n,z[,R]]',
##  and to *print* them as calls to 'SingerAlgebra'.
##
InstallMethod( ViewString,
    [ "IsSingerAlgebra" ],
    function( A )
    local paras;

    paras:= ParametersOfSingerAlgebra( A );
    if LeftActingDomain( A ) = Rationals then
      return Concatenation( "A[",
                 JoinStringsWithSeparator( List( paras, String ), "," ),
                 "]" );
    else
      return Concatenation( "A[",
                 JoinStringsWithSeparator( List( paras, String ), "," ),
                 ",", String( LeftActingDomain( A ) ), "]" );
    fi;
    end );

InstallMethod( String,
    [ "IsSingerAlgebra" ],
    function( A )
    local paras;

    paras:= ParametersOfSingerAlgebra( A );
    if LeftActingDomain( A ) = Rationals then
      return Concatenation( "SingerAlgebra( ",
                 JoinStringsWithSeparator( List( paras, String ), ", " ),
                 " )" );
    else
      return Concatenation( "SingerAlgebra( ",
                 JoinStringsWithSeparator( List( paras, String ), ", " ),
                 ", ", String( LeftActingDomain( A ) ), " )" );
    fi;
    end );

#T Currently the above installations are unfortunately *not* enough,
#T since some 'ViewObj' and 'PrintObj' methods for algebras interfere,
#T as well as a 'PrintString' method for magmas.

InstallMethod( ViewObj,
    [ "IsSingerAlgebra" ],
    function( A ) Print( ViewString( A ) ); end );

InstallMethod( PrintObj,
    [ "IsSingerAlgebra" ],
    function( A ) Print( PrintString( A ) ); end );

InstallMethod( PrintString,
    [ "IsSingerAlgebra" ],
    String );

#############################################################################
##
#M  CanonicalBasis( <A> )
##
##  This method provides the internal data for treating the Singer algebra
##  <A> as a structure constants algebra in GAP.
##
##  Formally, we require those filters that are required also by the
##  GAP library method for full s.c. algebras,
##  in order to guarantee a higher rank for our method;
##  these filters are set in algebras returned by 'SingerAlgebra'.
##
InstallMethod( CanonicalBasis,
    [ Concatenation( "IsSingerAlgebra and IsFreeLeftModule and ",
                     "IsSCAlgebraObjCollection and IsFullSCAlgebra" ) ],
    function( A )
    local paras, q, n, z, dim, coeffs, T, R, zero, one, empty, nonempty,
          i, j, Fam, gens;

    # Create the structure constants table.
    paras:= ParametersOfSingerAlgebra( A );
    q:= paras[1];
    n:= paras[2];
    z:= paras[3];
    dim:= Dimension( A );
    coeffs:= List( [ 0 .. z ],
                   k -> CoefficientsQadicReversed( k, z, q, n ) );
    T:= [];
    R:= LeftActingDomain( A );
    zero:= Zero( R );
    one:= One( R );
    empty:= MakeImmutable( [ [], [] ] );
    nonempty:= List( [ 1 .. dim ],
                     i -> MakeImmutable( [ [ i ], [ one ] ] ) );
    for i in [ 1 .. dim ] do
      T[i]:= [];
      for j in [ 1 .. i-1 ] do
        T[i][j]:= T[j][i];
      od;
      for j in [ i .. dim-i+1 ] do
        if q <= MaximumList( coeffs[i] + coeffs[j], 0 ) then
          T[i][j]:= empty;
        else
          T[i][j]:= nonempty[ i+j-1 ];
        fi;
      od;
      for j in [ dim-i+2 .. dim ] do
        T[i][j]:= empty;
      od;
    od;
    T[ dim+1 ]:= 1;  # commutativity flag
    T[ dim+2 ]:= zero;

    # Set the necessary entries in the family.
    Fam:= ElementsFamily( FamilyObj( A ) );
    Fam!.sctable:= T;
    Fam!.names:= MakeImmutable( List( [ 0 .. dim-1 ],
                     i -> Concatenation( "b", String( i ) ) ) );
    Fam!.zerocoeff:= zero;
    Fam!.defaultTypeDenseCoeffVectorRep :=
        NewType( Fam, IsSCAlgebraObj and IsDenseCoeffVectorRep );
    SetZero( Fam, ObjByExtRep( Fam, ListWithIdenticalEntries( dim, zero ) ) );

    # Set the algebra generators.
    gens:= MakeImmutable( List( IdentityMat( dim, R ),
               x -> ObjByExtRep( Fam, x ) ) );
    SetGeneratorsOfAlgebra( A, gens );
    SetGeneratorsOfAlgebraWithOne( A, gens );
    SetOne( A, gens[1] );
    Fam!.basisVectors:= gens;

    # Delegate to the library method for full s.c. algebras,
    # which has a lower rank.
    TryNextMethod();
    end );

#############################################################################
##
#M  Representative( <A> )
#M  GeneratorsOfAlgebra( <A> )
#M  GeneratorsOfAlgebraWithOne( <A> )
##
##  Note that we cannot use 'RedispatchOnCondition' here,
##  because we have only the attribute tester 'HasCanonicalBasis'
##  that may be missing,
##  whereas 'RedispatchOnCondition' checks for missing property values
##  plus their testers.
##
InstallMethod( Representative,
    [ "IsSingerAlgebra" ],
    function( A )
    if HasCanonicalBasis( A ) then
      TryNextMethod();
    fi;

    # Set the necessary data and redispatch.
    CanonicalBasis( A );
    return Representative( A );
    end );

InstallMethod( GeneratorsOfAlgebra,
    [ "IsSingerAlgebra" ],
    function( A )
    if HasCanonicalBasis( A ) then
      TryNextMethod();
    fi;

    # Set the necessary data and redispatch.
    CanonicalBasis( A );
    return GeneratorsOfAlgebra( A );
    end );

InstallMethod( GeneratorsOfAlgebraWithOne,
    [ "IsSingerAlgebra" ],
    function( A )
    if HasCanonicalBasis( A ) then
      TryNextMethod();
    fi;

    # Set the necessary data and redispatch.
    CanonicalBasis( A );
    return GeneratorsOfAlgebraWithOne( A );
    end );

#############################################################################
##
#M  GeneratingSubsetOfCanonicalBasisOfSingerAlgebra( <triv> )
##
##  Provide a method for the trivial subspace/subalgebra of a Singer algebra.
##  (We cannot be sure that this attribute gets set on creation.)
##
InstallMethod( GeneratingSubsetOfCanonicalBasisOfSingerAlgebra,
    [ "IsVectorSpace and HasParent" ],
    function( V )
    if IsSingerAlgebra( Parent( V ) ) and IsTrivial( V ) then
      return [];
    fi;
    TryNextMethod();
    end );

#############################################################################
##
#M  LoewyLengthGAP( <A> )
#M  LoewyLengthGAP( <q>, <z> )
#M  LoewyLengthGAP( <q>, <n>, <z>, <m> )
##
##  Use stored 'LoewyStructureInfoGAP' values,
##  use the cheap criteria when the upper bound is attained (using <m>),
##  compute 'LoewyStructureInfoGAP' only if necessary.
##
InstallMethod( LoewyLengthGAP,
    [ "IsSingerAlgebra and HasLoewyStructureInfoGAP" ],
    A -> LoewyStructureInfoGAP( A ).LL );

InstallMethod( LoewyLengthGAP,
    [ "IsSingerAlgebra" ],
    function( A )
    local m, paras;

    # We need m(q,e) anyhow in order to compute the upper bound
    # to which the cheap criteria refer.
    m:= MinimalDegreeOfSingerAlgebraGAP( A );
    if HasLoewyStructureInfoGAP( A ) then
      # Perhaps the computation of 'm' has triggered this.
      return LoewyStructureInfoGAP( A ).LL;
    else
      paras:= ParametersOfSingerAlgebra( A );
      return LoewyLengthGAP( paras[1], paras[2], paras[3], m );
    fi;
    end );

InstallMethod( LoewyLengthGAP,
    [ "IsPosInt", "IsPosInt" ],
    function( q, z )
    local n;

    n:= OrderMod( q, z );
    return LoewyLengthGAP( q, n, z,
               MinimalDegreeOfSingerAlgebraGAP( q, SingerAlgE( q, n, z ) ) );
    end );

InstallMethod( LoewyLengthGAP,
    [ "IsPosInt", "IsPosInt", "IsPosInt", "IsPosInt" ],
    function( q, n, z, m )
    local l;

    l:= LoewyLengthCheapGAP( q, n, z, m );
    if l = 0 then
      l:= LoewyLengthHardGAP( q, n, z, m );
    fi;
    return l;
    end );

#############################################################################
##
#F  OrderModExt( <n>, <m>[, <bound>] )
##
##  The optional argument <bound> is not yet supported for 'OrderMod'
##  in GAP 4.11.0, but it is useful in the context of this package.
##
InstallGlobalFunction( OrderModExt, function( n, m, bound... )
    local  x, o, d;

    # check the arguments and reduce $n$ into the range $0..m-1$
    if m <= 0  then Error("<m> must be positive");  fi;
    if n < 0   then n := n mod m + m;  fi;
    if m <= n  then n := n mod m;      fi;

    # return 0 if $m$ is not coprime to $n$
    if GcdInt(m,n) <> 1  then
        o := 0;

    # compute the order simply by iterated multiplying, $x= n^o$ mod $m$
    elif m < 100  then
        x := n;  o := 1;
        while x > 1  do
            x := x * n mod m;  o := o + 1;
        od;

    # otherwise try the divisors of $\lambda(m)$ and their divisors, etc.
    else
        if Length( bound ) = 1 then
            # We know a multiple of the desired order.
            o := bound[1];
        else
            # The default a priori known multiple is 'Lambda( m )'.
            o := Lambda( m );
        fi;
        for d in PrimeDivisors( o ) do
            while o mod d = 0  and PowerModInt(n,o/d,m) = 1  do
                o := o / d;
            od;
        od;

    fi;

    return o;
    end );

#############################################################################
##
#F  SufficientCriterionForLoewyBoundAttained( <q>, <n>, <z>, <m> )
##
InstallGlobalFunction( SufficientCriterionForLoewyBoundAttained,
    function( q, n, z, m )
    local qm, e, r, p, sum, d, phi6;

    if n <= 3 then
      return "Cor. I.7.1 (n <= 3)";
    elif z < 70 then
      # The bound is attained, by the classif. of small dim. Singer alg.
      return "z < 70";
    elif m = 2 then
      return "La. I.6.3";
    elif ( q - 1 ) mod m = 0 then
      return "Thm. I.7.1";
    elif n * (q-1) < 3 * m then
      return "La. I.7.1 (iii)";
    fi;

    qm:= q mod m;
    if 1 < qm and Int( n * ( qm-1 ) / m ) = 1 then
      return "Prop. II.3.15";
    fi;

    e:= ( q^n - 1 ) / z;

    if e <= 32 then
      return "Prop. II.6.1 (e <= 32)";
    elif e <= 3 * m then
      return "Prop. II.6.4";
    elif n <= 5 and ((q^n-1)/(q-1)) mod e = 0 then
      return "Prop. II.5.3, II.5.6 (e | (q^n-1)/(q-1), n <= 5)";
    elif n mod ( m / Gcd( m, q-1 ) * OrderModExt( q, e, n ) ) = 0 then
      return "La. II.5.2 (ii)";
    fi;

    if ( (q^n-1)/(q-1) ) mod e = 0 then
      # we know m <= n
      if m >= n-1 then
        return "Prop. II.5.1 (ii)";
      elif  m = n-2 and ( q mod m ) in [ 1 .. Int( (m+1)/2 ) ] then
        return "Prop. II.5.1 (iii)";
      elif Int( (n-m) * (q-1) / m ) mod ( n-m ) = 0 then
        return "Prop. II.5.1 (i)";
      fi;
    fi;

    if ( n mod m = 0 ) then
      if Sum( List( [ 0 .. m-1 ], i -> q^(i*n/m) ) ) mod e = 0 then
        return "La. II.4.1 for r = 1";
      fi;
      for r in DivisorsInt( n ) do
        if (n/r) mod m = 0 and
           Sum( [ 0 .. m-1 ], i -> q^(n*i/(m*r)) ) mod e = 0 then
          return "La. II.4.1";
        fi;
      od;
    fi;

    p:= SmallestRootInt( e );

    if p = 2 then
      return "Thm. II.4.3 (iii)";
    elif IsPrimeInt( p ) then
      # (Here we know that 'q mod p <> 1' holds.)
      # Check whether 'p' is a Pierpont prime.
      p:= p - 1;
      while p mod 2 = 0 do
        p:= p/2;
      od;
      while p mod 3 = 0 do
        p:= p/3;
      od;
      if p = 1 then
        return "Thm. II.4.3 (ii)";
      fi;
    fi;

    # Concerning Remark II.5.4,
    # we need not check \Phi_2 (since m = e if e divides q-1,
    # and we have checked whether m divides q-1),
    # \Phi_4, \Phi_8 (because then m = 2).
    phi6:= q^2 - q + 1;    # this is \Phi_6(q)
    if phi6 mod e = 0 or
       ( q^3*(q^3+1) + 1 ) mod e = 0 or     # this is \Phi_9(q)
       ( q^3*(q-1) + phi6 ) mod e = 0 then  # this is \Phi_{10}(q)
      return "Rem. II.5.4";
    fi;

    # We give up.
    return "";
    end );

InstallGlobalFunction( LoewyLengthCheapGAP, function( q, n, z, m )
    if SufficientCriterionForLoewyBoundAttained( q, n, z, m ) <> "" then
      return Int( n * ( q - 1 ) / m ) + 1;
    else
      return 0;
    fi;
    end );

InstallGlobalFunction( LoewyLengthHardGAP,
    { q, n, z, m } -> LoewyStructureInfoGAP( q, n, z ).LL );

#############################################################################
##
#M  LoewyLengthJulia( <A> )
#M  LoewyLengthJulia( <q>, <z> )
#M  LoewyLengthJulia( <q>, <n>, <z>, <m> )
##
##  Use the same criteria as for the &GAP; variant.
##  These methods are available only if &Julia; is available.
##
if IsPackageMarkedForLoading( "JuliaInterface", "" ) then

InstallMethod( LoewyLengthJulia,
    [ "IsSingerAlgebra and HasLoewyStructureInfoJulia" ],
    A -> JuliaToGAP( IsInt,
             Julia.Base.get( LoewyStructureInfoJulia( A ),
                             JuliaSymbol( "LL" ), 0 )  ));

InstallMethod( LoewyLengthJulia,
    [ "IsSingerAlgebra" ],
    function( A )
    local paras;

    paras:= ParametersOfSingerAlgebra( A );

    # Leave the computation of m(q,e) and of the Loewy length to Julia,
    # let Julia set also the 'MinimalDegreeOfSingerAlgebraJulia' value in A.
    return JuliaToGAP( IsInt,
               Julia.SingerAlg.LoewyLength( paras[1], paras[2], paras[3],
                   A ) );
    end );

InstallMethod( LoewyLengthJulia,
    [ "IsPosInt", "IsPosInt" ],
    function( q, z )
    return JuliaToGAP( IsInt, Julia.SingerAlg.LoewyLength( q, z ) );
    end );

InstallMethod( LoewyLengthJulia,
    [ "IsPosInt", "IsPosInt", "IsPosInt", "IsPosInt" ],
    function( q, n, z, m )
    return JuliaToGAP( IsInt,
               Julia.SingerAlg.LoewyLength( q, n, z, m ) );
    end );

fi;

#############################################################################
##
#M  RadicalOfAlgebra( <A> )
##
##  The Jacobson radical of a Singer algebra of dimension <M>z+1</M>
##  has the basis <M>[ b_1, b_2, \ldots, b_z ]</M>.
##
InstallMethod( RadicalOfAlgebra,
    [ "IsSingerAlgebra" ],
    function( A )
    local B, range, S;

    B:= BasisVectors( CanonicalBasis( A ) );
    range:= [ 2 .. Length( B ) ];
    S:= SubalgebraNC( A, B{ range }, "basis" );
    SetGeneratingSubsetOfCanonicalBasisOfSingerAlgebra( S, range );

    return S;
    end );

#############################################################################
##
#M  RadicalSeriesOfAlgebra( <A> )
##
InstallMethod( RadicalSeriesOfAlgebra,
    [ "IsSingerAlgebra" ],
    function( A )
    local result, data, B, i, S;

    result:= [ A, RadicalOfAlgebra( A ) ];
    data:= LoewyStructureInfoGAP( A );
    data:= SingerAlg.BasesOfRadicalSeries( data );
    B:= BasisVectors( CanonicalBasis( A ) );
    for i in [ 2 .. Length( data ) ] do
      S:= SubalgebraNC( A, B{ data[i] }, "basis" );
      SetGeneratingSubsetOfCanonicalBasisOfSingerAlgebra( S, data[i] );
      result[i+1]:= S;
    od;
    Add( result, TrivialSubalgebra( A ) );

    return result;
    end );

#############################################################################
##
#M  SocleSeriesOfAlgebra( <A> )
##
InstallMethod( SocleSeriesOfAlgebra,
    [ "IsSingerAlgebra" ],
    function( A )
    local result, data, B, i, S;

    result:= [ TrivialSubalgebra( A ) ];
    data:= LoewyStructureInfoGAP( A );
    data:= SingerAlg.BasesOfSocleSeries( data );
    B:= BasisVectors( CanonicalBasis( A ) );
    for i in [ 1 .. Length( data ) ] do
      S:= SubalgebraNC( A, B{ data[i] }, "basis" );
      SetGeneratingSubsetOfCanonicalBasisOfSingerAlgebra( S, data[i] );
      result[i+1]:= S;
    od;
    Add( result, A );

    return result;
    end );

#############################################################################
##
#M  Intersection2( <A>, <B> )
##
InstallMethod( Intersection2,
    IsIdenticalObj,
    [ "IsVectorSpace and HasGeneratingSubsetOfCanonicalBasisOfSingerAlgebra",
      "IsVectorSpace and HasGeneratingSubsetOfCanonicalBasisOfSingerAlgebra" ],
    100,  # beat generic methods for two 'IsFLMLORWithOne'
    function( A, B )
    local P, basis, data, S;

    P:= Parent( A );
    if not IsIdenticalObj( P, Parent( B ) ) then
      TryNextMethod();
    fi;
    basis:= BasisVectors( CanonicalBasis( P ) );
    data:= Intersection2(
               GeneratingSubsetOfCanonicalBasisOfSingerAlgebra( A ),
               GeneratingSubsetOfCanonicalBasisOfSingerAlgebra( B ) );
    S:= SubspaceNC( P, basis{ data }, "basis" );
    SetGeneratingSubsetOfCanonicalBasisOfSingerAlgebra( S, data );

    return S;
    end );

#############################################################################
##
#M  \+( <A>, <B> )
##
InstallOtherMethod( \+,
    IsIdenticalObj,
    [ "IsVectorSpace and HasGeneratingSubsetOfCanonicalBasisOfSingerAlgebra",
      "IsVectorSpace and HasGeneratingSubsetOfCanonicalBasisOfSingerAlgebra" ],
    100,  # beat generic methods for two rings, two left modules
    function( A, B )
    local P, basis, data, S;

    P:= Parent( A );
    if not IsIdenticalObj( P, Parent( B ) ) then
      TryNextMethod();
    fi;
    basis:= BasisVectors( CanonicalBasis( P ) );
    data:= Union(
               GeneratingSubsetOfCanonicalBasisOfSingerAlgebra( A ),
               GeneratingSubsetOfCanonicalBasisOfSingerAlgebra( B ) );
    S:= SubspaceNC( P, basis{ data }, "basis" );
    SetGeneratingSubsetOfCanonicalBasisOfSingerAlgebra( S, data );

    return S;
    end );

#############################################################################
##
#M  ProductSpace( <A>, <B> )
##
InstallMethod( ProductSpace,
    IsIdenticalObj,
    [ "IsVectorSpace and HasGeneratingSubsetOfCanonicalBasisOfSingerAlgebra",
      "IsVectorSpace and HasGeneratingSubsetOfCanonicalBasisOfSingerAlgebra" ],
    100,  # beat generic methods for two algebras, two left modules
    function( A, B )
    local P, basis, data, S;

    P:= Parent( A );
    if not IsIdenticalObj( P, Parent( B ) ) then
      TryNextMethod();
    fi;
    basis:= BasisVectors( CanonicalBasis( P ) );
    data:= LoewyStructureInfoGAP( P );
    data:= SingerAlg.BasisOfProductSpace( data,
               GeneratingSubsetOfCanonicalBasisOfSingerAlgebra( A ),
               GeneratingSubsetOfCanonicalBasisOfSingerAlgebra( B ) );
    S:= SubspaceNC( P, basis{ data }, "basis" );
    SetGeneratingSubsetOfCanonicalBasisOfSingerAlgebra( S, data );

    return S;
    end );

#############################################################################
##
#F  CoefficientsQadicReversed( <k>, <z>, <q>, <n> )
##
InstallGlobalFunction( CoefficientsQadicReversed, function( k, z, q, n )
    local v, i, r, rq, d;

    v:= [];

    if k = z then
      for i in [ 1 .. n ] do
        v[i]:= q-1;
      od;
    else
      r:= k;
      for i in [ 1 .. n ] do
        rq:= r * q;
        d:= QuoInt( rq, z );
        r:= rq - d * z;
        v[i]:= d;
      od;
    fi;

    return v;
    end );

#############################################################################
##
#F  SingerAlg.MultTable( <data> )
##
##  <#GAPDoc Label="SingerAlg.MultTable">
##  <ManSection>
##  <Func Name="SingerAlg.MultTable" Arg='data'/>
##
##  <Description>
##  Let <A>data</A> be the record returned by
##  <Ref Oper="LoewyStructureInfoGAP" Label="for parameters"/>
##  that describes a Singer algebra <M>A[q,z]</M>,
##  and let <M>B = B(A)</M>.
##  <Ref Func="SingerAlg.MultTable"/> returns the
##  <M>(z+1) \times (z+1)</M> matrix that contains at the position
##  <M>(i,j)</M> the value <M>i+j-1</M> if the product
##  <M>B_i \cdot B_j</M> is nonzero
##  (hence equal to <M>B_{{i+j-1}}</M>, the <M>i+j-1</M>-th basis vector),
##  and <M>0</M> otherwise.
##  <P/>
##  <Example><![CDATA[
##  gap> data:= LoewyStructureInfoGAP( 2, 3, 7 );;
##  gap> Display( SingerAlg.MultTable( data ) );
##  [ [  1,  2,  3,  4,  5,  6,  7,  8 ],
##    [  2,  0,  4,  0,  6,  0,  8,  0 ],
##    [  3,  4,  0,  0,  7,  8,  0,  0 ],
##    [  4,  0,  0,  0,  8,  0,  0,  0 ],
##    [  5,  6,  7,  8,  0,  0,  0,  0 ],
##    [  6,  0,  8,  0,  0,  0,  0,  0 ],
##    [  7,  8,  0,  0,  0,  0,  0,  0 ],
##    [  8,  0,  0,  0,  0,  0,  0,  0 ] ]
##  ]]></Example>
##  <P/>
##  The multiplication table of a Singer algebra of dimension <M>N</M>
##  has the value <M>N</M> on the antidiagonal (<M>i+j = N+1</M>),
##  is zero below the antidiagonal,
##  and contains only <M>N-i</M> and zero on the <M>i</M> parallel above the
##  antidiagonal.
##  </Description>
##  </ManSection>
##  <#/GAPDoc>
##
SingerAlg.MultTable:= function( data )
    local monoms, n, nn, mat, q, qhalf, i, j;

    if not IsBound( data.multtable ) then
      monoms:= data.monomials;
      n:= Length( monoms ); # z+1
      nn:= QuoInt( n+1, 2 );
      mat:= NullMat( n, n );
      q:= data.parameters[1];
      if q mod 2 = 1 then
        qhalf:= QuoInt( q, 2 ) + 1;
      else
        qhalf:= QuoInt( q, 2 );
      fi;
      for i in [ 1 .. nn ] do
        for j in [ 1 .. i-1 ] do
          if ForAll( monoms[i] + monoms[j], x -> x < q ) then
            mat[i,j]:= i+j-1;
            mat[j,i]:= i+j-1;
          fi;
        od;
        if ForAll( monoms[i], x -> x < qhalf ) then
          mat[i,i]:= i+i-1;
        fi;
      od;
      for i in [ nn+1 .. n ] do
        # The [i,j] entry can be nonzero only if i+j <= z+2 holds.
        for j in [ 1 .. n+1-i ] do
          if ForAll( monoms[i] + monoms[j], x -> x < q ) then
            mat[i,j]:= i+j-1;
            mat[j,i]:= i+j-1;
          fi;
        od;
      od;
      data.multtable:= mat;
    fi;

    return data.multtable;
    end;

#############################################################################
##
#F  SingerAlg.BasisOfSum( <data>, <I>, <J> )
#F  SingerAlg.BasisOfIntersection( <data>, <I>, <J> )
##
##  <#GAPDoc Label="SingerAlg.BasisOfSum">
##  <ManSection>
##  <Heading>SingerAlg.BasisOfSum and SingerAlg.BasisOfIntersection</Heading>
##  <Func Name="SingerAlg.BasisOfSum" Arg='data, I, J'/>
##  <Func Name="SingerAlg.BasisOfIntersection" Arg='data, I, J'/>
##
##  <Description>
##  For two subsets <A>I</A>, <A>J</A> of <M>\{ 1, 2, \ldots, z+1 \}</M>,
##  these functions just return the union and the intersection,
##  respectively, of <A>I</A> and <A>J</A>.
##  <P/>
##  <A>I</A> and <A>J</A> describe subsets of a basis, which generate the
##  spaces <M>U</M> and <M>V</M>, say, then the result describes the subset
##  of this basis that generates the sum and the intersection, respectively,
##  of <M>U</M> and <M>V</M>.
##  <P/>
##  <Example><![CDATA[
##  gap> SingerAlg.BasisOfSum( data, [ 1, 2, 3 ], [ 2, 4, 6 ] );
##  [ 1, 2, 3, 4, 6 ]
##  gap> SingerAlg.BasisOfIntersection( data, [ 1, 2, 3 ], [ 2, 4, 6 ] );
##  [ 2 ]
##  ]]></Example>
##  </Description>
##  </ManSection>
##  <#/GAPDoc>
##
SingerAlg.BasisOfSum:= { data, I, J } -> Union( I, J );

SingerAlg.BasisOfIntersection:= { data, I, J } -> Intersection( I, J );

#############################################################################
##
#F  SingerAlg.BasisOfProductSpace( <data>, <I>, <J> )
##
##  <#GAPDoc Label="SingerAlg.BasisOfProductSpace">
##  <ManSection>
##  <Func Name="SingerAlg.BasisOfProductSpace" Arg='data, I, J'/>
##
##  <Description>
##  Let <A>data</A> be the record returned by
##  <Ref Oper="LoewyStructureInfoGAP" Label="for parameters"/>
##  that describes a Singer algebra <M>A = A[q,z]</M>,
##  let <M>B = B(A)</M>,
##  and let <A>I</A>, <A>J</A> be subsets of
##  <M>\{ 1, 2, \ldots, z+1 \}</M>,
##  describing subspaces <M>U</M>, <M>V</M> of <M>A</M> with bases
##  <M>( B_i; i \in I )</M> and <M>( B_i; i \in J )</M>, respectively.
##  <Ref Func="SingerAlg.BasisOfProductSpace"/> returns the subset <M>K</M>
##  of <M>\{ 1, 2, \ldots, z+1 \}</M> such that
##  <M>( B_i; i \in K )</M> is a basis of the product space <M>U \cdot V</M>.
##  <P/>
##  <Example><![CDATA[
##  gap> data:= LoewyStructureInfoGAP( 2, 7 );;
##  gap> radser:= SingerAlg.BasesOfRadicalSeries( data );
##  [ [ 2 .. 8 ], [ 4, 6, 7, 8 ], [ 8 ] ]
##  gap> SingerAlg.BasisOfProductSpace( data, radser[1], radser[1] );
##  [ 4, 6, 7, 8 ]
##  gap> SingerAlg.BasisOfProductSpace( data, radser[1], radser[2] );
##  [ 8 ]
##  gap> SingerAlg.BasisOfProductSpace( data, radser[2], radser[2] );
##  [  ]
##  ]]></Example>
##  </Description>
##  </ManSection>
##  <#/GAPDoc>
##
SingerAlg.BasisOfProductSpace:= function( data, I, J )
    local mat, basis, i, j;

    mat:= SingerAlg.MultTable( data );
    basis:= [];

    for i in I do
      for j in J do
        if mat[i,j] <> 0 then
          AddSet( basis, i+j-1 );
        fi;
      od;
    od;

    return basis;
    end;

#############################################################################
##
#F  SingerAlg.BasisOfIdeal( <data>, <I> )
##
##  <#GAPDoc Label="SingerAlg.BasisOfIdeal">
##  <ManSection>
##  <Func Name="SingerAlg.BasisOfIdeal" Arg='data, I'/>
##
##  <Description>
##  Let <A>data</A> be the record returned by
##  <Ref Oper="LoewyStructureInfoGAP" Label="for parameters"/>
##  that describes a Singer algebra <M>A = A[q,z]</M>,
##  let <M>B = B(A)</M>,
##  and let <A>I</A> be a subset of <M>\{ 1, 2, \ldots, z+1 \}</M>,
##  describing a subspace <M>U</M> of <M>A</M> with basis
##  <M>( B_i; i \in I )</M>.
##  <Ref Func="SingerAlg.BasisOfIdeal"/> returns the subset <M>J</M>
##  of <M>\{ 1, 2, \ldots, z+1 \}</M> such that
##  <M>( B_i; i \in J )</M> is a basis of the ideal <M>U \cdot A</M>.
##  <P/>
##  <Example><![CDATA[
##  gap> data:= LoewyStructureInfoGAP( 2, 7 );;
##  gap> SingerAlg.BasisOfIdeal( data, [ 4 ] );
##  [ 4, 8 ]
##  ]]></Example>
##  </Description>
##  </ManSection>
##  <#/GAPDoc>
##
SingerAlg.BasisOfIdeal:= { data, I } -> SingerAlg.BasisOfProductSpace( data,
              I, [ 1 .. Length( data.monomials ) ] );

#############################################################################
##
#F  SingerAlg.BasisOfAnnihilator( <data>, <I> )
##
##  <#GAPDoc Label="SingerAlg.BasisOfAnnihilator">
##  <ManSection>
##  <Func Name="SingerAlg.BasisOfAnnihilator" Arg='data, I'/>
##
##  <Description>
##  Let <A>data</A> be the record returned by
##  <Ref Oper="LoewyStructureInfoGAP" Label="for parameters"/>
##  that describes a Singer algebra <M>A = A[q,z]</M>,
##  let <M>B = B(A)</M>,
##  and let <A>I</A> be a subset of <M>\{ 1, 2, \ldots, z+1 \}</M>,
##  describing a subspace <M>U</M> of <M>A</M> with basis
##  <M>( B_i; i \in I )</M>.
##  <Ref Func="SingerAlg.BasisOfAnnihilator"/> returns the subset <M>J</M>
##  of <M>\{ 1, 2, \ldots, z+1 \}</M> such that
##  <M>( B_i; i \in J )</M> is a basis of the annihilator
##  <M>\{ x \in A; x \cdot U = 0 \}</M> of <M>U</M> in <M>A</M>.
##  <P/>
##  <Example><![CDATA[
##  gap> data:= LoewyStructureInfoGAP( 2, 7 );;
##  gap> radser:= SingerAlg.BasesOfRadicalSeries( data );
##  [ [ 2 .. 8 ], [ 4, 6, 7, 8 ], [ 8 ] ]
##  gap> List( radser, I -> SingerAlg.BasisOfAnnihilator( data, I ) );
##  [ [ 8 ], [ 4, 6, 7, 8 ], [ 2, 3, 4, 5, 6, 7, 8 ] ]
##  ]]></Example>
##  </Description>
##  </ManSection>
##  <#/GAPDoc>
##
SingerAlg.BasisOfAnnihilator:= function( data, list )
    local mat;

    mat:= SingerAlg.MultTable( data );
    return Filtered( [ 1 .. Length( mat ) ],
                     i -> ForAll( list, j -> mat[i,j] = 0 ) );
    end;

#############################################################################
##
#F  SingerAlg.BasesOfRadicalSeries( <data> )
##
##  <#GAPDoc Label="SingerAlg.BasesOfRadicalSeries">
##  <ManSection>
##  <Func Name="SingerAlg.BasesOfRadicalSeries" Arg='data'/>
##
##  <Description>
##  Let <A>data</A> be the record returned by
##  <Ref Oper="LoewyStructureInfoGAP" Label="for parameters"/>
##  that describes a Singer algebra <M>A = A[q,z]</M>,
##  and let <M>B = B(A)</M>.
##  <Ref Func="SingerAlg.BasesOfRadicalSeries"/> returns the list
##  <M>[ I_1, I_2, \ldots, I_l ]</M> of subsets of
##  <M>\{ 1, 2, \ldots, z+1 \}</M> such that
##  <M>( B_i; i \in I_j )</M> is a basis of the <M>j</M>-th power of the
##  Jacobson radical <M>J</M> of <M>A</M>,
##  and such that <M>J^l</M> is nonzero and <M>J^{{l+1}}</M> is zero.
##  <P/>
##  <Example><![CDATA[
##  gap> data:= LoewyStructureInfoGAP( 2, 7 );;
##  gap> radser:= SingerAlg.BasesOfRadicalSeries( data );
##  [ [ 2 .. 8 ], [ 4, 6, 7, 8 ], [ 8 ] ]
##  ]]></Example>
##  </Description>
##  </ManSection>
##  <#/GAPDoc>
##
SingerAlg.BasesOfRadicalSeries:= function( data )
    local lay, ll, jbylayer, J, i;

    lay:= data.layers;
    ll:= data.LL;
    jbylayer:= List( [ 1 .. ll-1 ], i -> Positions( lay, i ) );
    J:= [ jbylayer[ ll-1 ] ];
    for i in [ 2 .. ll-1 ] do
      J[i]:= Union( J[ i-1 ], jbylayer[ ll - i ] );
    od;

    return Reversed( J );
    end;

#############################################################################
##
#F  SingerAlg.BasesOfSocleSeries( <data> )
##
##  <#GAPDoc Label="SingerAlg.BasesOfSocleSeries">
##  <ManSection>
##  <Func Name="SingerAlg.BasesOfSocleSeries" Arg='data'/>
##
##  <Description>
##  Let <A>data</A> be the record returned by
##  <Ref Oper="LoewyStructureInfoGAP" Label="for parameters"/>
##  that describes a Singer algebra <M>A = A[q,z]</M>,
##  and let <M>B = B(A)</M>.
##  <Ref Func="SingerAlg.BasesOfSocleSeries"/> returns the list
##  <M>[ I_1, I_2, \ldots, I_l ]</M> of subsets of
##  <M>\{ 1, 2, \ldots, z+1 \}</M> such that
##  <M>( B_i; i \in I_j )</M> is a basis of <M>S_j</M>,
##  where <M>S_1</M> is the socle of <M>A</M>,
##  <M>S_{{j+1}}/S_j</M> is the socle of <A>A</A><M>/S_j</M>,
##  and <A>A</A><M>/S_l</M> is nonzero and its own socle.
##  <P/>
##  <Example><![CDATA[
##  gap> socser:= SingerAlg.BasesOfSocleSeries( data );
##  [ [ 8 ], [ 4, 6, 7, 8 ], [ 2 .. 8 ] ]
##  ]]></Example>
##  </Description>
##  </ManSection>
##  <#/GAPDoc>
##
SingerAlg.BasesOfSocleSeries:= function( data )
    local lay, ll, sbylayer, S, i;

    lay:= Reversed( data.layers );
    ll:= data.LL;
    sbylayer:= List( [ 0 .. ll-2 ], i -> Positions( lay, i ) );
    S:= [ sbylayer[ 1 ] ];
    for i in [ 2 .. ll-1 ] do
      S[i]:= Union( S[ i-1 ], sbylayer[i] );
    od;

    return S;
    end;

#############################################################################
##
#F  SingerAlg.BasisOfPowers( <data>, <I>, <p>, <m> )
##
##  <#GAPDoc Label="SingerAlg.BasisOfPowers">
##  <ManSection>
##  <Func Name="SingerAlg.BasisOfPowers" Arg='data, I, p, m'/>
##
##  <Description>
##  Let <A>p</A> be a prime integer,
##  let <A>data</A> be the record returned by
##  <Ref Oper="LoewyStructureInfoGAP" Label="for parameters"/>
##  that describes a Singer algebra <M>A = A[q,z]</M>,
##  let <M>B = B(A_p)</M>,
##  let <A>I</A> be a subset of <M>\{ 1, 2, \ldots, z+1 \}</M>,
##  describing a subspace <M>U</M> of <M>A_p</M> with basis
##  <M>( B_i; i \in I )</M>,
##  and let <A>m</A> be a positive integer.
##  <Ref Func="SingerAlg.BasisOfPowers"/> returns the subset <M>J</M>
##  of <M>\{ 1, 2, \ldots, z+1 \}</M> such that
##  <M>( B_i; i \in J )</M> is a basis of the subspace
##  <M>\{ x^{{p^m}}; x \in U \}</M> of <M>A_p</M>.
##  <P/>
##  <Example><![CDATA[
##  gap> data:= LoewyStructureInfoGAP( 3, 8 );;
##  gap> SingerAlg.BasisOfPowers( data, [ 1 .. 9 ], 2, 1 );
##  [ 1, 3, 7, 9 ]
##  gap> SingerAlg.BasisOfPowers( data, [ 1 .. 9 ], 2, 2 );
##  [ 1 ]
##  gap> SingerAlg.BasisOfPowers( data, [ 1 .. 9 ], 3, 1 );
##  [ 1 ]
##  ]]></Example>
##  </Description>
##  </ManSection>
##  <#/GAPDoc>
##
SingerAlg.BasisOfPowers:= function( data, I, p, m )
    local result, q, pm, quo, k;

    result:= [];
    q:= data.parameters[1];
    pm:= p^m;
    quo:= QuoInt( q-1, pm );
    for k in I do
      if Maximum( data.monomials[k] ) <= quo then
        Add( result, pm * ( k-1 ) + 1 );
      fi;
    od;

    return result;
    end;

#############################################################################
##
#F  SingerAlg.BasisOfPMRoots( <data>, <I>, <p>, <m> )
##
##  <#GAPDoc Label="SingerAlg.BasisOfPMRoots">
##  <ManSection>
##  <Func Name="SingerAlg.BasisOfPMRoots" Arg='data, I, p, m'/>
##
##  <Description>
##  Let <A>p</A> be a prime integer,
##  let <A>data</A> be the record returned by
##  <Ref Oper="LoewyStructureInfoGAP" Label="for parameters"/>
##  that describes a Singer algebra <M>A = A[q,z]</M>,
##  let <M>B = B(A_p)</M>,
##  let <A>I</A> be a subset of <M>\{ 1, 2, \ldots, z+1 \}</M>,
##  describing a subspace <M>U</M> of <M>A_p</M> with basis
##  <M>( B_i; i \in I )</M>,
##  and let <A>m</A> be a positive integer.
##  (See Section <Ref Sect="sect:Singer algebras"/> for the definition of
##  <M>A_p</M>.)
##  <P/>
##  <Ref Func="SingerAlg.BasisOfPMRoots"/> returns the subset <M>J</M>
##  of <M>\{ 1, 2, \ldots, z+1 \}</M> such that
##  <M>( B_i; i \in J )</M> is a basis of the subspace
##  <M>\{ x \in A_p; x^{{p^m}} \in U \}</M> of <M>A_p</M>.
##  <P/>
##  <Example><![CDATA[
##  gap> data:= LoewyStructureInfoGAP( 3, 8 );;
##  gap> SingerAlg.BasisOfPMRoots( data, [], 2, 1 );
##  [ 3, 6, 7, 8, 9 ]
##  gap> SingerAlg.BasisOfPMRoots( data, [], 2, 2 );
##  [ 2, 3, 4, 5, 6, 7, 8, 9 ]
##  gap> SingerAlg.BasisOfPMRoots( data, [ 3 ], 2, 1 );
##  [ 2, 3, 6, 7, 8, 9 ]
##  ]]></Example>
##  </Description>
##  </ManSection>
##  <#/GAPDoc>
##
SingerAlg.BasisOfPMRoots:= function( data, I, p, m )
    local q, pm, quo;

    q:= data.parameters[1];
    pm:= p^m;
    quo:= QuoInt( q-1, pm );

    # The first condition means that B_k^{p^m} is zero
    # and hence inside the given subspace.
    # For the second condition, we may assume that B_k^{p^m}
    # is nonzero and hence equal to B_{p^m (k-1)+1}.
    return Filtered( [ 1 .. Length( data.monomials ) ],
                     k ->    Maximum( data.monomials[k] ) > quo
                          or pm * (k-1) + 1 in I );
    end;

#############################################################################
##
#F  SingerAlg.BasisOfPC( <data>, <I>, <J> )
##
##  <#GAPDoc Label="SingerAlg.BasisOfPC">
##  <ManSection>
##  <Func Name="SingerAlg.BasisOfPC" Arg='data, I, J'/>
##
##  <Description>
##  Let <A>data</A> be the record returned by
##  <Ref Oper="LoewyStructureInfoGAP" Label="for parameters"/>
##  that describes a Singer algebra <M>A = A[q,z]</M>,
##  let <M>B = B(A)</M>,
##  let <A>I</A> and <M>J</M> be subsets of <M>\{ 1, 2, \ldots, z+1 \}</M>,
##  describing subspaces <M>U</M> and <M>V</M> of <M>A</M> with bases
##  <M>( B_i; i \in I )</M> and <M>( B_i; i \in J )</M>,
##  respectively.
##  <Ref Func="SingerAlg.BasisOfPC"/> returns the subset <M>K</M>
##  of <M>\{ 2, 3, \ldots, z+1 \}</M> such that
##  <M>( B_i; i \in K )</M> is a basis of the subspace
##  <M>\{ x \in J(A); x \cdot U \subseteq V \}</M> of <M>J(A)</M>.
##  <P/>
##  (The perhaps strange name <Q>BasisOfPC</Q> was chosen because the
##  result contains the indices of those basis vectors such that the
##  <E>P</E>roduct with the space <M>U</M> is <E>C</E>ontained in the
##  space <M>V</M>.)
##  <P/>
##  <Example><![CDATA[
##  gap> data:= LoewyStructureInfoGAP( 23, 585 );;
##  gap> soc:= SingerAlg.BasesOfSocleSeries( data );;
##  gap> rad:= SingerAlg.BasesOfRadicalSeries( data );;
##  gap> I1:= SingerAlg.BasisOfPC( data, soc[3], rad[3] );;
##  gap> Length( I1 );
##  581
##  gap> data:= LoewyStructureInfoGAP( 212, 585 );;
##  gap> soc:= SingerAlg.BasesOfSocleSeries( data );;
##  gap> rad:= SingerAlg.BasesOfRadicalSeries( data );;
##  gap> I2:= SingerAlg.BasisOfPC( data, soc[3], rad[3] );;
##  gap> Length( I2 );
##  545
##  ]]></Example>
##  </Description>
##  </ManSection>
##  <#/GAPDoc>
##
SingerAlg.BasisOfPC:= { data, I, J } -> Filtered(
    # Note that we are interested in subspaces of the radical.
    [ 2 .. Length( data.monomials ) ],
    k -> IsSubset( J, SingerAlg.BasisOfProductSpace( data, [ k ], I ) ) );

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

max_stars_repo_name: oscar-system/SingerAlg

max_stars_repo_licenses: Naumen Condor-1.1 MS-PL

max_issues_repo_name: oscar-system/SingerAlg

max_issues_repo_licenses: Naumen Condor-1.1 MS-PL

max_forks_repo_name: oscar-system/SingerAlg

max_forks_repo_licenses: Naumen Condor-1.1 MS-PL


text:
#
# Read("~/Workspace/groupsSB/epi/G2/handle.gi");
#

type:="G";
rank:=2;
nr_pos_roots:=6;

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

#
# c_2=0 c_3=0
#
#
#
handle23:=function()
	local ua,ub,uab,uabb,u,e,nn,vals,rels,i,v;
	ua:=root_group(1,1);
	uabb:=root_group(4,1);
	u:=ua*uabb;
	id_mat:=u^0;
	e:=u-id_mat;

	nn:=List([1..2*nr_pos_roots+rank],i->xvars[i]);

	vals:=[];
	Append(vals,[[[xvars[7]],[Zero(APR)]]]);
	Append(vals,[[[xvars[8]],[Zero(APR)]]]);
	Append(vals,[[[xvars[5]],[Zero(APR)]]]);
	Append(vals,[[[xvars[10]],[Zero(APR)]]]); # for p>2
	Append(vals,[[[xvars[9]],[Zero(APR)]]]); # for p>2
	Append(vals,[[[xvars[2]],[-xvars[6]]]]);
	# ###
	# Append(vals,[[[xvars[1]],[Zero(APR)]]]);
	# Append(vals,[[[xvars[3]],[Zero(APR)]]]);
	# Append(vals,[[[xvars[4]],[Zero(APR)]]]);
	# Append(vals,[[[xvars[8]],[Zero(APR)]]]);
	# Append(vals,[[[xvars[9]],[Zero(APR)]]]);

	#
	# atentie la ordinea transpunerilor
	#
	rels:=nn*TransposedMat(u)-nn;
	rels:=evaluate_rels(rels,vals);
	nn:=evaluate_rels(nn,vals);

	return [u,nn,rels,vals];
end;

ext: gi

max_stars_repo_name: iuliansimion/groupsSB

max_issues_repo_name: iuliansimion/groupsSB

max_forks_repo_name: iuliansimion/groupsSB


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

########################################
# These functions work on formula level
########################################

MarkBlock := function(sums) sums.isBlock := true; return sums; end;

IsMarkedBlock := o -> IsRec(o) and IsBound(o.isBlock) and o.isBlock;

DoNotMark := s -> (IsBound(s.doNotMarkBB) and s.doNotMarkBB = true) or
                  (IsIterative(s) and IsSymbolic(s.domain));

IsBlockTransitive := o -> IsBound(o.isBlockTransitive) and o.isBlockTransitive;

Compose.isBlockTransitive := true;
SUM.isBlockTransitive := true;
SUMAcc.isBlockTransitive := true;
Data.isBlockTransitive := true;

markBlocks := function(bk, sums, criterion_func)
    local doNotMark, c;
    if not IsSPL(sums) then return sums; fi;

    #Error(">>>>>>>------", ObjId(sums), "------>>>>>>>>>");
    doNotMark := false;

    for c in sums.rChildren() do
        markBlocks(bk, c, criterion_func);
        # if any of our children cannot be marked BB, we cannot be marked BB
        if IsSPL(c) and DoNotMark(c) then
          doNotMark := true;

          sums.doNotMarkBB := true; # To communicate this up to parent.
          # This is a slight HACK. Maybe we should mark these with a different
          # label (like sums.childIsNoMarkBB) and check for this, so we don't confuse
          # these with preset .doNotMarkBB

          #Error("Not marking: ", sums.__name__, " because of child: ", c.__name__);
        fi;
    od;

    if DoNotMark(sums) then
      doNotMark := true;
      #Error("Not marking ", sums.__name__, " because of doNotMarkBB");
    fi;

    #Error("<<<<<<<------", ObjId(sums), "------<<<<<<<<<");

    if not doNotMark then
        if IsBlockTransitive(ObjId(sums)) then
            if ForAll(sums.children(), IsMarkedBlock) then
               MarkBlock(sums);
            fi;
        elif criterion_func(sums, bk) then
            MarkBlock(sums);
        fi;
    fi;

    return sums;

end;

standard_criterion_func := (sums, bk) -> not IsBound(sums._perm)
    and not (IsIterative(sums) and IsSymbolic(sums.domain))
    and not ForAll(Flat(sums.dims()), IsSymbolic)
    and     ForAny(Flat(sums.dims()), d -> d<=bk);

MarkBlocksDimSums := (bk,sums) -> SubstBottomUp(sums, @,
    e -> Cond(not IsSPL(e) or (IsBound(e.doNotMarkBB) and e.doNotMarkBB), e,
          not IsBound(e._perm) and not (IsIterative(e) and IsSymbolic(e.domain)) and
             (((not IsSymbolic(Cols(e))) and Cols(e)<=bk) or # finish this
              ((not IsSymbolic(Rows(e))) and Rows(e)<=bk)),    MarkBlock(e),
          IsBlockTransitive(ObjId(e)) and ForAll(e.children(), IsMarkedBlock), MarkBlock(e),
          e));

MarkBlocksAreaSums := (bk,sums) -> SubstBottomUpRules(sums, rec(
    blockAreaSums := Rule(@,
        e -> Cond(not IsSPL(e) or not IsBound(e.area),    e,
              not IsSymbolic(e.area()) and e.area() <= bk,  MarkBlock(e),
              IsBlockTransitive(ObjId(e)) and ForAll(e.children(), IsMarkedBlock), MarkBlock(e),
              e))));

MarkBlocksOps := (bk,sums) -> SubstBottomUp(sums, @,
    e -> Cond(not IsSPL(e), e,
          not IsBound(e._perm) and e.numops()<=bk,    MarkBlock(e),e));

##
#MarkBlocksSums := MarkBlocksDimSums;
MarkBlocksSums := (bk, sums) -> markBlocks(bk, sums, standard_criterion_func);
##

Class(PostProcessBBs, RuleSet);
RewriteRules(PostProcessBBs, rec(
    mergeBBsinsideCompose := ARule(Compose, [@(1, BB), @(2, BB)], e -> [BB(Compose(@(1).val.rChildren()[1], @(2).val.rChildren()[1]))]),
    kickOutNestedBBs := Rule( @@(1, BB, (x, cx) -> cx.isInside(BB)), e -> e.child(1)),
));

#F BlockSums(<bksize>, <sums>)
#F     Uses MarkBlocksSums with given block size <bksize> to partition the
#F     Sigma-SPL formula <sums> into unrolled blocks.
#F
#F     MarkBlocksSums = MarkBlocksDimSums || MarkBlocksAreaSums.
#F     ..DimSums partitions by the smallest dimensions
#F     ..AreaSums partitions by the total dense area of a (potentially sparse) matrix
#F
BlockSums := function(bk, sums)
    sums := MarkBlocksSums(bk, sums);
    sums := SubstTopDownNR(sums, @.cond(IsMarkedBlock), e->When(ObjId(e)=BB, e, BB(e)));
    return sums;
end;

BlockSumsOpts := function(sums, opts)
    if IsBound(opts.markBlock) then
        sums := opts.markBlock(opts.globalUnrolling, sums);
    else
        sums := MarkBlocksSums(opts.globalUnrolling, sums);
    fi;
    sums := SubstTopDownNR(sums, @.cond(IsMarkedBlock), e->When(ObjId(e)=BB, e, BB(e)));
    return PostProcessBBs(sums);
end;

# Alternative implementation can use rules
#RuleSet(RulesBlockSums);
#Rule(@.cond(IsMarkedBlock), e -> BB(e));
#ARule(Compose, [..., @(1, BB), @(2, BB), ...], e -> BB(Compose(@(1).val, @(2).val));
#ARule(SUM,     [..., @(1, BB), @(2, BB), ...], e -> BB(SUM    (@(1).val, @(2).val)));
#ARule(SUMAcc,  [..., @(1, BB), @(2, BB), ...], e -> BB(SUMAcc (@(1).val, @(2).val)));

####################################
# These functions work on code level
####################################

#F MarkForUnrolling(<bbnum>, <code>, [<opts>])
MarkForUnrolling := (arg) -> Cond(
    # if has options and no SimFlexKernelFlag then go without kern()
    Length(arg)=3 and not IsBound(arg[3].SimFlexKernelFlag),
       unroll_cmd(arg[2]),
    # else
       kern(arg[1], unroll_cmd(arg[2])));

IsMarkedForUnrolling := code -> ObjId(code) = unroll_cmd;

#F BlockCollect(<code>)
#F
#F Returns a list of blocks marked for unrolling (field doUnroll:=true must be set)
#F
BlockCollect := c -> Collect(c, @.cond(IsMarkedForUnrolling));

Declare(_BlockUnroll);

LAST_BLOCKUNROLL := 0;
LAST_BLOCKUNROLLOPTS := 0;
#F BlockUnroll(<code>, <opts>)
#F
#F Generates unrolled code for marked blocks (field doUnroll:=true must be set)
#F
BlockUnroll := function(c, opts)
    LAST_BLOCKUNROLL := Copy(c);
    LAST_BLOCKUNROLLOPTS := opts;
    return _BlockUnroll(c, opts, 1);
end;

_BlockUnroll := function(c, opts, ldepth)
    local i, last, ch, newch, tocompile;

    if IsMarkedForUnrolling(c) then return Compile(c.cmds[1], opts); # strip unroll_cmd
    elif not IsCommand(c) then return c;
    else
    # a chain with some subcommands to be unrolled
    if ObjId(c) = chain and ForAny(c.cmds, IsMarkedForUnrolling) then
        newch := [];
        i := 1;
        # group each sequence of marked subcommands and unroll together
        # to facilitate copy propagation in-between
        while i <= Length(c.cmds) do
            ch := c.cmds[i];
           # recurse on blocks that aren't marked
            if not IsMarkedForUnrolling(ch) then
                Add(newch, _BlockUnroll(ch, opts, ldepth));
                # otherwise find maximal sequence of marked blocks and unroll them together
            else
                last := i;
                tocompile := [ch.cmds[1]];
                while last < Length(c.cmds) and IsMarkedForUnrolling(c.cmds[last+1]) do
                    last := last+1;
                    ch := c.cmds[last];
                    Add(tocompile, ch.cmds[1]); # strip unroll_cmd
                od;
                Add(newch, Compile(chain(tocompile), opts));
                i := last;
            fi;
            i := i+1;
        od;
        return chain(newch);
    # for an unmarked command (that is not a chain) simply recurse on subcommands
    else
        # mark loop variables with their loop nesting depth
        if ObjId(c) in [loop, loopn] then
            c.var.ldepth := ldepth;
            ldepth := ldepth + 1;
        fi;

        i := 1;
        for ch in c.rChildren() do
            c.rSetChild(i, _BlockUnroll(ch, opts, ldepth));
            i := i+1;
        od;
        return c;
    fi;
    fi;
end;

# enumerate the blocks. useful for when you want to turn blocks on and off at
# the code level.

EnumBlocks := function(sums, opts)
    local blocks, numbered, notnumb, max, i;
    blocks := Collect(sums, BB);

    numbered := Filtered(blocks, e -> IsBound(e.bbnum));
    notnumb := Filtered(blocks, e -> not IsBound(e.bbnum));

    if numbered <> [] then
        max := Maximum(List(numbered, e -> e.bbnum));
    else
        max := 0;
    fi;

    for i in notnumb do
        max := max + 1;
        i.bbnum := max;
    od;

    return sums;
end;


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

################################################################################
##
#F DwG_simples(<group>,<int>,<function>) . . . . . compute simples in D^w(G)
##
InstallGlobalFunction( DwG_simples, function(G,n,w)
  local x,Zx,wx,CC,SimpleObjects,projective_characters,chi;
  CC:=List(ConjugacyClasses(G),Representative);;
  SimpleObjects:=[];
  for x in CC do
    Zx:=Centralizer(G,x);
    wx:=function(g,h) return w(g,h,x)*w(x,g,h)*(w(g,x,h)^-1); end;;
    projective_characters:=projective_reps(Zx,wx,n);
    if not(Sum(List(projective_characters,chi->chi[1]^2))=Order(Zx)) then
      Display("ERROR");
    fi;
    for chi in projective_characters do
      Append(SimpleObjects,[[x,chi,ConjugacyClass(G,x),wx,Zx]]);
    od;
  od;
  return SimpleObjects;
end );

################################################################################
##
#F DwG_S_and_T(<list>,<group>) . . . . . . . . . .compute S and T for D^w(G)
##
InstallGlobalFunction(DwG_S_and_T, function(Simples,G)
  local r,S,T,i,simple1,x,chi,Kx,betax,Zx,x_idx,Zx_Transversal,j,simple2,y,chi_prime,Ky,betay,Zy,Zy_Transversal,a,g,Zg,b,g_prime,h_prime,h,h_idx,h_prime_idx;

  r:=Length(Simples);
  S:=NullMat(r,r);;
  T:=List([1..r],i->0);
  for i in [1..r] do
    simple1:=Simples[i];
    x:=simple1[1];
    chi:=simple1[2];
    Kx:=simple1[3];
    betax:=simple1[4];
    Zx:=simple1[5];
    #x_idx:=Position(GSet,x);
    x_idx:=Position(AsSet(Zx),x);
    T[i]:=chi[x_idx]/chi[Position(AsSet(Zx),Identity(Zx))];
    Zx_Transversal:=RightTransversal(G,Zx);
    for j in [i..r] do
      simple2:=Simples[j];
      y:=simple2[1];
      chi_prime:=simple2[2];
      Ky:=simple2[3];
      betay:=simple2[4];
      Zy:=simple2[5];
      #y_idx:=Position(GSet,y);
      Zy_Transversal:=RightTransversal(G,Zy);
      for a in Zx_Transversal do
        g:=(a^-1)*x*a;
        if not(g in Kx) then
          continue;
        fi;
        Zg:=Centralizer(G,g);
        for b in Zy_Transversal do
          g_prime:=(b^-1)*y*b;
          if not((g_prime in Ky) and (g_prime in Zg)) then
            continue;
          fi;
          h_prime:=b*g*(b^-1);
          h:=a*g_prime*(a^-1);
          h_idx:=Position(AsSet(Zx),h);
          h_prime_idx:=Position(AsSet(Zy),h_prime);
          S[i][j]:=S[i][j]+betax(a,g_prime)*betax(a*g_prime,a^-1)*betay(b,g)*betay(b*g,b^-1)*(betax(a,a^-1)^-1)*(betay(b,b^-1)^-1)*chi[h_idx]*chi_prime[h_prime_idx];
        od;
      od;
      S[j][i]:=S[i][j];
    od;
  od;

  return [S,T];
end);

################################################################################
##
#F DwG_data(<group>,<function>,<int>) . . . . . compute data of D^w(G)
InstallGlobalFunction(DwG_data,function(G,w,n)
  local Simples, S_and_T, pseudounitary_data, S, T, N, FPdimC, FPdim,sorted_data;
  Simples := DwG_simples(G,n,w);;
  S_and_T := DwG_S_and_T(Simples,G);;
  S := S_and_T[1];
  T := S_and_T[2];
  pseudounitary_data := pseudounitary_sort(S,T,Simples);
  Simples := pseudounitary_data[1];
  S := pseudounitary_data[2];
  T := pseudounitary_data[3];
  FPdim := pseudounitary_data[4];
  FPdimC := pseudounitary_data[5];
  sorted_data := dimension_sort(S,T,Simples,FPdim);
  Simples := sorted_data[1];
  S := sorted_data[2];
  T := sorted_data[3];
  FPdim := sorted_data[4];
  N:=compute_fusion_rules(S,FPdimC);;
  return [Simples, S, T, N, FPdim, FPdimC];;
end);

#E DwG_data.gi . . . . . . . . . . . . . . . . . . . . . . . . . . . . ends here

max_stars_repo_name: pnnl/GroupTheoretical

max_issues_repo_name: pnnl/GroupTheoretical

max_forks_repo_name: pnnl/GroupTheoretical


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

Class(UTwid, DiagFunc, rec(
    def := (N, n, k, a, b, i) -> rec(size:=n),
    lambda := self >> let(N:=self.params[1], n:=self.size, k:=self.params[3],
	a := self.params[4], b := self.params[5], i := self.params[6],
	na := Numerator(a),   nb := Numerator(b),
	da := Denominator(a), db := Denominator(b),
	w := E(N*da*db)^k, j := Ind(n),
	Lambda(j, cond(eq(j,n/2), omega(da*db,k), omega(N*da*db, k*(i*da+na)*(j*db+nb)))))
));
Class(UTwid1, DiagFunc, rec(
    def := (N, n, k, a, b, i) -> rec(size:=n),
    lambda := self >> let(N:=self.params[1], n:=self.size, k:=self.params[3],
	a := self.params[4], b := self.params[5], i := self.params[6],
	na := Numerator(a),   nb := Numerator(b),
	da := Denominator(a), db := Denominator(b),
	w := E(N*da*db)^k, j := Ind(n),
	Lambda(j, cond(eq(j,n/2), 1, omega(N*da*db, k*(i*da+na)*(j*db+nb)))))
));

Class(UDFT_NonTerm, DFT_NonTerm, rec(
    isReal := False,
));

Class(UDFT, UDFT_NonTerm, rec(
   abbrevs := [
    (n)    -> Checked(IsPosIntSym(n), AnySyms(n) or (n mod 2=0), [n, 1]),
    (n,k)  -> Checked(IsPosIntSym(n), AnySyms(n) or (n mod 2=0),
	              IsIntSym(k), AnySyms(n,k) or Gcd(n,k)=1,
		      [n, When(AnySyms(n,k), k, k mod n)]) ],
   omega4pow := (r,c)->4*r*c,
   terminate := self >> let(
	n:=self.params[1], k:=self.params[2], m:=MatSPL(DFT(2))^-1,
	DirectSum(Mat(m), I(n-2))^L(n,n/2) *
	DFT(n,k).terminate()),

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

Class(UDFT2, UDFT_NonTerm, rec(
   abbrevs := [
    (n)    -> Checked(IsPosIntSym(n), AnySyms(n) or (n mod 2=0), [n, 1]),
    (n,k)  -> Checked(IsPosIntSym(n), AnySyms(n) or (n mod 2=0),
	              IsIntSym(k), AnySyms(n,k) or Gcd(n,k)=1,
		      [n, When(AnySyms(n,k), k, k mod (2*n))]) ],
    omega4pow := (r,c)->2*r*(2*c+1),
    terminate := self >> let(
	n:=self.params[1], k:=self.params[2], m:=MatSPL(DFT2(2,k))^-1,
	DirectSum(Mat(m), I(n-2))^L(n,n/2) *
	DFT2(n,k).terminate()),
));

Class(UUDFT, UDFT_NonTerm, rec(
   abbrevs := [
    (n)    -> Checked(IsPosIntSym(n), AnySyms(n) or (n mod 4=0), [n, 1]),
    (n,k)  -> Checked(IsPosIntSym(n), AnySyms(n) or (n mod 4=0),
	              IsIntSym(k), AnySyms(n,k) or Gcd(n,k)=1,
		      [n, When(AnySyms(n,k), k, k mod n)]) ],
    omega4pow := (r,c)->4*r*c,
    terminate := self >> let(
	n:=self.params[1], k:=self.params[2], m:=MatSPL(DFT(2))^-1,
	DirectSum(Mat(m), I((n-4)/2), Mat(m), I((n-4)/2))^L(n,n/2) *
	DFT(n,k).terminate()),
));
UUDFT1 := UUDFT;

Class(UUDFT2, UDFT_NonTerm, rec(
   abbrevs := [
    (n)    -> Checked(IsPosIntSym(n), AnySyms(n) or (n mod 4=0), [n, 1]),
    (n,k)  -> Checked(IsPosIntSym(n), AnySyms(n) or (n mod 4=0),
	              IsIntSym(k), AnySyms(n,k) or Gcd(n,k)=1,
		      [n, When(AnySyms(n,k), k, k mod (2*n))])  ],
    omega4pow := (r,c)->2*r*(2*c+1),
    terminate := self >> let(
	n:=self.params[1], k:=self.params[2], m:=MatSPL(DFT2(2,k))^-1,
	DirectSum(Mat(m), I((n-2)/2), Mat(m), I((n-2)/2))^L(n,n/2) *
	DFT2(n,k).terminate()),
));


text:
Filtered([2008 .. 2121], y -> WeekDay([25, 12, y]) = "Sun");
# [ 2011, 2016, 2022, 2033, 2039, 2044, 2050, 2061, 2067, 2072, 2078, 2089, 2095, 2101, 2107, 2112, 2118 ]

# A possible implementation of WeekDayAlt

days := ["Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun"];;

WeekDayAlt := function(args)
   local d, m, y, k;
   d := args[1];
   m := args[2];
   y := args[3];
   if m < 3 then
      m := m + 12;
      y := y - 1;
   fi;
   k := 1 + RemInt(d + QuoInt((m + 1)*26, 10) + y + QuoInt(y, 4)
          + 6*QuoInt(y, 100) + QuoInt(y, 400) + 5, 7);
   return days[k];
end;

Filtered([2008 .. 2121], y -> WeekDayAlt([25, 12, y]) = "Sun");
# [ 2011, 2016, 2022, 2033, 2039, 2044, 2050, 2061, 2067, 2072, 2078, 2089, 2095, 2101, 2107, 2112, 2118 ]

ext: gap

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


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

###############################################################################
##
#R  IsTreeAutomorphismRep
##
## XXX remove it, use IsTreeHomomorphismRep
DeclareRepresentation("IsTreeAutomorphismRep",
                      IsComponentObjectRep and IsAttributeStoringRep,
                      ["states", "perm", "deg"]);

###############################################################################
##
#R  IsTreeAutomorphismFamilyRep
##
DeclareRepresentation("IsTreeAutomorphismFamilyRep",
                      IsComponentObjectRep and IsAttributeStoringRep,
                      ["spher_index"]);

###############################################################################
##
##  AG_CreatedTreeAutomorphismFamilies
##
BindGlobal("AG_CreatedTreeAutomorphismFamilies", []);

###############################################################################
##
#M  TreeAutomorphismFamily(<sph_ind>)
##
InstallMethod(TreeAutomorphismFamily, [IsRecord],
function(sph_ind)
  local p, red_ind, fam;
  red_ind := AG_ReducedSphericalIndex(sph_ind);
  for p in AG_CreatedTreeAutomorphismFamilies do
    if p[1] = red_ind then
      return p[2]; fi;
  od;
  fam := NewFamily( Concatenation("Automorphisms of ",
                      red_ind.start, ", (", red_ind.period, ")-tree"),
                    IsTreeAutomorphism,
                    IsTreeAutomorphism,
                    IsTreeAutomorphismFamily and IsTreeAutomorphismFamilyRep);
  fam!.spher_index := red_ind;
  AddSet(AG_CreatedTreeAutomorphismFamilies, [red_ind, fam]);
  return fam;
end);

BindGlobal("AG_TreeAutomorphism",
function(list_states, permutation)
  local top_deg, bot_deg, ind, fam, s, Orbs, orb;

  if Length(list_states)=0 then
    Error("The list of states can not be empty");
  fi;

  top_deg := Length(list_states);
  if not IsOne(permutation) and
      top_deg < Maximum(MovedPoints(permutation)) then
    Error("The root permutation ", permutation, " must move only points from 1 to the degree ", top_deg, " of the tree");
  fi;

  Orbs := OrbitsPerms([permutation], [1..Length(list_states)]);

  for orb in Orbs do
    for s in [2..Length(orb)] do
      if list_states[orb[s]]!.deg<>list_states[orb[1]]!.deg then
        Error("Sections in one orbit are acting on different trees");
      fi;
    od;
  od;

  bot_deg := DegreeOfTree(list_states[1]);
  ind := rec(start := [top_deg], period := [bot_deg]);
  fam := TreeAutomorphismFamily(ind);

  return Objectify( NewType(fam, IsTreeAutomorphism and IsTreeAutomorphismRep),
                    rec(states := list_states,
                        perm := permutation,
                        deg := top_deg));
end);

###############################################################################
##
#M  TreeAutomorphism(<states_list>, <perm>)
##
InstallMethod(TreeAutomorphism, "for [IsList, IsPerm]", [IsList, IsPerm],
function(states, perm)
  local autom, nstates, s;

  autom := fail;

  for s in states do
    if IsTreeAutomorphism(s) then
      autom := s;
    elif not IsOne(s) then
      Error("Invalid state `", s, "'");
    fi;
  od;

  if autom = fail then
    # XXX homogenous tree, stupid!
    Error("Can't create an automaton with all trivial states ",
          "without information about the tree");
  fi;

  nstates := List(states, function(x)
                            if IsInt(x) and IsOne(x) then
                              return One(autom);
                            else
                              return x;
                            fi;
                          end);

  return AG_TreeAutomorphism(nstates, perm);
end);

###############################################################################
##
#M  TreeAutomorphism(<state_1>, <state_2>, ..., <state_n>, <perm>)
##
InstallMethod(TreeAutomorphism, [IsObject, IsObject, IsPerm],
  function(a1, a2, perm) return TreeAutomorphism([a1, a2], perm); end);
InstallMethod(TreeAutomorphism, [IsObject, IsObject, IsObject, IsPerm],
  function(a1, a2, a3, perm) return TreeAutomorphism([a1, a2, a3], perm); end);
InstallMethod(TreeAutomorphism, [IsObject, IsObject, IsObject, IsObject, IsPerm],
  function(a1, a2, a3, a4, perm) return TreeAutomorphism([a1, a2, a3, a4], perm); end);
InstallMethod(TreeAutomorphism, [IsObject, IsObject, IsObject, IsObject, IsObject, IsPerm],
  function(a1, a2, a3, a4, a5, perm) return TreeAutomorphism([a1, a2, a3, a4, a5], perm); end);
#InstallMethod(TreeAutomorphism, [IsObject, IsObject, IsObject, IsObject, IsObject, IsObject, IsPerm],
#  function(a1, a2, a3, a4, a5, a6, perm) return TreeAutomorphism([a1, a2, a3, a4, a5, a6], perm); end);

###############################################################################
##
#M  ViewObj(<a>)
##
InstallMethod(ViewObj, [IsTreeAutomorphism],
function (a)
    local deg, printword, i;

    deg := AG_TopDegreeInSphericalIndex(FamilyObj(a)!.spher_index);
    Print("(");
    for i in [1..deg] do
        View(a!.states[i]);
        if i <> deg then Print(", "); fi;
    od;
    Print(")");
    if not IsOne(a!.perm) then Print(a!.perm); fi;
end);

###############################################################################
##
#M  PrintObj(<a>)
##
InstallMethod(PrintObj, "for [IsTreeAutomorphism and IsTreeAutomorphismRep]",
                             [IsTreeAutomorphism and IsTreeAutomorphismRep],
function (a)
    local deg, printword, i;

    deg := AG_TopDegreeInSphericalIndex(FamilyObj(a)!.spher_index);
    Print("(");
    for i in [1..deg] do
      if IsAutom(a!.states[i]) then
        View(a!.states[i]);
      else
        Print(a!.states[i]);
      fi;
      if i <> deg then Print(", "); fi;
    od;
    Print(")");
    if not IsOne(a!.perm) then Print(a!.perm); fi;
end);

###############################################################################
##
#M  String(<a>)
##
InstallMethod(String, "for [IsTreeAutomorphism]", [IsTreeAutomorphism],
function (a)
    local deg, printword, i, perm, states, str;

    states := Sections(a);
    deg := Length(states);
    perm := PermOnLevel(a, 1);
    str:= "(";

    for i in [1..deg] do
        Append(str, String(states[i]));
        if i <> deg then Append(str, ", "); fi;
    od;
    Append(str, ")");
    if not IsOne(perm) then
      Append(str, AG_TransformationString(perm));
    fi;
    return str;
end);

###############################################################################
##
#M  SphericalIndex(<a>)
##
InstallMethod(SphericalIndex, [IsTreeAutomorphism and IsTreeAutomorphismRep],
function(a)
  return FamilyObj(a)!.spher_index;
end);

###############################################################################
##
#M  Perm(<a>)
##
InstallMethod(Perm, "for [IsTreeAutomorphism and IsTreeAutomorphismRep]",
              [IsTreeAutomorphism and IsTreeAutomorphismRep],
function(a)
  return a!.perm;
end);

###############################################################################
##
#M  Perm (<a>, <k>)
##
InstallMethod(Perm, "for [IsTreeAutomorphism, IsPosInt]",
              [IsTreeAutomorphism, IsPosInt],
function(a, k)
  return PermOnLevel(a, k);
end);

###############################################################################
##
#M  PermOnLevel (<a>, <k>)
##
InstallMethod(PermOnLevelOp, "for [IsTreeAutomorphism, IsPosInt]",
              [IsTreeAutomorphism, IsPosInt],
function(a, k)
  local states, top, first_level, i, j, d1, d2, permuted, p;

  if k = 1 then
    return Perm(a);
  fi;

  # XXX test this function

  # TODO: it goes through all vertices of the second level, it may be
  # unnecessary for sparse actions
  d1 := a!.deg;
  d2 := 1;
  for i in [2 .. k] do
    d2 := d2 * DegreeOfLevel(a, i);
  od;
  states := Sections(a);
  top := Perm(a);
  first_level := List(states, s -> PermOnLevel(s, k-1));
  permuted := [];
  for i in [1..d1] do
    for j in [1..d2] do
      permuted[d2*(i-1) + j] := d2*(i^top - 1) + j^first_level[i];
    od;
  od;

  p := PermList(permuted);
  if p<>fail then
    return p;
  else
    Error("An element ",a," does not induce a permutation on the level ",k);
  fi;
end);

###############################################################################
##
#M  k ^ a
##
InstallMethod(\^, "for [IsPosInt, IsTreeAutomorphism]",
                   [IsPosInt, IsTreeAutomorphism],
function(k, a)
    return k ^ Perm(a);
end);

###############################################################################
##
#M  seq ^ a
##
InstallMethod(\^, "for [IsList, IsTreeAutomorphism]",
                   [IsList, IsTreeAutomorphism],
function(seq, a)
  if Length(seq) = 0 then return []; fi;
  if Length(seq) = 1 then return [seq[1]^Perm(a)]; fi;
  return Concatenation([seq[1]^Perm(a)], seq{[2..Length(seq)]}^Section(a, seq[1]));
end);

###############################################################################
##
#M  FixesLevel(<a>, <k>)
##
InstallMethod(FixesLevel, "for [IsTreeAutomorphism, IsPosInt]",
              [IsTreeAutomorphism, IsPosInt],
function(a, k)
  if HasIsSphericallyTransitive(a) then
    if IsSphericallyTransitive(a) then
      return false; fi; fi;

  if IsOne(PermOnLevel(a, k)) then
    Info(InfoAutomGrp, 3, "IsSphericallyTransitive(a): false");
    Info(InfoAutomGrp, 3, "  a is not transitive on level", k);
    Info(InfoAutomGrp, 3, "  a = ", a);
    SetIsSphericallyTransitive(a, false);
    return true;
  else
    return false;
  fi;
end);

###############################################################################
##
#M  FixesVertex(<a>, <v>)
##
InstallMethod(FixesVertex,  "for [IsTreeAutomorphism, IsObject]",
                   [IsTreeAutomorphism, IsObject],
function(a, v)
  if HasIsSphericallyTransitive(a) then
    if IsSphericallyTransitive(a) then
      Info(InfoAutomGrp, 3, "FixesVertex(a, v): false");
      Info(InfoAutomGrp, 3, "  IsSphericallyTransitive(a)");
      Info(InfoAutomGrp, 3, "  a = ", a);
      return false;
    fi;
  fi;

  if v^a = v then
    Info(InfoAutomGrp, 3, "IsSphericallyTransitive(a): false");
    Info(InfoAutomGrp, 3, "  a fixes vertex ", v);
    Info(InfoAutomGrp, 3, "  a = ", a);
    SetIsSphericallyTransitive(a, false);
    return true;
  else
    return false;
  fi;
end);

###############################################################################
##
#M  IsSphericallyTransitive (<a>)
##
InstallMethod(IsSphericallyTransitive,
              "for [IsTreeAutomorphism and IsActingOnBinaryTree]",
              [IsTreeAutomorphism and IsActingOnBinaryTree],
function(a)
  local ab;
  Info(InfoAutomGrp, 4, "IsSphericallyTransitive(a): using AbelImage");
  Info(InfoAutomGrp, 4, "  a = ", a);
  ab := AbelImage(a);
  return ab = One(ab)/(One(ab)+IndeterminateOfUnivariateRationalFunction(ab));
end);

RedispatchOnCondition(IsSphericallyTransitive, true, [IsTreeAutomorphism],
                      [IsTreeAutomorphism and IsActingOnBinaryTree], 0);

InstallMethod(IsSphericallyTransitive, "for [IsTreeAutomorphism]",
              [IsTreeAutomorphism],
function(a)
  if not IsTransitive(Group(PermOnLevel(a, 1)), [1..Degree(a)]) then
    Info(InfoAutomGrp, 4, "IsSphericallyTransitive(a): false");
    Info(InfoAutomGrp, 4, "  PermOnLevel(a, 1) isn't transitive");
    Info(InfoAutomGrp, 4, "  a = ", a);
    return false;
  fi;
  TryNextMethod();
end);

InstallMethod(IsSphericallyTransitive,
              "for [IsTreeAutomorphism and HasOrder]",
              [IsTreeAutomorphism and HasOrder],
function(a)
  if Order(a) < infinity then
    Info(InfoAutomGrp, 4, "IsSphericallyTransitive(a): false");
    Info(InfoAutomGrp, 4, "  Order(a) < infinity");
    Info(InfoAutomGrp, 4, "  a = ", a);
    return false;
  fi;
  TryNextMethod();
end);

###############################################################################
##
#M  Order (<a>)
##
InstallImmediateMethod(Order, IsTreeAutomorphism and HasIsSphericallyTransitive, 0,
function(a)
  if IsSphericallyTransitive(a) then return infinity; fi;
  TryNextMethod();
end);

###############################################################################
##
#M  Order(<a>)
##
InstallMethod(Order, "for [IsTreeAutomorphism]", [IsTreeAutomorphism],
function(a)
  local i, perm, stab, stab_order, ord, exp, states;

  perm := Perm(a);
  if IsOne(perm) then
    exp := 1;
    stab := a;
  else
    exp := Order(perm);
    stab := a^exp;
  fi;

  if IsOne(stab) then
    return exp;
  fi;

  states := Sections(stab);
  stab_order := 1;

  for i in [1..Length(states)] do
    ord := Order(states[i]);
    if ord = infinity then
      return infinity;
    else
      stab_order := Lcm(stab_order, ord);
    fi;
  od;

  return exp * stab_order;
end);

###############################################################################
##
#M  Section(<a>, <k>)
##

InstallMethod(Section, [IsTreeAutomorphism and IsTreeAutomorphismRep, IsPosInt],
function(a, k)
  return a!.states[k];
end);

###############################################################################
##
#M  Sections(<a>)
##
InstallMethod(Sections, [IsTreeAutomorphism and IsTreeAutomorphismRep],
function(a)
  return a!.states;
end);

###############################################################################
##
#M  Decompose(<a>, <k>)
##
InstallMethod(Decompose, "for [IsTreeAutomorphism, IsPosInt]",
              [IsTreeAutomorphism, IsPosInt],
function(a, level)
  return TreeAutomorphism(Sections(a, level), PermOnLevel(a, level));
end);

###############################################################################
##
#M  Decompose(<a>)
##
InstallMethod(Decompose, "for [IsTreeAutomorphism]", [IsTreeAutomorphism],
function(a)
  return Decompose(a, 1);
end);

###############################################################################
##
#M  IsOne(<a>)
##
InstallMethod(IsOne, "for [IsTreeAutomorphism and IsTreeAutomorphismRep]",
                   [IsTreeAutomorphism and IsTreeAutomorphismRep],
function(a)
  local i;
  if a!.perm <> () then return false; fi;
  for i in [1..a!.deg] do
    if not IsOne(a!.states[i]) then return false; fi;
  od;
  return true;
end);

###############################################################################
##
#M  IsOne(<a>)
##
InstallMethod(IsOne, [IsTreeAutomorphism],
function(a)
  local i;
  if not IsOne(Perm(a)) then return false; fi;
  for i in [1..TopDegreeOfTree(a)] do
    if not IsOne(Section(a, i)) then return false; fi;
  od;
  return true;
end);

###############################################################################
##
#M  \=(<a1>, <a2>)
##
# TODO: can lead to infinite recursion
InstallMethod(\=, "for [IsTreeAutomorphism, IsTreeAutomorphism]", ReturnTrue,
              [IsTreeAutomorphism, IsTreeAutomorphism],
function(a1, a2)
  return Perm(a1) = Perm(a2) and Sections(a1) = Sections(a2);
end);

###############################################################################
##
#M  \<(<a1>, <a2>)
##
InstallMethod(\<, [IsTreeAutomorphism and IsTreeAutomorphismRep,
                   IsTreeAutomorphism and IsTreeAutomorphismRep],
function(a1, a2)
  return AG_TreeHomomorphismCmp(a1, a2) < 0;
end);

###############################################################################
##
#M  OneOp(<a>)
##
InstallMethod(OneOp, [IsTreeAutomorphism and IsTreeAutomorphismRep],
function(a)
  return Objectify( NewType(FamilyObj(a), IsTreeAutomorphism and IsTreeAutomorphismRep),
                    rec(states := List([1..a!.deg], i -> One(a!.states[1])),
                        perm := (),
                        deg := a!.deg));
end);

###############################################################################
##
#M  \*(<a1>, <a2>)
##
InstallMethod(\*, [IsTreeAutomorphism and IsTreeAutomorphismRep,
                   IsTreeAutomorphism and IsTreeAutomorphismRep],
function(a1, a2)
  local a;
  a := Objectify(NewType(FamilyObj(a1), IsTreeAutomorphism and IsTreeAutomorphismRep),
        rec(states := List([1..a1!.deg], i -> a1!.states[i] * a2!.states[i^(a1!.perm)]),
            perm := a1!.perm * a2!.perm,
            deg := a1!.deg));
  SetIsActingOnBinaryTree(a, IsActingOnBinaryTree(a1));
  SetIsActingOnRegularTree(a, IsActingOnRegularTree(a1));
  return a;
end);

###############################################################################
##
#M  \*(<a1>, <a2>)
##
InstallMethod(\*, [IsTreeAutomorphism, IsTreeAutomorphism],
function(a1, a2)
  local s1, s2, p1, p2, states, perm, d, a;
  s1 := Sections(a1); p1 := Perm(a1);
  s2 := Sections(a2); p2 := Perm(a2);
  states := List([1..Length(s1)], i -> s1[i] * s2[i^p1]);
  return TreeAutomorphism(states, p1*p2);
end);

###############################################################################
##
#M  \[\](<a1>, <a2>)
##
InstallOtherMethod(\[\], [IsTreeAutomorphism, IsPosInt],
function(a, k)
  return Section(a, k);
end);

###############################################################################
##
#M  InverseOp(<a>)
##
InstallMethod(InverseOp, "for [IsTreeAutomorphism and IsTreeAutomorphismRep]",
              [IsTreeAutomorphism and IsTreeAutomorphismRep],
function(a)
  local inv;
  inv := Objectify(NewType(FamilyObj(a), IsTreeAutomorphism and IsTreeAutomorphismRep),
            rec(states := List([1..a!.deg], i -> a!.states[i^(a!.perm^-1)]^-1),
                perm := a!.perm ^ -1,
                deg := a!.deg) );
  SetIsActingOnBinaryTree(inv, IsActingOnBinaryTree(a));
  SetIsActingOnRegularTree(inv, IsActingOnRegularTree(a));
  return inv;
end);

InstallMethod(InverseOp, "for [IsTreeAutomorphism]", [IsTreeAutomorphism],
function(a)
  local states, inv_states, perm;
  states := Sections(a);
  perm := Inverse(Perm(a));
  inv_states := List([1..Length(states)], i -> Inverse(states[i^perm]));
  return TreeAutomorphism(inv_states, perm);
end);

###############################################################################
##
#M  AbelImage(<a>)
##
##  XXX  Works for IsAutom or IsSelfSim only !!!!
##
InstallMethod(AbelImage, "for [IsTreeAutomorphism]",
              [IsTreeAutomorphism],
function(a)
  local abels, w, i;
  w := LetterRepAssocWord(Word(a));
  for i in [1..Length(w)] do
    if w[i] < 0 then w[i] := -w[i]+FamilyObj(a)!.numstates; fi;
  od;
  abels := AG_AbelImagesGenerators(FamilyObj(a));
  if not IsEmpty(w) then
    return Sum(List(w, x -> abels[x]));
  else
    return Zero(abels[1]);
  fi;
end);

#E

max_stars_repo_name: gap-packages/automgrp

max_issues_repo_name: gap-packages/automgrp

max_forks_repo_name: gap-packages/automgrp


text:
InstallMethod( JupyterRender, [ IsRecord ],
               r -> Objectify( JupyterRenderableType
                             , rec( data := rec( text\/plain := String(r) )
                                   , metadata := rec() ) ) );

# This is still an ugly hack, but its already much better than before!
BindGlobal("JupyterSplashDot",
function(dot)
    local fn, fd, r;

    fn := TmpName();
    fd := IO_File(fn, "w");
    IO_Write(fd, dot);
    IO_Close(fd);

    fd := IO_Popen(IO_FindExecutable("dot"), ["-Tsvg", fn], "r");
    r := IO_ReadUntilEOF(fd);
    IO_close(fd);
    IO_unlink(fn);

    return JupyterRenderable( rec( ("image/svg+xml") := r )
                            , rec( ("image/svg+xml") := rec( width := 500, height := 500 ) ) );
end);

# Splash the subgroup lattice of a group
BindGlobal("JupyterSplashSubgroupLattice",
function(group)
    local fn, fd, r, L, dot;

    fn := TmpName();

    L := LatticeSubgroups(group);
    DotFileLatticeSubgroups(L, fn);

    fd := IO_Popen(IO_FindExecutable("dot"), ["-Tsvg", fn], "r");
    r := IO_ReadUntilEOF(fd);
    IO_close(fd);
    IO_unlink(fn);

    return JupyterRenderable( rec( ("image/svg+xml") := r )
                            , rec( ("image/svg+xml") := rec( width := 500, height := 500 ) ) ) ;

end);

# To show TikZ in a GAP jupyter notebook
BindGlobal("JupyterSplashTikZ",
function(tikz)
    local tmpdir, fn, header, ltx, svgfile, stream, svgdata, tojupyter;

    header:=Concatenation( "\\documentclass[crop,tikz]{standalone}\n",
                    "\\usepackage{pgfplots}",
                    "\\makeatletter\n",
                    "\\batchmode\n",
                    "\\nonstopmode\n",
                    "\\begin{document}",
                    "\\begin{tikzpicture}");
    header:=Concatenation(header, tikz);
    header:=Concatenation(header,"\\end{tikzpicture}\n\\end{document}");

    tmpdir := DirectoryTemporary();
    fn := Filename( tmpdir, "svg_get" );

    PrintTo( Concatenation( fn, ".tex" ), header );

    ltx := Concatenation( "pdflatex -shell-escape --output-directory ",
                   Filename( tmpdir, "" ), " ",
                   Concatenation( fn, ".tex" ), " > ", Concatenation( fn, ".log2" ) );
    Exec( ltx );

    if not( IsExistingFile( Concatenation(fn, ".pdf") ) ) then
        tojupyter := rec( json := true, name := "stdout",
                          data := "No pdf was created; pdflatex is installed in your system?" );
    else
        svgfile := Concatenation( fn, ".svg" );
        ltx := Concatenation( "pdf2svg ", Concatenation( fn, ".pdf" ), " ",
                       svgfile, " >> ", Concatenation( fn, ".log2" ) );
        Exec( ltx );

        if not( IsExistingFile( svgfile ) ) then
            tojupyter := rec( json := true, name := "stdout",
                              data := "No svg was created; pdf2svg is installed in your system?" );
        else
            stream := InputTextFile( svgfile );
            if stream <> fail then
                svgdata := ReadAll( stream );
                tojupyter := rec( json := true, source := "gap",
                                  data := rec( ( "image/svg+xml" ) := svgdata ),
                                  metadata := rec( ( "image/svg+xml" ) := rec( width := 500, height := 500 ) ) );
                CloseStream( stream );
            else
                tojupyter := rec( json := true, name := "stdout",
                                  data := Concatenation( "Unable to render ", tikz ), metadata := rec() );
            fi;
        fi;
    fi;

    return JupyterRenderable(tojupyter.data, tojupyter.metadata);
end);

# This is really not what I should be doing here...
InstallGlobalFunction(ISO8601Stamp,
function()
    local tz, gm, pad;

    tz := IO_gettimeofday();
    pad := function(i, l, c)
        local s;
        s := String(i);
        if Length(s) < l then
            return Concatenation(RepeatedString(c, l - Length(s)), s);
        else
            return s;
        fi;
    end;

    gm := IO_gmtime(tz.tv_sec);
    return STRINGIFY( 1900 + gm.tm_year, "-"
                      , pad(gm.tm_mon + 1, 2, '0'), "-"
                      , pad(gm.tm_mday, 2, '0'), "T"
                      , pad(gm.tm_hour, 2, '0'), ":"
                      , pad(gm.tm_min, 2, '0'), ":"
                      , pad(gm.tm_sec, 2, '0'), "."
                      , pad(tz.tv_usec, 6, '0') );
end);

max_stars_repo_name: ZachNewbery/JupyterKernel

max_issues_repo_name: ZachNewbery/JupyterKernel

max_forks_repo_name: ZachNewbery/JupyterKernel