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

## Hackity.
#Class(fcall_addr, fcall, rec(
#        computeType := self >> self.args[3].t
#));

# A dist_loop is a parallel loop with domain=#spus

Class(dist_loop, loop_base, rec(
   __call__ := meth(self, P, loopvar, range, cmd)
       local result;
       Constraint(IsVar(loopvar));
       Constraint(IsCommand(cmd));
       range := toRange(range);

       loopvar.setRange(range);
       #loopvar.isLoopIndex := true;
       return WithBases(self, rec(
           operations := CmdOps,
           P := P,
           var := loopvar,
           cmd := cmd,
           range := listRange(range)
       ));
   end,

   rChildren := self >> [self.P, self.var, self.range, self.cmd],
   rSetChild := rSetChildFields("P", "var", "range", "cmd"),

   #rChildren := self >> [self.P, self.var, self.cmd],
   #rSetChild := rSetChildFields("P", "var", "cmd", "range"),
   #from_rChildren := (self, rch) >> ApplyFunc(ObjId(self), [rch[1], rch[2], self.range, rch[3]]),

   print := (self, i, is) >> Print(self.name, "(", self.P, ", ", self.var, ", ",
       Blanks(i+is),
       self.cmd.print(i+is, is),
       Print("\n", Blanks(i), ")")),

   #NOTE: This should go in loop_bases, which should free self.vars only if it's bound.
   free := meth(self) local c;
       c := self.cmd.free();
       return c;
   end

));

Class(DistCodegen, VectorCodegen, rec(
    # Wed 16 Jul 2008 07:17:33 PM EDT
    #NOTE: commenting this out! Check to see if this is now equivalent to DefaultCodegen.Formula.
    # This had BlockSums() commented out. Why?

    #NOTE: Need to find a way (hack?) to callback DefaultCodegen's formula with one extra stage, etc.
    #Formula := meth(self, o, y, x, opts)
    #    local icode, datas, prog, params, sub, initsub, io;
    #    if IsBound(opts.XType) then x.t := TPtr(opts.XType); fi;
    #    if IsBound(opts.YType) then y.t := TPtr(opts.YType); fi;

    #    o := o.child(1);
    #    params := Set(Collect(o, param));

    #    datas := Collect(o, FDataOfs);
    #    #o := BlockSums(opts.globalUnrolling, o);
    #    icode := self(o, y, x, opts);
    #    icode := RemoveAssignAcc(icode);
    #    icode := BlockUnroll(icode, opts);
    #    # icode := PowerOpt(icode);
    #    icode := DeclareHidden(icode);
    #    # icode := InsertBarriers(icode);
    #    if IsBound(opts.isFixedPoint) and opts.isFixedPoint then
    #        icode := FixedPointCode(icode, opts.bits, opts.fracbits);
    #    fi;

    #    io := When(X=Y, [X], [Y, X]);
    #    #io := When(IsBound(opts.multibuffer_its) and opts.multibuffer_its > 1,
    #    #                Concatenation(io, [Yprev, Xnext]),
    #    #                io);
    #    sub := Cond(IsBound(opts.subName), opts.subName, "transform");
    #    initsub := Cond(IsBound(opts.subName), Concat("init_", opts.subName), "init");
    #    icode := func(TVoid, sub, Concatenation(io, params), icode);

    #    if IsBound(opts.generateInitFunc) and opts.generateInitFunc then
    #        prog := program(
    #            decl(List(datas, x->x.var),
    #                chain(
    #                    func(TVoid, initsub, params, chain(List(datas, x -> SReduce(x.var.init, opts)))),
    #                    icode
    #                )));
    #    else
    #        prog := program( func(TVoid, initsub, params, chain()), icode);
    #    fi;
    #    prog.dimensions := o.dimensions;
    #    return prog;
    #end,

    # ---------------
    # Container
    # ---------------

#    DContainer := meth(self, o, y, x, opts)
#         local spuid;
#         spuid := var.fresh_t("spuid", TInt);
#         return( decl(
#             spuid,
#             self(o.child(1), y, x, opts)
#             ));
#      end,

    DContainer := (self, o, y, x, opts) >> self(o.child(1), y, x, opts),

    # ---------------
    # Gather
    # ---------------

    #GathDist := (self, o, y, x, opts) >> self(Gath(fId( (o.N/o.P)*o.pkSize )), y, x, opts),
    GathDist := (self, o, y, x, opts) >> call(var("// GathDist"), y, x),

    # Pass the base pointer of the data to be transferred as a param to the DMA_GET

    GathRecv := meth(self, o, y, x, opts)
        local i, func, rfunc, ix;
        i := Ind(); func := o.func.lambda();

        # Standard Gather
        #return loop(i, o.func.domain()/o.P, assign(nth(y, i), nth(x, func.at(i))));

        # We shouldn't use func.at(i), since we've done ScatSend-side normalizing already
        #return chain(
        #    #dist_barrier(),
        #    # Moving dist_barrier to scatsend for now
        #    call(var("// BLOCK_ON_READ")),
        #    loop(i, o.pkSize*o.func.domain()/o.P, assign(nth(y, i), nth(x, i)))
        #);

        return call(var("// Gath_Recv"), y, x);
    end,

    # ---------------
    # Scatter
    # ---------------

    #ScatDist := (self, o, y, x, opts) >> self(Scat(fId( (o.N/o.P)*o.pkSize )), y, x, opts),
    ScatDist := (self, o, y, x, opts) >> call(var("// ScatDist"), y, x),

    #NOTE: Assuming that Length(Y) is the full size of the transform.

                   #idiv(func.at( (i*pkSize)+(chunkSize*o.i) ), chunkSize), # This is the SPU#
    ScatSend := meth(self, o, y, x, opts)
        local i, func, rfunc, ix, pkSize, numPktsPerSPU;
        func          := o.func.lambda();
        pkSize        := o.pkSize;
        numPktsPerSPU := o.func.domain()/o.P;
        i             := Ind(numPktsPerSPU);

        return chain(loop(i, numPktsPerSPU, chain(                            # Loop over number of packets
           call(var("SCATSEND_PUT"),
               x,                                                       # Base address of X
               i*pkSize,                                                # Offset of X in elements
               fcall(var("ADDR", TFunc(TInt, TInt, TULongLong)),
                   idiv(func.at((numPktsPerSPU*o.i)+i), numPktsPerSPU), # This is the SPU#
                   y                                                    # Base address of Y
               ),
               mul(
                imod(func.at((numPktsPerSPU*o.i)+i), numPktsPerSPU),    # Offset of Y in elements
                pkSize),
               pkSize                                                   # Size of DMA transfer in elements
           )
        )),
        dist_barrier());

        # Standard Scatter

        # return loop(i, numPktsPerSPU,                                    # Loop over number of packets * pkSize
        #     assign(nth(y, func.at(i)), nth(x, i))
        # );
    end,

    # ---------------
    # Sum
    # ---------------

    # For now, assume loop range = # of spus

    DistSum := (self, o, y, x, opts) >> let(
        dist_loop(o.P, o.var, o.domain, self(o.child(1), y, x, opts))
        ),

    DistSumLoop := (self, o, y, x, opts) >> chain(
        dist_loop(o.P, o.var, o.domain, self(o.child(1), y, x, opts)),
        dist_barrier()
        ),

# Doesn't work for unrolled code.
#   DistSum := (self, o, y, x, opts) >> let(
#       Constraint(o.domain = opts.spus),
#       # If we're GathRecv'ing inside, we must sync before getting into the loop.
#       When(o._children[1]._children[Length(o._children[1]._children)].name = "GathRecv",
#         chain(dist_barrier(), dist_loop(o.P, o.var, o.domain, self(o.child(1), y, x, opts))),
#         dist_loop(o.P, o.var, o.domain, self(o.child(1), y, x, opts))
#       )
#    )

    # PTensor and Comm_Cell are used by ParCellDMP only.
    PTensor := (self, o, y, x, opts) >> self(o.L, y, x, opts),

    # Looped version (NOTE: works, except since we'll never do more than 8 SPUs for now, we always want a fully unrolled version

    Comm_Cell := meth(self, o, y, x, opts)
        local i;
        i := Ind(o.P);

        return(chain(loop(i, o.P,             # Loop over P
            call(var("SCATSEND_PUT"),
                x,                      # Base address of X
                i*o.pkSize,             # Offset of X in elements
                fcall(var("ADDR"),
                    i,                  # This is the SPU#
                    y                   # Base address of Y
                ),
                fcall(var("SPUID_TIMES"), o.pkSize),             # Offset of Y in elements (MUST BE spuid*pkSize!)
                o.pkSize                # Size of DMA transfer in elements
                )
        ),
        dist_barrier()
        ));
    end

));


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

_symGreater1 := o -> IsSymbolic(o) or o > 1;

