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

# Higher-order functions support
#

# Functional expression base class
Class(FuncExp, Exp, rec(
    isFuncExp := true,

    at := arg >> let(
	self := arg[1],
	Cond(Length(arg)=2 and IsList(arg[2]),
	     ApplyFunc(fcall, [self] :: arg[2]),
	     ApplyFunc(fcall, arg))),

    mkVars := self >>
        List(DropLast(self.computeType().params, 1), x->var.fresh_t("q", x)),

    getRankType := (self, r) >> let(
	t := self.computeType(),
	rank := Length(t.params)-2,
	Cond(r > rank,
	         TDummy,
	     t.params[ Length(t.params) - 1 - r ])
    )
));

Function.mkVars      := FuncExp.mkVars;
Function.getRankType := FuncExp.getRankType;

Symbolic.mkVars      := FuncExp.mkVars;
Symbolic.getRankType := FuncExp.getRankType;

IsFuncExp := o -> IsRec(o) and
    ((IsBound(o.isFuncExp) and o.isFuncExp) or (IsBound(o.t) and ObjId(o.t)=TFunc));

# Same as IsFuncExp(<o>) but returns false for param's with t = TFunc
# (these are considered 'trivial')
IsNonTrivFuncExp := o -> IsRec(o) and (IsBound(o.isFuncExp) and o.isFuncExp);

# !!! Below is no longer needed/used, autolib handles this via Parametrizer
# These are hacks for dealing with types which are integers like 2
# these really mean the range of an integer variable
# _convType converts 2 to TInt not to confuse other parts of the system
#

# !!! Below is no longer needed/used, autolib handles this via Parametrizer
#_convType := t -> Cond(IsSymbolic(t) or IsInt(t) or IsValue(t), TInt,
#    SubstBottomUp(Copy(t), TFunc,
#        e -> ApplyFunc(TFunc, List(e.params, p -> When(IsInt(p) or IsSymbolic(p) or IsValue(p), TInt, p)))));
_convType := t->t;

# !!! Below is no longer needed/used, autolib handles this via Parametrizer
#_convTypeExp := t ->
#    SubstBottomUp(Copy(t), TFunc,
#        e -> ApplyFunc(TFunc, List(e.params, p -> When(IsInt(p) or IsSymbolic(p) or IsValue(p), TInt, p))));
_convTypeExp := t->t;

#F fcurry(<func>, <pos>, <arg>) - symbolic representation of function currying
#F     Given a function with n arguments returns a functino with n-1 arguments
#F     where <arg> is plugged into position <pos>
#F
#F Example: p := param(TFunc(TInt, 4, TComplex), "p");  # ie. rank-1 4-elt diagonal func
#F          i1 := Ind(4);
#F          fcurry(p, 1, i1).t;
#F            => TFunc(4, TComplex);
#F
Class(fcurry, FuncExp, rec(
    computeType := self >> let(
        n := self.args[2],
        ct := self.args[1].computeType(),
        Checked(IsValue(n), ObjId(ct)=TFunc, n.v < Length(ct.params),
            ApplyFunc(TFunc, ListWithout(ct.params, n.v)))),

    eval := self >> let(
        torig := self.args[1].t.params,
        tcurried := self.t.params,
        n := self.args[2].v,
        vars := List([1..Length(tcurried)-1], i -> var.fresh_t("q", tcurried[i])),
        plugin := List([1..Length(torig)-1], i -> Cond(i<n, vars[i], i=n, self.args[3],vars[i-1])),
        Lambda(vars, ApplyFunc(fcall, Concatenation([self.args[1].eval()], plugin)))),

    at := (self, i) >> self.eval().at(i)
));

#F flift(<func>, <t>) - symbolic representation of function "lifting"
#F
#F  Given a function with n arguments returns a function with n+1 arguments,
#F  extra argument is of type <t> and is the "one but last" argument of new function.
#F  It is the ignored when computing the value of the function.
#F
#F  The purpose of flift is to increase the rank of functions. It implicitly creates
#F  a new inner loop, which loop variable is ignored. This is needed for rewrite rules
#F  like GT(T, ..) * Gath(f) -> GT(T*Gath(f), ...), where f is pulled into the loop.
#F
Class(flift, FuncExp, rec(
    computeType := self >> let(ct := self.args[1].computeType(), Checked(
        ObjId(ct)=TFunc, IsType(self.args[2]),
        let(tt := ct.params, n := Length(tt),
        ApplyFunc(TFunc, Concatenation(tt{[1..n-2]}, [self.args[2]], tt{[n-1..n]}))))),

    # NOTE: define what exactly .lambda() does in these cases
    lambda := self >> self.eval(),

    eval := self >> let(vars := List(DropLast(self.t.params,1), x->var.fresh_t("q", x)), n := Length(vars),
        Lambda(vars, ApplyFunc(fcall, Concatenation([self.args[1].eval()], ListWithout(vars, n-1)))))
));

#F fsplit(<func>, <loopid>, <inner_its>, <outer_its>)
Class(fsplit, FuncExp, rec(
    computeType := self >> let(ct := self.args[1].computeType(), Checked(
        ObjId(ct)=TFunc, IsValue(self.args[2]), self.args[2].t=TInt,
        let(tt := ct.params, loopid := self.args[2].v, pos := Length(tt)-1-loopid,
            ApplyFunc(TFunc, Concatenation(tt{[1..pos-1]}, [tt[pos], tt[pos]], tt{[pos+1..Length(tt)]}))))),

    # NOTE: define what exactly .lambda() does in these cases
    lambda := self >> self.eval(),

    eval := self >> let(
        vars      := List(DropLast(self.computeType().params, 1), x->var.fresh_t("q", x)),
        loopid    := self.args[2].v,
        inner_its := self.args[3],
        pos       := Length(vars)-loopid,
        callargs  := vars{[1..pos-2]} ::
	             Cond(vars[pos].t=TDummy, [0], [inner_its * vars[pos-1] + vars[pos]]) ::
                     vars{[pos+1..Length(vars)]},
        Lambda(vars, ApplyFunc(fcall, Concatenation([self.args[1].eval()], callargs))))
));

#F frotate(<func>, <n>)
#F   switch the order of loops (=ranks), by making <n>-th loop innermost
Class(frotate, FuncExp, rec(
    computeType := self >> let(ct := self.args[1].computeType(), nn := self.args[2], Checked(
        ObjId(ct)=TFunc, IsValue(nn), nn.t=TInt,
	let(tt := ct.params, rank := Length(tt)-2, n := nn.v, pos := rank+1-n,
            Cond(rank=0 or (rank=1 and n=1), ct,
		 n > rank,
		            ApplyFunc(TFunc, tt{[1..rank]} :: [TDummy] :: [tt[rank+1], tt[rank+2]]),
		 # else
                            ApplyFunc(TFunc, tt{[1..pos-1]} :: tt{[pos+1..rank]} :: [tt[pos], tt[rank+1], tt[rank+2]]))))),

    # NOTE: define what exactly .lambda() does in these cases
    lambda := self >> self.eval(),

    eval := self >> let(f := self.args[1], rank := f.rank(), n := self.args[2].v, pos := rank+1-n,
	Cond(rank=0 or (rank=1 and n=1), f,
	     n > rank,
             let(vars := List(DropLast(self.computeType().params, 1), x->var.fresh_t("q", x)),
		 Lambda(vars, ApplyFunc(fcall, [f] :: vars{[Length(vars)-1-rank..Length(vars)-2]} :: [Last(vars)]))),
	     # else
             let(vars := List(DropLast(self.computeType().params, 1), x->var.fresh_t("q", x)),
		 Lambda(vars, ApplyFunc(fcall, [f] :: vars{[1..pos-1]} :: [vars[rank]] :: vars{[pos..rank-1]} :: [vars[rank+1]])))
	))
));

## Ranked functions support
## Ranked functions == functions with implicit dependencies on loop variables
##
## These are used in paradigms.common.GT and autolib.*
##
## NOTE: get rid of _, these functions are not private, but public exports
##
_rankManip := (obj, newfunc) >>
    SubstTopDownNR(Copy(obj), @.cond(e->IsFunction(e) or IsFuncExp(e)), e -> newfunc(e));

_rank     := o -> Cond(
    IsList(o), Maximum0(List(o, _rank)),
    not IsRec(o) or not IsBound(o.rank), 0, o.rank());
_upRank   := o -> Cond(
    IsList(o), List(o, _upRank),
    not IsRec(o) or not IsBound(o.upRank), o, o.upRank());
_upRankNeq0   := o -> Cond(
    IsList(o), List(o, _upRankNeq0),
    not IsRec(o) or not IsBound(o.upRank) or o.rank()=0, o, o.upRank());
_upRankBy := (o,n) -> Cond(IsList(o), # NOTE: =0 ??
    List(o, x->_upRankBy(x,n)),
    not IsRec(o) or not IsBound(o.upRankBy) or o.rank()=0, o, o.upRankBy(n));
_downRank := (o,loopid,ind) -> Cond(
    IsList(o), List(o, x->_downRank(x, loopid, ind)),
    not IsRec(o) or not IsBound(o.downRank), o, o.downRank(loopid, ind));
_downRankFull := (o,inds) -> Cond(
    IsList(o), List(o, x->_downRankFull(x, inds)),
    not IsRec(o) or not IsBound(o.downRankFull), o, o.downRankFull(inds));
_split := (o, loopid, iits, oits) -> Cond(
    IsList(o), List(o, x->_split(x, loopid, iits, oits)),
    not IsRec(o) or not IsBound(o.split), o, o.split(loopid, iits, oits));
_rotate := (o, n) -> Cond(
    IsList(o), List(o, x->_rotate(x, n)),
    not IsRec(o) or not IsBound(o.rotate), o, o.rotate(n));

_rch_rank      := self >> Maximum0(List(self.rChildren(), _rank));
_rch_upRank    := self >> self.from_rChildren(List(self.rChildren(), _upRank));
_rch_upRankBy  := (self, n) >>  self.from_rChildren(List(self.rChildren(), c->_upRankBy(c,n)));
_rch_downRank  := (self, loopid, ind) >>
    self.from_rChildren(List(self.rChildren(), c->_downRank(c,loopid,ind)));
_rch_downRankFull := (self, inds) >>
    self.from_rChildren(List(self.rChildren(), c->_downRankFull(c,inds)));
_rch_split     := (self, loopid, iits, oits) >>
    self.from_rChildren(List(self.rChildren(), c->_split(c,loopid, iits, oits)));
_rch_rotate    := (self, n) >>
    self.from_rChildren(List(self.rChildren(), c->_rotate(c, n)));

Symbolic.domain := self >> Checked(ObjId(self.t)=TFunc, Length(self.t.params) > 1,
    self.computeType().params[Length(self.t.params)-1]);

Symbolic.range := self >> Checked(ObjId(self.t)=TFunc, Length(self.t.params) > 1,
    self.computeType().params[Length(self.t.params)]);

Symbolic.rank := self >> Cond(ObjId(self.t)=TFunc and Length(self.t.params) > 1,
    Length(self.t.params) - 2,
    _rch_rank(self));

Symbolic.at := (self, vars) >> Checked(ObjId(self.t)=TFunc, self.lambda().at(vars));

Symbolic.lambda := self >> Checked(ObjId(self.t)=TFunc, let(
    selft := self.computeType(),
    vars := List([1..self.rank()+1], x -> let(t:=selft.params[x],
            Cond(IsType(t), var.fresh_t("w", t), var.fresh("w", TInt, t)))),
    Lambda(vars, ApplyFunc(fcall, Concatenation([self], vars)))));

Symbolic.upRank := self >> Cond(ObjId(self.t)<>TFunc, _rch_upRank(self), flift(self, TDummy));

Symbolic.upRankBy := (self, n) >> Cond(
    ObjId(self.t)<>TFunc, _rch_upRankBy(self,n),
    Checked(IsPosInt0(n), FoldL([1..n], (f, i) -> f.upRank(), self)));

Symbolic.downRankFull := (self, inds) >> Cond(
    ObjId(self.t)<>TFunc, _rch_downRankFull(self, inds),
    FoldL(Reversed([1..Minimum(self.rank(), Length(inds))]), (f, i) -> f.downRank(i, inds[i]), self));

Symbolic.downRank := (self, loopid, ind) >> let(rank := self.rank(), Cond(
    loopid > rank, self,
    ObjId(self.t)<>TFunc, _rch_downRank(self, loopid, ind),
    fcurry(self, rank+1-loopid, ind)));  # !! loop variables are ordered with decreasing rank, highest-rank (outermost) is first var

Symbolic.split := (self, loopid, inner_its, outer_its) >> Cond(
    loopid > self.rank(), self,
    ObjId(self.t)<>TFunc, _rch_split(self, loopid, inner_its, outer_its),
    fsplit(self, loopid, inner_its, outer_its));

Symbolic.rotate := (self, n) >> Cond(
    self.rank() <= 1 and n <= 1, self,
    ObjId(self.t)<>TFunc, _rch_rotate(self, n),
    frotate(self, n));

# NOTE: below is really a hack, since the first argument might also have ranked
# things, although this is not obvious at first glance, an example of such first
# argument would be fcurry(func, lambdaWrap(rank-n function))

fcall.rank := self >> Cond(ObjId(self.t)=TFunc and Length(self.t.params) > 1,
    Length(self.t.params) - 2,
    _rank(Drop(self.args, 1)));

fcall.upRank := self >> Cond(ObjId(self.t)<>TFunc,
    ApplyFunc(fcall, [self.args[1]] :: _upRank(Drop(self.args, 1))),
    flift(self, TDummy));

fcall.upRankBy := (self, n) >> Cond(ObjId(self.t)<>TFunc,
    ApplyFunc(fcall, [self.args[1]] :: _upRankBy(Drop(self.args, 1), n)),
    Checked(IsPosInt0(n), FoldL([1..n], (f, i) -> f.upRank(), self)));

fcall.downRankFull := (self, inds) >> Cond(ObjId(self.t)<>TFunc,
    ApplyFunc(fcall, [self.args[1]] :: _downRankFull(Drop(self.args, 1), inds)),
    FoldL(Reversed([1..Minimum(self.rank(), Length(inds))]), (f, i) -> f.downRank(i, inds[i]), self));

fcall.downRank := (self, loopid, ind) >> let(rank := self.rank(), Cond(
    loopid > rank, self,
    ObjId(self.t) <> TFunc,
        ApplyFunc(fcall, [self.args[1]] :: _downRank(Drop(self.args, 1), loopid, ind)),
    # else
    fcurry(self, rank+1-loopid, ind)));  # !! loop variables are ordered with decreasing rank, highest-rank (outermost) is first var

fcall.split := (self, loopid, inner_its, outer_its) >> Cond(
    loopid > self.rank(), self,
    ObjId(self.t) <> TFunc,
        ApplyFunc(fcall, [self.args[1]] :: _split(Drop(self.args, 1), loopid, inner_its, outer_its)),
    fsplit(self, loopid, inner_its, outer_its));

fcall.rotate := (self, n) >> Cond(
    self.rank() <= 1 and n <= 1, self,
    ObjId(self.t) <> TFunc,
        ApplyFunc(fcall, [self.args[1]] :: _rotate(Drop(self.args, 1), n)),
    frotate(self, n));

Function.rank      := _rch_rank;
Function.upRank    := _rch_upRank;
Function.upRankBy  := _rch_upRankBy;
Function.downRank  := _rch_downRank;
Function.downRankFull := _rch_downRankFull;
Function.split     := _rch_split;
Function.rotate    := _rch_rotate;
Function.computeType := self >> self.lambda().t;

Lambda.rank      := Symbolic.rank;
Lambda.upRank    := Symbolic.upRank;
Lambda.upRankBy  := Symbolic.upRankBy;
#Lambda.downRank  := Symbolic.downRank;  downRank now defined in lambda.gi
Lambda.downRankFull := Symbolic.downRankFull;
Lambda.split     := Symbolic.split;
Lambda.rotate    := Symbolic.rotate;

# ind(<range>, <n>) - "nameless" reference to a loop index of n-th inner most loop
#            (eg. ind(1) is inner most, ind(n) is outermost in n-loop nest)
#            loop counter runs from 0..range-1
Class(ind, Loc, rec(
    __call__ := (self, range, n) >> Checked(IsInt(n),
	WithBases(self, rec(operations := ExpOps, range:=range, n:=n))),
    print := self >> Print(self.name, "(", self.range, ", ", self.n, ")"),
    rChildren := self >> [self.range, self.n],
    rSetChild := rSetChildFields("range", "n"),
    t := TInt,
    eval := self >> self,
    can_fold := False,

    upRank := self >> ObjId(self)(self.range, self.n+1),

    split := (self, loopid, inner_its, outer_its) >>
        Cond( loopid > self.n, self,
              loopid < self.n, ObjId(self)(self.range, self.n+1),
              ObjId(self)(inner_its, self.n) + inner_its*ObjId(self)(outer_its, self.n+1)),
    rotate := (self, n) >>
        Cond( n < self.n, self,
              n > self.n, ObjId(self)(self.range, self.n+1),
              ObjId(self)(self.range, 1)),

));

# NOTE: ind.downRank might be a hack
ind.downRankFull := (self, inds) >> inds[self.n];
ind.downRank := (self, loopid, ind) >> Cond(loopid=self.n, ind, self);
ind.rank := self >> self.n;

_hofnew := true;
# ExpMarkActiveRank(<s>)
#   performs a recursive walk over <s> and sets ._expMarkActiveRank attribute to the "active rank"
#   of each subexpression.
#
#   Active rank of an expression denotes the maximum implicit loop id that the expression refers to.
#   Implicit loop id's are introduced by objects such as GT and Lambda.
#
ExpMarkActiveRank := function(s)
    local c, rch, rank;
    rch := Cond(IsRec(s) and IsBound(s.rChildren), s.rChildren(), IsList(s) and not IsString(s), s, []);
    if ObjId(s) = ind then
        rank := s.n;
    else
	if _hofnew then
	    rank := _rank(s);
	    DoForAll(rch, ExpMarkActiveRank);
	else
# this was invalid (fcurry, fcall, etc) ->
	    rank := Maximum(_rank(s), Maximum0(List(rch, ExpMarkActiveRank)));
	fi;
    fi;

    if IsRec(s) and IsSymbolic(s) then s._expMarkActiveRank := rank; fi;
    return rank;
end;

_ExpMarkPassiveRank := function(s, parent_rank)
    local c, rch, rank, my_rank;
    rch := Cond(IsRec(s) and IsBound(s.rChildren), s.rChildren(), IsList(s) and not IsString(s), s, []);
    my_rank := _rank(s);

    # NOTE: this is a terrible hack! WHAT ABOUT Lambda?
    if spiral.paradigms.common.IsGT(s)        then rank := my_rank + parent_rank;
    elif ObjId(s)=ind then rank := Maximum(parent_rank, s.n);
    else                   rank := Maximum(parent_rank, my_rank);
    fi;
    if IsRec(s) and (IsSymbolic(s) or IsFuncExp(s)) then s._expMarkPassiveRank := rank; fi;

    DoForAll(rch, x -> _ExpMarkPassiveRank(x, rank));
    return rank;
end;
# ExpMarkPassiveRank(<s>)
#   performs a recursive walk over <s> and sets ._expMarkPassiveRank attribute to the "passive rank"
#   of each subexpression.
#
#   Passive rank of an expression denotes the maximum implicit loop id that is defined in the expression.
#   Regardless of whether expression refers to it or now. In contracst, an "active rank" is the loop id
#   that is actually referred to.
#
#   Implicit loop id's are introduced by objects such as GT and Lambda.
#
ExpMarkPassiveRank := s -> _ExpMarkPassiveRank(s, 0);

Class(lambdaWrap, Exp, rec(
    computeType := self >> Checked(ObjId(self.args[1].t) = TFunc, Last(self.args[1].t.params))
));


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

RulesFor(DFT, rec(
    #F DFT_DFT1and3 : DFT1_2n -> L (DFT1_n dirsum DFT3_n) L (I2 tensor F2) L
    #F Derived using polynomial factorization:
    #F    x^2n -> (x^n-1) (x^n+1)
    #F
    DFT_DFT1and3 := rec(
        isApplicable := P -> P[1] > 2 and P[1] mod 2 = 0,
        allChildren := P -> [[ DFT1(P[1]/2, P[2]), DFT3(P[1]/2, P[2]) ]],
        rule := (P,C) -> let(n := P[1]/2,
            L(2*n, n) * DirectSum(C[1], C[2]) * L(2*n, 2) *
            Tensor(I(n), F(2)) *
            L(2*n, n))
    )
));

lperm := p -> Z(p, (p-1)/2);
rperm := p -> Z(p, (p+1)/2);
pdiag := (p,k) -> When(IsEvenInt(k),
    I(p),
    Diag(List([0..p-1], i -> (-1)^i)));

RulesFor(DFT2, rec(
    DFT2_CT := Inherit(DFT_CT, rec(
	switch := false,
	allChildren  := P -> Map2(DivisorPairs(P[1]),
	    (m,n) -> [ DFT2(m, P[2]), DFT(n, P[2]) ]),

	rule := (P,C) -> let(mn := P[1], m := Rows(C[1]), n := Rows(C[2]),
	    Tensor(C[1], I(n)) * Diag(fPrecompute(Tw2(mn, n, P[2]))) *
	    Tensor(I(m), C[2]) * L(mn, m))
    )),
));

NewRulesFor(DFT2, rec(
    DFT2_PF := rec(
    	info         := "DFT2_n*k -> diag * perm * (DFT2_n_a tensor DFT2_k_b) * perm",
    	forTransposition := true,
    	applicable     := nt -> nt.params[1] > 2 and DivisorPairsRP(nt.params[1]) <> [] and not nt.hasTags(),
    	allChildren := nt -> let(N := nt.params[1], k:= nt.params[2], Map2(DivisorPairsRP(N),
    		(r,s) -> [ DFT2(r, 1/s mod r), DFT2(s, 1/r mod s) ])),

        apply := (nt, C, cnt) -> let(
    	    N := nt.params[1],
    	    r := Rows(C[1]),
    	    s := Rows(C[2]),
    	    alpha := 1 / s mod r,
    	    beta  := 1 / r mod s,
    	    i := Ind(r*s),

    	    CRT(r,s).transpose() *
    	    Diag(Lambda(i, -sign(imod(s*alpha*idiv(i,s) + r*beta*imod(i,s), 2*N) - N))) *
    	    Tensor(C[1], C[2]) * CRT(r,s)
    	)
#D    	isApplicable     := P -> P[1] > 2 and DivisorPairsRP(P[1]) <> [] and Length(P[3]) = 0,
#D    	allChildren := P -> let(N := P[1], k:= P[2], Map2(DivisorPairsRP(N),
#D    		(r,s) -> [ DFT2(r, 1/s mod r), DFT2(s, 1/r mod s) ])),
#D
#D    	rule := (P,C) -> let(
#D    	    N := P[1],
#D    	    r := Rows(C[1]),
#D    	    s := Rows(C[2]),
#D    	    alpha := 1 / s mod r,
#D    	    beta  := 1 / r mod s,
#D    	    i := Ind(r*s),
#D
#D    	    CRT(r,s).transpose() *
#D    	    Diag(Lambda(i, -sign(imod(s*alpha*idiv(i,s) + r*beta*imod(i,s), 2*N) - N))) *
#D    	    Tensor(C[1], C[2]) * CRT(r,s)
#D    	)
    )
));

RulesFor(DFT3, rec(
    DFT3_Base := BaseRule(DFT3, [2, ...]),

    DFT3_CT := Inherit(DFT_CT, rec(
	switch := false,
	allChildren  := P -> Map2(DivisorPairs(P[1]),
	    (m,n) -> [ DFT(m, P[2]), DFT3(n, P[2]) ]),

	rule := (P,C) -> let(mn := P[1], m := Rows(C[1]), n := Rows(C[2]),
	    Tensor(C[1], I(n)) * Diag(fPrecompute(Tw3(mn, n, P[2]))) *
	    Tensor(I(m), C[2]) * L(mn, m))
    )),

    DFT3_CT_Radix2 := rec(
	isApplicable := P -> P[1] > 2 and (P[1] mod 2) = 0,
	allChildren  := P -> [[ DFT(P[1]/2, P[2]) ]],
	rule := (P,C) -> let(n := P[1]/2,
	    Tensor(C[1], I(2)) * Diag(fPrecompute(Tw3(2*n, 2, P[2]))) *
	    Tensor(I(n), F(2)*Diag(1,E(4)^P[2])) * L(2*n, n))
    ),

    DFT3_OddToDFT1 := rec(
	isApplicable := P -> IsOddInt(P[1]),
	allChildren := P -> [[ DFT(P[1], P[2]) ]],
	rule := (P, C) -> rperm(P[1]) * C[1] * pdiag(P[1], P[2])
    ),

    DFT3_2xOddToDFT1 := rec(
	isApplicable := P -> P[1] mod 4 = 2 and P[1] > 2,
	allChildren := P -> [[ DFT(P[1], P[2]) ]],
	rule := (P, C) -> let(n:=P[1], k:=P[2],
	    Z(n,-(n-2)/4) * C[1] * Diag(List([0..n-1], x->E(4)^(k*x))))
    )

));