Class(RulesRankedIndf, RuleSet);
RewriteRules(RulesRankedIndf, rec(
   fCompose_flatten := ARule(fCompose, [@(1,fCompose)], e->@(1).val.children()),

   # Degenerate HH rules (HH -> fAdd)
   #HH_to_fAdd := Rule([HH, @, @, @(1).cond(e->e[1]=1 and Sum(e)=1)], e -> fAdd(e.N, e.n, 0)),

   fAddElim := ARule(fCompose,
       [@(1), @(2, fAdd, e -> e.params[3] = 0 and range(e) = range(@(1).val) and domain(e)=domain(@(1).val))],
        e -> [ @(1).val ]),

   HHZ_fAdd := ARule(fCompose, [@(1, HHZ), @(2, [fAdd])], e -> When(@(2).val.params[3]=0, [ @(1).val ],
	[ HHZ(@(1).val.params[1], @(2).val.params[2], @(1).val.params[3] + @(2).val.params[3] * @(1).val.params[4][1], @(1).val.params[4])])),

   HH_fAdd := ARule(fCompose, [@(1, HH), @(2, fAdd)], e ->
	[ HH(@(1).val.params[1], @(2).val.params[2], @(1).val.params[3] + @(2).val.params[3] * @(1).val.params[4][1],
	                           @(1).val.params[4]) ]),

   fAdd_HH := ARule(fCompose, [@(2, fAdd), @(1, HH)], e ->
	[ HH(range(@(2).val), domain(@(1).val), @(1).val.params[3] + @(2).val.params[3], @(1).val.params[4]) ]),

   # Maude:
   #  var b1 b2 m n MM NN : NatExp . var d : Diagf . var i : Indf .
   #  var w ww : NatExpList . var v vv : NeNatExpList . vars gg gf : Genf .
   #  eq  H (    b1, n w) o H (    b2, m ww)  =  H (    b1 + n * b2, (n * m)  (w + lmul(n, ww))) .
   #  eq  HZ(NN, b1, n w) o H (    b2, m ww)  =  HZ(NN, b1 + n * b2, (n * m)  (w + lmul(n, ww))) .
   #  eq  HZ(NN, b1, n w) o HZ(MM, b2, m ww)  =  HZ(NN, b1 + n * b2, (n * m)  (w + lmul(n, ww))) .
   #
   HHoHH := ARule(fCompose, [@(1, HH), @(2, HH)], e -> [ let(
       s1 := @(1).val.params[4], s2 := @(2).val.params[4],
       b1 := @(1).val.params[3], b2 := @(2).val.params[3],
       n := s1[1],  m := s2[1],  w := Drop(s1, 1),  ww := List(Drop(s2, 1), x -> n*x),
       HH(@(1).val.params[1], @(2).val.params[2], b1 + n*b2, Concatenation([n*m], ListAddZP(w, ww)))) ]),

   # This rule is invalid when HH stride is odd, so in general it should not be used
   # NOTE: this rule is needed for DCTs, it seems that in that case stride is always even, so its ok to use it
   #        how to guard this for libraries??
   BHHoHH := ARule(fCompose, [@(1, BHH), @(2, HH, e->_symGreater1(e.params[4][1]))], e -> [ let(
       s1 := @(1).val.params[4], s2 := @(2).val.params[4],
       b1 := @(1).val.params[3], b2 := @(2).val.params[3],
       n := s1[1],  m := s2[1],  w := Drop(s1, 1),  ww := List(Drop(s2, 1), x -> n*x),
       BHH(@(1).val.params[1], @(2).val.params[2], b1 + n*b2, Concatenation([n*m], ListAddZP(w, ww)),
           @(1).val.params[5])) ]),

   #BHH_pull_out_base := Rule([BHH, @, @, @.cond(e->e<>0), @, @], e ->
   #    fCompose(HH(e.params[1], e.params[1], e.params[3], [1]), BHH(e.params[1], e.params[2], 0, e.params[4], e.params[5]-2*e.params[3]))),

   HH_1_ftensor_pull_out := Rule([fTensor, [@(1,fCompose), [@(2,HH), @, @, @, [ListClass, _1]], ...], @(3, fId)], e->
       let(h := @(2).val, id := @(3).val.params[1],
           fCompose(HH(id*h.params[1], id*h.params[2], id*h.params[3], [1]),
                    fTensor(fCompose(Drop(@(1).val.rChildren(), 1)), @(3).val)))),

   # HHxI2_base_pull_out is for closure reduction in autolib RDFT.
   # Commented out to get smaller closure in DFT and it conflicts with HHxfId_base_pull_in rule below
   #HHxI2_base_pull_out := Rule([fTensor, [@(2,HH), @, @, @.cond(e->e<>0), @], @(3, fId)], e->
   #    let(hp := @(2).val.params, id := @(3).val.params[1],
   #        fCompose(HH(id*hp[1], id*hp[1], id*hp[3], [1]),
   #                 fTensor(HH(hp[1], hp[2], 0, hp[4]), @(3).val)))),

   # Pull in HH into fTensor.
   #
   HHxfId_base_pull_in := ARule(fCompose, [[HH, @(1), @(2), @(3), [ListClass, _1]],
       [fTensor, [@(4,HH), @, @, @, @], @(5, fId, x -> ForAll([@(1).val, @(2).val, @(3).val], e -> _divides(x.params[1], e)))]],
       e -> let( t := @(5).val.params[1], hh := @(4).val.params,
                 [fTensor( HH( div(@(1).val, t), hh[2], hh[3] + div(@(3).val, t), hh[4] ), @(5).val )])),

   HHoBHH := ARule(fCompose, [@(1, HH, e->e.params[3]=0 and _symGreater1(e.params[4][1])), @(2, BHH)], e -> [ let(
       s1 := @(1).val.params[4], s2 := @(2).val.params[4],
       b1 := @(1).val.params[3], b2 := @(2).val.params[3],
       n := s1[1],  m := s2[1],  w := Drop(s1, 1),  ww := List(Drop(s2, 1), x -> n*x),
       BHH(@(1).val.params[1], @(2).val.params[2], b1 + n*b2, Concatenation([n*m], ListAddZP(w, ww)),
           n*@(2).val.params[5])) ]),

   HHZoHH := ARule(fCompose, [@(1, HHZ), @(2, HH)], e -> [ let(
	       s1 := @(1).val.params[4],
	       s2 := @(2).val.params[4],
	       b := @(1).val.params[3], bb := @(2).val.params[3],
	       n := s1[1], m := s2[1], w := Drop(s1, 1),  ww := List(Drop(s2, 1), x->n*x),
	       HHZ(@(1).val.params[1], @(2).val.params[2], b + n*bb, Concatenation([n*m], ListAddZP(w, ww)))) ]),

   HHZoHHZ := ARule(fCompose, [@(1, HHZ), @(2, HHZ)], e -> [ let(
	       s1 := @(1).val.params[4],
	       s2 := @(2).val.params[4],
	       b := @(1).val.params[3], bb := @(2).val.params[3],
	       n := s1[1], m := s2[1], w := Drop(s1, 1),  ww := List(Drop(s2, 1), x->n*x),
	       HHZ(@(1).val.params[1], @(2).val.params[2], b + n*bb, Concatenation([n*m], ListAddZP(w, ww)))) ]),

   KHoKH := ARule(fCompose, [@(1, KH),
                            [@(2, KH), @, @, _0, [ListClass, @, _1], [ListClass, _0, _0]]],
       e -> [ let(
	       s1 := @(1).val.params[4],
	       s2 := @(2).val.params[4],
               corr := @(1).val.params[5],
	       b := @(1).val.params[3], bb := @(2).val.params[3],
	       n := s1[1], m := s2[1], w := Drop(s1, 1),  ww := List(Drop(s2, 1), x->n*x),
	       KH(@(1).val.params[1], @(2).val.params[2], b,
                  Concatenation([n*m], ListAddZP(w, ww)), corr)) ]),

   # HHxIoHH assumes @(3) cannot be pulled into fTensor
   HHxIoHH := ARule(fCompose, [[fTensor, [@(1, HH), @, @, @, [ListClass, _1]], @(2, fId)], @(3, HH)],
       e -> let( p := @(1).val.params, n := @(2).val.params[1],
                 [ HH(n*p[1], n*p[2], n*p[3], [1]), @(3).val ] )),

   # eq tr(n, m) o H(0, 1 n) = H(0, m 1) .
   # H.domain must be less or equal to n
   TroHH := ARule(fCompose,
       [@(1, Tr), @(2, HH, e -> let(n := @(1).val.params[1], Cond(AnySyms(n, e.params[2]), n=e.params[2], n>=e.params[2]) and e.params[4][1]=1 and ForAll([e.params[3]]::Drop(e.params[4],1), x -> _divides(n, x))))],
           e -> [ let(
	       n := @(1).val.params[1], m := @(1).val.params[2], b := @(2).val.params[3], str := Drop(@(2).val.params[4], 1),
	       HH(@(2).val.params[1], @(2).val.params[2], b/n, [m] :: List(str, x -> x/n))) ]),

   Refl0_u_HH0 := ARule(fCompose, [@(1, Refl0_u), @(2, HH, e -> e.params[3]=0 and e.params[4]=[1])], e -> [ let(
	       k := @(1).val.params[1],
	       HH(@(1).val.range(), @(2).val.domain(), 0, [k])) ]),

   Refl0_u_HHrest := ARule(fCompose, [@(1, Refl0_u), @(2, HH, e -> e.params[3]=e.params[2] and e.params[4]=[1, @(1).val.params[2]])], e -> [ let(
	       k := @(1).val.params[1],
	       BHH(@(1).val.range(), @(2).val.domain(), 1, [2*k,1], 2*@(1).val.range())) ]),

   Refl0_odd_HH0 := ARule(fCompose, [@(1, Refl0_odd), @(2, HH, e -> e.params[3]=0 and e.params[4]=[1])], e -> [ let(
	       k := @(1).val.params[1],
	       HH(@(1).val.range(), @(2).val.domain(), 0, [2*k+1])) ]),

   Refl0_odd_HHrest := ARule(fCompose, [@(1, Refl0_odd), @(2, HH, e ->  e.params[4]=[1, @(1).val.params[2]])], e -> [ let(
	       k := @(1).val.params[1],
	       BHH(@(1).val.range(), @(2).val.domain(), 1, [2*k+1,1], 2*@(1).val.range())) ]),

   Refl1_HH := ARule(fCompose, [@(1, Refl1), @(2, HH, e -> e.params[3]=0 and e.params[4]=[1, @(1).val.params[2]])], e -> [ let(
	       k := @(1).val.params[1],
	       BHH(@(1).val.range(), @(2).val.domain(), 0, [2*k,1], 2*@(1).val.range() - 1)) ]),

   MMoHH := ARule(fCompose, [@(1, MM), @(2, HH, e -> e.params[3]=0 and e.params[4]=[1, @(1).val.params[1]])], e -> [ let(
	       m := @(1).val.params[2],
	       KH(@(2).val.params[1], @(2).val.params[2], 0, [m, 1], [0,0])) ]),

   fTensor_AxI_HH_stride_1 := ARule(fCompose, [[@(1, fTensor), ..., [fId, @(2)]], # NOTE: assumes d | base
                                        [@(3,HH), @,
                                                  @.cond(e->_dividesUnsafe(@(2).val,e)),
                                                  @.cond(e->_dividesUnsafe(@(2).val,e)),
                                                  [ListClass, _1, ...]]], e ->
       let(p := @(3).val.params, d := @(2).val, vs := Drop(p[4],1), base := p[3],
           [ fTensor( fCompose( fTensor(DropLast(@(1).val.children(), 1)), HH(div(p[1],d), div(p[2],d), div(base,d),
                                  Concatenation([1], List(vs, v -> div(v,d))))), fId(d))])),

   fTensor_AxI_HH_stride_n := ARule(fCompose, # NOTE: check assumptions
       [[@(1,fTensor), ..., [fId, @(2)]],
        [@(3,HH), @, @, _0, [ListClass, @.cond(e->e=@(2).val), _0, _1]]],

       e -> let(
	   d := @(2).val,
	   p := @(3).val.params,
	   ch := @(1).val.children(),
	   A := fTensor(DropLast(ch, 1)),
           [ fCompose(HH(range(A)*d, range(A), 0, [d, 0, 1]), A) ])
   ),

   fTensor_HH_domain1_fId := ARule(fTensor, [ [@(1, HH), @, _1, @, @], @(2, fId) ],
       e -> let(
	   p := @(1).val.params,
	   n := @(2).val.params[1],
	   [ HH(n*p[1], n, n*p[3], [1] :: List(Drop(p[4], 1), x->n*x)) ])
   ),

   HH_domain1_s1_to_fId1 := Rule( [HH, _1, _1, _0, [ListClass, _1]], e -> fId(1) ),

   fTensor_fId_X_HH := ARule(fCompose, [[@(4,fTensor), @(1, fId), @(2)],
                                        [@(3,HH), @, @.cond(e->e=@(2).val.domain()), _0, [ListClass, _1, @.cond(e->e=@(2).val.domain())]]],
       e ->let(p := @(3).val.params, vs := p[4][2],
           [ HH(@(4).val.range(), @(2).val.domain(), 0, [1, vs]), @(2).val ])),

   # This rule is a hack because it changing advdims of SPL. Range of the outer loop is also
   # assumed to be less than @(5).val (though if it's not then there is overlap in the
   # inner loop.
   f2DTensor_HH_HH_HH :=ARule( fCompose, [
       [f2DTensor, [HH, @(1), @(2), @(3), [ListClass, _1]],
                   [HH, @(4), @(5), @(6), [ListClass, _1]]],
       [HH, @(7).cond(x->x=@(2).val*@(5).val or x=@(5).val*@(2).val),
            @(8).cond(x->x=@(2).val),
            _0,
            [ListClass, @(9).cond(x->x=@(5).val), _1]]],
       e -> [HH( @(1).val*@(4).val, @(2).val, @(6).val+@(4).val*@(3).val, [@(4).val, 1] )]),

   # simplified merge_tensors rule
   fTensor_fTensor := ARule(fCompose,
       [[fTensor, @(1), @(2)], [fTensor, @(3), @(4).cond(
                   e->e.range()=@(2).val.domain() or @(3).val.range() = @(1).val.domain())]],
       e -> [ fTensor(fCompose(@(1).val, @(3).val), fCompose(@(2).val, @(4).val)) ]),

   # eq crt(n, m) o H(0, 1 m) = HZ(m * n, 0, m n) .
   # eq crt(n, m) o H(0, m 1) = HZ(m * n, 0, n m) .
   CRToHH := ARule(fCompose, [ [CRT, @(1), @(2), 1, 1],
	                       @(3,HH,e -> e.params[3]=0 and e.params[4] in [[1,@(2).val], [@(2).val, 1]]) ], e -> [ let(
		n := @(1).val,
                m := @(2).val,
		s1 := @(3).val.params[4][1],
		HHZ(@(3).val.params[1], @(3).val.params[2], 0, When(s1=1, [m, n], [n,m]))) ]),

   ## Scat * ISumn,  ISumn * Gath
   Scat_ISumn := ARule(Compose,  [ @(1, Scat), @(2, ISumn) ],
     e -> [ ISumn(@(2).val.domain,
	          Scat(@(1).val.func.upRank()) * @(2).val.child(1)) ]),

   ISumn_Gath  := ARule(Compose, [ @(1, ISumn), @(2, Gath) ],
     e -> [ ISumn(@(1).val.domain,
	          @(1).val.child(1) * Gath(@(2).val.func.upRank())) ]),

   ## Prm * ISumn,  ISumn * Prm
   Prm_ISumn := ARule(Compose,  [ @(1, Prm), @(2, ISumn) ],
     e -> [ ISumn(@(2).val.domain,
	          Prm(@(1).val.func.upRank()) * @(2).val.child(1)) ]),

   ISumn_Prm  := ARule(Compose, [ @(1, ISumn), @(2, Prm) ],
     e -> [ ISumn(@(1).val.domain,
	          @(1).val.child(1) * Prm(@(2).val.func.upRank())) ]),

   ## Diag * ISumn,  ISumn * Diag
   Diag_ISumn := ARule(Compose,  [ @(1, Diag), @(2, ISumn) ],
     e -> [ ISumn(@(2).val.domain,
	          Diag(@(1).val.element.upRank()) * @(2).val.child(1)) ]),

   ISumn_Diag  := ARule(Compose, [ @(1, ISumn), @(2, Diag) ],
     e -> [ ISumn(@(1).val.domain,
	          @(1).val.child(1) * Diag(@(2).val.element.upRank())) ]),

   ##UU
   #Implementing non-zero base would probably require to have a 2D base
   #Implicit requirements 2 has base=0 and both leading dims are the same.
    UU_UU := ARule(fCompose, [@(1, UU), @(2, UU)],
        e-> let(b1X := @(1).val.params[3], b2X := @(2).val.params[3], b1Y := @(1).val.params[4], b2Y := @(2).val.params[4], s1 := @(1).val.params[5], s2 := @(2).val.params[5],
            c2 := @(2).val.params[6], ld := @(1).val.params[7], ss1 := @(1).val.params[8],
            ss2 := @(2).val.params[8], ss:=ListAddLP(When(Length(ss2)>0,Transposed([List(Transposed(ss2)[1],x->x*s1)
                    ,Transposed(ss2)[2]]),[]), ss1),
            [UU(@(1).val.params[1], @(2).val.params[2], b1X + s1*b2X, b1Y+b2Y, s1*s2, c2, ld, ss)])),

#    UU_H := ARule(fCompose, [@(1, UU), @(2, H)],
#        e-> let(b1X := @(1).val.params[3], b2X := @(2).val.params[3], b1Y := @(1).val.params[4], s1 := @(1).val.params[5], s2 := @(2).val.params[4],
#            ld := @(1).val.params[7], ss1 := @(1).val.params[8], c1 := @(1).val.params[6],
#            [UU(@(1).val.params[1], @(2).val.params[2], b1X + s1*b2X, b1Y, s1*s2, c1, ld, ss1)])),

));