RulesFor(DFT4, rec(
    DFT4_Base := rec(
	isApplicable := P -> P[1]=2,
	forTransposition := false,
	rule := (P,C) -> E(8)*Mat((1/E(8))*MatSPL(ApplyFunc(DFT4, P)))
    ),

    DFT4_CT := Inherit(DFT_CT, rec(
	allChildren  := P -> Map2(DivisorPairs(P[1]),
	    (m,n) -> [ DFT2(m, P[2]), DFT3(n, P[2]) ]),

	rule := (P,C) -> let(mn := P[1], m := Rows(C[1]), n := Rows(C[2]),
	    Tensor(C[1], I(n)) * Diag(fPrecompute(Tw4(mn, n, P[2]))) *
	    Tensor(I(m), C[2]) * L(mn, m))
    )),
));

NewRulesFor(DFT4, rec(
    DFT4_PF := rec(
    	info         := "DFT4_n*k -> diag * perm * (DFT4_n_a tensor DFT4_k_b) * perm",
    	forTransposition := true,
    	applicable     := nt -> nt.params[1] > 2 and DivisorPairsRP(nt.params[1]) <> [] and not nt.hasTags(),
        children := nt -> let(N := nt.params[1], k:= nt.params[2], Map2(DivisorPairsRP(N),
    		(r,s) -> [ DFT4(r, 1/s mod r), DFT4(s, 1/r mod s) ])),

        apply := (nt, C, cnt) -> let(
    	    N := nt.params[1],
    	    r := Rows(C[1]),
    	    s := Rows(C[2]),
    	    alpha := 1 / s mod r,
    	    beta  := 1 / r mod s,
    	    i := Ind(r*s),
    	    d := Diag(Lambda(i, -sign(imod(s*alpha*idiv(i,s) + r*beta*imod(i,s), 2*N) - N))),
    	    CRT(r,s).transpose() *
    	    d * (-E(4))*
    	    Tensor(C[1], C[2]) *
    	    d *
    	    CRT(r,s)
    	)
#D    	isApplicable     := P -> P[1] > 2 and DivisorPairsRP(P[1]) <> [] and Length(P[3]) = 0,
#D    	allChildren := P -> let(N := P[1], k:= P[2], Map2(DivisorPairsRP(N),
#D    		(r,s) -> [ DFT4(r, 1/s mod r), DFT4(s, 1/r mod s) ])),
#D
#D    	rule := (P,C) -> let(
#D    	    N := P[1],
#D    	    r := Rows(C[1]),
#D    	    s := Rows(C[2]),
#D    	    alpha := 1 / s mod r,
#D    	    beta  := 1 / r mod s,
#D    	    i := Ind(r*s),
#D    	    d := Diag(Lambda(i, -sign(imod(s*alpha*idiv(i,s) + r*beta*imod(i,s), 2*N) - N))),
#D    	    CRT(r,s).transpose() *
#D    	    d * (-E(4))*
#D    	    Tensor(C[1], C[2]) *
#D    	    d *
#D    	    CRT(r,s)
#D    	)
    )
));

ext: gi


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

LRB_16x32f.setRules(rec(
  binrules := [ rec(
          perm := rec(
              spl := L(32, 16),
              l := 1,
              N := 32,
              n := 16,
              r := 1 ),
          instr := [ rec(
                  instr := vunpacklo_16x32f,
                  p := [  ] ), rec(
                  instr := vunpackhi_16x32f,
                  p := [  ] ) ],
          v := 16,
          vperm := VPerm(L(32, 16), (y, x) -> chain(assign(vref(y, 0, 16), vunpacklo_16x32f(vref(x, 0, 16), vref(x, 16, 16), [  ])), assign(vref(y, 16, 16), vunpackhi_16x32f(vref(x, 0, 16), vref(x, 16, 16), [  ]))), 16, 2) ), rec(
          perm := rec(
              spl := Tensor(
                  L(16, 8),
                  I(2)
                ),
              l := 1,
              N := 16,
              n := 8,
              r := 2 ),
          instr := [ rec(
                  instr := vunpacklo2_16x32f,
                  p := [  ] ), rec(
                  instr := vunpackhi2_16x32f,
                  p := [  ] ) ],
          v := 16,
          vperm := VPerm(Tensor(
                L(16, 8),
                I(2)
              ), (y, x) -> chain(assign(vref(y, 0, 16), vunpacklo2_16x32f(vref(x, 0, 16), vref(x, 16, 16), [  ])), assign(vref(y, 16, 16), vunpackhi2_16x32f(vref(x, 0, 16), vref(x, 16, 16), [  ]))), 16, 2) ) ],
  unrules := [  ],
  x_I_vby2 := [ [ [  ], [  ] ], [ [  ], [  ] ] ] ));


text:
############################################################################
##
##  loading.gi                    IRREDSOL                  Burkhard Höfling
##
##  Copyright © 2003–2016 Burkhard Höfling
##

############################################################################
##
#F  IsAvailableAbsolutelyIrreducibleSolubleGroupData(<n>, <q>)
##
##  see IRREDSOL documentation
##
InstallGlobalFunction(IsAvailableAbsolutelyIrreducibleSolubleGroupData,
    function(n, q)

        if not IsPosInt(n) or not IsPosInt(q) or not IsPPowerInt(q)  then
            Error("n and q must be positive integers and q must be a prime power");
        fi;
        return TryLoadAbsolutelyIrreducibleSolubleGroupData(n, q);
    end);

############################################################################
##
#F  TryLoadAbsolutelyIrreducibleSolubleGroupData(<n>, <q>)
##
InstallGlobalFunction(TryLoadAbsolutelyIrreducibleSolubleGroupData,
    function(n, q)

        local filename, pathname, dirs, dir, inds, desc, guardianData, guardianDesc, data,
            i, j, p, d, e, orders, divs, t, info, q0, m, G, H, hom, inv,
            pos, gdPos, maxes, cpcgscode, pcgs;

        if not IsPosInt(n) or not IsPosInt(q) or not IsPPowerInt(q)  then
            Error("n and q must be positive integers and q must be a prime power");
        fi;

        if not IsBound(IRREDSOL_DATA.GAL_PERM[n]) then
            IRREDSOL_DATA.GAL_PERM[n] := [];
        fi;
        if not IsBound(IRREDSOL_DATA.GUARDIANS[n]) then
            IRREDSOL_DATA.GUARDIANS[n] := [];
        fi;

        if n = 1 then
            if not IsBound(IRREDSOL_DATA.GROUPS_DIM1[q]) then
                Info(InfoIrredsol, 2, "Computing irreducible soluble group data for ",
                  "GL(", n, ", ", q, ")");
                p := SmallestRootInt(q);
                e := LogInt(q, p);
                orders := ShallowCopy(DivisorsInt(q-1));
                # remove groups which are over a proper subfield
                if e > 1 then
                  for d in Set(Factors(e)) do
                     SubtractSet(orders, DivisorsInt(p^(e/d) - 1));
                  od;
                fi;

                if p < q then
                  divs := Reversed(Difference(DivisorsInt(e), [1]));
                  info := List(orders, o -> [o, []]);
                  for j in [1..Length(divs)] do
                     d := divs[j];
                     q0 := p^(e/d);

                     for i in info do

                        for t in DivisorsInt(GcdInt(i[1], d)) do
                            m := d/t;
                            if (q0^m -1) mod (i[1]/t) = 0 then
                              i[2][j] := m;
                            fi;
                        od;
                     od;
                  od;
                else
                  info := List(orders, o -> [o]);
                fi;

                IRREDSOL_DATA.GUARDIANS[1][q] := [CyclicGroup(IsPcGroup, q-1)];
                IRREDSOL_DATA.GROUPS_DIM1[q] := info;
                IRREDSOL_DATA.GAL_PERM[1][q] := ();
            fi;
            return true;
        fi;

       if not IsBound(IRREDSOL_DATA.GROUPS[n]) then
            IRREDSOL_DATA.GROUPS[n] := [];
        fi;
        if not IsBound(IRREDSOL_DATA.MAX[n]) then
            IRREDSOL_DATA.MAX[n] := [];
        fi;
        if not IsBound(IRREDSOL_DATA.GROUPS_LOADED[n]) then
            IRREDSOL_DATA.GROUPS_LOADED[n] := [];
        fi;

        if not IsBound(IRREDSOL_DATA.GUARDIANS[n][q]) or not IsBound(IRREDSOL_DATA.GROUPS[n][q])
                or not IsBound(IRREDSOL_DATA.MAX[n][q]) or not IsBound(IRREDSOL_DATA.GAL_PERM[n][q])
                or not IsBound(IRREDSOL_DATA.GROUPS_LOADED[n][q]) then

            Unbind(IRREDSOL_DATA.GROUPS_LOADED[n][q]); # if anything fails during loading
                # we won't have inconsistencies

            pathname := Concatenation("data/gl_", String(n), "_",String(q),".grp");

            Info(InfoIrredsol, 2, "Reading data file ", pathname);
            if not ReadPackage("irredsol", pathname) then
                return false;
            fi;

            if not IsBound(IRREDSOL_DATA.GUARDIANS[n][q]) or not IsBound(IRREDSOL_DATA.GROUPS[n][q])
                    or not IsBound(IRREDSOL_DATA.GAL_PERM[n][q]) then
                Error("Panic: reading data file didn't define required data");
            fi;

            # convert guardian data

            guardianData := IRREDSOL_DATA.GUARDIANS[n][q];
            maxes := [];

            for gdPos in [1..Length(guardianData)]  do
                guardianDesc := guardianData[gdPos];
                guardianDesc[1] := List(guardianDesc[1], m -> FFMatrixByNumber(m, n, q));

                G := PcGroupCode(guardianDesc[3], guardianDesc[2]);
                H := Group(guardianDesc[1]);
                SetSize(H, Size(G));
                hom := GroupGeneralMappingByImagesNC(G, H, FamilyPcgs(G), guardianDesc[1]);
                SetIsGroupHomomorphism(hom, true);
                SetIsBijective(hom, true);
                guardianDesc[3] := hom;
                if TestFlag(guardianDesc[5], 0) then
                    cpcgscode := 2^Length(FamilyPcgs(G)) - 1;
                    pos := PositionProperty(IRREDSOL_DATA.GROUPS[n][q],
                        desc -> desc[2] = cpcgscode and desc[1] = gdPos);
                    if pos = fail then
                        Error("panic: did not find guardian in list of groups");
                    fi;
                    AddSet(maxes, pos);
                fi;
            od;

            # compupte group orders

            data := IRREDSOL_DATA.GROUPS[n][q];
            for d in data do
                pcgs := FamilyPcgs(Source(guardianData[d[1]][3])); # get pcgs of guardian
                d[4] := OrderGroupByCanonicalPcgsByNumber(pcgs, d[2]);
            od;

            IRREDSOL_DATA.GUARDIANS[n][q] := guardianData; # now attach converted data
            IRREDSOL_DATA.MAX[n][q] := maxes;
            MakeImmutable(IRREDSOL_DATA.GUARDIANS[n][q]);
            MakeImmutable(IRREDSOL_DATA.GROUPS[n][q]);
            MakeImmutable(IRREDSOL_DATA.MAX[n][q]);
            IRREDSOL_DATA.GROUPS_LOADED[n][q] := true;
            Info(InfoIrredsol, 2, "irreducible soluble group data for GL(", n, ", ", q, ") loaded");
        fi;
        return true;
    end);

###########################################################################
##
#F  LoadAbsolutelyIrreducibleSolubleGroupData(<n>, <q>)
##
##  see IRREDSOL documentation
##
InstallGlobalFunction(LoadAbsolutelyIrreducibleSolubleGroupData,
    function(n, q)

        if not IsPosInt(n) or not IsPosInt(q) or not IsPPowerInt(q)  then
            Error("n and q must be positive integers and q must be a prime power");
        fi;

        if not TryLoadAbsolutelyIrreducibleSolubleGroupData(n, q) then
            Error("Panic: missing data file for GL(",n,", ", q, ")");
        fi;
    end);