ext: gi


text:
local reps,tori,labels,I,II,ords;

mizuno:=[[1,2,3,4,5,6,7,
          9,10,11,12,
          13,14,16,17,
          19,18,20,21,
          23,24,25,27,
          26,28,30,31,
          32,33,35,34,
          37,38,39,40,
          41,44,45,46,
          48,49,51,52,
          53,55,57,58,
          59,61,63,64,
          66,69,70,71,
          75,76,80,82,
          85,89,93,97],
         [8,
          15,22,29,36,
          43,42,47,50,
          54,56,60,62,
          65,67,68,72,
          73,74,77,78,
          79,81,83,84,
          86,87,88,90,
          91,92,94,95,
          96,98,99,100,
          101,102,103,104,
          105,106,107,108,
          109,110,111,112,
          114,113,115,116,
          117,118,119,120]];

reps:=[[[1,1],[2,1],[3,1],[4,1],[5,1],[6,1],[7,1],[8,1]],
       [[1,1],[2,1],[10,1],[11,1],[5,1],[6,1],[7,1],[8,1]],
       [[1,1],[3,1],[2,1],[10,1],[12,1],[13,1],[14,1],[8,1]],
    [[8,1],[18,1],[7,1],[6,1],[19,1],[1,1],[17,1],[4,1]],
    [[18,1],[6,1],[19,1],[1,1],[17,1],[21,1],[4,1],[8,1]],
    [[1,1],[2,1],[3,1],[4,1],[5,1],[6,1],[7,1]],
    [[8,1],[7,1],[14,1],[26,1],[19,1],[1,1],[17,1],[4,1]],
    [[18,1],[22,1],[19,1],[1,1],[17,1],[20,1],[21,1],[4,1],[6,1]],
    [[1,1],[2,1],[11,1],[10,1],[5,1],[6,1],[7,1]],
    [[18,1],[8,1],[31,1],[7,1],[33,1],[4,1],[19,1],[23,1]],
    [[8,1],[7,1],[6,1],[5,1],[4,1],[2,1],[3,1]],
    [[18,1],[31,1],[7,1],[33,1],[4,1],[29,1],[23,1],[19,1]],
    [[1,1],[2,1],[3,1],[10,1],[12,1],[13,1],[14,1]],
    [[1,1],[2,1],[3,1],[4,1],[5,1],[6,1],[8,1]],

    [[8,1],[7,1],[6,1],[5,1],[10,1],[11,-1],[2,1],[3,1]],
    [[31,1],[33,1],[18,1],[4,1],[39,-1],[41,1],[8,1],[19,1],[23,1]],
    [[17,1],[1,1],[19,1],[6,1],[7,1],[18,1],[4,1]],
    [[3,1],[5,1],[1,1],[6,1],[11,1],[12,1],[2,1],[8,1]],
       [[1,1],[3,1],[4,1],[5,1],[6,1],[7,1],[8,1]],
       [[8,1],[7,1],[13,1],[10,1],[11,-1],[12,1],[2,1],[3,1]],
       [[1,1],[2,1],[3,1],[4,1],[5,1],[6,1]],
       [[7,1],[6,1],[5,1],[4,1],[2,1],[3,1]],
       [[1,1],[3,1],[4,1],[2,1],[5,1],[7,1],[8,1]],
       [[3,1],[5,1],[1,1],[6,1],[11,1],[12,1],[2,1]],
       [[17,1],[1,1],[19,1],[14,1],[26,1],[7,1],[18,1],[4,1]],
       [[2,1],[4,1],[5,1],[6,1],[7,1],[8,1],[1,1]],
       [[7,1],[6,1],[5,1],[10,1],[11,-1],[2,1],[3,1]],
       [[2,1],[4,1],[5,1],[6,1],[7,1],[8,1]],
       [[47,1],[49,1],[5,1],[45,1],[44,1],[42,1],[27,1],[46,1]],

       [[1,1],[3,1],[4,1],[2,1],[5,1],[7,1]],
       [[24,1],[26,1],[7,1],[27,1],[23,1],[4,1],[25,1]],
       [[17,1],[1,1],[19,1],[6,1],[18,1],[4,1],[8,1]],
       [[7,1],[13,1],[10,1],[11,-1],[12,1],[2,1],[3,1]],
       [[1,1],[3,1],[10,1],[12,1],[2,1],[5,1],[7,1],[8,1]],
       [[2,1],[4,1],[5,1],[6,1],[7,1],[1,1]],
       [[1,1],[3,1],[4,1],[2,1],[6,1],[7,1],[8,1]],
       [[1,1],[3,1],[4,1],[2,1],[5,1]],
       [[17,1],[1,1],[19,1],[6,1],[18,1],[4,1]],
       [[3,1],[4,1],[2,1],[5,1],[7,1],[8,1]],
       [[5,1],[6,1],[7,1],[8,1],[1,1],[3,1],[2,1]],
       [[1,1],[3,1],[10,1],[12,1],[2,1],[5,1],[7,1]],
       [[1,1],[3,1],[4,1],[5,1],[6,1]],
       [[1,1],[3,1],[4,1],[2,1],[6,1],[7,1]],
       [[1,1],[3,1],[4,1],[2,1],[6,1],[8,1]],
       [[1,1],[3,1],[10,1],[12,1],[2,1],[5,1]],

       [[1,1],[3,1],[4,1],[6,1],[7,1],[8,1]],
       [[1,1],[3,1],[4,1],[2,1],[6,1]],
       [[3,1],[10,1],[12,1],[2,1],[5,1],[7,1],[8,1]],
       [[3,1],[4,1],[2,1],[5,1],[7,1]],
       [[5,1],[6,1],[7,1],[1,1],[3,1],[2,1]],
       [[1,1],[3,1],[4,1],[2,1]],
       [[3,1],[4,1],[2,1],[6,1],[7,1]],
       [[3,1],[10,1],[12,1],[2,1],[5,1],[7,1]],
       [[5,1],[6,1],[7,1],[2,1],[3,1]],
       [[1,1],[3,1],[5,1],[6,1],[2,1],[8,1]],
       [[3,1],[4,1],[2,1],[5,1]],
       [[3,1],[10,1],[12,1],[2,1],[5,1]],
       [[1,1],[3,1],[4,1],[6,1]],
       [[1,1],[3,1],[5,1],[6,1],[2,1]],
       [[1,1],[3,1],[5,1],[6,1]],
       [[1,1],[3,1],[2,1],[5,1],[7,1]],
       [[1,1],[3,1],[4,1]],
       [[2,1],[4,1],[1,1],[6,1]],
       [[1,1],[3,1],[2,1]],
       [[2,1],[3,1],[5,1],[7,1]],
       [[1,1],[3,1]],
       [[1,1],[4,1],[6,1]],
       [[1,1],[6,1]],
       [[1,1]]
       ];

I:=[[],
    [4], #??? E_8(a_1)
    [4,6], #??? E_8(a_2)
    [2,3,5], # E_8(a_3)
    [2,3,5,7], # E_8(a_4)
    [],
    [2,3,5,6], # E_8(b_4)
    [2,3,5,6,8], # E_8(a_5)
    [4], # E_7(a_1)
    [1,2,3,5,6], # E_8(b_5)
    [],
    [1,2,3,5,6,8], # E_8(a_6)
    [4,6], #??? E_7(a_2)
    [],
    [4], # D_7(a_1)
    [1,2,3,5,6,7], # E_8(b_6)
    [2,3,5], # E_7(a_3)
    [4], # E_6(a_1)A_1
    [],
    [4,5], #??? D_7(a_2)

    [],[],[],
    [4], # E_6(a_1)
    [2,3,5,6], # E_7(a_4)
    [],
    [4],# D_6(a_1)
    [],
    [1,2,3,4,6,7,8], # E_8(a_7)
    [],
    [1,2,3,5,6],# E_7(a_5)
    [2,3,5],# E_6(a_3)A_1
    [4,5],#??? D_6(a_2)
    [4],# D_5(a_1)A_2
    [],[],[],
    [2,3,5],# E_6(a_3)
    [],[],
    [4],# D_5(a_1)A_1
    [],[],[],

    [4],# D_5(a_1) check this also
    [],[],
    [4],# D_4(a_1)A_2
    [],[],[],[],
    [4],# D_4(a_1)A_1
    [],[],[],
    [4],# D_4(a_1)
    [],[],[],[],[],[],[],[],[],[],[],[]
    ];

#
# The parabolic where the representative is distinguished
#
II:=[[1,2,3,4,5,6,7,8],
     [1,2,3,4,5,6,7,8],
     [1,2,3,4,5,6,7,8],
     [1,2,3,4,5,6,7,8],
     [1,2,3,4,5,6,7,8],
     [1,2,3,4,5,6,7],
     [1,2,3,4,5,6,7,8],
     [1,2,3,4,5,6,7,8],
     [1,2,3,4,5,6,7],
     [1,2,3,4,5,6,7,8],
     [2,3,4,5,6,7,8],
     [1,2,3,4,5,6,7,8],
     [1,2,3,4,5,6,7],
     [1,2,3,4,5,6,8],

     [2,3,4,5,6,7,8],
     [1,2,3,4,5,6,7,8],
     [1,2,3,4,5,6,7],
     [1,2,3,4,5,6,8],
     [1,3,4,5,6,7,8],
     [2,3,4,5,6,7,8],
     [1,2,3,4,5,6],
     [2,3,4,5,6,7],
     [1,2,3,4,5,7,8],
     [1,2,3,4,5,6],
     [1,2,3,4,5,6,7],
     [1,2,4,5,6,7,8],
     [2,3,4,5,6,7],
     [2,4,5,6,7,8],
     [1,2,3,4,5,6,7,8],

     [1,2,3,4,5,7],
     [1,2,3,4,5,6,7],
     [1,2,3,4,5,6,8],
     [2,3,4,5,6,7],
     [1,2,3,4,5,7,8],
     [1,2,4,5,6,7],
     [1,2,3,4,6,7,8],
     [1,2,3,4,5],
     [1,2,3,4,5,6],
     [2,3,4,5,7,8],
     [1,2,3,5,6,7,8],
     [1,2,3,4,5,7],
     [1,3,4,5,6],
     [1,2,3,4,6,7],
     [1,2,3,4,6,8],
     [1,2,3,4,5],

     [1,3,4,6,7,8],
     [1,2,3,4,6],
     [2,3,4,5,7,8],
     [2,3,4,5,7],
     [1,2,3,5,6,7],
     [1,2,3,4],
     [2,3,4,6,7],
     [2,3,4,5,7],
     [2,3,5,6,7],
     [1,2,3,5,6,8],
     [2,3,4,5],
     [2,3,4,5],
     [1,3,4,6],
     [1,2,3,5,6],
     [1,3,5,6],
     [1,2,3,5,7],
     [1,3,4],
     [1,2,4,6],
     [1,2,3],
     [2,3,5,7],
     [1,3],
     [1,4,6],
     [1,6],
     [1]
     ];

labels:=["E_8","E_8(a_1)","E_8(a_2)","E_8(a_3)","E_8(a_4)","E_7","E_8(b_4)","E_8(a_5)","E_7(a_1)","E_8(b_5)","D_7","E_8(a_6)","E_7(a_2)","E_6A_1",
         "D_7(a_1)","E_8(b_6)","E_7(a_3)","E_6(a_1)A_1","A_7","D_7(a_2)","E_6","D_6","D_5A_2","E_6(a_1)","E_7(a_4)","A_6A_1","D_6(a_1)","A_6","E_8(a_7)",
         "D_5A_1","E_7(a_5)","E_6(a_3)A_1","D_6(a_2)","D_5(a_1)A_2","A_5A_2","A_4A_3","D_5","E_6(a_3)","D_4A_2","A_4A_2A_1","D_5(a_1)A_1","A_5","A_4A_2","A_4A_1^2","D_5(a_1)",
         "A_3^2","A_4A_1","D_4(a_1)A_2","D_4A_1","A_3A_2A_1","A_4","A_3A_2","D_4(a_1)A_1","A_3A_1^2","A_2^2A_1^2","D_4","D_4(a_1)","A_3A_1","A_2^2A_1","A_2^2","A_2A_1^3","A_3","A_2A_1^2","A_2A_1","A_1^4","A_2","A_1^3","A_1^2","A_1"];

tori:=[
    [2,2,2,2,2,2,2,2],
    [2,2,2,0,2,2,2,2],
    [2,2,2,0,2,0,2,2],
    [2,0,0,2,0,2,2,2],
    [2,0,0,2,0,2,0,2],
    [2,2,2,2,2,2,2,-27],
    [2,0,0,2,0,0,2,2],
    [2,0,0,2,0,0,2,0],
    [2,2,2,0,2,2,2,-21],
    [0,0,0,2,0,0,2,2],
    [-21,2,2,2,2,2,2,2],
    [0,0,0,2,0,0,2,0],
    [2,2,2,0,2,0,2,-17],
    [2,2,2,2,2,2,-17,2],

    [-16,2,0,2,2,2,2,2],
    [0,0,0,2,0,0,0,2],
    [2,0,0,2,0,2,2,-15],
    [2,2,2,0,2,2,-13,2],
    [2,-15,2,2,2,2,2,2],
    [-13,2,2,0,2,0,2,2],
    [2,2,2,2,2,2,-16,0],
    [-15,2,2,2,2,2,2,-10],
    [2,2,2,2,2,-12,2,2],
    [2,2,2,0,2,2,-12,0],
    [2,0,0,2,0,0,2,-11],
    [2,2,-11,2,2,2,2,2],
    [-11,2,2,0,2,2,2,-8],
    [0,2,-10,2,2,2,2,2],
    [0,0,0,0,2,0,0,0],

    [2,2,2,2,2,-11,2,-1],
    [0,0,0,2,0,0,2,-9],
    [2,0,0,2,0,2,-9,2],
    [-9,2,2,0,2,0,2,-6],
    [2,2,2,0,2,-9,2,2],
    [2,2,-9,2,2,2,2,-5],
    [2,2,2,2,-9,2,2,2],
    [2,2,2,2,2,-10,0,0],
    [2,0,0,2,0,2,-8,0],
    [-6,2,2,2,2,-8,2,2],
    [2,2,2,-7,2,2,2,2],
    [2,2,2,0,2,-8,2,-1],
    [2,-9,2,2,2,2,-5,0],
    [2,2,2,2,-8,2,2,-2],
    [2,2,2,2,-7,2,-2,2],
    [2,2,2,0,2,7,0,0],

    [2,-3,2,2,-6,2,2,2],
    [2,2,2,2,-7,2,-1,0],
    [-4,2,2,0,2,-6,2,2],
    [-6,2,2,2,2,-7,2,-1],
    [2,2,2,-6,2,2,2,-3],
    [2,2,2,2,-6,0,0,0],
    [-3,2,2,2,-6,2,2,-2],
    [-4,2,2,0,2,-5,2,-1],
    [-1,2,2,-5,2,2,2,-3],
    [2,2,2,-5,2,2,-3,2],
    [-6,2,2,2,2,-6,0,0],
    [-4,2,2,0,2,-4,0,0],
    [2,-3,2,2,-4,2,-1,0],
    [2,2,2,-5,2,2,-2,0],
    [2,0,2,-4,2,2,-2,0],
    [2,2,2,-4,2,-2,2,-1],
    [2,-3,2,2,-3,0,0,0],
    [2,2,-3,2,-3,2,-1,0],
    [2,2,2,-3,0,0,0,0],
    [-1,2,2,-3,2,-2,2,-1],
    [2,2,0,-2,0,0,0,0],
    [2,-1,-2,2,-2,2,-1,0],
    [2,0,-1,0,-1,2,-1,0],
    [2,0,-1,0,0,0,0,0]
    ];

#ords:=[];
ords:=[ 125, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 5, 25, 25, 25, 25, 25, 25, 5, 25, 25, 25, 5, 25, 25, 5, 5, 25, 5, 5, 5, 25, 5, 5, 5, 5, 5, 5, 25, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5 ];

return [reps,tori,labels,I,II,ords];


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

#############################################################################
##
## The following functions can be used to determine Betti-numbers of a
## torsion-free polycyclic group given by a pcp presentation. All
## Betti-numbers can be obtained if G has Hirsch length at most 5.
##
## The Betti-numbers B(G,m) are defined as the ranks of H_m(G,Z) for the
## trivial G-module Z. If M is the orientation G-module, then we can also
## characterise B(G,m) for n >= m >= n-2 as the ranks of H^n-m(G,M).
## Further, the alternating sum of all Betti-numbers is 0 using the
## Euler characteristic.
##

#############################################################################
##
#F OrientationModule( G )
##
InstallMethod( OrientationModule, "for pcp groups", true, [IsPcpGroup], 0,
function( G )
    local pcps, gens, mats, acts, dets, i, pcp;
    pcps := PcpsOfEfaSeries( G );
    pcps := Filtered( pcps, x -> RelativeOrdersOfPcp(x)[1] = 0 );
    gens := Igs(G);
    mats := List( gens, x -> IdentityMat( 1 ) );
    for pcp in pcps do
        acts := LinearActionOnPcp( gens, pcp );
        dets := List( acts, x -> Determinant( x ) );
        for i in [1..Length(mats)] do
            mats[i] := dets[i] * mats[i];
        od;
    od;
    return mats;
end );

#############################################################################
##
#F IsOrientedMatGroup( G )
##
IsOrientedMatGroup := function( G )
    return ForAll( GeneratorsOfGroup(G), x -> Determinant(x) = 1 );
end;