############################################################################
##
#F  LoadedAbsolutelyIrreducibleSolubleGroupData()
##
##  see IRREDSOL documentation
##
InstallGlobalFunction(LoadedAbsolutelyIrreducibleSolubleGroupData,
    function()

        local n, p, data, fields;

        data := [];
        for n in [1..Length(IRREDSOL_DATA.GROUPS)] do
            if IsBound(IRREDSOL_DATA.GROUPS[n]) then
                fields := [];
                for p in [1..Length(IRREDSOL_DATA.GROUPS[n])] do
                    if IsBound(IRREDSOL_DATA.GROUPS[n][p]) then
                        Add(fields, p);
                    fi;
                od;
                if not IsEmpty(fields) then
                    Add(data, [n, fields]);
                fi;
            fi;
        od;
        return data;
    end);

############################################################################
##
#F  UnloadAbsolutelyIrreducibleSolubleGroupData([<n>[, <q>]])
##
##  see IRREDSOL documentation
##
InstallGlobalFunction(UnloadAbsolutelyIrreducibleSolubleGroupData,
    function(arg)

        local UnbindIfBound;

        UnbindIfBound := function(arg)
            local data, i;

            data := arg[1];
            i := 2;
            while i < Length(arg)  do
                if not IsBound(data[arg[i]]) then
                    return false;
                fi;
                data := data[arg[i]];
                i := i + 1;
            od;
            if IsBound(data[arg[i]]) then
                Unbind(data[arg[i]]);
                return true;
            else
                return false;
            fi;
        end;

        if Length(arg) = 0 then
            IRREDSOL_DATA.GUARDIANS := [];
            IRREDSOL_DATA.GROUPS := [];
            IRREDSOL_DATA.GROUPS_LOADED := [];
            IRREDSOL_DATA.GAL_PERM := [];
            IRREDSOL_DATA.MAX := [];
            IRREDSOL_DATA.GROUPS_DIM1 := [];
            IRREDSOL_DATA.PRIM_GUARDIANS := [];
        elif IsPosInt(arg[1]) and Length(arg) = 1 then
            UnbindIfBound(IRREDSOL_DATA.GAL_PERM, arg[1]);
            UnbindIfBound(IRREDSOL_DATA.GUARDIANS, arg[1]);
            UnbindIfBound(IRREDSOL_DATA.PRIM_GUARDIANS, arg[1]);
            if arg[1] = 1 then
                UnbindGlobal("IRREDSOL_DATA.GROUPS_DIM1");
                BindGlobal("IRREDSOL_DATA.GROUPS_DIM1", []);
            else
                UnbindIfBound(IRREDSOL_DATA.GROUPS, arg[1]);
                UnbindIfBound(IRREDSOL_DATA.MAX, arg[1]);
            fi;
        elif Length(arg) = 2 and IsPosInt(arg[1]) and IsPPowerInt(arg[2]) then
            UnbindIfBound(IRREDSOL_DATA.GUARDIANS, arg[1], arg[2]);
            UnbindIfBound(IRREDSOL_DATA.PRIM_GUARDIANS, arg[1], arg[2]);
            UnbindIfBound(IRREDSOL_DATA.GAL_PERM, arg[1], arg[2]);
            if arg[1] = 1 then
                UnbindIfBound(IRREDSOL_DATA.GROUPS_DIM1, arg[2]);
            else
                UnbindIfBound(IRREDSOL_DATA.GROUPS, arg[1], arg[2]);
                UnbindIfBound(IRREDSOL_DATA.MAX, arg[1], arg[2]);
            fi;
        else
            Error("Usage: `UnloadAbsolutelyIrreducibleSolubleGroupData( [n [, q]] )'");
        fi;
    end);

############################################################################
##
#F  IsAvailableIrreducibleSolubleGroupData(<n>, <q>)
##
##  see IRREDSOL documentation
##
InstallGlobalFunction(IsAvailableIrreducibleSolubleGroupData,
    function(n, q)

        local d;

        if not IsPosInt(n) or not IsPosInt(q) or not IsPPowerInt(q)  then
            Error("n and q must be positive integers and q must be a prime power");
        fi;
        for d in DivisorsInt(n) do
            if not IsAvailableAbsolutelyIrreducibleSolubleGroupData(n/d, q^d) then
                return false;
            fi;
        od;
        return true;
    end);

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

max_stars_repo_name: fingolfin/irredsol

max_issues_repo_name: fingolfin/irredsol

max_forks_repo_name: fingolfin/irredsol


text:
Bottles := function(n)
	local line, i, j, u;
	line := function(n)
		s := String(n);
		if n < 2 then
			return Concatenation(String(n), " bottle of beer");
		else
			return Concatenation(String(n), " bottles of beer");
		fi;
	end;
	for i in [1 .. n] do
		j := n - i + 1;
		u := line(j);
		Display(Concatenation(u, " on the wall"));
		Display(u);
		Display("Take one down, pass it around");
		Display(Concatenation(line(j - 1), " on the wall"));
		if i <> n then
			Display("");
		fi;
	od;
end;

ext: gap

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text:
#######################################################################
# Description: A function to initialize a generator for random        #
#              permutations, described in section 4.1.                #
# Input: - list of generators S                                       #
# Output: - list of permutations X including repetitions of generators#
#           from S                                                    #
#######################################################################

InitGenerator := function(S)
    local i, X, k;
    X := List(S); k := Length(X);
    for i in [k+1..Maximum(11, Length(S))] do
        Add(X, X[i-k]);
    od;
    if IsReadOnlyGlobal("X") = true then
        MakeReadWriteGlobal("X");
    fi;
    return X;
end;

#######################################################################
#######################################################################

#######################################################################
# Description: A function to calculate a random permutation,          #
#              described in section 4.1.                              #
# Input: - list of generators S, generator X returned by InitGenerator#
# Output: - random permutation a as cycle                             #
#######################################################################

RandomPerm := function(S, X)
    local Randomize, a, r, i;

    r := Maximum(11, Length(S));
    Randomize := function()
        local s, t, e, range;

        # indices s and t specify which permutations X[s], X[t]
        # will be multiplied with each other
        # e specifies whether X[t] or its inverse will be used
        range := [1..r];  s := Random(range);  e := Random([-1,1]);
        t := Random(range);  Remove(range, s);

        # alternate the order of multiplication
        # of X[t] and X[s] randomly
        # update generator X

        if Random([1,2]) = 1 then
           X[s] := X[s]*X[t]^e;
           a := a*X[s];
        else
           X[s] := X[t]^e*X[s];
           a := X[s]*a;
        fi;
        return a;
    end;

    # calculate random permutation a as
    # product of random permutations from X
    a := ();
    for i in [1..50] do
        Randomize();
    od;
    return a;
end;

max_stars_repo_name: AlexanderKlemps/parallel-schreier-sims

max_issues_repo_name: AlexanderKlemps/parallel-schreier-sims

max_forks_repo_name: AlexanderKlemps/parallel-schreier-sims


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

# NOTE: Support odd sizes ie with (I_4 dirsum J_3)
condIJ := (N, m) -> let(k := N/m, Checked(IsSymbolic(k) or IsEvenInt(k),
    fCond(diagTensor(II(k/2), II(2,0,1), II(m)), fTensor(fId(k), fId(m)),
                                                 fTensor(fId(k), J(m)))
));

# NOTE: Support odd sizes ie with (I_4 dirsum OS_3)
condIOS := (N, m) -> let(k := N/m, Checked(IsSymbolic(k) or IsEvenInt(k),
    fCond(diagTensor(II(k/2), II(2,0,1), II(m)), fTensor(fId(k), fId(m)),
                                                 fTensor(fId(k), OS(m, -1)))
));

condM := (N, m) -> L(N, m) * Prm(condIJ(N, m));
condK := (N, m) -> Prm(condIJ(N, N/m)) * L(N, m);

condMp := (N, m) -> L(N, m) * Prm(condIOS(N, m));
condKp := (N, m) -> Prm(condIOS(N, N/m)) * L(N, m);

condB := (N, m) -> Checked(IsEvenInt(m), let(k := N/m,
    condK(N, m) * Tensor(I(N/m), condM(m,m/2))));

condBp := (N, m) -> Checked(IsEvenInt(m), let(k := N/m,
    condKp(N, m) * DirectSum(L(m,m/2), Tensor(I(k-1), condM(m,m/2)))));

condBpu := (N, m) -> Checked(IsEvenInt(m), let(k := N/m,
    condKp(N, m) * DirectSum(I(m), Tensor(I(k-1), condM(m,m/2)))));

# condBpu(k*m,m) == Scat(Refl0_u(k,m)) = Refl0_u(k,m).transpose()
# u = stands for URDFT. i.e. this is for URDFT based algorithms
# NOTE: not valid for odd <k> (but will need new ASPF rule anyway with an extra child)
# GOOD:  valid for odd <m>.
Class(Refl0_u, PermClass, rec(
    def := (k, m) -> rec(),
    domain := self >> self.params[1] * self.params[2],
    range := self >> self.params[1] * self.params[2],
    lambda := self >> let(k := self.params[1], m := self.params[2], i := Ind(k*m),
        Lambda(i, cond(leq(i, m-1), k*i,
                                    cond(leq(imod(i, m), idiv(m-1,2)), 2*k*imod(i,m) + idiv(i,m),
                                                                      2*k*m - 2*k*imod(i,m) - idiv(i,m)))))
));

# Valid for odd <k>, even <m>
# <omit> is how many initial points we want to omit from the permutation
# the only sensible values for it are 0 and 1
Class(Refl0_odd, PermClass, rec(
    def := (k, m, omit) -> rec(), #size := k*m+idiv(m+1,2)-omit),
    range := self >> self.domain(),
    domain := self >> let(k:=self.params[1], m:=self.params[2], omit:=self.params[3],
        k*m + idiv(m+1, 2) - omit),

    lambda := self >> let(
        k := self.params[1],   m := self.params[2], mc := idiv(m+1, 2),
        N := self.domain(), omit := self.params[3], mf := idiv(m, 2),
        N2:= 2*k*m+m,         ii := Ind(N-omit),    i  := ii+omit,

        Lambda(ii, cond(leq(i, mc-1), (2*k+1)*i,
                                     cond(leq(imod(i-mc, m), mc-1),       (2*k+1)*imod(i-mc,m) + idiv(i+mf,m),
                                                                    N2  - (2*k+1)*imod(i-mc,m) - idiv(i+mf,m)))
            -omit))
));

# condB(k*m,m) == Scat(Refl1(k,m)) = Refl1(k,m).transpose()
# NOTE: not valid for odd <k> (but will need new ASPF rule anyway with an extra child)
# GOOD:  valid for odd <m>.
Class(Refl1, PermClass, rec(
    def := (k, m) -> rec(), #size := k*m),
    domain := self >> self.params[1]*self.params[2],
    range := self >> self.params[1]*self.params[2],
    lambda := self >> let(k := self.params[1], m := self.params[2], i := Ind(k*m),
        Lambda(i, cond(leq(imod(i, m), idiv(m-1,2)), 2*k*imod(i,m) + idiv(i,m),
                                       2*k*m - 1 - 2*k*imod(i,m) - idiv(i,m))))
));

Declare(KK, MM, KKp, MMp);

Class(KK, PermClass, rec(
    def := (k, m) -> rec(), #size := k*m),
    domain := self >> self.params[1]*self.params[2],
    range := self >> self.params[1]*self.params[2],
    transpose := self >> MM(self.params[2], self.params[1]),
    lambda := self >> let(k := self.params[1], m := self.params[2],
        fCompose(Tr(k,m), condIJ(k*m, k)).lambda()),
));

Class(MM, PermClass, rec(
    def := (k, m) -> rec(), #size := k*m),
    domain := self >> self.params[1]*self.params[2],
    range := self >> self.params[1]*self.params[2],
    transpose := self >> KK(self.params[2], self.params[1]),
    lambda := self >> let(k := self.params[1], m := self.params[2],
        fCompose(condIJ(k*m, m), Tr(k, m)).lambda()),
));

Class(KKp, PermClass, rec(
    def := (k, m) -> rec(), #size := k*m),
    domain := self >> self.params[1]*self.params[2],
    range := self >> self.params[1]*self.params[2],
    transpose := self >> MMp(self.params[2], self.params[1]),
    lambda := self >> let(k := self.params[1], m := self.params[2],
        fCompose(Tr(k,m), condIOS(k*m, k)).lambda()),
));

Class(MMp, PermClass, rec(
    def := (k, m) -> rec(), #size := k*m),
    domain := self >> self.params[1]*self.params[2],
    range := self >> self.params[1]*self.params[2],
    transpose := self >> KKp(self.params[2], self.params[1]),
    lambda := self >> let(k := self.params[1], m := self.params[2],
        fCompose(condIOS(k*m, m), Tr(k, m)).lambda()),
));

Declare(IP);

libB := (k, m) -> Checked(IsSymbolic(m) or IsEvenInt(m),
    KK(k, m) * Tensor(I(k), MM(2,m/2)));

libBp := (k, m) -> Checked(IsSymbolic(m) or IsEvenInt(m),
    KKp(k, m) * DirectSum(Tr(2,m/2), Tensor(I(k-1), MM(2,m/2))));

libBpu := (k, m) -> Checked(IsSymbolic(m) or IsEvenInt(m),
    KKp(k, m) * DirectSum(I(m), Tensor(I(k-1), MM(2,m/2))));

#F IJ(<size>, <n>) - I_n dirsum J_n dirsum I_n dirsum ... (repeats N/n times, n|N)
#F
#
# Subst(i -> BinXor(i, BinAnd(i,$n) - BinShr(BinAnd(i,$n), $(LogInt(n,2))))),
# Class(IJ, PermClass, rec(
#   def := (N, n) -> Checked(IsInt(N), IsInt(n), N mod n = 0,
#       rec(
#       size := N,
#       direct  := Subst(x -> x + (QuoInt(x,$n) mod 2) * ($(n-1) - 2*(x mod $n))),
#       inverse := Subst(x -> x + (QuoInt(x,$n) mod 2) * ($(n-1) - 2*(x mod $n)))
#       )),
#   transpose := self >> self
#);
IJ := (N, n) -> IP(N, J(n));

#F IP(<N>, <splperm>) - I_n dirsum P_n dirsum I_n dirsum ... (repeats N/n times, n|N)
#F
#F    P_n given by <splperm>, which is a permutation function
#F
Class(IP, PermClass, rec(
    def := (N, splperm) -> Checked(IsInt(N), IsSPL(splperm), N mod Rows(splperm) = 0, rec()),
    range := self >> self.params[1],
    domain := self >> self.params[1],

    lambda := self >> let(
	N := self.params[1],
	P := self.params[2],
	Pl := P.lambda(),
	n := P.domain(),
	i := Ind(N),
	Lambda(i, n*idiv(i,n) + imod(1+idiv(i,n), 2)*imod(i, n) +
                                imod(  idiv(i,n), 2)*Pl.at(imod(i,n)))
    ),

    transpose := self >> IP(self.params[1], self.params[2].transpose())
));

# -----------------------------------------------------------------------------#
#F LIJ(<N>) - L(N, N/2) * (I(N/2) dirsum J(N/2))
##   works also for odd sizes
Class(LIJ, PermClass, rec(
    def := N -> Checked(IsInt(N), rec()),
    range := self >> self.params[1],
    domain := self >> self.params[1],
    lambda := self >> let(
	N := self.params[1],
	N1 := When(IsOddInt(N), N, N-1),
	i := Ind(N),
	When(not self.transposed,
            Lambda(i, QuoInt(N,2) + idiv(N-i, 2)*cond(neq(imod(N1-i,2),0), -1, 1)),
            Lambda(i, 2*i + idiv(i, QuoInt(N+1,2)) * (2*N - 1 - 4*i))
        ))
));

# -- Originals
# K ( <n>, <m> )
#   returns the equivalent of the stride permutation for Cooley-Tukey FFT
#   for the DCT3 and 4. <m> has to divide <n>
#     (I_n/m dirsum J_n/m dirsum I_n/m ....) * L^n_m  (direct sum alternately)
#   This permutation is to be used if it occurs on the left side in the rule.
#   It is K(n, m)^-1 =  M(n, n/m).
K := (mn, m) -> Checked(IsPosInt(mn), IsPosInt(m), mn mod m = 0,
    IJ(mn, mn/m) * L(mn, m));

M := (mn, m) -> Checked(IsPosInt(mn), IsPosInt(m), mn mod m = 0,
    L(mn, m) * IJ(mn, m));
# -------------

########################
# Generating functions
########################

L_or_OS := (n, str) -> When((n mod str) <> 0, OS(n, str), L(n, str));

Kp := (NN, m) -> let(n := NN / m,
    Prm(IP(NN, OS(n, -1)))* L(NN, m));

Mp := (NN, m) -> Kp(NN, NN/m).transpose();

# Type 0 symmetries  (ce0, co0) aka whole-point
# BH0(<Nreal>, <n>, <base>, <stride>), Nreal = number of output real pts = Rows(PRDFT(.))
# BH0 := (Nreal,n,base,stride) -> BH(Int(Nreal/2)+1, Nreal, n, base, stride);

# Type 1 symmetries  (ce1, co1) aka half-point
# BH1(<Nreal>, <n>, <base>, <stride>), Nreal = number of output real pts = Rows(PRDFT3(.))
# BH1 := (Nreal,n,base,stride) -> BH(Int((Nreal+1)/2), Nreal-1, n, base, stride);

# BH(<N>, <reflect>, <n>, <b>, <s>) - Reflective stride function
#   Defined as  i -> When(i < ceil(n/2)-1, si + b, reflect - si - b)
#
Class(BH, FuncClass, rec(
    def := (N,reflect,n,base,stride) -> rec(), #N := N, n := n),
    domain := self >> self.params[3],
    range  := self >> self.params[1],
    lambda := self >> let(
	N:=self.params[1], reflect:=self.params[2], n:=self.params[3], b:=self.params[4],
	s:=self.params[5], i:=Ind(n),
	Lambda(i, cond(leq(i, idiv(n-1, 2)), s*i + b, reflect - s*i - b)))
));

Declare(Refl);
# Refl(<N>, <reflect>, <n>, <func>) - Reflective arbitrary function
#   Defined as  i -> When(func(i) < N, func(i), reflect - func(i))
#   with domain of <n> points (usually smaller than domain(func))
#
Class(Refl, FuncClass, rec(
    def := (N,reflect,n,func) -> rec(), #N := N, n := n),
    domain := self >> self.params[3],
    range  := self >> self.params[1],
    transpose := self >> Refl(self.params[1], self.params[2], self.params[3], self.params[4].transpose()),
    lambda := self >> let(
	func:=self.params[4].lambda(), reflect:=self.params[2], N:=self.params[1],
	i := Ind(self.params[3]),
	Lambda(i, cond(leq(func.at(i), N-1), func.at(i), reflect-func.at(i))))
));


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

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

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[19,17,37]);
    uu:=Unipotent(chevalleyAdj(cent),[[3,1]]);
    tmp2:=Conj(tmp2,uu);
    #tmp2:=ConjugateByTori(tmp2,[3],[-1]);

    return [tmp1,tmp2];