#############################################################################
##
#F BettiNumber( G, m )
##
BettiNumberPcpGroup := function(G,m)
    local n, pcp, mats, CR, one, two;

    if not IsTorsionFree( G ) then
        Print("the input group must be torsion-free \n");
        return fail;
    fi;

    # catch the trivial case
    if IsFinite(G) then
        if m = 0 then
            return 1;
        else
            return 0;
        fi;
    fi;

    # the hirsch length
    n := HirschLength( G );

    if m < 0 or m > n then return 0; fi;

    if m = 0 then return 1; fi;

    if m = 1 then
        pcp := Pcp( G, DerivedSubgroup(G) );
        return Length( Filtered( RelativeOrdersOfPcp( pcp ),x -> x=0 ));
    fi;

    if m = n then
        mats := OrientationModule( G );
        if ForAny( mats, x -> x[1][1] = -1 ) then
            return 0;
        else
            return 1;
        fi;
    fi;

    if m = 2 then
        mats := List( Pcp(G), x -> IdentityMat(1) );
        CR := CRRecordByMats( G, mats );
        two := TwoCohomologyCR( CR ).factor.rels;
        return Length( Filtered( two, x -> x = 0 ) );
    fi;

    if m = n-1 then
        mats := OrientationModule( G );
        CR := CRRecordByMats( G, mats );
        one := OneCohomologyCR( CR ).factor.rels;
        return Length( Filtered( one, x -> x = 0 ) );
    fi;

    if m = n-2 then
        mats := OrientationModule( G );
        CR := CRRecordByMats( G, mats );
        two := TwoCohomologyCR( CR ).factor.rels;
        return Length( Filtered( two, x -> x = 0 ) );
    fi;

    Print("Betti-number is out of range for our methods \n");
    return fail;
end;

InstallMethod( BettiNumber, "for torsion-free pcp groups", true,
   [IsPcpGroup, IsInt], 0,
function(G, m)
    if not IsTorsionFree(G) then TryNextMethod(); fi;
    if m in [3..HirschLength(G)-3] then TryNextMethod(); fi;
    return BettiNumberPcpGroup(G,m);
end);