end;

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

    rep:=Representative(o);
    #pr[53]=[1,1,1,2,1,1,1,0]
    #pr[52]=[1,1,1,2,2,1,0,0]
    #pr[93]=[1,2,3,4,3,2,1,0]
    tmp1:=ApplyRootsReflections(rep,[53,52,93]);
    uu:=Unipotent(chevalleyAdj(rep),[[1,-1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);
    tmp1:=ConjugateByTori(tmp1,[1],[-1]);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[53,52,93]);
    uu:=Unipotent(chevalleyAdj(cent),[[1,-1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);
    tmp2:=ConjugateByTori(tmp2,[1],[-1]);

    return [tmp1,tmp2];
end;

#
# cls 64: A_2^2
#
handle64char2:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[64];
    info:=infos[64];

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

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[43,41,79]);
    uu:=Unipotent(chevalleyAdj(cent),[[1,1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);
    tmp2:=ConjugateByTori(tmp2,[3],[-1]);

    return [tmp1,tmp2];
end;

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

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

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

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

    tmp4:=0;
    tmp4:=ApplyRootsReflections(cent,[26,27]);
    uu:=Unipotent(chevalleyAdj(cent),[[3,1],[4,1],[5,-1],[11,1],[12,-1]],Ordering(cent));
    tmp4:=Conj(tmp4,uu);
    tmp4:=ConjugateByTori(tmp4,[2,4],[-1,-1,-1,-1]);

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

#
# NOT FINISHED .. Should be an S_3
#
# cls 57: D_4(a_1)A_1 ([2,5,9,10,7] in E_7) [2,5,10,11,7] in E_8  ?CHANGE THE REPRESENTATIVE IN CHARACTERISTIC 2?
#
handle57char2:=function(orbs,infos)
    local tmp1,tmp2,tmp3,tmp4,o,info,rep,cent,uu;
    o:=Classes(orbs)[57];
    info:=infos[57];

    rep:=Representative(o);

    #pr[23]=[1,1,1,1,0,0,0,0]
    #pr[24]=[1,0,1,1,1,0,0,0]
    tmp1:=ApplyRootsReflections(rep,[23,24]);
    tmp1:=ConjNegUa(tmp1,4,1);
    uu:=Unipotent(chevalleyAdj(rep),[[2,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);

    #pr[65]=[1,1,2,2,1,1,1,1]
    #pr[67]=[1,1,1,2,2,1,1,1]
    tmp3:=ApplyRootsReflections(rep,[65,67,4]);
    uu:=Unipotent(chevalleyAdj(rep),[[12,1],[11,1]],Ordering(rep));
    tmp3:=Conj(tmp3,uu);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[23,24]);
    tmp2:=ConjNegUa(tmp2,4,1);
    uu:=Unipotent(chevalleyAdj(cent),[[2,1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);

    tmp4:=0;
    tmp4:=ApplyRootsReflections(cent,[65,67,4]);
    uu:=Unipotent(chevalleyAdj(cent),[[12,1],[11,1]],Ordering(cent));
    tmp4:=Conj(tmp4,uu);

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

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

    rep:=Representative(o);

    #pr[33]=[0,1,1,1,1,1,0,0]
    #pr[31]=[1,0,1,1,1,1,0,0]
    #pr[64]=[1,1,2,2,2,1,1,0]
    #pr[89]=[1,2,2,4,3,2,1,0]
    tmp1:=ApplyRootsReflections(rep,[33,31,64,89]);
    uu:=Unipotent(chevalleyAdj(rep),[[2,1],[3,1],[9,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[33,31,64,89]);
    uu:=Unipotent(chevalleyAdj(cent),[[2,1],[3,1],[9,1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);

    return [tmp1,tmp2];
end;

#
# cls 51: D_4(a_1)A_2 [2,5,10,11,7,8] in E_8
#
handle51char2:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[51];
    info:=infos[51];

    rep:=Representative(o);

    #pr[23]=[1,1,1,1,0,0,0,0]  asta merge si cu rep=[2,3,10,12,7,8] urmat de ConjNegUa(tmp1,4,1)
    #pr[24]=[1,0,1,1,1,0,0,0]
    tmp1:=ApplyRootsReflections(rep,[23,24]);
    tmp1:=ConjNegUa(tmp1,4,1);
    uu:=Unipotent(chevalleyAdj(rep),[[5,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[23,24]);
    tmp2:=ConjNegUa(tmp2,4,1);
    uu:=Unipotent(chevalleyAdj(cent),[[5,1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);

    return [tmp1,tmp2];
end;

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

    rep:=Representative(o);

    #pr[41]=[0,1,1,1,1,1,1,0]
    #pr[39]=[1,0,1,1,1,1,1,0]
    #pr[57]=[1,1,2,2,2,1,0,0]
    #pr[89]=[1,2,2,4,3,2,1,0]
    tmp1:=ApplyRootsReflections(rep,[41,39,57,89]);
    uu:=Unipotent(chevalleyAdj(rep),[[2,1],[3,1],[9,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[41,39,57,89]);
    uu:=Unipotent(chevalleyAdj(cent),[[2,1],[3,1],[9,1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);

    return [tmp1,tmp2];
end;

#
# cls 48: D_5(a_1)
#
handle48char2:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[48];
    info:=infos[48];

    rep:=Representative(o);

    #pr[82]=[1,1,2,3,3,2,1,0]
    #pr[80]=[1,2,2,3,2,2,1,0]
    tmp1:=ApplyRootsReflections(rep,[82,80]);
    uu:=Unipotent(chevalleyAdj(rep),[[4,1],[12,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);
    tmp1:=ConjugateByTori(tmp1,[2,4],[-1,-1,-1,-1]);

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

    return [tmp1,tmp2];
end;

#
# cls 47: A_4A_1^2
#
handle47char2:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[47];
    info:=infos[47];

    rep:=Representative(o);

    #pr[71]=[1,1,2,2,2,2,1,0]
    #pr[72]=[1,1,2,2,2,1,1,1]
    #pr[112]=[1,3,3,5,4,3,2,1]
    #pr[111]=[2,2,3,5,4,3,2,1]
    tmp1:=ApplyRootsReflections(rep,[71,72,112,111]);
    uu:=Unipotent(chevalleyAdj(rep),[[2,1],[3,1],[9,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[71,72,112,111]);
    uu:=Unipotent(chevalleyAdj(cent),[[2,1],[3,1],[9,1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);

    return [tmp1,tmp2];
end;

# THE COMPONENT GROUP IS GENERATED BY u IN CHARACTERISTIC 2
#
# cls 40: D_4A_2
#
#handle40char3:=function(orbs,infos)
#    local tmp1,tmp2,o,info,rep,cent,uu;
#    o:=Classes(orbs)[40];
#    #info:=infos[40];
#
#    rep:=Representative(o);
#
#    #pr[34]=[0,1,0,1,1,1,1,0]
#    #pr[35]=[0,0,1,1,1,1,1,0]
#    #pr[68]=[0,1,1,2,2,2,1,1]
#    tmp1:=ApplyRootsReflections(rep,[34,35,68]);
#    uu:=Unipotent(chevalleyAdj(rep),[[8,-1]],Ordering(rep));
#    tmp1:=Conj(tmp1,uu);
#    tmp1:=ConjugateByTori(tmp1,[4,6,7,8],[-1,-1,-1,-1]);
#
#    tmp2:=0;
#    #cent:=FromPositiveBorel(o,info[2]);
#    #tmp2:=ApplyRootsReflections(cent,[36,34,46,61]);
#    #uu:=Unipotent(chevalleyAdj(cent),[[2,-1],[4,1]],Ordering(cent));
#    #tmp2:=Conj(tmp2,uu);
#    #tmp2:=ConjugateByTori(tmp2,[2,3,4],[-1,-1,-1,-1]);
#
#    return [tmp1,tmp2];
#end;

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

    rep:=Representative(o);
    tmp1:=ConjNegUa(rep,3,1);
    tmp1:=ConjNegUa(tmp1,2,1);
    uu:=Unipotent(chevalleyAdj(rep),[[5,1],[6,1],[9,1],[10,1],[17,1],[18,1],[25,1],[26,1],[33,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ConjNegUa(cent,3,1);
    tmp2:=ConjNegUa(tmp2,2,1);
    uu:=Unipotent(chevalleyAdj(cent),[[5,1],[6,1],[9,1],[10,1],[17,1],[18,1],[25,1],[26,1],[33,1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);

    return [tmp1,tmp2];
end;

#
# cls 35: D_6(a_2)
#
handle35char2:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[35];
    info:=infos[35];

    rep:=Representative(o);

    #pr[42]=[0,1,0,1,1,1,1,0]
    #pr[43]=[0,0,1,1,1,1,1,0]
    tmp1:=ApplyRootsReflections(rep,[42,43]);
    uu:=Unipotent(chevalleyAdj(rep),[[4,1],[6,1],[10,1],[13,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[42,43]);
    uu:=Unipotent(chevalleyAdj(cent),[[4,1],[6,1],[10,1],[13,1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);

    return [tmp1,tmp2];
end;

#
# cls 34: E_6(a_3)A_1
#
handle34char2:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[34];
    info:=infos[34];

    rep:=Representative(o);
    tmp1:=ConjNegUa(rep,3,1);
    tmp1:=ConjNegUa(tmp1,2,1);
    uu:=Unipotent(chevalleyAdj(rep),[[5,1],[6,1],[9,1],[10,1],[17,1],[18,1],[25,1],[26,1],[33,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ConjNegUa(cent,3,1);
    tmp2:=ConjNegUa(tmp2,2,1);
    uu:=Unipotent(chevalleyAdj(cent),[[5,1],[6,1],[9,1],[10,1],[17,1],[18,1],[25,1],[26,1],[33,1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);

    return [tmp1,tmp2];
end;

#
# cls 33: E_7(a_5)
#
handle33char2:=function(orbs,infos)
    local tmp1,tmp2,tmp3,tmp4,o,info,rep,cent,uu;
    o:=Classes(orbs)[33];
    info:=infos[33];

    rep:=Representative(o);
    tmp1:=ApplyRootsReflections(rep,[2,3,5]);
    uu:=Unipotent(chevalleyAdj(rep),[[4,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);

    tmp3:=ConjugateByTori(rep,[2,3,5],[aaa,aaa,aaa]);

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

    tmp4:=0;
    tmp4:=ConjugateByTori(cent,[2,3,5],[aaa,aaa,aaa]);

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

#
# cls 31: E_8(a_7)       NOT FINISHED IN CHARACTERISTIC 2 !!!!
#
handle31char2:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent;
    o:=Classes(orbs)[31];
    #info:=infos[31];

    rep:=Representative(o);

    tmp2:=0;
    #cent:=FromPositiveBorel(o,info[2]);

    return [tmp1,tmp2];
end;

#
# cls 30: A_6       NOT FINISHED IN CHARACTERISTIC 2 !!!!
#
handle30char2:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[30];
    info:=infos[30];

    rep:=Representative(o);
    tmp1:=rep^-1;              # longest words in E7 and E6: n_0*w_0=(17)(36)(45)
    #uu:=Unipotent(chevalleyAdj(rep),[[3,1],[5,1],[7,1],[10,1],[11,1],[24,1],[25,1]],Ordering(rep)); in E_7
    uu:=Unipotent(chevalleyAdj(rep),[[3,1],[5,1],[7,1],[11,1],[12,1],[27,1],[28,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);# =[27,28,24,39]
    tmp2:=coefficients(cent);# =[27,28,24,39]
    tmp2:=Unipotent(chevalleyAdj(cent),[tmp2[1],tmp2[3],tmp2[2],tmp2[4]],Ordering(cent));# =[27,24,28,39]
    uu:=Unipotent(chevalleyAdj(cent),[[3,1],[5,1],[7,1],[11,1],[12,1],[27,1],[28,1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);

    return [tmp1,tmp2];
end;

#
# cls 29: D_6(a_1)
#
handle29char2:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[29];
    info:=infos[29];

    rep:=Representative(o);

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

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

    return [tmp1,tmp2];
end;

# NO COMPONENT GROUP IN CHARACTERISTIC 2 FOR THIS CLASS
#
# cls 27: E_7(a_4)
#
#handle27char2:=function(orbs,infos)
#    local tmp1,tmp2,o,info,rep,cent;
#    o:=Classes(orbs)[27];
#    #info:=infos[27];
#
#    rep:=Representative(o);
#    tmp1:=ConjugateByTorus(rep,4,-1);# h_1 --> h_4
#
#    tmp2:=0;
#    #cent:=FromPositiveBorel(o,info[2]);
#    #tmp2:=ConjugateByTorus(cent,4,-1);# h_1 --> h_4
#
#    return [tmp1,tmp2];
#end;

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

    #pr[61]=[0,1,1,2,2,2,1]
    #pr[59]=[1,1,1,2,2,1,1]
    #pr[58]=[1,1,2,2,1,1,1]

    rep:=Representative(o);
    tmp1:=ApplyRootsReflections(rep,[58,59,61]);
    tmp1:=ConjNegUa(tmp1,4,1);
    uu:=Unipotent(chevalleyAdj(rep),[[2,1],[6,1],[10,1],[11,1],[16,1],[17,1],[23,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[58,59,61]);
    tmp2:=ConjNegUa(tmp2,4,1);
    uu:=Unipotent(chevalleyAdj(cent),[[2,1],[6,1],[10,1],[11,1],[16,1],[17,1],[23,1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);

    return [tmp1,tmp2,uu];
end;

# THE COMPONENT GROUP IS GENERATED BY u IN CHARACTERISTIC 2
#
# cls 24: D_5A_2
#
#handle24char3:=function(orbs,infos)
#    local tmp1,tmp2,o,info,rep,cent,uu;
#    o:=Classes(orbs)[24];
#    #info:=infos[24];
#
#    rep:=Representative(o);
#
#    #pr[87]=[1,1,2,3,3,2,1,1]
#    #pr[86]=[1,2,2,3,2,2,1,1]
#    #pr[119]=[2,3,4,6,5,4,3,1]
#    tmp1:=ApplyRootsReflections(rep,[87,86,119]);
#    uu:=Unipotent(chevalleyAdj(rep),[[8,-1]],Ordering(rep));
#    tmp1:=Conj(tmp1,uu);
#    tmp1:=ConjugateByTori(tmp1,[1,4,7],[-1,-1,-1,-1,-1]);
#
#    tmp2:=0;
#    #cent:=FromPositiveBorel(o,info[2]);
#    #tmp2:=ApplyRootsReflections(cent,[36,34,46,61]);
#    #uu:=Unipotent(chevalleyAdj(cent),[[2,-1],[4,1]],Ordering(cent));
#    #tmp2:=Conj(tmp2,uu);
#    #tmp2:=ConjugateByTori(tmp2,[2,3,4],[-1,-1,-1,-1]);
#
#    return [tmp1,tmp2];
#end;

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

    rep:=Representative(o);

    #pr[113]=[2,3,3,5,4,3,2,1]
    #pr[114]=[0,0,1,1,1,1,1,0]
    tmp1:=ApplyRootsReflections(rep,[113,114]);
    uu:=Unipotent(chevalleyAdj(rep),[[4,1],[6,1],[11,1],[13,1],[20,1],[27,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[113,114]);
    uu:=Unipotent(chevalleyAdj(cent),[[4,1],[6,1],[11,1],[13,1],[20,1],[27,1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);

    return [tmp1,tmp2];
end;
#
# THIS IS FOR THE D_7(a_2) class in E_8 characteristic 2
#
#set:=Filtered(subs,i->i<=63);
#Sort(set);
#uu:=Unipotent(alg,List(set,i->[i,avars[i]]),Ordering(rep));
#vals:=ShallowCopy(avars);
#vals[3]:=-vals[2];
#vals[4]:=-One(APR);
#vals[5]:=0;
#vals[6]:=One(APR);
#vals[7]:=0;
#vals[8]:=-vals[3];
#vals[11]:=-(vals[2]+vals[10]+One(APR));
#vals[12]:=-vals[3];
#vals[13]:=-(vals[2]+One(APR));
#vals[14]:=-vals[8];
#vals[15]:=-(vals[2]^2+vals[2]);
#vals[18]:=-(vals[2]^2+vals[2]+vals[17]);
#vals[19]:=-(vals[2]^2+vals[2]+vals[18]);
#vals[20]:=-(vals[2]+vals[10]+vals[13]);
#vals[21]:=-vals[15];
#vals[22]:=-(vals[2]^2+vals[19]);
#vals[26]:=-(vals[2]^2+vals[2]*vals[10]+vals[10]^2+vals[2]+vals[10]+vals[17]);
#vals[27]:=-(vals[2]^2+vals[2]*vals[13]+vals[10]*vals[13]+vals[11]*vals[13]+vals[2]+vals[10]+vals[11]+vals[13]+vals[17]+vals[20]+vals[26]+One(APR));
#vals[28]:=-vals[19];
#vals[29]:=-(vals[2]^3+vals[2]^2+vals[17]+vals[25]);

#vals[33]:=-(vals[2]^3+vals[2]^2*vals[10]+vals[2]*vals[10]+vals[2]+vals[17]+vals[32]);
#vals[34]:=-(vals[2]^3+vals[2]^2+vals[25]+vals[32]);
#vals[35]:=-(vals[2]^3+vals[2]^2+vals[17]+vals[25]+vals[34]);
#vals[36]:=-(vals[2]^2*vals[10]+vals[2]^2+vals[2]*vals[10]+vals[2]+vals[29]+vals[34]);
#vals[41]:=-(vals[2]^2*vals[17]+vals[2]*vals[17]+vals[17]^2+vals[25]);
#vals[43]:=-(vals[2]^4+vals[2]^2*vals[17]+vals[2]^2+vals[2]*vals[17]+vals[17]+vals[25]+vals[34]+vals[35]+vals[36]+vals[41]+vals[42]);
#vals[48]:=-(vals[2]^3*vals[10]+vals[2]^3+vals[2]*vals[10]+vals[10]*vals[17]+vals[2]+vals[17]+vals[25]+vals[32]+vals[42]);
#vals[49]:=-(vals[2]^4+vals[2]^3*vals[10]+vals[2]^2*vals[10]^2+vals[2]^3+vals[2]*vals[10]^2+vals[2]*vals[10]+vals[2]*vals[17]+vals[25]+vals[32]+vals[34]+vals[41]+vals[42]);

#vals[55]:=-(vals[2]^5+vals[2]^4*vals[10]+vals[2]^4+vals[2]^3+vals[2]^2*vals[10]+vals[2]^2+vals[17]^2+vals[17]*vals[25]+vals[25]+vals[50]);
#vals[61]:=-(vals[2]^5+vals[2]^4*vals[10]+vals[2]^4+vals[2]^3*vals[10]+vals[2]^3*vals[17]+vals[2]^2*vals[10]^2+vals[2]^2*vals[10]*vals[17]+vals[2]^2*vals[17]+vals[2]*vals[10]^2+vals[2]*vals[10]*vals[17]+vals[2]^2+vals[2]*vals[17]+vals[10]*vals[17]+vals[2]+vals[17]+vals[32]+vals[48]+vals[56]);
#vals[62]:=-(vals[2]^6+vals[2]^4*vals[10]^2+vals[2]^5+vals[2]^4*vals[10]+vals[2]^4*vals[17]+vals[2]^4+vals[2]^2*vals[10]^2+vals[2]^2*vals[10]*vals[17]+vals[2]^2*vals[10]*vals[25]+vals[2]^2*vals[17]^2+vals[2]^3+vals[2]^2*vals[10]+vals[2]^2*vals[25]+vals[2]*vals[10]*vals[17]+vals[2]*vals[10]*vals[25]+vals[17]^3+vals[2]*vals[17]+vals[2]*vals[25]+vals[17]^2+vals[17]*vals[25]+vals[32]^2+vals[25]+vals[50]);

#vv:=Value(uu,avars,vals);
#vvv:=Value(vv,avars,List(avars,i->0));

#
# cls 19: E_6(a_1)A_1
#
handle19char2:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[19];
    info:=infos[19];

    #pr[106]=[1,2,2,4,4,3,2,1]
    #pr[105]=[1,2,3,4,3,3,2,1]
    #pr[104]=[2,2,3,4,3,2,2,1]

    rep:=Representative(o);
    tmp1:=ApplyRootsReflections(rep,[106,105,104]);
    tmp1:=ConjNegUa(tmp1,4,1);
    uu:=Unipotent(chevalleyAdj(rep),[[2,1],[11,1],[10,1],[6,1],[16,1],[17,1],[23,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[106,105,104]);
    tmp2:=ConjNegUa(tmp2,4,1);
    uu:=Unipotent(chevalleyAdj(cent),[[2,1],[11,1],[10,1],[6,1],[16,1],[17,1],[23,1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);

    return [tmp1,tmp2,uu];
end;

#
# cls 18: E_7(a_3)
#
handle18char2:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[18];
    info:=infos[18];

    rep:=Representative(o);
    tmp1:=ConjNegUa(rep,3,1);
    tmp1:=ConjNegUa(tmp1,2,1);
    uu:=Unipotent(chevalleyAdj(rep),[[5,1],[18,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ConjNegUa(cent,3,1);
    tmp2:=ConjNegUa(tmp2,2,1);
    uu:=Unipotent(chevalleyAdj(cent),[[5,1],[18,1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);

    return [tmp1,tmp2];
end;

#
# cls 17: E_8(b_6)
#
handle17char2:=function(orbs,infos)
    local tmp1,tmp2,tmp3,tmp4,o,info,rep,cent,uu;
    o:=Classes(orbs)[17];
    info:=infos[17];

    rep:=Representative(o);
    tmp1:=ApplyRootsReflections(rep,[1,2,5]);
    uu:=Unipotent(chevalleyAdj(rep),[[7,1],[8,1],[18,1],[23,1],[31,1],[39,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);

    tmp3:=ConjugateByTori(rep,[1,2,5],[aaa,aaa,aaa^2]);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[1,2,5]);
    uu:=Unipotent(chevalleyAdj(cent),[[7,1],[8,1],[18,1],[23,1],[31,1],[39,1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);

    tmp4:=0;
    tmp4:=ConjugateByTori(cent,[1,2,5],[aaa,aaa,aaa^2]);

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

# THIS WORKS BUT HERE THE COMPONENT GROUP IS GENERATED BY U
#
# cls 16: D_7(a_1)
#
handle16char2:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[16];
    info:=infos[16];

    rep:=Representative(o);

    #pr[113]=[2,3,3,5,4,3,2,1]
    #pr[114]=[0,0,1,1,1,1,1,0]
    tmp1:=ApplyRootsReflections(rep,[113,114]);
    #uu:=Unipotent(chevalleyAdj(rep),[[4,1],[10,-1]],Ordering(rep)); si asta merge
    uu:=Unipotent(chevalleyAdj(rep),[[4,1],[11,-1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);
    #tmp1:=ConjugateByTori(tmp1,[2,3,4],[aaa,aaa,-1]); IN CHAR 2 I DON'T NEED THIS

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[113,114]);
    uu:=Unipotent(chevalleyAdj(cent),[[4,1],[11,-1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);
    #tmp2:=ConjugateByTori(tmp2,[2,3,4],[-1,-1,-1,-1]);

    return [tmp1,tmp2];
end;

#
# cls 12: E_8(a_6)
#
handle12char2:=function(orbs,infos)
    local tmp1,tmp2,tmp3,tmp4,o,info,rep,cent,uu;
    o:=Classes(orbs)[12];
    info:=infos[12];

    rep:=Representative(o);
    tmp1:=ApplyRootsReflections(rep,[2,3,8,13]);
    uu:=Unipotent(chevalleyAdj(rep),[[18,1],[29,1],[31,1],[33,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);

    tmp3:=ConjugateByTori(rep,[2,3,5,6,8],[aaa,aaa,aaa,aaa,aaa^2]);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[2,3,8,13]);
    uu:=Unipotent(chevalleyAdj(cent),[[18,1],[29,1],[31,1],[33,1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);

    tmp4:=0;
    tmp4:=ConjugateByTori(cent,[2,3,5,6,8],[aaa,aaa,aaa,aaa,aaa^2]);

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

#
# cls 10: E_8(b_5)
#
handle10char2:=function(orbs,infos)
    local tmp1,tmp2,tmp3,tmp4,o,info,rep,cent,uu;
    o:=Classes(orbs)[10];
    info:=infos[10];

    rep:=Representative(o);
    tmp1:=ApplyRootsReflections(rep,[1,2,5]);
    uu:=Unipotent(chevalleyAdj(rep),[[19,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);

    tmp3:=ConjugateByTori(rep,[1,2,5],[aaa,aaa,aaa^2]);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[1,2,5]);
    uu:=Unipotent(chevalleyAdj(cent),[[19,1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);

    tmp4:=0;
    tmp4:=ConjugateByTori(cent,[1,2,5],[aaa,aaa,aaa^2]);

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

#
# cls 8: E_8(a_5)       NOT FINISHED IN CHARACTERISTIC 2 !!!!
#
handle8char2:=function(orbs,infos)
    local tmp1,tmp2,tmp3,tmp4,o,info,rep,cent,uu;
    o:=Classes(orbs)[8];
    info:=infos[8];

    rep:=Representative(o);
    uu:=Unipotent(chevalleyAdj(rep),[ [ 2, 1 ], [ 3, 1 ], [ 6, 1 ], [ 8, 1 ], [ 21, 1 ], [ 25, 1 ] ],Ordering(rep));
    tmp1:=Conj(rep,uu);

    uu:=Unipotent(chevalleyAdj(rep),[ [ 1, 1 ], [ 7, 1 ], [ 10, 1 ], [ 12, 1 ], [ 14, 1 ], [ 16, 1 ], [ 19, 1 ], [ 21, 1 ], [ 23, 1 ], [ 25, 1 ], [ 28, 1 ], [ 37, 1 ], [ 44, 1 ], [ 46, 1 ], [ 53, 1 ], [ 55, 1 ], [ 59, 1 ], [ 66, 1 ], [ 67, 1 ], [ 69, 1 ], [ 71, 1 ], [ 74, 1 ], [ 76, 1 ], [ 77, 1 ], [ 81, 1 ], [ 107, 1 ], [ 110, 1 ], [ 114, 1 ], [ 115, 1 ] ],Ordering(rep));
    tmp3:=Conj(rep,uu);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    uu:=Unipotent(chevalleyAdj(cent),[ [ 2, 1 ], [ 3, 1 ], [ 6, 1 ], [ 8, 1 ], [ 21, 1 ], [ 25, 1 ] ],Ordering(cent));
    tmp2:=Conj(cent,uu);

    tmp4:=0;
    uu:=Unipotent(chevalleyAdj(cent),[ [ 1, 1 ], [ 7, 1 ], [ 10, 1 ], [ 12, 1 ], [ 14, 1 ], [ 16, 1 ], [ 19, 1 ], [ 21, 1 ], [ 23, 1 ], [ 25, 1 ], [ 28, 1 ], [ 37, 1 ], [ 44, 1 ], [ 46, 1 ], [ 53, 1 ], [ 55, 1 ], [ 59, 1 ], [ 66, 1 ], [ 67, 1 ], [ 69, 1 ], [ 71, 1 ], [ 74, 1 ], [ 76, 1 ], [ 77, 1 ], [ 81, 1 ], [ 107, 1 ], [ 110, 1 ], [ 114, 1 ], [ 115, 1 ] ],Ordering(cent));
    tmp4:=Conj(cent,uu);

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

# THE COMPONENT GROUP IS GENERATED BY u IN CHARACTERISTIC 2
#
# cls 7: E_8(b_4) [1,4,7,8,14,17,19,26]  this doesn't work in char2
#
#handle7char3:=function(orbs,infos)
#    local tmp1,tmp2,o,info,rep,cent;
#    o:=Classes(orbs)[7];
#    #info:=infos[7];
#
#    rep:=Representative(o);
#    tmp1:=ConjugateByTorus(rep,4,-1);# h_1 --> h_4
#
#    tmp2:=0;
#    #cent:=FromPositiveBorel(o,info[2]);
#    #tmp2:=ConjugateByTorus(cent,4,-1);# h_1 --> h_4
#
#    return [tmp1,tmp2];
#end;

#
# cls 5: E_8(a_4)       NOT FINISHED IN CHARACTERISTIC 2 !!!!
#
handle5char2:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[5];
    info:=infos[5];

    rep:=Representative(o);
    uu:=Unipotent(chevalleyAdj(rep),[[2,1],[3,1],[5,1],[7,1]],Ordering(rep));
    tmp1:=Conj(rep,uu);
    uu:=Unipotent(chevalleyAdj(rep),[ [ 1, 1 ], [ 4, 1 ], [ 6, 1 ], [ 8, 1 ], [ 11, 1 ], [ 12, 1 ], [ 14, 1 ], [ 16, 1 ], [ 19, 1 ], [ 25, 1 ], [ 27, 1 ], [ 28, 1 ], [ 30, 1 ], [ 31, 1 ], [ 32, 1 ], [ 33, 1 ], [ 35, 1 ], [ 37, 1 ], [ 38, 1 ], [ 39, 1 ], [ 40, 1 ], [ 41, 1 ], [ 45, 1 ], [ 46, 1 ], [ 48, 1 ], [ 49, 1 ], [ 52, 1 ], [ 53, 1 ], [ 59, 1 ], [ 64, 1 ], [ 66, 1 ], [ 67, 1 ], [ 72, 1 ], [ 77, 1 ], [ 80, 1 ], [ 81, 1 ], [ 85, 1 ], [ 98, 1 ], [ 99, 1 ], [ 100, 1 ], [ 102, 1 ], [ 106, 1 ], [ 107, 1 ], [ 110, 1 ], [ 112, 1 ], [ 114, 1 ], [ 115, 1 ], [ 116, 1 ] ],Ordering(rep));
    tmp1:=Conj(tmp1,uu);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    uu:=Unipotent(chevalleyAdj(cent),[[2,1],[3,1],[5,1],[7,1]],Ordering(cent));
    tmp2:=Conj(cent,uu);
    uu:=Unipotent(chevalleyAdj(cent),[ [ 1, 1 ], [ 4, 1 ], [ 6, 1 ], [ 8, 1 ], [ 11, 1 ], [ 12, 1 ], [ 14, 1 ], [ 16, 1 ], [ 19, 1 ], [ 25, 1 ], [ 27, 1 ], [ 28, 1 ], [ 30, 1 ], [ 31, 1 ], [ 32, 1 ], [ 33, 1 ], [ 35, 1 ], [ 37, 1 ], [ 38, 1 ], [ 39, 1 ], [ 40, 1 ], [ 41, 1 ], [ 45, 1 ], [ 46, 1 ], [ 48, 1 ], [ 49, 1 ], [ 52, 1 ], [ 53, 1 ], [ 59, 1 ], [ 64, 1 ], [ 66, 1 ], [ 67, 1 ], [ 72, 1 ], [ 77, 1 ], [ 80, 1 ], [ 81, 1 ], [ 85, 1 ], [ 98, 1 ], [ 99, 1 ], [ 100, 1 ], [ 102, 1 ], [ 106, 1 ], [ 107, 1 ], [ 110, 1 ], [ 112, 1 ], [ 114, 1 ], [ 115, 1 ], [ 116, 1 ] ],Ordering(cent));
    tmp2:=Conj(tmp2,uu);

    return [tmp1,tmp2,uu];
end;

#
# cls 4: E_8(a_3)       NOT FINISHED IN CHARACTERISTIC 2 !!!!
#
handle4char2:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[4];
    info:=infos[4];

    rep:=Representative(o);
    uu:=Unipotent(chevalleyAdj(rep),[[2,1],[3,1],[5,1],[19,1]],Ordering(rep));
    tmp1:=Conj(rep,uu);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    uu:=Unipotent(chevalleyAdj(cent),[[2,1],[3,1],[5,1],[19,1]],Ordering(cent));
    tmp2:=Conj(cent,uu);

    return [tmp1,tmp2];
end;


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

mizuno:=[[ 1, 2, 3, 4, 5, 6, 7,
           8, 9, 10, 11,
           12, 13, 14, 15,
           17, 16, 18, 19,
           20, 21, 22, 24,
           23, 25, 26, 27,
           28, 29, 31, 30,
           32, 33, 34, 35,
           36, 37, 38, 39,
           40, 41, 42, 43,
           44, 45, 46, 47,
           48, 49, 50, 51,
           52, 53, 54, 55,
           56, 57, 58, 59,
           60, 61, 62, 63 ],
          []];

reps:=[[[1,1],[2,1],[3,1],[4,1],[5,1],[6,1],[7,1]],            # E_7
       [[1,1],[3,1],[10,1],[9,1],[5,1],[6,1],[7,1]],           # E_7(a_1)
       [[1,1],[2,1],[3,1],[9,1],[11,1],[12,1],[13,1]],         # E_7(a_2)      ? 3/8
       [[1,1],[9,1],[10,1],[17,1],[22,1],[12,1],[7,1]],        # D_6+A_1
       [[1,1],[3,1],[9,1],[11,1],[12,1],[13,1]],               # E_6           ? 3/8
       [[8,1],[9,1],[11,1],[16,1],[12,1],[13,1]],              # E_6(a_1)
       [[1,1],[15,1],[16,1],[17,1],[6,1],[7,1]],               # D_6
       [[1,1],[15,1],[16,1],[17,1],[18,1],[23,1],[13,1]],      # D_6(a_1)+A_1
       [[14,1],[15,1],[16,1],[17,1],[12,1],[13,1]],            # A_6
       [[1,1],[15,1],[17,1],[23,1],[24,1],[7,1]],              # D_6(a_1)
       [[1,1],[15,1],[16,1],[17,1],[18,1],[19,1]],             # D_5+A_1
       [[14,1],[15,1],[27,1],[17,1],[18,1],[19,1],[23,1]],     # D_6(a_2)+A_1
       [[1,1],[15,1],[16,1],[18,1],[19,1]],                    # D_5
       [[8,1],[15,1],[16,1],[20,1],[24,1],[25,1]],             # (A_5+A_1)'
       [[15,1],[20,1],[21,1],[23,1],[13,1],[17,1]],            # D_6(a_2)
       [[20,1],[21,1],[23,1],[24,1],[28,1],[7,1]],             # (A_5+A_1)''
       [[28,1],[12,1],[13,1],[20,1],[21,1]],                   # A_5'
       [[14,1],[22,1],[23,1],[24,1],[25,1],[19,1]],            # D_5(a_1)+A_1
       [[20,1],[21,1],[23,1],[24,1],[7,1]],                    # A_5''
       [[20,1],[21,1],[22,1],[23,1],[24,1],[25,1]],            # A_4+A_2
       [[8,1],[28,1],[24,1],[25,1],[30,1]],                    # D_5(a_1)
       [[20,1],[21,1],[24,1],[30,1]],                          # A_4
       [[20,1],[21,1],[24,1],[28,1],[30,1]],                   # A_4+A_1
       [[1,1],[28,1],[29,1],[30,1],[31,1]],                    # D_4+A_1
       [[26,1],[27,1],[28,1],[29,1],[30,1],[31,1]],            # A_3+A_2+A_1
       [[1,1],[28,1],[29,1],[30,1]],                           # D_4
       [[27,1],[29,1],[30,1],[31,1],[32,1],[37,1]],            # A_3+A_2
#       [[27,1],[29,1],[30,1],[31,1],[32,1]],                   # (A_3+A_2)_2
       [[28,1],[30,1],[31,1],[33,1],[34,1]],                   # D_4(a_1)+A_1
       [[28,1],[30,1],[33,1],[34,1]],                          # D_4(a_1)
       [[30,1],[31,1],[32,1],[33,1],[40,1]],                   # A_3+2A_1
       [[26,1],[27,1],[40,1],[41,1]],                          # (A_3+A_1)'
       [[30,1],[31,1],[32,1],[33,1]],                          # (A_3+A_1)''
       [[20,1],[21,1],[49,1]],                                 # A_3
       [[34,1],[36,1],[37,1],[38,1],[40,1]],                   # 2A_2+A_1
       [[34,1],[36,1],[38,1],[40,1]],                          # 2A_2
       [[37,1],[38,1],[39,1],[40,1],[41,1]],                   # A_2+3A_1
       [[42,1],[43,1],[44,1],[45,1]],                          # A_2+2A_1
       [[44,1],[46,1],[49,1]],                                 # A_2+A_1
       [[44,1],[46,1]],                                        # A_2
       [[47,1],[48,1],[49,1],[53,1]],                          # 4A_1
       [[53,1],[54,1],[55,1]],                                 # (3A_1)'
       [[47,1],[48,1],[49,1]],                                 # (3A_1)''
       [[58,1],[59,1]],                                        # 2A_1
       [[63,1]]                                                # A_1
       ];

labels:=["E_7",
         "E_7(a_1)",
         "E_7(a_2)",
         "D_6+A_1",
         "E_6",
         "E_6(a_1)",
         "D_6",
         "D_6(a_1)+A_1",
         "A_6",
         "D_6(a_1)",
         "D_5+A_1",
         "D_6(a_2)+A_1",
         "D_5",
         "(A_5+A_1)'",
         "D_6(a_2)",
         "(A_5+A_1)''",
         "A_5'",
         "D_5(a_1)+A_1",
         "A_5''",
         "A_4+A_2",
         "D_5(a_1)",
         "A_4",
         "A_4+A_1",
         "D_4+A_1",
         "A_3+A_2+A_1",
         "D_4",
         "A_3+A_2",
#         "(A_3+A_2)_2"
         "D_4(a_1)+A_1",
         "D_4(a_1)",
         "A_3+2A_1",
         "(A_3+A_1)'",
         "(A_3+A_1)''",
         "A_3",
         "2A_2+A_1",
         "2A_2",
         "A_2+3A_1",
         "A_2+2A_1",
         "A_2+A_1",
         "A_2",
         "4A_1",
         "(3A_1)'",
         "(3A_1)''",
         "2A_1",
         "A_1",
         ];

I:=[ [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ] ];#, [  ] ];

II:=[ [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ] ];#, [  ] ];

tori:=[ [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ] ];#, [  ] ];

#blocks:=[];
blocks:=[[[27, 3], [19, 1], [15, 1], [9, 2]], [[27, 1], [19, 1], [18, 2], [15, 1], [9, 4]], [[19, 1], [18, 2], [15, 1], [9, 7]], [[19, 1], [17, 1], [15, 2], [11, 2], [9, 3], [6, 2], [3, 2]], [[19, 1], [15, 4], [9, 5], [3, 3]], [[9, 14], [7, 1]], [[19, 1], [16, 2], [15, 1], [11, 1], [10, 2], [9, 2], [6, 2], [3, 1], [1, 3]], [[9, 12], [7, 1], [6, 2], [3, 2]], [[9, 11], [7, 4], [3, 1], [1, 3]], [[9, 11], [8, 2], [6, 2], [3, 1], [1, 3]], [[9, 11], [8, 2], [6, 2], [3, 1], [1, 3]], [[9, 9], [7, 1], [6, 2], [5, 3], [3, 6]], [[9, 11], [7, 4], [1, 6]], [[9, 6], [7, 1], [6, 8], [3, 8]], [[9, 9], [7, 1], [6, 2], [5, 1], [4, 4], [3, 3], [1, 3]], [[9, 9], [3, 16], [2, 2]], [[9, 5], [8, 2], [6, 8], [3, 7], [1, 3]], [[9, 5], [7, 4], [6, 6], [3, 7], [1, 3]], [[9, 9], [3, 15], [1, 7]], [[9, 3], [7, 1], [6, 8], [5, 3], [3, 12]], [[9, 3], [8, 4], [7, 2], [6, 6], [3, 4], [2, 4], [1, 4]], [[9, 1], [7, 7], [5, 9], [3, 7], [1, 9]], [[9, 1], [8, 2], [7, 3], [6, 4], [5, 5], [4, 2], [3, 8], [2, 2], [1, 2]], [[9, 2], [8, 4], [7, 5], [6, 4], [3, 2], [2, 4], [1, 10]], [[7, 1], [6, 8], [5, 3], [3, 21]], [[9, 2], [7, 13], [3, 1], [1, 21]], [[7, 1], [6, 8], [5, 1], [4, 4], [3, 18], [1, 3]], [[7, 1], [6, 6], [5, 4], [4, 4], [3, 16], [1, 6]], [[7, 1], [6, 2], [5, 12], [3, 15], [1, 9]], [[7, 1], [6, 4], [5, 7], [4, 4], [3, 11], [2, 6], [1, 6]], [[7, 1], [6, 2], [5, 7], [4, 10], [3, 6], [2, 6], [1, 9]], [[7, 1], [5, 15], [3, 10], [1, 21]], [[7, 1], [5, 7], [4, 16], [3, 1], [1, 24]], [[3, 43], [2, 2]], [[3, 42], [1, 7]], [[3, 42], [1, 7]], [[3, 37], [2, 8], [1, 6]], [[3, 34], [2, 8], [1, 15]], [[3, 33], [1, 34]], [[3, 28], [2, 14], [1, 21]], [[3, 19], [2, 26], [1, 24]], [[3, 27], [1, 52]], [[3, 10], [2, 32], [1, 39]], [[3, 1], [2, 32], [1, 66]]];

#ords:=[];
ords:=[ 27, 27, 27, 27, 27, 9, 27, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3 ];

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


text:
Load(graph);
Import(graph);
conf := SparseDefaultConf;
opts := conf.getOpts(V(1));
opts.vec := true;
i := Ind(opts.symbol[1]);
front := var.fresh_t("frontier",TArray(TSemiring_Boolean(TInt), opts.symbol[1]));
fdataofs := FDataOfs(front, opts.symbol[1], 0);
md := MakeDiag(fdataofs);
eye := I(i);
sid := SUM(eye, md);
matmul := MatMul(opts.symbol[1],opts.symbol[1]);
spmv := matmul * sid;
spmv2 := Rewrite(spmv, RulesIDiag, opts);
spmv3 := Rewrite(spmv2, RulesMatMul, opts);
rt := RandomRuleTree(last, opts);
srt := SumsRuleTree(rt, opts);
push := Rewrite(Copy(srt), RulesSigSPMV2, opts);
cs := CodeSums(BFSKernel(PushBFSKernel(push)), opts);
PrintCode("spmv", cs, opts);


text:
for i in [1 .. 5] do
    for j in [1 .. i] do
        Print("*");
    od;
    Print("\n");
od;

# *
# **
# ***
# ****
# *****