#############################################################################
##
#F BettiNumbers( G )
##
InstallMethod( BettiNumbers, "for torsion-free pcp groups", true,
    [IsPcpGroup], 0,
function( G )
    local n, betti;

    n := HirschLength( G );
    if not IsTorsionFree( G ) or n > 6 then TryNextMethod(); fi;

    # set up the Betti-numbers
    betti := [1];
    if n = 0 then return betti; fi;
    betti[2] := BettiNumber( G, 1 );
    if n = 1 then return betti; fi;
    betti[3] := BettiNumber( G, 2 );
    if n = 2 then return betti; fi;
    betti[4] := betti[1] - betti[2] + betti[3];
    if n = 3 then return betti; fi;
    if n > 3 then
        betti[5] := BettiNumber( G, 4 );
        betti[4] := betti[4] + betti[5];
    fi;
    if n > 4 then
        betti[6] := BettiNumber( G, 5 );
        betti[4] := betti[4] - betti[6];
    fi;
    if n > 5 then
        betti[7] := BettiNumber( G, 6 );
        betti[4] := betti[4] + betti[7];
    fi;
    return betti;
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

DropMat := function(fmat)
    local mat;
    mat := Copy(fmat);
    if IsBound(mat.mat) then Unbind(mat.mat); fi;
    return mat;
end;

#   All size x size digit swap permutations that could be implemented in register for v
AllDigitSwaps := function(size, v)
    local ll, l, Nl, N, nl, n, r, mats, L1, L2, spl;

    mats := [];
    ll := DropLast(DivisorsInt(v), 1);
    for l in ll do
        Nl := Filtered(DivisorsInt(size/l), i -> i >= 4);
        for N in Nl do
            r:= size/(N*l);
            nl := DivisorsIntDrop(N);
            for n in nl do
                if l=1 then L1 := []; else L1 := [I(l)]; fi;
                if r=1 then L2 := []; else L2 := [I(r)]; fi;
                spl := Tensor(Concat(L1, [L(N,n)], L2));
                Add(mats, rec(spl:=spl, mat:=MatSPL(spl), l:=l, N:=N, n:=n, r:=r));
            od;
        od;
    od;

    return mats;
end;

#   NOTE: Better algorithm: outer loop over AllDigitSwaps

BuildCandidates := function(ops, v)
    local mats, op, p, v;
    mats := [];
    for op in ops do
        for p in op.permparams() do
            Add(mats,
                rec(instr := op,
                    p := p,
#   NOTE - how to get op.v correctly? maybe introduce op.getVLen() method
                    mat := let(ov := When(IsInt(op.v), op.v, op.v()), When(op.isUnop(),
                        TransposedMat(List([0..ov-1], i->op.semantic(BasisVec(ov, i), p))),
                        TransposedMat(Concat(List([0..ov-1], i->op.semantic(BasisVec(ov, i), Replicate(v, 0), p)),
                            List([0..ov-1], i->op.semantic(Replicate(v, 0), BasisVec(ov, i), p))))
                        ))
                ));
        od;
    od;
    return mats;
end;

#   Build rules from binary operations
#
BuildBinRules := function(arch)
    local archrules, binops, p, op1, op2, i1, i2, v, binmats, binperms, dsperms, opmat, fmat;

    v := arch.v;
    archrules := [];
    binops := Filtered(arch.instr, i->i.isBinop());

    #   find al candidate instructions to build permutations
    binmats := BuildCandidates(binops, v);

    # find all pairs of binops to build valid swap perms
    dsperms := AllDigitSwaps(2*v, v);
    for op1 in binmats do
        for op2 in binmats do
            opmat := Concat(op1.mat, op2.mat);
            if IsPermMat(AMatMat(opmat)) then
                fmat := Filtered(dsperms, i->i.mat=opmat);
                if Length(fmat)=1 then
                    i1 := op1.instr; i2 := op2.instr;
                    Add(archrules, rec(perm := DropMat(fmat[1]), instr:=[DropMat(op1), DropMat(op2)], v:=v, # NOTE: use full list!
                                vperm := VPerm(fmat[1].spl,
                                Subst((y,x) -> chain(
                                assign(vref(y, 0,$v), $i1(vref(x,0,$v), vref(x,$v,$v), $(op1.p))),
                                assign(vref(y,$v,$v), $i2(vref(x,0,$v), vref(x,$v,$v), $(op2.p))))),
                            v, op1.instr.vcost + op2.instr.vcost)
                            ));
                    dsperms := Filtered(dsperms, i->opmat<>i.mat);
                fi;
            fi;
        od;
    od;
    return archrules;
end;

#   Build rules from binary operations - first shuffle Inputs by LxI
#
BuildBinRulesLxI := function(arch, binrules, left, right)
    local archrules, binops, p, op1, op2, i1, i2, v, binmats, binperms, dsperms, opmat, fmat, lm, lmleft, lmright,
        lmc, lmvleft, lmvright, opmatconcat, i2v, splleft, splright, lval, rval;

    v := arch.v;
    archrules := [];
    binops := Filtered(arch.instr, i->i.isBinop());

    lm := MatSPL(Tensor(L(4,2), I(v/2)));
    lmleft := When(left, lm, MatSPL(I(2*v)));
    splleft := When(left, Tensor(L(4,2), I(v/2)), I(2*v));
    lmright := When(right, lm, MatSPL(I(2*v)));
    splright := When(right, Tensor(L(4,2), I(v/2)), I(2*v));
    i2v := MatSPL(I(2*v));

    lmc := Filtered(binrules, i->(i.perm.l=1 and i.perm.r=v/2 and i.perm.N=4 and i.perm.n=2));

    if Length(lmc) = 0 then return []; fi;

    lmvleft := When(left, lmc[1].vperm, VTensor(I(2), v));
    lmvright := When(right, lmc[1].vperm, VTensor(I(2), v));

    lval := When(left, lmc[1].vperm._vcost + 0.2, 0);
    rval := When(right, lmc[1].vperm._vcost + 0.2, 0);

    #   find al candidate instructions to build permutations
    binmats := BuildCandidates(binops, v);

    # find all pairs of binops to build valid swap perms
    dsperms := AllDigitSwaps(2*v, v);
    for op1 in binmats do
        for op2 in binmats do
            opmatconcat := Concat(op1.mat, op2.mat);
            if opmatconcat  <> i2v then
                opmat := lmleft * opmatconcat * lmright;
                if IsPermMat(AMatMat(opmat)) then
                    fmat := Filtered(dsperms, i->i.mat=opmat);
                    if Length(fmat)=1 then
                        i1 := op1.instr; i2 := op2.instr;
                        Add(archrules, rec(perm := DropMat(fmat[1]), instr:=[DropMat(op1), DropMat(op2)], v:=v, # NOTE: use full list!
                                    vperm := lmvleft * VPerm(splleft * fmat[1].spl * splright,
                                    Subst((y,x) -> chain(
                                    assign(vref(y, 0,$v), $i1(vref(x,0,$v), vref(x,$v,$v), $(op1.p))),
                                    assign(vref(y,$v,$v), $i2(vref(x,0,$v), vref(x,$v,$v), $(op2.p))))),
                                v, 2 + lval + rval) * lmvright
                                ));
                        dsperms := Filtered(dsperms, i->opmat<>i.mat);
                    fi;
                fi;
            fi;
        od;
    od;
    return archrules;
end;

#   Build rules from unary operations
#
BuildUnRules := function(arch)
    local archrules, unops, unmats, unperms, unperms2, unmats2, p, op, op1, op2, v, dsperms, opmat, fmat, i1, i2;

    v := arch.v;
    archrules := [];

    #   unops -> perms
    unops := Filtered(arch.instr, i->i.isUnop());
    unmats := BuildCandidates(unops, v);

    # find all unops to build valid swap perms
    dsperms := AllDigitSwaps(v, v);
    unmats2 := [];
    for op in unmats do
        opmat := op.mat;
        if IsPermMat(AMatMat(opmat)) then
            fmat := Filtered(dsperms, i->i.mat=opmat);

            if Length(fmat)>=1 then
        i1 := op.instr;
        Add(archrules, rec(perm := DropMat(fmat[1]), instr:=DropMat(op), v:=v, # NOTE: use full list!
                    vperm := VPerm(fmat[1].spl,
                                   Subst((y,x) -> assign(vref(y, 0, $v), $i1(vref(x,0,$v), $(op.p)))),
                   v, 0.9*op.instr.vcost) # NOTE: avoid unops as they may not exist
        ));
                dsperms := Filtered(dsperms, i->opmat<>i.mat);
            else
                Add(unmats2, op);
            fi;
        fi;
    od;

    #   product of unops -> perms
    for op1 in unmats2 do
        for op2 in unmats2 do
            opmat := op2.mat * op1.mat;
            if IsPermMat(AMatMat(opmat)) then
                fmat := Filtered(dsperms, i->opmat=i.mat);

                if Length(fmat)=1 then
                    i1 := op1.instr;
                    i2 := op2.instr;
                    Add(archrules, rec(
                        perm  := DropMat(fmat[1]),
                        instr := [DropMat(op1), DropMat(op2)],
                        v     := v, # NOTE: use full list!
                        vperm := VPerm( fmat[1].spl, Subst((y,x) -> let(t := TempVec(TArray(x.t.t, $v)),
                                        decl(t, chain(assign(vref(t,0,$v), $i1(vref(x,0,$v), $(op1.p))),
                                                      assign(vref(y,0,$v), $i2(vref(t,0,$v), $(op2.p))))))),
                                        v, 0.9*(op1.instr.vcost+op2.instr.vcost))
                    ));
                    dsperms := Filtered(dsperms, i->opmat<>i.mat);
                fi;
            fi;
        od;
    od;

    return archrules;
end;

Find_x_I_vby2 := function(arch, i, j)
    local gmat, v, binops, binmats, fmat, op;

    v := arch.v;
    binops := Filtered(arch.instr, i->i.isBinop());

    #   find al candidate instructions to build permutations
    binmats := BuildCandidates(binops, v);
    gmat := MatSPL(Gath(fTensor(fDirsum(fBase(2,i), fBase(2,j)), fId(v/2))));

    fmat := Filtered(binmats, i->i.mat=gmat);
    return When(Length(fmat) > 0, DropMat(fmat[1]), []);
end;

Build_x_I_vby2 := arch -> List([0,1], i-> List([0,1], j-> Find_x_I_vby2(arch,i,j)));


text:
BindGlobal( "MenuOpsForFiniteGroups",
[ rec( name := "All Subgroups",
       func := function(poset, G)
            local L, cls, len, sz, graphHasse, max, rep, z, t, i, j, k, nodes, last, knownArgs, m, cb;

            L:=LatticeSubgroups(G);;
            cls:=ConjugacyClassesSubgroups(L);
            len:=[];
            sz:=[];
            for i in cls do
              Add(len,Size(i));
              AddSet(sz,Size(Representative(i)));
            od;

            graphHasse := poset!.graph;
            graphHasse!.nodes := rec(); # initialize nodes
            graphHasse!.links := rec(); # initialize links

            nodes := [];
            sz:=Reversed(sz);

            # subgroup nodes, also acccording to size
            for i in [1..Length(cls)] do
              nodes[i] := [];
              for j in [1..len[i]] do
                if len[i]=1 then
                    nodes[i][j] := Shape(ShapeType.DIAMOND, String(i));
                else
                    nodes[i][j] := Shape(ShapeType.CIRCLE, String(i));
                fi;
                SetLayer(nodes[i][j], -Size(Representative(cls[i])));
                for m in poset!.latticeType.contextMenus do
                  knownArgs := [poset, Representative(cls[i])];
                  if m.group = true then
                    Add(knownArgs, G);
                  fi;
                  cb := Callback(m.func, knownArgs);
                  Add(nodes[i][j], Menu(m.name, cb));
                od;
                if i = Length(cls) and j = len[i] then
                  SetTitle(nodes[i][j], "G"); # set G
                fi;
                Add(graphHasse, nodes[i][j]);
              od;
            od;
            # uniform layers hack
            last:=rec(o:=0,n:=0);
            for i in [1..Length(cls)] do
              for j in [1..len[i]] do
                if Layer(nodes[i][j]) <> last.o then
                  last.o := Layer(nodes[i][j]);
                  last.n := last.n - 2;
                fi;
                SetLayer(nodes[i][j], last.n);
              od;
            od;

            max:=MaximalSubgroupsLattice(L);
            for i in [1..Length(cls)] do
              for j in max[i] do
                rep:=ClassElementLattice(cls[i],1);
                for k in [1..len[i]] do
                  if k=1 then
                    z:=j[2];
                  else
                     t:=cls[i]!.normalizerTransversal[k];
                    z:=ClassElementLattice(cls[j[1]],1); # force computation of transv.
                    z:=cls[j[1]]!.normalizerTransversal[j[2]]*t;
                    z:=PositionCanonical(cls[j[1]]!.normalizerTransversal,z);
                  fi;
                  Add(graphHasse, Link(nodes[j[1]][z], nodes[i][k]));
                od;
              od;
            od;
            return Draw(poset);
        end, multiple := false, group := true ),
  #rec( name := "Centralizers", func := Centralizer ),
  #rec( name := "Centres", func := Centre ),
  #rec( name := "Closures", func := ClosureGroup ),
  #rec( name := "Commutator Subgroups", func := CommutatorSubgroup ),
  #rec( name := "Conjugate Subgroups",
  #     func := {G, H} -> AsList(ConjugacyClassSubgroups(G, H))),
  #rec( name := "DerivedSeries", func := DerivedSeriesOfGroup ),
  #rec( name := "Fitting Subgroups", func := FittingSubgroup ),
  #rec( name := "Normal Closures", func := NormalClosure ),
  #rec( name := "Normal Subgroups", func := NormalSubgroups ),
  #rec( name := "Sylow Subgroups", func := GGLSylowSubgroup ),
] );

BindGlobal( "MenuOpsCommon",
[ #rec( name := "DerivedSubgroups", func := DerivedSubgroup ),
  #rec( name := "Cores", func := Core ),
  #rec( name := "Closure", func := GGLClosureGroup ),
  #rec( name := "Intermediate Subgroups", func := IntermediateSubgroups ),
  #rec( name := "Intersection", func := Intersection ),
  #rec( name := "Intersections", func := Intersection ),
  #rec( name := "Normalizers", func := Normalizer ),
] );

BindGlobal( "MenuOpsForFpGroups",
[ #rec( name := "Abelian Prime Quotient", func := GGLAbelianPQuotient ),
  #rec( name := "All Overgroups", func := IntermediateSubgroups ),
  #rec( name := "Compare Subgroups", func := GGLCompareSubgroups ),
  #rec( name := "Conjugacy Class",
  #     func := {G, H} -> Filtered(AsList(ConjugacyClassSubgroups(G,H)), h -> h <> H)),
  #rec( name := "Epimorphisms (GQuotients)", func := GGLEpimorphisms ),
  #rec( name := "Low Index Subgroups", func := GGLLowIndexSubgroups ),
  #rec( name := "Prime Quotient", func := GGLPrimeQuotient ),
  #rec( name := "Test Conjugacy", func := GGLTestConjugacy )
] );

BindGlobal( "ContextMenusForFiniteGroups",
  [ rec( name := "Size",
         func := function(poset, n)
              local message;
              message := FrancyMessage(FrancyMessageType.INFO, "Size",  String(Size(n)));
              Add(poset, message);
              return Draw(poset);
         end, group := false ),
    #rec( name := "Index",
    #     func := function(poset, n)
    #          local message;
    #          message := FrancyMessage(FrancyMessageType.INFO, "Index",  String(Index(n)));
    #          Add(poset, message);
    #          return Draw(poset);
    #     end, group := true ),
    rec( name := "IsAbelian",
         func := function(poset, n)
              local message;
              message := FrancyMessage(FrancyMessageType.INFO, "Is commutative?",  String(IsCommutative(n)));
              Add(poset, message);
              return Draw(poset);
         end, group := false ),
    #rec( name := "IsCentral",
    #     func := function(poset, n)
    #          local message;
    #          message := FrancyMessage(FrancyMessageType.INFO, "Is central?",  String(IsCentral(n)));
    #          Add(poset, message);
    #          return Draw(poset);
    #     end, group := false ),
    rec( name := "IsCyclic",
         func := function(poset, n)
              local message;
              message := FrancyMessage(FrancyMessageType.INFO, "Is cyclic?",  String(IsCyclic(n)));
              Add(poset, message);
              return Draw(poset);
         end, group := false ),
    rec( name := "IsNilpotent",
         func := function(poset, n)
              local message;
              message := FrancyMessage(FrancyMessageType.INFO, "Is NilpotentGroup?",  String(IsNilpotentGroup(n)));
              Add(poset, message);
              return Draw(poset);
         end, group := false ),
    rec( name := "IsNormal",
         func := function(poset, n)
              local message;
              message := FrancyMessage(FrancyMessageType.INFO, "Is normal?",  String(IsNormal(n)));
              Add(poset, message);
              return Draw(poset);
         end, group := true ),
    rec( name := "IsPerfect",
         func := function(poset, n)
              local message;
              message := FrancyMessage(FrancyMessageType.INFO, "Is perfect group?",  String(IsPerfectGroup(n)));
              Add(poset, message);
              return Draw(poset);
         end, group := false ),
    rec( name := "IsSimple",
         func := function(poset, n)
              local message;
              message := FrancyMessage(FrancyMessageType.INFO, "Is simple group?",  String(IsSimpleGroup(n)));
              Add(poset, message);
              return Draw(poset);
         end, group := false ),
    rec( name := "IsSolvable",
         func := function(poset, n)
              local message;
              message := FrancyMessage(FrancyMessageType.INFO, "Is solvable group?",  String(IsSolvableGroup(n)));
              Add(poset, message);
              return Draw(poset);
         end, group := false ),
    rec( name := "Isomorphism",
         func := function(poset, n)
              local message;
              message := FrancyMessage(FrancyMessageType.INFO, "Isomorphism",  String(IdGroup(n)));
              Add(poset, message);
              return Draw(poset);
         end, group := false ),
  ] );

BindGlobal( "ContextMenusForFpGroups",
  [ rec( name := "Index", func := Index ),
    rec( name := "IsNormal", func := IsNormal ),
    rec( name := "IsFpGroup", func := IsFpGroup ),
    rec( name := "IsSubgroupFpGroup", func := IsSubgroupFpGroup ),
    rec( name := "Abelian Invariants", func := AbelianInvariants ),
    rec( name := "CosetTable", func := CosetTableInWholeGroup ),
    rec( name := "IsomorphismFpGroup", func := IsomorphismFpGroup ),
    rec( name := "FactorGroup",
         func := function(G, N);
          if IsNormal(G, N) then
            return FactorGroup(G, N);
          else
            return "subgroup is not normal";
          fi;
         end )
  ] );

max_stars_repo_name: mcmartins/subgroup-lattice

max_issues_repo_name: mcmartins/subgroup-lattice

max_forks_repo_name: mcmartins/subgroup-lattice


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

#expand FData to report range as size
Class(FDataIndirect,FData,rec (
    range := self >> self.var.t.size,
));

Class(ICFDataIndirect,FDataIndirect,rec (
    range := self >> 2*self.var.t.size,
    domain := self >> 2*self.var.t.size,
));

Class(ExpandMult,NonTerminal,rec (
        abbrevs := [(height,rows,isSplit) -> [height,rows,isSplit,[]]],
        dmn := self >> let(x:=When(self.params[3],2,1),[TArray(TReal,x*self.params[1] * self.params[2]),TArray(TReal,x*self.params[2])]),
        rng := self >> let(x:=When(self.params[3],2,1),[TArray(TReal,x*self.params[1]* self.params[2])]),
        isReal := self >> true,
        tagpos := 4,
        transpose := self >> Copy(self)));

NewRulesFor(ExpandMult,rec(

       ExpandMult_Base:= rec (
        forTransposition := false,
        applicable := (self,nt) >> nt.params[1]=1,
        apply := (nt,c,cnt) -> When(nt.params[3],ICMultiplication(2,2*nt.params[2],true),ICMultiplication(2,nt.params[2]))),

       ExpandMult_One := rec (
        forTransposition := false,
        applicable := (self,nt) >> nt.params[1] > 1,
        freedoms := nt -> [[1]],
        child := (self,nt,freedoms) >> [TTensorI_OL(ExpandMult(1,nt.params[2],nt.params[3]),[APar,AOne],[[1,0]],
                        [nt.params[1]/freedoms[1],1])],
        apply := (nt,c,cnt) -> c[1]
       )));

Class(SAR_Interpolation,NonTerminal,rec(
        abbrevs := [(iters,height,rows,isSplit) -> Checked(IsPosInt(rows) and IsPosInt(iters) and
                     IsPosInt(height),[iters,height,rows,isSplit,[]]),
                    (iters,height,rows,isSplit,t) -> Checked(IsPosInt(rows) and IsPosInt(iters) and
                        IsPosInt(height), [iters,height,rows,isSplit,t])],
        dmn := self >> [TArray(TReal,2*self.params[1] * self.params[2] * self.params[3]),
                        TArray(TReal,2*self.params[1] * self.params[3])],
        rng := self >> [TArray(TReal,2*self.params[1] * self.params[3])],
        isReal := self >> true,
        isSplit := self >> self.params[4],
        numops := self >> self.params[1] * self.params[2] * self.params[3] * self.params[1],
        tagpos := 5,
        transpose := self >> Copy(self)));

NewRulesFor(SAR_Interpolation,rec(

     SAR_Interpolation_One := rec (
        forTransposition := false,
        applicable := (self,nt) >> nt.params[1] > 1 and nt.params[2] > 1,
        freedoms := nt -> [[1]],
        child := function(nt,freedoms)
		local size,tt,fd;
                size:=(nt.params[1]*nt.params[2]*nt.params[3]);
		fd:=When(nt.params[4],
			ICFDataIndirect(param(TArray(TDouble,size),"indirect")),FDataIndirect(param(TArray(TDouble,size),"indirect")));
                tt:=When(nt.params[4],ExpandMult(nt.params[2],nt.params[3],true),
                        ExpandMult(nt.params[2],nt.params[3],false));
                nt.params[5] := When(nt.params[4],nt.params[1],nt.params[2]);
		return [ICScatAcc(fPrecompute(fd))*
			   TTensorI_OL(tt,[APar,AMul(1,APar)],[[1,0]],
                        [nt.params[1]/freedoms[1],nt.params[5]/freedoms[1]])];
                end,
        apply := (nt,c,cnt) ->  let(len1:=2*nt.params[1]*nt.params[2]*nt.params[3],
			            len2:=2*nt.params[1]*nt.params[3],
						When(nt.isSplit(),Split(len1,2)*Prm(L(len1,2))*c[1]*
                                			Cross(Prm(L(len1,len1/2)),Prm(L(len2,len2/2)))*
                                			Cross(Glue(2,len1/2),Glue(2,len2/2)),c[1])),

     )
   )
);

Class(SAR2,TaggedNonTerminal,rec(
       abbrevs := [(n,s) -> Checked(IsPosInt(n),[n,s])],
       dmn := self >> [TArray(TReal,When(self.params[1]=1,2,self.params[1])),TArray(TReal,When(self.params[1]=1,2,self.params[1]))],
       rng := self >> [TArray(TReal,When(self.params[1]=1,2,self.params[1]))],
       isSplit := self >> self.params[2],
       isReal := self >> true,
       transpose := self >> Copy(self)));

NewRulesFor(SAR2,rec(

      SAR_Base := rec (
      	forTransposition := false,
      	applicable := (self,nt) >> (nt.params[1]=1),
      	apply := (nt,c,cnt) -> ICMultiplication(2,2)),

      SAR_MatchFilter := rec (
	forTransposition := false,
        applicable := (self,nt) >> nt.params[1] >= 2,
        freedoms := nt -> [[2]],
        child := (self,nt,freedoms) >> [TTensorI_OL(SAR2(1,nt.params[2]),[APar,AMul(1,APar)],[[1,0]],
                        [nt.params[1]/freedoms[1],nt.params[1]/freedoms[1]])],
        apply := (nt,c,cnt) -> When(nt.isSplit(),Split(nt.params[1],2)*Prm(L(nt.params[1],2))*c[1]*
				Cross(Prm(L(nt.params[1],nt.params[1]/2)),Prm(L(nt.params[1],nt.params[1]/2)))*
				Cross(Glue(2,nt.params[1]/2),Glue(2,nt.params[1]/2)),c[1]),

      ),

      SAR_MatchFilterSMP := rec (
      	forTransposition := false,
 	#applicable := (self,nt) >> nt.params[1] >= 1 and nt.hasTags() and nt.isTag(1,AParSMP),
        applicable := (self,nt) >>  false,
        freedoms := nt -> Cond(nt.isTag(1, AParSMP) and
               nt.params[1] mod nt.firstTag().params[1]=0,
               [[nt.params[1]/nt.firstTag().params[1]]],
               [[1]]),
        child := (self,nt,freedoms) >> [TTensorI_OL(SAR2(freedoms[1],nt.params[2]),[APar,AMul(1,APar)],[[1,0]],
			[nt.params[1]/freedoms[1],nt.params[1]/freedoms[1]]).withTags(nt.getTags())],
        apply := (nt,c,cnt) -> c[1]
     ),

     SAR_MatchFilterVec := rec (
         forTransposition := false,
         #applicable := (self,nt) >> nt.hasTags() and nt.isTag(1,AVecReg) and nt.params[1] > nt.firstTag().v,
         applicable := (self,nt) >> false,
         freedoms := nt -> Cond(nt.isTag(1,AVecReg) and Mod(nt.params[1],nt.firstTag().v) = 0,
			[[nt.firstTag().v]],[[1]]),
         child := (self,nt,freedoms) >> [TTensorI_OL(SAR2(nt.firstTag().v,nt.firstTag()),[AVec,AMul(1,AVec)],[[0,1]],
			 [nt.params[1]/freedoms[1],nt.params[1]/freedoms[1]]).withTags(nt.getTags())],
         apply := (nt,c,cnt) -> c[1]
     ),

    SAR_MatchFilterVecBase := rec (
        forTransposition := false,
        #applicable := (self,nt) >> nt.hasTags() and nt.isTag(1,AVecReg) and nt.params[1] = nt.firstTag().v,
        applicable := (self,nt) >> false,
        apply := (nt,c,cnt) ->
                        VTensor_OL(Multiplication(2,1),nt.firstTag().v)

        )

    )

);


text:
# Codes related to Rationality of Algebraic Tori.
# This file contains all codes written by myself.

Sub_Matrix:= function(G,L,M)
#Inputs: G as a finite group, L list with 2 entries showing the interval of desired rows, M list with 2
# entries showing the interval of desired columns. The output is a list containing extracted submatrices
# from generators of group G.
local g,H;
H:= [];
	for g in GeneratorsOfGroup(G) do
		Append(H,[g{[L[1]..L[2]]}{[M[1]..M[2]]}]);
	od;
return(H);
end;

GDual:= function(G)
# this function returns the dual G-lattice, by transposing the generators of G.

	local i,L,n,T,Gp;

	T:= GeneratorsOfGroup(G);

	n:= Size(T);

	L:=[];
	for i in [1..n] do
		L[i]:= TransposedMat(T[i]);
	od;

	Gp:= GroupByGenerators(L);
	return(Gp);
end;

Rank1PermQuot:= function(G)

local   Denom,		# List of all denominators of entries of a vector.
	SDenom,		# Size of Denom.
	gens,		# A generating list for the rank m-1 subgroup (sublattice
			# of the dual lattice).
	H,		# The Corresponding group of the dual lattice.
	HH,		# A list of Generators of H.
	I,		# Identity matrix of size equal to rank of G
	K,		# An appropriate conjugate of the dual group in which,
			# the invariant sublattice of rank m-1 is readable.
	L,		# The record containing the vectors which extend
			# the sign vector to a basis fo lattice.
	LCM,		# The least common multiple of Denom.
	m,		# The rank of group
	W,		# A list of vectors which is formed by bases of eigenspaces.
	R,		# A mutable copy of extended basis.
	r,		# A mutable copy of the basis for the invariant space.
	S,		# GAP ID of the rank n-1 invariant for the original lattice.
	G1,		# The group asocciated to an invariant subspace of rank m-1 in the original lattice.
	Res;		# Record containing all the information of the reduction process.
	H:= GDual(G);
	HH:= GeneratorsOfGroup(H);
	m:= Size(HH[1]);
	I:= HH[1]^0;
	W:= BaseFixedSpace(HH);

	if Size(W)<>0 then
		if ForAll(W[1], IsInt) = false then
			Denom := List(W[1], w-> DenominatorRat(w));
			LCM:= Denom[1];
			SDenom:= Size(Denom);
			for i in [2..SDenom] do
				LCM:= LcmInt(LCM,Denom[i]);
			od;
			W[1]:= LCM*W[1];
		else
			W[1]:= W[1];
		fi;
		L:= ComplementIntMat(I,[W[1]]);
		R:= ShallowCopy(L.complement);
		r:= ShallowCopy(L.sub);
		Append(R,r);
		if IsInt(Determinant(R)) = true and AbsInt(Determinant(R))=1 then
			K:= ConjugateGroup(H,R^(-1));
			gens:= Sub_Matrix(K,[1,m-1],[1,m-1]);
			G1:= GDual(GroupByGenerators(gens));
			S:= CrystCatZClass(G1);
			Res:= rec(Structure:= StructureDescription(G), DualGens:= HH, Trans:= R^(-1), 		 				NewG:=GeneratorsOfGroup(GDual(K)),SubLatGens:= GeneratorsOfGroup(G1), ZClassSubLat:= S);
			return(S);
		else
			return(fail);
		fi;
	else
		#Print("No Permutation lattice of Rank 1\n");
		return(fail);

	fi;
end;

###### For a given lattice (represented by finite number of unimodular matrices)
# The function calculates intersections of Eigenspaces for eigenvalues one and
# negative one. In Case that there is a nonzero intersection, it returns the basis.
# Input:  a finite subgroup of GL(n,Z).
# Output: The conjugation element for the dual lattice which is a matrix in GL(n,Z).
#	  An invariant sublattice of the corresponding lattice to G, of rank 4 where
#	  the quotient of the original lattice and the output is a sign permutation.

Rank1SgnPermQuot:= function(G)

local
	BEigSP,		# Duplicate free list of bases of intersection of eigenspaces.
	EigSp,		# List of all eigenspaces of all generators of G.
	EigSP,		# List of Cartesian product of EigSp.
	Denom,		# List of all denominators of entries of a vector.
	gens,		# A generating list for the rank 4 subgroup (invariant sublattice
			# of the dual lattice).
	H,		# The Corresponding group of the dual lattice.
	HH,		# A list of Generators of H.
	I,		# Identity Matrix of rank 5.
	InterEigSP,	# A list of Intersections of eigenspaces.
	K,		# An appropriate conjugate of the dual group in which,
			# the invariant sublattice of rank 4 is easily readable.
	l,		# rank of lattice
	L,		# The record containing the vectors which extend
			# the sign vector to a basis fo lattice.
	LCM,		# The least common multiple of Denom.
	W,		# A list of vectors which is formed by bases of eigenspaces.
	R,		# A mutable copy of extended basis.
	Rec,
	r,		# A mutable copy of the basis for the sign invariant space.
	S,		# GAP ID of the rank n-1 invariant for the original lattice.
	Sgen,		# The size of a generating set of H.
	G1;		# An invariant subspace of rank 4 in the original lattice.

	H:= GDual(G);
	HH:= GeneratorsOfGroup(H);
	I:= HH[1]^0;
	l:= Size(I);
	Sgen:= Size(HH);
	EigSp:= List([1..Sgen],x->[]);
	for i in [1..Sgen] do
		Append(EigSp[i],Eigenspaces(Rationals,HH[i]));
	od;
	EigSP:= Cartesian(EigSp);
	InterEigSP:= [];
	for a in EigSP do
		n:= Size(a);
		S:= a[1];
		for i in [2..n] do
			S:= Intersection(S,a[i]);
		od;
		Append(InterEigSP,[S]);
	od;
	BEigSP:= Set(List(InterEigSP,t->GeneratorsOfVectorSpace(t)));
	Remove(BEigSP,Position(BEigSP,[]));
	if BEigSP <> [] then
		W:= ShallowCopy(BEigSP[1]);

		if Size(W)<>[] then
			if ForAll(W[1], IsInt) = false then
				Denom := List(W[1], w-> DenominatorRat(w));
				LCM:= Denom[1];
				for i in [2..Size(Denom)] do
					LCM:= LcmInt(LCM,Denom[i]);
				od;
				W[1]:= LCM*W[1];
			else
				W[1]:= W[1];
			fi;
			L:= ComplementIntMat(I,W);
			R:= ShallowCopy(L.complement);
			r:= ShallowCopy(L.sub);
			Append(R,r);

				K:= ConjugateGroup(H,R^(-1));
				gens:= Sub_Matrix(K,[1,l-1],[1,l-1]);
				G1:= GDual(GroupByGenerators(gens));
				S:= CrystCatZClass(G1);
				Print(R^(-1),"\n");
				Rec:= rec(DualGens:= HH, Trans:= R^(-1), NewG:= GeneratorsOfGroup(GDual(K)),
                                         SubLatGens:= GeneratorsOfGroup(G1), ZClassSubLat:= S);
					return(S);
		else
			Print("No sign permutation lattice of rank 1\n");

		fi;
	else
		Print("No sign permutation lattice of rank 1\n");
	fi;
end;

ext: gap

max_stars_repo_name: ajamshid/Algebraic-Tori

max_issues_repo_name: armin-jamshidpey/Algebraic-Tori

max_forks_repo_name: armin-jamshidpey/Algebraic-Tori


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

# Configuration of system dependent parameters
# ============================================

#F Global code generation options record
#F
SpiralDefaults := rec();

SPL_DEFAULTS := SpiralDefaults;


text:
# Creating a random source
rs := RandomSource(IsMersenneTwister);

# Generate a random number between 1 and 10
Random(rs, 1, 10);

# Same with default random source
Random(1, 10);

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


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

#F ResolveMemos(<o>)
#F   This function removes all precomputed variables (with Data(v, expr, ...))
#F   by inserting the corresponding expressions
#F
ResolveMemos := o ->
   SubstTopDown(o, @(1,var,e->IsBound(e.mapping)), e->e.mapping);

_GenerateData := function(map, free, lambda)
    local newlambda, v, val, d, free1;
    if free = [] then
	newlambda := Copy(lambda);
	# NOTE: hardcoded [var, param] is ugly
	newlambda.expr := SubstParamsCustom(newlambda.expr, map, [var, param]);
	return newlambda.tolist();
    else
	v := free[1];
	d := [];
	free1 := Drop(free, 1);
	for val in [0..v.range-1] do
	    map.(v.id) := V(val);
	    Append(d, _GenerateData(map, free1, lambda));
	od;
	return d;
    fi;
end;

#Declare(fComputeOnline);

# PrecalculatedData - accumulating precalculated data variables so
# same constant data table (variable) could be reused in different places in code.
# Beware: search for table by key is inefficient here.

Class(PrecalculatedData, rec(

     # _tables is a list of rec(key, data, opts) items.
    _tables := [],

    __call__ := (self) >> self._tables,

    # keyForFunc - returns search key for function 'func' that has 'var' variables.
    # The key is just the same function with unnamed variables.

    keyForFunc := meth( self, func, vars)
        local i, j, map;
        map := rec();
        i := 1; j := Length(vars);
        while j>0 do
            map.(vars[i].id) := ind(vars[i].range, j);
            i := i+1; j := j-1;
        od;
	# NOTE: hardcoded [var, param] is ugly
        return SubstParamsCustom( Copy(func), map, [var, param]);
    end,

    # find(key, opts) - searching if data table already exists.
    find := meth( self, key, opts)
        local item;
        for item in self._tables do
            #NOTE: opts.dataType - just to distinguish single precision from double
            if item.key=key and item.opts.dataType=opts.dataType then
                return item.data;
            fi;
        od;
        return false;
    end,

    # add - remember data table assosiated with key, data is FData object.
    add := meth( self, key, data, opts)
        Add(self._tables, rec( key := key, data := data, opts := opts));
    end,

    clear := meth(self)
        self._tables := [];
    end,

));

GenerateData := function ( func, opts )
    local f, n, free, lambda, data, t, olvars, key, map, i;

    if ObjId(func)=FDataOfs then return func; fi;

    olvars := [];
    map := rec();

    func := SubstTopDown(Copy(func), fComputeOnline,
        function(s)
            local v;
            v := var.fresh_t("o", s._children[1].t);
            v.range := s._children[1].range();
            Add(olvars, v);
            map.(v.id) := s._children[1];
            return v;
        end
    );

    lambda := func.lambda();
    # NOTE: Collect is a hack
    # NOTE: hardcoded 'param' is ugly
    free := Difference(lambda.expr.free() :: Collect(lambda.expr, param), lambda.vars);
    free := Filtered(free, IsLoopIndex);

    for i in free do map.(i.id) := fBase(i); od;

    Append(free, olvars);

    Constraint(ForAll(free, v->IsBound(v.range)));
    Sort(free, (x,y) -> DoubleString(Drop(x.id,1)) <= DoubleString(Drop(y.id,1)));

    if IsBound(opts.dataSharing) and opts.dataSharing.enabled then
        key := PrecalculatedData.keyForFunc(func, free);
        data := PrecalculatedData.find(key, opts);
        if ObjId(data)=false then
            data := FData(_GenerateData(tab(), free, lambda));
            PrecalculatedData.add(key, data, opts);
        fi;
    else data := FData(_GenerateData(tab(), free, lambda));
    fi;

    n := lambda.domain();

    if (free=[] and olvars=[]) then return data.part(n,0);
    else return data.part(n, n * fTensor(List(free, i->map.(i.id))).at(0));
    fi;
end;

_GenerateInitCode := function(data, N, idx, free, lambda)
    local i, v, res;
    if free = [] then
	i := lambda.vars[1];
	res := loop(i, lambda.domain(),
		assign(nth(data, idx+i),
		       When(data.t.t = lambda.expr.t,
			        lambda.expr,
			    # else
				tcast(data.t.t, lambda.expr))));

        # ObjId<>loop, when i.range was 1, loop just returned its body
	if i.range <= 16 and ObjId(res)=loop then return SReduceSimple(res.unroll());
	else return SReduceSimple(res);
	fi;
    else
	v := free[1];
        # careful here, since N and v.range can be symbols we have to use div ,
        # since N/v.range will become an fdiv (floating point division) in the code,
        # div = integer division with divisible arguments (allows normal simplification as in fdiv)
	N := div(N, v.range);
	return loop(v, v.range,
	    _GenerateInitCode(data, N, idx+v*N, Drop(free,1), lambda));
    fi;
end;

GenerateInitCode := function ( func, data_container, opts )
    local f, N, n, free, lambda, data, map, initcode, res, range, key;
    if ObjId(func)=FDataOfs then return func; fi;
    lambda := func.lambda();

    # NOTE: Collect is a hack
    # NOTE: hardcoded 'param' is ugly
    free := Difference(lambda.expr.free() :: Collect(lambda.expr, param), lambda.vars);
    free := Filtered(free, IsLoopIndex);

    Constraint(ForAll(free, v->IsBound(v.range)));
    Sort(free, (x,y) -> DoubleString(Drop(x.id,1)) <= DoubleString(Drop(y.id,1)));

    n := func.domain();
    N := n * Product(free, x->x.range);

    range := func.range();

    if IsBound(opts.dataSharing) and opts.dataSharing.enabled then
        key := PrecalculatedData.keyForFunc(func, free);
        data := PrecalculatedData.find(key, opts);
    else
        data := false;
    fi;

    if data=false then
        if IsInt(range) or IsSymbolic(range) then
            data := data_container(TInt, N); # index mapping function, not diag
        else
	    if opts.unifyStoredDataType = "input" then
		range := UnifyTypes([range, opts.XType.t]);
	    elif opts.unifyStoredDataType = "output" then
		range := UnifyTypes([range, opts.YType.t]);
	    elif IsType(opts.unifyStoredDataType) then
		range := UnifyTypes([range, opts.unifyStoredDataType]);
	    fi;
            data := data_container(range, N);
        fi;
        data.init := _GenerateInitCode(data, N, 0, free, lambda);

        if IsBound(opts.dataSharing) and opts.dataSharing.enabled then
            PrecalculatedData.add(key, data, opts);
        fi;
    fi;

    if free=[] then
        res := FData(data).part(n,0);
    else
        res := FData(data).part(n, n * fTensor(List(free, fBase)).at(0));
    fi;
    res := res.setRange(range);
    res.init := res.var.init;

    return res;
end;

_Process_fPrecompute := (o, precompute_func, precond) -> SubstBottomUpRules(o, [
    [ [fPrecompute, @(1)],
      (e,cx) -> When(precond(e,cx) and e.rank()=0, precompute_func(ResolveMemos(e)), e),
      "ContextDependentLocalize_fPrecompute"],
]);

# Process_fPrecompute(<sums>, <opts>)
#
Process_fPrecompute := function(sums, opts)
    local bb, b, has_free_vars, precompute_func;

    # Autolib needs ability to disable fPrecompute when it uses strategies from LocalConfig.bench
    if IsBound(opts.disablePrecompute) and opts.disablePrecompute then
        return sums;
    fi;

    if IsBound(opts.dataSharing) and opts.dataSharing.enabled then
    	sums := ApplyStrategy(sums, opts.dataSharing.precomputeStrategy, UntilDone, opts);
    fi;

    bb := Collect(sums, @(1, [BB, RecursStep]));
    for b in bb do
        b.precomputed_free := b.free();
    od;

    has_free_vars := (e, cxentry) ->
        (CollectNR(e, param) <> []) or
        (Intersection(e.free(), cxentry[1].precomputed_free)<>[]);

    precompute_func := When(opts.generateInitFunc,
        e -> GenerateInitCode(e, Dat1d, opts),
        e -> GenerateData(e, opts));

    sums := _Process_fPrecompute(sums, precompute_func,
           (e, cx) -> (
                        ((not IsBound(cx.BB) or cx.BB=[]) and
                        (not IsBound(cx.RecursStep) or cx.RecursStep=[]) and
                        (not IsBound(cx.RS) or cx.RS=[]) and
                        (not IsBound(cx.RSBase) or cx.RSBase=[])) or
                        (IsBound(cx.BB) and cx.BB<>[] and has_free_vars(e, cx.BB)) or
	                    (IsBound(cx.RecursStep) and cx.RecursStep<>[] and has_free_vars(e, cx.RecursStep))
                    )
                    );
    return sums;
end;