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

# equality of scalars
scalareq := (a,b) -> AbsFloat(a-b) < 1e-10;
# equality of lists
lseq := (a,b) -> Length(a)=Length(b) and ForAll([1..Length(a)], scalareq(a[i], b[i]));

lequiv := (a,b) -> Cond(
    lseq(a,b),   "eq",
    lseq(a, -b), "eqneg",
    lseq(a, Reversed(b)), "rev",
    lseq(a, -Reversed(b)), "revneg",
    "none");

tovec := lst -> TransposedMat([lst]);
tolst := vec -> TransposedMat(vec)[1];
rpart := vec -> List(vec, x->ReComplex(Complex(x)));
ipart := vec -> List(vec, x->ImComplex(Complex(x)));

rand := arg -> RandomList([-100..100]);
#rand := x->x;

##
## Constructors for symmetric sequences
##
real := n -> List([1..n], rand);
cplx := (rsym, isym) -> (n -> rsym(n) + E(4)*isym(n));

even0 := n -> let(len := Int((n-1)/2), lst := List([2..len+1], rand),
    Flat([rand(1), lst, When(IsEvenInt(n),rand(len+2),[]), Reversed(lst)]));

even00 := n -> let(len := Int((n-1)/2), lst := List([1..len], rand),
    Flat([0, lst, When(IsEvenInt(n),rand(len+1),[]), Reversed(lst)]));
odd0 := n -> let(len := Int((n-1)/2), lst := List([2..len+1], rand),
    Flat([rand(1), lst, When(IsEvenInt(n),0,[]), -Reversed(lst)]));
odd00 := n -> let(len := Int((n-1)/2), lst := List([1..len], rand),
    Flat([0, lst, When(IsEvenInt(n),0,[]), -Reversed(lst)]));

even1 := n -> let(len := Int(n/2), lst := List([1..len], rand),
    Flat([lst, When(IsOddInt(n),rand(len+1),[]), Reversed(lst)]));
odd1 := n -> let(len := Int(n/2), lst := List([1..len], rand),
    Flat([lst, When(IsOddInt(n),0,[]), -Reversed(lst)]));

jeven1 := n -> Checked(IsEvenInt(n), let(lst:=List([1..n/2], rand), Flat([lst,lst])));
jodd1  := n -> Checked(IsEvenInt(n), let(lst:=List([1..n/2], rand), Flat([lst,-lst])));
jeven0 := n -> Checked(IsOddInt(n), Concatenation([rand(0)], jeven1(n-1)));
jeven00 := n -> Checked(IsOddInt(n), Concatenation([0], jeven1(n-1)));
jodd0  := n -> Checked(IsOddInt(n), Concatenation([rand(0)], jodd1(n-1)));
jodd00 := n -> Checked(IsOddInt(n), Concatenation([0], jodd1(n-1)));

ce0 := cplx(even0, odd00);
co0 := cplx(odd0, even00);
ce1 := cplx(even1, odd1);
co1 := cplx(odd1, even1);

jce0 := cplx(jeven0, jodd00);
jco0 := cplx(jodd0, jeven00);
jce1 := cplx(jeven1, jodd1);
jco1 := cplx(jodd1, jeven1);

raderce0 := n -> Checked(IsPrime(n+1), Drop(
    tolst(MatSPL(DFT_Rader.raderMid(n+1,1,PrimitiveRootMod(n+1))) *
	  tovec(Concat([rand(0)], ce0(n)))), 1));

upsample0 := n -> Checked(IsEvenInt(n),  List([1..n], i -> When(IsEvenInt(i),0,rand())));

upsample1 := n -> Checked(IsEvenInt(n),  List([1..n], i -> When(IsOddInt(i),0,rand())));

_even := n -> ((n+1) mod 2);
_odd := n -> (n mod 2);

_minsetPoints := function(points, n)
    local uniq, cur, p;
    cur := []; uniq := Set([]);
    for p in points do
        if not (p in uniq) and p < n then
	    Add(cur, p);
	    AddSet(uniq, p);
	fi;
    od;
    return cur;
end;

_CommuteExt := function(perm, sym, refl_diff, n, drop0)
    local mat, N, ext, points, minset, pmat;
    mat:=MatSPL(Prm(perm));
    N := perm.domain();
    ext := sym.ext(N,false);
    points := List(Refl(n, N-refl_diff, N, perm).tolist(), x->x.ev());
    minset := When(not drop0, _minsetPoints(points, n), Drop(_minsetPoints(points,n), 1)-1);
    pmat := MatSPL(Scat(FList(minset).setRange(Cols(ext))));
    return
       mat * MatSPL(ext) * pmat;
end;

# CommuteExt(spl, sym, scaled) - find X, such that spl * sym.ext(n, scaled) = X * pspl
#  X = spl * sym.ext(n, scaled) * pspl^-1
#

CommuteExt1 := (perm, sym) -> _CommuteExt(perm, sym, -1, Cols(sym.ext(perm.domain(),false)), false);
CommuteExt0 := (perm, sym) -> _CommuteExt(perm, sym, 0,  Cols(sym.ext(perm.domain(),false)), false);
CommuteExt00 := (perm, sym) -> _CommuteExt(perm, sym, 0, 1+Cols(sym.ext(perm.domain(),false)), true);

##
## Constructors for symmetric sequences
##
Class(Symmetry, rec(
    isSymmetry := true,
    __call__ := (self, n) >> self.seq(n),
    red := (self, n, scaled) >> self.ext(n, scaled).transpose(),
    verify := (self, n) >> MatSPL(self.red(n,false)*self.ext(n, true))
));

IsSymmetry := x -> IsRec(x) and IsBound(x.isSymmetry) and x.isSymmetry;

_gathI := (nn, top_hole, bot_hole) -> let(n:=Int(nn), Gath(fAdd(n, n-top_hole-bot_hole, top_hole)));
_gathI2 := (nn, top_hole, bot_hole) -> let(n:=Int(nn), Gath(H(n, Int((n-top_hole-bot_hole+1)/2), top_hole, 2)));
_gathJ := (nn, top_hole, bot_hole) -> let(n:=Int(nn), Gath(fCompose(fAdd(n, n-top_hole-bot_hole, top_hole), J(n-top_hole-bot_hole))));

Class(GathExtend, Sym, rec(def := (n, symmetry) -> symmetry.ext(n, true)));
Class(GathExtendU, Sym, rec(def := (n, symmetry) -> symmetry.ext(n, false)));
Class(GathExtendZ, Sym, rec(def := (n, symmetry) -> DirectSum(I(1),-J(n-1))*symmetry.ext(n, true)));
Class(GathExtendZU, Sym, rec(def := (n, symmetry) -> DirectSum(I(1),-J(n-1))*symmetry.ext(n, false)));
Class(ScatReduce, Sym, rec(def := (n, symmetry) -> symmetry.red(n, true)));
Class(ScatReduceU, Sym, rec(def := (n, symmetry) -> symmetry.red(n, false)));

Class(EOSymmetry, Symmetry, rec(tsize := "cols"));

Class(Upsample0, EOSymmetry, rec(
    seq:=  n -> Checked(IsEvenInt(n),  List([1..n], i -> When(IsEvenInt(i),0,rand()))),
    ext := (n, scaled) -> Checked(IsEvenInt(n), Tensor(I(n/2),Mat([[1],[0]]))),
    red := (n, scaled) -> Checked(IsEvenInt(n), Tensor(I(n/2),Mat([[1,0]]))),
));

Class(Upsample1, EOSymmetry, rec(
    seq:=  n -> Checked(IsEvenInt(n),  List([1..n], i -> When(IsEvenInt(i),0,rand()))),
    ext := (n, scaled) -> Checked(IsEvenInt(n), Tensor(I(n/2),Mat([[0],[1]]))),
    red := (n, scaled) -> Checked(IsEvenInt(n), Tensor(I(n/2),Mat([[0,1]]))),
));

Class(Even0, EOSymmetry, rec(
    seq := n -> let(len := Int((n-1)/2), lst := List([2..len+1], rand),
	Flat([rand(1), lst, When(IsEvenInt(n),rand(len+2),[]), Reversed(lst)])),
    ext := (n, scaled) -> let(hf := When(scaled, 1/2, 1),
	DirectSum(I(1),
	            VStack(
			DirectSum(hf*I(Int((n-1)/2)), I(_even(n))),
			hf*_gathJ(n/2, 0, _even(n))))),
#    red := n -> Gath(fAdd(n, Int((n+2)/2), 0))
));

Class(Even0star, EOSymmetry, rec(
    seq := n -> let(len := Int((n-1)/2), lst := List([2..len+1], rand),
	Flat([rand(1), lst, When(IsEvenInt(n),0,[]), Reversed(lst)])),
    ext := (n, scaled) -> let(hf := When(scaled, 1/2, 1),
	DirectSum(I(1),
	            VStack(
			hf*I(Int((n-1)/2)),
                        When(IsEvenInt(n), O(1, Int((n-1)/2)), []),
			hf*J(Int((n-1)/2))))),
#    red := n -> Gath(fAdd(n, Int((n+2)/2), 0))
));

Class(Even00, EOSymmetry, rec(
    seq := n -> let(len := Int((n-1)/2), lst := List([2..len+1], rand),
	Flat([0, lst, When(IsEvenInt(n),rand(len+2),[]), Reversed(lst)])),
    ext := (n, scaled) -> let(hf := When(scaled, 1/2, 1),
	VStack(O(1,Int(n/2)),
	       DirectSum(hf*I(Int((n-1)/2)), I(_even(n))),
	       hf*_gathJ(n/2, 0, _even(n)))),
#    red := n -> Gath(fAdd(n, Int(n/2), 1))
));

Class(Even00star, EOSymmetry, rec(
    seq := n -> let(len := Int((n-1)/2), lst := List([2..len+1], rand),
	Flat([0, lst, When(IsEvenInt(n),0,[]), Reversed(lst)])),
    ext := (n, scaled) -> let(hf := When(scaled, 1/2, 1),
	VStack(O(1,Int(n/2)),
	       DirectSum(hf*I(Int((n-1)/2)), O(_even(n))),
	       hf*_gathJ(n/2, 0, _even(n)))),
#    red := n -> Gath(fAdd(n, Int(n/2), 1))
));

Class(Even1, EOSymmetry, rec(
    seq := n -> let(len := Int(n/2), lst := List([1..len], rand),
	Flat([lst, When(IsOddInt(n),rand(len+1),[]), Reversed(lst)])),
    ext := (n, scaled) -> let(hf := When(scaled, 1/2, 1),
	When(IsEvenInt(n),
	    hf * VStack(I(n/2), J(n/2)),
	    VStack(DirectSum(hf*I((n-1)/2), I(1)), hf*_gathJ((n+1)/2, 0, 1))))
#    red := n -> Gath(fAdd(n, Int((n+1)/2), 0))
));

Class(Odd0, EOSymmetry, rec(
    seq := n -> let(len := Int((n-1)/2), lst := List([2..len+1], rand),
	Flat([rand(1), lst, When(IsEvenInt(n),0,[]), -Reversed(lst)])),
    ext := (n, scaled) -> let(hf := When(scaled, 1/2, 1),
	DirectSum(I(1),
	         hf*VStack(I(Int((n-1)/2)),
		            When(IsEvenInt(n), O(1, Int((n-1)/2)), []),
			    -J(Int((n-1)/2)))))
#    red := n -> Gath(fAdd(n, Int((n+2)/2), 0))
));

Class(Odd00, EOSymmetry, rec(
    seq := n -> let(len := Int((n-1)/2), lst := List([2..len+1], rand),
	Flat([0, lst, When(IsEvenInt(n),0,[]), -Reversed(lst)])),
    ext := (n, scaled) -> let(hf := When(scaled, 1/2, 1),
	        hf * VStack(O(1,Int((n-1)/2)),
	                     I(Int((n-1)/2)),
			     When(IsEvenInt(n), O(1,Int((n-1)/2)), []),
	                     -J(Int((n-1)/2)))),
#    red := n -> Gath(fAdd(n, Int((n-1)/2), 1))
));

Class(Odd00star, EOSymmetry, rec(
    seq := n -> let(len := Int((n-1)/2), lst := List([2..len+1], rand),
	Flat([0, lst, When(IsEvenInt(n),rand(len+2),[]), -Reversed(lst)])),
    ext := (n, scaled) -> let(hf := When(scaled, 1/2, 1),
	        hf * VStack(O(1,Int((n-1)/2)),
	                    I(Int((n-1)/2)),
                            When(IsEvenInt(n), RowVec( Replicate(Int((n-3)/2),0)::[1]  ), []),
                            -J(Int((n-1)/2)))),
#    red := n -> Gath(fAdd(n, Int((n-1)/2), 1))
));

Class(Odd1, EOSymmetry, rec(
    seq := n -> let(len := Int(n/2), lst := List([1..len], rand),
	Flat([lst, When(IsOddInt(n),0,[]), -Reversed(lst)])),
    ext := (n, scaled) -> let(hf := When(scaled, 1/2, 1),
	     hf * VStack(I(Int(n/2)),
		         When(IsOddInt(n), O(1, Int(n/2)),[]),
			 -J(Int(n/2)))),
#    red := n -> Gath(fAdd(n, Int(n/2), 0))
));

Class(Real, EOSymmetry, rec(
    seq := n -> List([1..n], rand),
    ext := (n, scaled) -> I(n),
    tsize := "rows", # does not matter, since ext is square
));

# NOTE: red != ext^T, how to handle?
# - silently drop scaling?

# DirectSum((I tensor W), j * I(1))
#
Class(Conj_Base, Symmetry, rec(
    __call__ := (self, w, j) >> Checked(IsMat(w), DimensionsMat(w)=[2,2],
	WithBases(self, rec(
		w := w,
		j := j,
		operations := PrintOps
    ))),
    tsize := "rows",
    print := self >> Print(self.name, "(", self.w, ", ", self.j, ")"),
    W := self >> When(IsBound(self.isOdd) and self.isOdd, self.w * DiagonalMat([1, -1]), self.w),
    invertW := self >> let(base := self.__bases__[1], base(TransposedMat(self.w^-1), 1/self.j)),
    seq := (self, n) >> self.re.seq(n) + E(4) * self.im.seq(n),
    verify := (self, n) >> MatSPL(self.red(n,false)*self.ext(n, false))
));

Class(Conj0_Base, Conj_Base, rec(
    P := n -> DirectSum(I(1), L_or_OS(n-1, Int(n/2)) * DirectSum(I(Int(n/2)), J(Int((n-1)/2)))),
    ext := (self, n, scaled) >> self.invertW().red(n, scaled).transpose(),
    red := (self, n, scaled) >> let(
	W := Mat(When(scaled, self.W(), self.W())),
        self.PP(n) * DirectSum(I(1), Tensor(I(Int((n-1)/2)), W), self.j * I(_even(n))) * self.P(n)
    )
));

Class(Conj1_Base, Conj_Base, rec(
    P := n -> L_or_OS(n, Int((n+1)/2)) * DirectSum(I(Int((n+1)/2)), J(Int(n/2))),
    seq := (self, n) >> self.re.seq(n) + E(4) * self.im.seq(n),
    ext := (self, n, scaled) >> self.invertW().red(n, scaled).transpose(),
    red := (self, n, scaled) >> let(
	W := Mat(When(scaled, self.W(), self.W())),
        self.PP(n) * DirectSum(Tensor(I(Int(n/2)), W), self.j * I(_odd(n))) * self.P(n)
    )
));

Class(Conj0_R2R, Conj0_Base, rec(PP := (self, n) >> self.P(n).transpose()));
Class(Conj1_R2R, Conj1_Base, rec(PP := (self, n) >> self.P(n).transpose()));

Class(CE0_R2R, Conj0_R2R, rec(re := Even0, im := Odd00,  isOdd := false));
Class(CO0_R2R, Conj0_R2R, rec(re := Odd0,  im := Even00, isOdd := true, j := -E(4)));
Class(CE1_R2R, Conj1_R2R, rec(re := Even1, im := Odd1,   isOdd := false));
Class(CO1_R2R, Conj1_R2R, rec(re := Odd1,  im := Even1,  isOdd := true, j := -E(4)));

_m := n -> When(n=1, Mat([[1], [0]]), I(0));
_mj := n -> When(n=1, Mat([[0], [1]]), I(0));

Class(Conj0_R2Cpx, Conj0_Base, rec(PP := (self, n) >>
	DirectSum(_m(1), I(n-1-_even(n)), _m(_even(n)))));

Class(Conj1_R2Cpx, Conj1_Base, rec(PP := (self, n) >>
	DirectSum(I(n-_odd(n)), When(self.isOdd, _mj(_odd(n)), _m(_odd(n))))));

Class(CE0_R2Cpx, Conj0_R2Cpx, rec(re := Even0, im := Odd00,  isOdd := false));
Class(CO0_R2Cpx, Conj0_R2Cpx, rec(re := Odd0,  im := Even00, isOdd := true, j := -E(4)));
Class(CE1_R2Cpx, Conj1_R2Cpx, rec(re := Even1, im := Odd1,   isOdd := false));
Class(CO1_R2Cpx, Conj1_R2Cpx, rec(re := Odd1,  im := Even1,  isOdd := true, j := -E(4)));

W_RFT := 1/2*[[1, 1], [-E(4), E(4)]];
W_RFTT := 1/2*[[1, 1], [E(4), -E(4)]];
W_URFT := [[1, 1], [-E(4), E(4)]];
W_URFTT := [[1, 1], [E(4), -E(4)]];
W_DHT := 1/2*[[1-E(4), 1+E(4)], [1+E(4), 1-E(4)]];

rsym := function(x)
    local n, nc, nf, mid, uniq, zeros, fst_0, mid_0, fstmid;
    n := Length(x); nc := Int((n+1)/2); nf := Int(n/2); mid := nf;
    if n=0 then return "empty"; fi;

    fst_0 := scalareq(x[1], 0);
    mid_0 := scalareq(x[nf+1], 0);
    fstmid := [Cond(fst_0,0, 1), Cond(mid_0, 0, 1)];

    if ForAll(x, e->scalareq(e,0)) then return "zero";

    # normal symmetries
    elif ForAll([1..nc-1], i->scalareq(x[1+ i],  x[1+ (n-i) mod n])) then
        fstmid[2] := Cond(IsOddInt(n), 1, fstmid[2]);
        if   fstmid=[0,0] then return "even00*";
	elif fstmid=[0,1] then return "even00";
        elif fstmid=[1,0] then return "even0*";
        elif fstmid=[1,1] then return "even0";
        fi;
    elif ForAll([1..nc-1], i->scalareq(x[1+ i], -x[1+ (n-i) mod n])) then
        fstmid[2] := Cond(IsOddInt(n), 0, fstmid[2]);
        if   fstmid=[0,0] then return "odd00";
	elif fstmid=[0,1] then return "odd00*";
        elif fstmid=[1,0] then return "odd0";
        elif fstmid=[1,1] then return "odd0*";
        fi;
    elif ForAll([0..mid-1], i->scalareq(x[1+ i],  x[1+ (n-1-i)])) then
        if IsOddInt(n) and mid_0 then return "even1*";
        else return "even1";
        fi;
    elif ForAll([0..mid-1], i->scalareq(x[1+ i], -x[1+ (n-1-i)])) then
        if IsOddInt(n) and not mid_0 then return "odd1*";
        else return "odd1";
        fi;

    # symmetries with flipped bottom, i.e. abcabc
    elif n>2 and IsOddInt(n) and ForAll([1..nc-1], i->scalareq(x[1+ i],  x[1+ mid+i])) then
	if fst_0 then return "jeven00"; else return "jeven0"; fi;
    elif n>2 and IsOddInt(n) and ForAll([1..nc-1], i->scalareq(x[1+ i], -x[1+ mid+i])) then
	if fst_0 then return "jodd00"; else return "jodd0"; fi;
    elif n>2 and ForAll([0..mid-1], i->scalareq(x[1+ i],  x[1+ mid+i])) then return "jeven1";
    elif n>2 and ForAll([0..mid-1], i->scalareq(x[1+ i], -x[1+ mid+i])) then return "jodd1";

    # rader symmetries, similar to normal but every other sample in bottom half is negated
    elif  ForAll([1..nc-1], i->scalareq(x[1+i], (-1)^i * x[1+(n-i) mod n])) then
	if fst_0 then return "rader_even00"; else return "rader_even0"; fi;
    elif  ForAll([1..nc-1], i->scalareq(x[1+i], - (-1)^i * x[1+(n-i) mod n])) then
	if fst_0 then return "rader_odd00"; else return "rader_odd0"; fi;

    # symmetries with normal shape, but elements negated in some non-recognized way
    elif  ForAll([1..nc-1], i->scalareq(x[1+i], x[1+(n-i) mod n]) or scalareq(x[1+i], -x[1+(n-i) mod n])) then
	if fst_0 then return "eo00"; else return "eo0"; fi;
    elif  ForAll([0..mid-1], i->scalareq(x[1+i], x[1+(n-1-i)]) or scalareq(x[1+i], -x[1+(n-1-i)])) then
	return "eo1";
    else
	uniq := Length(Set(List(x, e->V(AbsFloat(e)))));
	zeros := Length(Filtered(x, e->scalareq(e,0)));
	if uniq <= Int((n+2)/2) then
            # no obvious symmetry, but only approx half unique elements
	    # so sequence is still redundant
	    return Concat("red_", String(uniq-zeros), "/", String(zeros), "/",String(Length(x)));
	else return "none";
	fi;
    fi;
end;

risym := function(x)
    local rs, is, re, im;
    re := List(x, e->ReComplex(Complex(e)));
    im := List(x, e->ImComplex(Complex(e)));
    rs := rsym(re);
    is := rsym(im);
    return [rs,is];
end;

csym := function(x)
   local s;
   s := sym(x);
   if s = ["even0", "odd00"] then return "ce0";
   elif s = ["odd0", "even00"] then return "co0";
   elif s = ["even1", "odd1"] then return "ce1";
   elif s = ["odd1", "even1"] then return "co1";
   elif s = ["jeven1", "jodd1"] then return "jce1";
   elif s = ["jodd1", "jeven1"] then return "jco1";
   elif s = ["rader_even0", "rader_odd00"] then return "rader_ce0";
   elif s = ["rader_odd0", "rader_even00"] then return "rader_co0";
   elif s = ["none", "zero"] then return "real";
   elif s = ["zero", "none"] then return "imag";
   elif s = ["zero", "zero"] then return "zero";
   elif s = ["none", "none"] then return "cplx";
   else return s;
   fi;
end;

part := (v, str, n) -> v{1 + [ 0 .. (Length(v) / str - 1) ] * str + n};

fxi_inputs := (inp, fsize) -> let(
    m := Length(inp)/fsize,
    List([0..m-1], x->part(inp,m,x)));

fxi_rsyms := (inp, fsize) -> List(fxi_inputs(inp, fsize), rsym);
fxi_syms := (inp, fsize) -> List(fxi_inputs(inp, fsize), risym);

transf_syms := (n, transforms, syms) ->
   List(syms, s ->
       List(transforms, transform ->
	   let(inpvec := TransposedMat([s(n)]),
	       outvec := TransposedMat(MatSPL(transform(n))*inpvec)[1],
	       outsym := sym(outvec),
	       Cond(outsym[2]="zero", outsym[1],
		    outsym[1]="zero", Concat("j*", outsym[2]),
		    "none"))));

ctransf_syms := (n, transforms, syms) ->
   List(syms, s ->
       List(transforms, transform ->
	   let(inpvec := TransposedMat([s(n)]),
	       outvec := TransposedMat(MatSPL(transform(n))*inpvec)[1],
	       risym(outvec))));

# transf_syms(18, [ DFT1, DFT2, DFT3, DFT4 ], [even00, even1, odd00, odd1]);

# fxi_syms(even0(36),6);
# fxi_syms(even00(36),6);
# fxi_syms(even1(36),6);
# fxi_syms(odd0(36),6);
# fxi_syms(odd00(36),6);
# fxi_syms(odd1(36),6);

#    if rs="none" and is="none" then return "none";
#    elif rs="none" and is="zero" then return "real";
#    elif rs="even" and is="odd" then return "ce";
#    elif rs="even" and is="odd" then return "ce";
#end;
twpow :=  (N,n) -> List([0..N-1], j -> (j mod n) * Int(j / n));

tw1 :=  (N,n) -> List([0..N-1], j -> E(4*N)^((2*(j mod n))   * (2*Int(j / n))));
tw2 :=  (N,n) -> List([0..N-1], j -> E(4*N)^((2*(j mod n))   * (2*Int(j / n)+1)));
tw3 :=  (N,n) -> List([0..N-1], j -> E(4*N)^((2*(j mod n)+1) * (2*Int(j / n))));
tw4 :=  (N,n) -> List([0..N-1], j -> E(4*N)^((2*(j mod n)+1) * (2*Int(j / n)+1)));
dct2diag := N -> List([0..N-1], i->E(2*N)^i);
dct4diag := N -> List([0..N-1], i->E(4*N)^(2*i+1));
idct2diag := N -> List([0..N-1], i->E(2*N)^-i);
listmul := (l1,l2) -> List([1..Length(l1)], i->l1[i]*l2[i]);

ComplexFFT2 := input -> TransposedMat(MatSPL(DFT2(Length(input))) * TransposedMat([input]))[1];
ComplexFFT3 := input -> TransposedMat(MatSPL(DFT3(Length(input))) * TransposedMat([input]))[1];
ComplexFFT4 := input -> TransposedMat(MatSPL(DFT4(Length(input))) * TransposedMat([input]))[1];

# Macro Extensions
#
MacroExt0 := (k,a,b,c,d) -> Checked(IsSymmetry(a), IsSPL(b), IsSPL(d), b.dims()=d.dims(),
    IsOddInt(k) or IsSymmetry(c),
    Cond(IsOddInt(k),
         let(kk := (k-1)/2, n := Cols(b),
             top := DirectSum(GathExtendU(n, a), Tensor(I(kk), b)),
             VStack(top,
                    Tensor(J(kk), d) * Gath(fAdd(Cols(top), kk*Cols(d), Cols(a.ext(n,true)))))),
         let(kk := (k-2)/2, n := Cols(b),
             top := DirectSum(GathExtendU(n, a), Tensor(I(kk), b), GathExtendU(n, c)),
             VStack(top,
                    Tensor(J(kk), d) * Gath(fAdd(Cols(top), kk*Cols(d), Cols(a.ext(n,true))))))
         ));

MacroExt1 := (k,b,c,d) -> Checked(IsSPL(b), IsSPL(d), b.dims()=d.dims(),
    IsEvenInt(k) or IsSymmetry(c),
    Cond(IsOddInt(k),
         let(kk := (k-1)/2, n := Cols(b),
             top := DirectSum(Tensor(I(kk), b), GathExtendU(n, c)),
             VStack(top,
                    Tensor(J(kk), d) * Gath(fAdd(Cols(top), kk*Cols(d), 0)))),
         let(kk := k/2,
             VStack(Tensor(I(kk), b),
                    Tensor(J(kk), d)))
         ));

Jp := x -> OS(x, -1);

listExtMat := mat -> List(mat, row -> let(n:=PositionProperty(row, x->x<>0),
    Cond(n=false, 0,
         row[n]<0, -n,
         +n)));

analyze := (stride, extmat) -> let(
    lst := listExtMat(MatSPL(extmat)),
    fsize := Length(lst)/stride,
    rec(inp := fxi_inputs(lst, fsize),
        sym := fxi_rsyms(lst, fsize)));

analyzeT := (stride, extmat, partsizes) -> let(
    fsize := Rows(extmat)/stride,
    lst := listExtMat(TransposedMat(MatSPL(L(fsize*stride, stride)*extmat))),
    part := Drop(ScanL(partsizes, (p, x) -> [Last(p)+1..Last(p)+x], [0]), 1),
    rec(out := List(part, p->lst{p})));

panalyze := function(stride, extmat)
    local r, zip, i, e, sym, n, mid, j, cur, absprev, prev, sign;
    r := analyze(stride, extmat);
    zip := Zip2(r.sym, r.inp);
    n := Length(zip);
    mid := Int((n+1)/2);
    for i in [1..Length(zip)] do
        e := zip[i];
        if e[1] = "none" then
            sym := "I";
            if i > mid then
                cur := List(e[2], AbsInt);
                for j in [1..mid] do
                   prev := r.inp[j];
                   absprev := List(r.inp[j], AbsInt);
                   if cur=Reversed(absprev) then
                       prev := Reversed(prev);
                       sign := Cond(e[2]=prev, "", e[2]=-prev, "-", "#");
                       sym := sign :: "J/" :: StringInt(j);
                   elif cur=[absprev[1]]::Reversed(Drop(absprev, 1)) then
                       prev := [prev[1]]::Reversed(Drop(prev, 1));
                       sign := Cond(e[2]=prev, "", e[2]=-prev, "-", "#");
                       sym := sign :: "J'/"::StringInt(j);
                   fi;
                od;
            fi;
        else
            sym := YellowStr(e[1]);
        fi;
        Print(sym, ": ", e[2], "\n");
    od;
end;


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

Import(fftx.platforms.cuda);
Import(fftx.codegen);
Import(simt);

Class(WarpXCUDADeviceOpts, FFTXCUDADeviceOpts, rec(
#    tags := [ASIMTKernelFlag(ASIMTGridDimY), ASIMTBlockDimX, ASIMTLoopDim() ],
#    tags := [ASIMTGridDimX, ASIMTBlockDimX, ASIMTKernelFlag(ASIMTLoopDim()) ],
#    tags := [ASIMTBlockDimZ, ASIMTKernelFlag(ASIMTBlockDimX), ASIMTKernelFlag(ASIMTBlockDimY) ],
    tags := [ASIMTKernelFlag(ASIMTGridDimY), ASIMTGridDimX, ASIMTBlockDimZ],

    operations := rec(Print := s -> Print("<FFTX WarpX CUDA Device options record>"))
));

warpXCUDADeviceOpts := function(arg) # specific to WarpX size 80 and 100... # the magic list [20, 80, 100] needs to be exposed from the API level
    local opts, rfs;

    rfs := Copy(RulesFuncSimp);
    Unbind(rfs.rules.TensorIdId);
    Unbind(rfs.rules.Const_fbase);

    opts := Copy(WarpXCUDADeviceOpts);
    opts.breakdownRules.Circulant := [Circulant_PRDFT_FDataNT];
    opts.breakdownRules.PRDFT := List([PRDFT1_Base1, PRDFT1_Base2, CopyFields(PRDFT1_CT,
            rec(allChildren := P ->Filtered(PRDFT1_CT.allChildren(P), i->When(P[1] in [20, 80, 100], Cols(i[1]) = 4, true)))),
        PRDFT_PD], _noT);
    opts.breakdownRules.IPRDFT := List([ IPRDFT1_Base1, IPRDFT1_Base2, CopyFields(IPRDFT1_CT,
            rec(allChildren := P ->Filtered(IPRDFT1_CT.allChildren(P), i->When(P[1] in [20, 80, 100], Cols(i[2]) = 4, true)))),
        IPRDFT_PD], _noT);
    opts.breakdownRules.IPRDFT2 := List([ IPRDFT2_Base1, IPRDFT2_Base2, IPRDFT2_CT ], _noT);
    opts.breakdownRules.PRDFT3 := List([ PRDFT3_Base1, PRDFT3_Base2, PRDFT3_CT ], _noT);
    opts.breakdownRules.DFT := [ DFT_Base,
        CopyFields(DFT_CT, rec(children := nt ->Filtered(DFT_CT.children(nt), i->When(nt.params[1] in [20, 80, 100], Cols(i[1]) = 4, true)))),
        DFT_PD ];

    opts.breakdownRules.TTensorInd := [TTensorInd_SIMT];
    opts.breakdownRules.TL := [L_SIMT];
    opts.breakdownRules.TIterHStack := [TIterHStack_SIMT];
    opts.breakdownRules.TIterVStack := [TIterVStack_SIMT];
    opts.breakdownRules.MDPRDFT := [MDPRDFT_3D_SIMT];
    opts.breakdownRules.IMDPRDFT := [IMDPRDFT_3D_SIMT];
    opts.breakdownRules.TRC := [TRC_SIMT];
    opts.breakdownRules.MDDFT := [ MDDFT_Base, CopyFields(MDDFT_tSPL_RowCol, rec(switch := true)), MDDFT_RowCol];
    opts.breakdownRules.TTensor := [ AxI_IxB, IxB_AxI ];
    opts.breakdownRules.TTensorI := [ IxA_SIMT,  AxI_SIMT ];
    opts.breakdownRules.TSparseMat := [CopyFields(TSparseMat_base, rec(max_rows := 1)), TSparseMat_VStack];

    opts.globalUnrolling := 12; #23; ## 12 for size 80, 23 for size 100--needs to be exposed at the API level
    opts.codegen.GathPtr := fftx.codegen.MultiPtrCodegenMixin.GathPtr;
    opts.codegen.ScatPtr := fftx.codegen.MultiPtrCodegenMixin.ScatPtr;
    opts.codegen.OO := (self, o, y, x, opts) >> skip();

    opts.sumsgen.IterHStack := MultiPtrSumsgenMixin.IterHStack;
    opts.preProcess := t -> ApplyStrategy(t,
                    [ RulesFFTXPromoteWarpX1,
                      MergedRuleSet(rfs, RulesSums, RulesFFTXPromoteWarpX2),
                      RulesFFTXPromoteNT,
                      MergedRuleSet(rfs, RulesSums, RulesFFTXPromoteWarpX3) ],
        BUA, opts);

    opts.postProcessSums := (s, opts) -> let(s1 := ApplyStrategy(s, [ MergedRuleSet(rfs, RulesSums, RulesSIMTFission) ], BUA, opts),
        FixUpCUDASigmaSPL(s1, opts));

#    opts.useDeref := false;
    # for debugging...
#    opts.generateInitFunc := false;
#    opts.codegen := DefaultCodegen;
#    opts.arrayBufModifier := "";
#    opts.arrayDataModifier := "static const";

    # lets see if that gets all the temps into device memory
    opts.max_shmem_size := 16384;
    opts.use_shmem := false;
    opts.max_threads := 1024;
    opts.fixUpZ := true;

    return opts;
end;

confWarpXCUDADevice := rec(
    defaultName := "confWarpXCUDADevice",
    defaultOpts := (arg) >> rec(useWarpXCUDADevice := true),
    confHandler := warpXCUDADeviceOpts
);

fftx.FFTXGlobals.registerConf(confWarpXCUDADevice);

ext: gi

max_stars_repo_name: franzfranchetti/spiral-package-fftx

max_issues_repo_name: franzfranchetti/spiral-package-fftx

max_forks_repo_name: franzfranchetti/spiral-package-fftx


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


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

Import(compiler, code, fpgen);

# temporary fix to declare variables used in _mm_loadl_xx/_mm_loadh_xx
CellCompileStrategyVector := Concatenation(
#    [ c -> vref.setNoScalar(c) ], # must be at the begin; otherwise subvector access breaks
    BaseCS,
    [
    BinSplit, CSE,
    # MarkDefUse, FFTWScheduleAssignments, CopyPropagate, <- NOTE: breaks
#    MarkDefUse, DFSChain, # currently not used
    CopyPropagate,
    Compile.declareVars,
    (c,opts) -> opts.vector.isa.fixProblems(c, opts),
    (c, opts) -> ESReduce(c, opts)
#    c -> DeadCodeElim(c),
#    c -> vref.resetNoScalar(c) # must be at the end right before declaring the missing vars!!
    #DeclareVars
]);

CellCompileStrategyVectorFP := (bits, fracbits) -> Concatenation(
    CellCompileStrategyVector,
    [ c -> FixedPointCode(c, bits, fracbits) ]
);


text:
# Built-in
Combinations([1, 2, 3]);
# [ [  ], [ 1 ], [ 1, 2 ], [ 1, 2, 3 ], [ 1, 3 ], [ 2 ], [ 2, 3 ], [ 3 ] ]

# Note that it handles duplicates
Combinations([1, 2, 3, 1]);
# [ [  ], [ 1 ], [ 1, 1 ], [ 1, 1, 2 ], [ 1, 1, 2, 3 ], [ 1, 1, 3 ], [ 1, 2 ], [ 1, 2, 3 ], [ 1, 3 ],
#   [ 2 ], [ 2, 3 ], [ 3 ] ]

ext: gap

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


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

# slab allocator.

_LuaInit := function(datas, slab, opts)
    local init;

    Error("haha");

end;

#
## _DefaultWrapInitCompute
#
# this is the first function which gets called on the code. It wraps
# the pure transform code in a function and also extracts any necessary
# init code.
#
# in addition, for the slab allocator, this code defines the slab struct.

Class(_SlabInitCompute, CodegenStrat, rec(
    __call__ := function(code, sums, opts)
        local x, y, io, params, datas, slab_struct,
            slab, sub, initsub, init, compute, stackvars,
            replacer,  type, length;

        # replacer function: replaces reference of a slab
        # variable in the code with a reference to the variable
        # in the slab. EG: X ==> slab->X
        replacer := (slab, code) ->
            DoForAll(slab.t.t.getVars(), e ->
                SubstTopDown(code, e, ee -> fld(e.t, slab, e.id))
            );

        type := (a) -> a;
        length := (a) -> a;

        # -------
        # here we build the slab structure
        # -------

        # get the standard input output
        [x, y] := _getXY(sums, opts);
        io := When(x=y, x, Concat(y, x));

        # any parameters in our sigma-spl expression
        params := Set(Collect(sums, param));

        # any data required by the transform
        datas := Collect(sums, FDataOfs);

        [stackvars, code] := Pull(code,

            # compare shape
            @(1, decl, e -> ForAny(e.vars, ee -> ObjId(ee.t) = TArray)),

            # subst shape: always drop array declarations, possibly drop decl.
            e -> let(nonarrayvars := Filtered(e.vars, ee -> ObjId(ee.t) <> TArray),
                When(nonarrayvars = [],
                    @(1).val.cmd,
                    decl(nonarrayvars, @(1).val.cmd)
                )
            ),

            # pull shape: pull out the array vars
            e -> Filtered(@(1).val.vars, ee -> ObjId(ee.t) = TArray)
        );

        stackvars := Flat(stackvars);

        # build the slab structure, and a variable for it
        slab_struct := T_Struct("slab_t", Concat(params, List(datas, e -> e.var), io, stackvars));
        slab := var("slab", TPtr(slab_struct));

        # names
        sub := Cond(IsBound(opts.subName), opts.subName, "compute");
        initsub := Cond(IsBound(opts.subName), Concat("init_", opts.subName), "initz");

        # generate the 'init' code
        if IsBound(opts.generateInitFunc) and opts.generateInitFunc then

            if IsBound(opts.luaInit) and opts.luaInit then
                init := _LuaInit(datas, slab, opts);
            else
                init := chain( List(datas, e -> SReduce(e.var.init, opts)) );
                replacer(slab, init);
            fi;

            init := func(TVoid, initsub, [slab], init);

        else
            init := func(TVoid, initsub, [slab], code);
        fi;

        # convert the direct array references to input, output, and data to
        # references to fields inside the slab alloced block.

        replacer(slab, code);

        compute := func(TVoid, sub, [slab], code);

        return program(
            define([slab_struct]),
            init,
            compute
        );

    end
));

Class(_SlabAlloc, CodegenStrat, rec(
    __call__ := function(code, sums, opts)
        local sub, initsub, funcs, slab_struct, slab, pslab, alloc, free, tmp;

        # names
        sub := Cond(IsBound(opts.subName), opts.subName, "compute");
        initsub := Cond(IsBound(opts.subName), Concat("init_", opts.subName), "initz");

        # we expect a program wrapper
        Constraint(ObjId(code) = program);

        # we expect a slab to be defined first.
        Constraint(ObjId(code.cmds[1]) = define);
        Constraint(ObjId(code.cmds[1].types[1]) = T_Struct);

        # extract functions for additional checks.
        funcs := Collect(code, func);

        # ensure ordering, code must already be wrapped by init/compute
        Constraint(2 = Length(Filtered(funcs, e -> e.id = sub or e.id = initsub)));

        slab_struct := code.cmds[1].types[1];
        slab := var("slab", TPtr(slab_struct));
        pslab := var("pslab", TPtr(TPtr(slab_struct)));
        tmp := var("tmp", slab_struct);
        tmp.size := 1;

        # allocation function.
        alloc := func(TVoid, "alloc", [pslab],
            allocate(deref(pslab), tmp)
        );

        free := func(TVoid, "dealloc", [slab],
            deallocate(slab, slab_struct)
        );

        Add(code.cmds, alloc);
        Add(code.cmds, free);

        return code;
    end
));

Class(_SlabTimer, CodegenStrat, rec(
    __call__ := function(code, sums, opts)
        local sub, initsub, funcs, i, numruns, slab_struct, slab, t, timer;

        # names
        sub := Cond(IsBound(opts.subName), opts.subName, "compute");
        initsub := Cond(IsBound(opts.subName), Concat("init_", opts.subName), "initz");

        # we expect a program wrapper
        Constraint(ObjId(code) = program);

        # we expect a slab to be defined first.
        Constraint(ObjId(code.cmds[1]) = define);
        Constraint(ObjId(code.cmds[1].types[1]) = T_Struct);

        # extract functions.
        funcs := Collect(code, func);

        # ensure ordering, make sure we have an init and compute
        Constraint(2 = Length(Filtered(funcs, e -> e.id = sub or e.id = initsub)));

        # ensure ordering, make sure we have the alloc/dealloc functions too!
        Constraint(2 = Length(Filtered(funcs, e -> e.id = "alloc" or e.id = "dealloc")));

        # these two are used interchangably, MUST be same type.
        i := var.fresh_t("i", TInt);
        numruns := var.fresh_t("numruns", TInt);

        slab_struct := code.cmds[1].types[1];
        slab := var("slab", TPtr(slab_struct));

        t := var("t", TPtr(TVoid));

        timer := func(TVoid, "timer", [t, numruns],
            decl([slab], chain(
                _fCall("alloc", [addrof(slab)]),
                _fCall(initsub, [slab]),
                # this is used by simics to switch from a fast functional
                # to a slow timed execute mode
                When(IsBound(opts.extraTimerCall) and opts.extraTimerCall,
                    _fCall("timer_start", [t]),
                    skip()
                ),
                When(IsBound(opts.coldcache) and opts.coldcache,
                    skip(),
                    _fCall(sub, [slab]) # warmup the cache by default.
                ),
                _fCall("timer_start", [t]),
                loop(i, numruns,
                    _fCall(sub, [slab])
                ),
                _fCall("timer_end", [t]),
                _fCall("dealloc", [slab])
            ))
        );

        Add(code.cmds, timer);

        return code;
    end,
));

Class(_SlabVerify, CodegenStrat, rec(
    __call__ := function(code, sums, opts)
        Add(code.cmds, func(TVoid, "verify", [], chain()));
        return code;
    end
));

#
## _SlabHackInit
#
# this object removes the sines/cosines from the init function. It was written to support
# the simplescalar backend. Simplescalar takes forever (and a day) to compute sin/cos, so
# this hack was necessary to improve the execution turnaround time.
#
# NOTE: It ONLY removes cos/sin calls from the init function, not any other.
#
Class(_SlabHackInit, CodegenStrat, rec(
    __call__ := function(code, sums, opts)
        local initsub, init;

        # get the init func name, and find it in the code.
        initsub := Cond(IsBound(opts.subName), Concat("init_", opts.subName), "initz");
        init := Collect(code, @1(1, func, e -> e.id = initsub));

        # replace all sin/cos calculations with constant 1.0
        SubstTopDown(init, sinpi, e -> V(1.0));
        SubstTopDown(init, cospi, e -> V(1.0));

        # replace the init function in the main code.
#        SubstTopDown(code, @(1, func, e -> e.id = initsub), e -> Error("a"));

        return code;
    end
));


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

Declare(VRC, VRCL, VRCR, VRCLR, VBlkInt);

##############################################################################
Class(VRC, RC, rec(
    toAMat := (self) >> AMatMat(RCMatCyc(MatSPL(self.child(1)))),
    new := meth(self, spl, v)
        local res;
        res := SPL(WithBases(self, rec(_children:=[spl], v:=v,
                                dimensions := spl.dimensions)));
        res.dimensions := res.dims();
        return res;
    end,
    print := (self, i, is) >> Print(self.__name__, "(\n", Blanks(i+is), self.child(1).print(i+is,is), ", ",
	#"\n", Blanks(i+is),
	self.v,
	#"\n", Blanks(i),
	")", self.printA()),

    unroll := self >> self,
    transpose := self >> VRC(self.child(1).conjTranspose(), self.v),
    vcost := self >> self.child(1).vcost(),

    from_rChildren := (self, rch) >> ObjId(self)(rch[1], self.v)
));

Class(VRCL, VRC, rec(
    toAMat := self >> AMatMat(MatSPL(Tensor(I(self.child(1).dims()[1]/self.v), L(2*self.v, 2)))) *
                      VRC(self.child(1), self.v).toAMat(),
    unroll := self >> self,
    transpose := self >> VRCR(self.child(1).conjTranspose(), self.v)
));

Class(VRCR, VRC, rec(
    toAMat := self >> VRC(self.child(1), self.v).toAMat() *
                      AMatMat(MatSPL(Tensor(I(self.child(1).dims()[2]/self.v), L(2*self.v, self.v)))),
    unroll := self >> self,
    transpose := self >> VRCL(self.child(1).conjTranspose(), self.v)
));

Class(VRCLR, VRC, rec(
    toAMat := self >> let(
    rows := Rows(self.child(1)), cols := Cols(self.child(1)), v := self.v,
    AMatMat(MatSPL(
        Tensor(I(rows/v), L(2*v, 2)) *
        RC(self.child(1)) *
        Tensor(I(cols/v), L(2*v, v))))
    ),
    unroll := self >> self,
    transpose := self >> VRCLR(self.child(1).conjTranspose(), self.v)
));

Class(vRC, RC, rec(
    print := (self, i, is) >> Print(self.name, "(", self.child(1).print(i+is,is), ")", self.printA()),
    vcost := self >> self.child(1).vcost()
));

Class(VBlkInt, BaseMat, SumsBase, rec(
    new := (self, spl, v) >> SPL(WithBases(self,
        rec(dimensions := spl.dims(), _children:=[spl], v := v))),
    #-----------------------------------------------------------------------
    dims := self >> self.dimensions,
    child := (self, i) >> self._children[i],
    rChildren := self >> [self._children[1], self.v],
    rSetChild := meth(self, n, what)
        if n=1 then self._children[1] := what;
        elif n=2 then self.v := what;
        else Error("<n> must be in [1..2]");
        fi;
    end,
    #-----------------------------------------------------------------------
    print := (self, i, is) >> Print(self.name, "(", self.child(1).print(i+is,is), ", ", self.v, ")", self.printA()),
    unroll := self >> self,
    transpose := self >> VBlkInt(self.child(1).transpose(), self.v),
    #-----------------------------------------------------------------------
    toAMat := self >> let(
        rows := Rows(self.child(1)), cols := Cols(self.child(1)), v := self.v,
        AMatMat(MatSPL(
            Tensor(I(rows/(2*v)), L(2*v, v)) *
            self.child(1) *
            Tensor(I(cols/(2*v)), L(2*v, 2))))
        )
));


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

# Declaration of a sorter of one dimension within tuples of size n
#n is the number of inputs in total. So, there are n/2 tuples
#w should be minimum 2
Class(SortVecBase, BaseMat, rec(
   abbrevs   := [()-> []],
   new       := (self) >> SPL( WithBases(self, rec()) ).setDims(),
   dims      := self >> [ 4, 4 ],
   isReal    := True,
   sums     := self >> self,
   rChildren := self >> [],
   rSetChild := rSetChildFields(),
   toAMat := self >> Error("not supported"),
   transpose := self >> self,
));

Class(SortVecConfigBase_w2, BaseMat, rec(
   abbrevs   := [(a)-> [a] , (b)-> [b]],
   new       := (self, a, b) >> SPL( WithBases(self, rec(dimensions:=[2,2], a := a, b:=b))),
   print := (self, i, is) >> Print(self.name, "(", self.a, ",", self.b, ")"),
   dims      := self >> [ 2, 2 ],
   isReal    := True,
   sums     := self >> self,
   rChildren := self >> [self.a, self.b],
   rSetChild := rSetChildFields("a","b"),
   toAMat := self >> Error("not supported"),
   transpose := self >> self,
));

# Declaration of SortConfigBase, a 2x2 Configurable sorter.
Class(SortVecConfigBase, BaseMat, rec(
   abbrevs   := [(a)-> [a]],
   new       := (self, a) >> SPL( WithBases(self, rec(dimensions:=[4,4], a := a))),
   print := (self, i, is) >> Print(self.name, "(", self.a, ")"),
   dims      := self >> [ 4, 4 ],
   isReal    := True,
   sums     := self >> self,
   rChildren := self >> [self.a],
   rSetChild := rSetChildFields("a"),
   toAMat := self >> Error("not supported"),
   transpose := self >> self,
));

Class(SortVecBase_w2, BaseMat, rec(
   abbrevs   := [(a)-> [a]],
   new       := (self, a) >> SPL( WithBases(self, rec(dimensions:=[2,2], a := a))),
   print := (self, i, is) >> Print(self.name, "(", self.a, ")"),
   dims      := self >> [ 2, 2 ],
   isReal    := True,
   sums     := self >> self,
   rChildren := self >> [self.a],
   rSetChild := rSetChildFields("a"),
   toAMat := self >> Error("not supported"),
   transpose := self >> self,
));

HDLCodegen.SortVecBase_w2 := (self, o, y, x, opts) >>
   let(
     t0 := TempVar(x.t.t),
     t1 := TempVar(x.t.t),
     t3 := TempVar(x.t.t),
     t5 := TempVar(x.t.t),
     t6 := TempVar(x.t.t),

     t20 := TempVar(x.t.t),
     t21 := TempVar(x.t.t),
     t23 := TempVar(x.t.t),
     t25 := TempVar(x.t.t),
     t26 := TempVar(x.t.t),

     chain(
        assign(t3, imod(o.a, 2)),
        regassign(t0, cond(t3, t0, nth(x,0))),
        assign(t5, cond(leq(t0, nth(x,0)), nth(x,0), t0)),
        assign(t6, cond(leq(t0, nth(x,0)), t0, nth(x,0))),
        regassign(t1, cond(t3, t5, t1)),
        assign(nth(y,0), cond(t3, t6, t1)),

        regassign(t20, cond(t3, t20, nth(x,1))),
        assign(t25, cond(leq(t0, nth(x,0)), nth(x,1), t20)),
        assign(t26, cond(leq(t0, nth(x,0)), t20, nth(x,1))),
        regassign(t21, cond(t3, t25, t21)),
        assign(nth(y,1), cond(t3, t26, t21))

     )
  );

HDLCodegen.SortVecConfigBase_w2 := (self, o, y, x, opts) >>
   let(
     t0 := TempVar(x.t.t),
     t1 := TempVar(x.t.t),
     t3 := TempVar(x.t.t),
     t5 := TempVar(x.t.t),
     t6 := TempVar(x.t.t),
     t7 := TempVar(x.t.t),
     t8 := TempVar(x.t.t),
     t9 := TempVar(x.t.t),
     t10 := TempVar(x.t.t),
     t2 := TempVar(x.t.t),

     t20 := TempVar(x.t.t),
     t21 := TempVar(x.t.t),
     t23 := TempVar(x.t.t),
     t25 := TempVar(x.t.t),
     t26 := TempVar(x.t.t),

     chain(
	assign(t3, imod(o.a, 2)),
	regassign(t0, cond(t3, t0, nth(x,0))),
	assign(t7, eq(o.b,0)),
	assign(t8, eq(o.b,1)),

	assign(t2, leq(t0, nth(x,0))),
	assign(t9, cond(t2, t0, nth(x,0))),
	assign(t10, cond(t2, nth(x,0), t0)),
	assign(t5, cond(t7, nth(x,0) , t8, t9, t10)),
	assign(t6, cond(t7, t0 , t8, t10, t9)),

	regassign(t1, cond(t3, t5, t1)),
	assign(nth(y,0), cond(t3, t6, t1)),

        regassign(t20, cond(t3, t20, nth(x,1))),
        assign(t25, cond(leq(t0, nth(x,0)), nth(x,1), t20)),
        assign(t26, cond(leq(t0, nth(x,0)), t20, nth(x,1))),
        regassign(t21, cond(t3, t25, t21)),
        assign(nth(y,1), cond(t3, t26, t21))

     )
    );

HDLCodegen.SortVecBase := (self, o, y, x, opts) >>
    chain(
        assign(nth(y,0), cond(leq(nth(x,0), nth(x,2)), nth(x,0), nth(x,2))),
        assign(nth(y,2), cond(leq(nth(x,0), nth(x,2)), nth(x,2), nth(x,0))),

        assign(nth(y,1), cond(leq(nth(x,0), nth(x,2)), nth(x,1), nth(x,3))),
        assign(nth(y,3), cond(leq(nth(x,0), nth(x,2)), nth(x,3), nth(x,1)))
    );

HDLCodegen.SortVecConfigBase := (self, o, y, x, opts) >>
    let(
	t0 := TempVar(x.t.t),
	t1 := TempVar(x.t.t),
	t2 := TempVar(x.t.t),
	t3 := TempVar(x.t.t),
	chain(
	    assign(t2, nth(x,0)),
	    assign(t3, nth(x,2)),
	    assign(t0, cond(leq(t2, t3), t2, t3)),
	    assign(t1, cond(leq(t2, t3), t3, t2)),
	    assign(nth(y,0), cond(eq(o.a,0), t2, eq(o.a,1), t1, t0)),
	    assign(nth(y,2), cond(eq(o.a,0), t3, eq(o.a,1), t0, t1)),

            assign(nth(y,1), cond(leq(nth(x,0), nth(x,2)), nth(x,1), nth(x,3))),
            assign(nth(y,3), cond(leq(nth(x,0), nth(x,2)), nth(x,3), nth(x,1)))

	)
    );

Class(SortVec, TaggedNonTerminal, rec(
    abbrevs := [
    (n)       -> Checked(IsPosIntSym(n), [_unwrap(n)]),
    ],

    hashAs := self >> ObjId(self)(self.params[1]).withTags(self.getTags()),

    dims := self >> [ self.params[1], self.params[1] ],

    terminate := self >> Error("not supported"), # we could probably support this
));

NewRulesFor(SortVec, rec(

    Sort_Stream_Vec := rec(
        info         := "Streaming sorting network",

        applicable   := nt -> Length(nt.params) = 1 and IsTwoPower(nt.params[1]),

        children := (self, nt) >> let(

	    tag_w_tmp := nt.tags[1],
	    tag_w := tag_w_tmp.bs,
	    t := Log2Int(nt.params[1]),
	    p := Ind(2^t),
	    #get_bb := w -> TTensorI(SortVecBase(), 2^(t-1), APar, APar),
	    get_bb := w -> Cond(w=2, TTensorInd(SortVecBase_w2(p), p, APar, APar),TTensorI(SortVecBase(), 2^(t-1), APar, APar)),

	    [[ TCompose(
	           [TCompose(List([1..t-1], i ->
	               TCompose([
		           #TTensorI(SortBase(), 2^(t-1), APar, APar),
		           get_bb(tag_w),
		           TCompose(List([2..(t-i+1)], j ->
		               TCompose([

			           # The one below *should* work, but it causes some rewriting problems.
			           #TPrm(Tensor(Tensor(I(2^(t-j)),(Tensor(I(2), L(2^(j-1), 2^(j-2))) * L(2^j,2))),I(2))),
				   # So, I'm going to pull the first I() out of the tensor product.  This
				   # changes it from:
				   #      TPrm(I x (IxL)*L x I)
				   # to:
				   #      I x (TPrm(IxL)*L x I)
				   TTensorI(TPrm(Tensor(Tensor(I(2), L(2^(j-1), 2^(j-2))) * L(2^j,2), I(2))), 2^(t-j), APar, APar),

				   # This one doesn't work because it has nested TPrms and it has TTensorI in a TPrm.
			           #TPrm(Tensor(TTensorI(TPrm(Tensor(I(2), L(2^(j-1), 2^(j-2))) * L(2^j,2)), 2^(t-j), APar, APar),I(2))),
			           get_bb(tag_w)
			           #TTensorI(SortBase(), 2^(t-1), APar, APar)
			       ])
		           )),

			   # TPrm(Tensor(TTensorI(TPrm(L(2^(t-i+1), 2^(t-i)) * SortIJPerm(2^(t-i+1))), 2^(i-1), APar, APar),I(2)))
			   # Like above, re-ordering these terms
			   # TPrm(Tensor(Tensor(I(2^(i-1)), L(2^(t-i+1), 2^(t-i)) * SortIJPerm(2^(t-i+1) ) ) ,I(2)))
			   TTensorI(TPrm(Tensor(L(2^(t-i+1), 2^(t-i)) * SortIJPerm(2^(t-i+1)), I(2))), 2^(i-1), APar, APar)

		       ])
	            )),
		    get_bb(tag_w)]
		    #TTensorI(SortBase(), 2^(t-1), APar, APar)]
                ).withTags(nt.getTags())
            ]]
	),

        apply        := (nt, c, cnt) -> c[1],

    ),

   Sort_Stream4_Vec := rec(
        info         := "",

	depth := 1,

        applicable   := (self, nt) >> Length(nt.params) = 1 and IsTwoPower(nt.params[1]) and
	                              let (d := self.depth, t := Log2Int(nt.params[1]), IsInt(t*t/self.depth) and
				      (IsInt(d/t) or IsInt(t/d))),

        children := (self, nt) >> let(
   	    t := Log2Int(nt.params[1]),
	    d := self.depth,
	    d1 := cond(leq(d,t), d, t).ev(),
	    d2 := cond(leq(d,t), 1, d/t).ev(),

	    k := Ind(2^(t-1)),
	    j := Ind(t),
	    l := Ind(t),
	    n := Ind(t/d2),
	    v := Ind(t/d1),

            d_tmp := cond(leq(d,t), t/d, 0).ev(),
            d2_tmp := cond(leq(t,d), d/t, 0).ev(),
	    s_tmp := ((t*t)/d),
            m2 := Ind(d_tmp),
	    s2 := Ind(s_tmp),
	    v2 := Ind(d2_tmp),
	    n2 := Ind(d),
	    n3 := Ind(d),
	    l2 := Ind(t),

	    tag_w_tmp := nt.tags[1],
	    tag_w := tag_w_tmp.bs,
	    p := Ind(2^t),

	    c1 := (lp, jp) >> lt((t-1), (lp+jp)),
	    z := (lp, jp) >> (t-1)-(lp+jp),
	    z_w1 := (lp, jp) >> (t-1)-(lp+jp)+1,
	    c2 := (lp, jp) >> logic_and(eq(bit_sel(k, z(lp, jp)), 1), neq(lp, 0)),
	    c2_w1 := (lp, jp) >> logic_and(eq(bit_sel(p, z_w1(lp, jp)), 1), neq(lp, 0)),

	    access_f := (lp, jp) >> cond(c1(lp, jp), 0, c2(lp, jp), 1, 2),
	    access_f_w1 := (lp, jp) >> cond(c1(lp, jp), 0, c2_w1(lp, jp), 1, 2),

	    get_bb := (lp,jp) -> Cond(tag_w=2, TTensorInd(SortVecConfigBase_w2(p,access_f_w1(lp,jp)), p, APar, APar),TTensorInd(SortVecConfigBase(access_f(lp, jp)), k, APar, APar)),

	    stage := (lp, jp) >> TCompose([
		       get_bb(lp, jp),
		       #TTensorInd(SortConfigBase(access_f(lp, jp)), k, APar, APar),
		       TPrm(Tensor(L(2^t, 2^(t-1)),I(2)))
	           ]),

	    full_stage := np_1 >> TCompose(List([0..t-1], m_1 -> TCompose(List([0..t-1], j_1 -> stage(m_1, j_1))))),
	    full_stage1 := np >> TCompose(List([0..d2-1], m -> TCompose(List([0..t-1], j -> stage(d2*np+m, j))))),
	    full_stage2 := vp >> TCompose(List([0..d1-1], s -> stage(l, vp+s))),
            full_stage1b := np2 >> TICompose(m2,d_tmp, TICompose(j,t, stage((t/d)*np2+m2, j))),

            # Old: problem is that it's assuming the l2 above, which is an unassigned iterator.
	    # There is also a problem with the vp2+s2 parameter: you need to multiply vp2 by the number of iterations.
            # full_stage2b := (vp2) >> TICompose(s2,s_tmp, stage(l2, vp2+s2)),
            full_stage2b := (vp2, l3) >> TICompose(s2,s_tmp, stage(l3, vp2*s_tmp+s2)),

  	    [[ Cond(d=t*t, full_stage(0).withTags(nt.getTags()),
		    d<t, TCompose(List([0..d-1], n2 -> full_stage1b(n2))).withTags(nt.getTags()),
                    d=t, TCompose(List([0..d-1], n3 -> full_stage1b(n3))).withTags(nt.getTags()),
                    d>t, TCompose(List([0..t-1], l3 -> TCompose(List([0..d2_tmp-1], v2 -> full_stage2b(v2, l3))))).withTags(nt.getTags())) #the problem seems to be the outer most TCompose works if it was TICompose
		    #d>t, TICompose(n, t/d2, full_stage1(n)).withTags(nt.getTags())) #old one but will leave it as new one does not work yet
	        ]]
	),

#

        apply        := (nt, c, cnt) -> c[1],

    ),

));


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

Class(Function, rec(
   isFunction := true,

   range := self >> self._range,
   tolist := self >> self.lambda().tolist(),

   lambda := self >> Error("Not implemented"),

   free   := self >> self.lambda().free(),
   domain := self >> self.lambda().domain(),
   at := (self, pos) >> self.lambda().at(pos),

   from_rChildren := (self, rch) >> ApplyFunc(ObjId(self), rch),

   _domain := false,
   _range := false,
   setRange := meth(self, range)
       self._range := range;
       return self;
   end,
   setDomain := meth(self, domain)
       self._domain := spiral.code.RulesStrengthReduce(domain);
       return self;
   end,

   # NOTE: below is a hack that already causes problems if we compute w/ symbolic higher order functions
   # 2 methods below allow to use constant functions (ie domain=1) as expressions

   # NOTE: don't remember where the hack above is used, so I removed it, now need to test what happens
   eval := self >> self, #When(self.domain() = 1 and self.rank() = 0, self.at(0), self),
   ev := self >> self,   #Checked(self.domain() = 1, self.at(0).ev()),
));

IsFunction := x -> IsRec(x) and IsBound(x.isFunction) and x.isFunction;

Declare(Lambda, LambdaOps, FList, FData);

#F Lambda(<vars>, <expr>)  - symbolic representation of a function mapping vars -> expr
#F
#F <vars> must be a list of variables or a single variable
#F <expr> is a symbolic expression using variables from <vars>
#F
Class(Lambda, Function, rec(
   rChildren := self >> [self.vars, self.expr, self._domain, self._range],
   rSetChild := rSetChildFields("vars", "expr", "_domain", "_range"),
   isLambda := true,
   lambda := self >> self,

   from_rChildren := (self, rch) >> ObjId(self)(rch[1], rch[2]).setDomain(rch[3]).setRange(rch[4]),

   domain := self >> When(self._domain <> false,
       self._domain,
       let(v := Last(self.vars), # ignore rank-associated vars
           When(IsBound(v.range) and not IsMeth(v.range), v.range, v.t))),

   # FF: need to be able to override range
   range := self >> When(self._range <> false, self._range, self.expr.t),

   tolist := self >> Checked(Length(self.vars)=1, let(dom := _unwrapV(self.domain()),
       Cond(IsType(dom), Error("Can't convert <self> to a list, ",
	       "because its .domain() is not a bounded interval"),
	    IsSymbolic(dom), Error("Can't convert <self> to a list, ",
	       "because its .domain() interval size is not known"),
	    List([0 .. dom-1], r -> self.at(r).eval())))),

   free := self >> Difference(self.expr.free(), self.vars),

   # FF: need to be able to override range
   print := self >> Print(self.name, "(", self.vars, ", ", self.expr, ")",
        When(self._range <> false, Print(".setRange(", self._range, ")"), ""),
        When(self._domain <> false, Print(".setDomain(", self._domain, ")"), "")
   ),

   # at(args), strict version, must have Length(args) = Length(vars)
   at := meth(arg)
       local self, args, i, bind;
       self := arg[1];
       args := Drop(arg, 1);
       # support variable-length call through lists; required by OL
       if IsList(args[1]) then args := args[1]; fi;

       Constraint(Length(args) = Length(self.vars));
       bind := tab();
       for i in [1..Length(self.vars)] do
           if IsBound(args[i]) then
               bind.(self.vars[i].id) := toExpArg(args[i]);
           fi;
       od;
       return SubstVars(Copy(self.expr), bind);
    end,

    # relaxed_at(lvars, v), relaxed version of at() for multi-ranks, allows lvars to be of any length
    # equivalent to at(Concatenation(lvars, [v])), if Length(lvars)=Length(self.vars)-1
    # v is the "standard" variable, everything in lvars is a loop index
    relaxed_at := meth(self, lvars, v)
       local i, bind, la, lv, args, vars;
       args := Concatenation(lvars, [v]);

       # for multi-rank functions loop indices can stay unbound
       vars := self.vars;
       [la, lv] := [Length(args), Length(vars)];
       if   la > lv then args := Drop(args, la - lv);
       elif lv > la then vars := Drop(vars, lv - la);
       fi;

       bind := tab();
       for i in [1..Length(vars)] do
           bind.(vars[i].id) := toExpArg(args[i]);
       od;
       return SubstVars(Copy(self.expr), bind);
    end,

#    computeType := self >> ApplyFunc(TFunc, Concatenation(
#           List(self.vars, x ->
#               When(IsBound(x.range) and not IsMeth(x.range), x.range, x.t)),
#           [self.expr.t])),

    computeType := self >> ApplyFunc(TFunc,
           List(self.vars, x ->
                Cond(IsBound(x.range) and not IsMeth(x.range), x.range,
                     Collect(self.expr, x)=[], TDummy,
                     x.t)
           ) :: [self.expr.t]),

    normalize := self >> let(
	t := self.computeType(),
	first_non_dummy := PositionProperty(DropLast(t.params,2), x->x<>TDummy),
	k := Cond(first_non_dummy=false, Length(self.vars), first_non_dummy),
	Lambda(self.vars{ [k .. Length(self.vars)]},
	       self.expr)),

    forCmp := meth(self)
        local map, i, t;
        map := tab();
        for i in [1..Length(self.vars)] do
            t := self.vars[i].t;
            map.(self.vars[i].id) := param(When(IsType(t), t, InferType(t)), "v"::StringInt(i));
        od;
        return SubstParamsCustom([self.vars, self.expr], map, [var]);
    end,

    __call__ := meth(self, vars, expr)
        local usage, res, rank, loopvars;
        usage := Concatenation("Usage: Lambda(<argList>, <expr>)\n",
                            "Usage: Lambda(<arg>, <expr>)");
            if not IsList(vars) then vars := [vars]; fi;
        if not ForAll(vars, x -> IsVar(x))
            then Error("all elements in <argList> must be variables"); fi;
        if Length(vars) <> Length(Set(vars))
            then Error("<argList> should not contain any duplicates"); fi;

        expr := toExpArg(expr);
            rank := _rank(expr);
            # if <expr> is ranked, add extra arguments to lambda
            if rank > 0 and Length(vars)=1 then
                loopvars := List([1..rank], i->var.fresh_t("q", TInt));
                vars := Concatenation(loopvars, vars);
                expr := _downRankFull(expr, loopvars);
            fi;
        res := WithBases(self, rec(expr := expr, vars := vars, operations := LambdaOps));
        res.t := res.computeType();
	# NOTE: the below line should not be necessary, not sure what the purpose was
        # if IsExp(Last(res.vars).range) then res:=res.setDomain(Last(res.vars).range); fi;
        return res;
    end,

    downRank := meth(self, loopid, ind)
        local v, vars, exp;
        if Length(self.vars) <= loopid then
            return self;
        else
            v    := self.vars[ Length(self.vars) - loopid ];
            vars := ListWithout(self.vars, Length(self.vars)-loopid);
            exp  := SubstVars(Copy(self.expr), tab((v.id) := ind));
            return ObjId(self)(vars, exp);
        fi;
    end,
));

IsLambda := obj -> IsRec(obj) and IsBound(obj.isLambda) and obj.isLambda = true;

CompatLambdas := (x, y) -> Length(x.vars) = Length(y.vars);

LambdaCompose := function(l1, l2)
    if ObjId(l1)<>Lambda then l1 := l1.lambda(); fi;
    if ObjId(l2)<>Lambda then l2 := l2.lambda(); fi;
    Constraint(Length(l1.vars)=1); Constraint(Length(l2.vars)=1);
    return let(v:=l2.vars[1].clone(), Lambda(v, l1.at(l2.at(v))));
end;

Class(LambdaOps, rec(
    Print := x->x.print(),
    \+ := (l1, l2) -> let(ll1:=l1.lambda(), ll2:=l2.lambda(), v:=l1.vars[1].clone(),
    Lambda(v, ll1.at(v) + ll2.at(v))),
    \- := (l1, l2) -> let(ll1:=l1.lambda(), ll2:=l2.lambda(), v:=l1.vars[1].clone(),
    Lambda(v, ll1.at(v) - ll2.at(v))),
    \* := (l1, l2) -> let(ll1:=l1.lambda(), ll2:=l2.lambda(), v:=l1.vars[1].clone(),
    Lambda(v, ll1.at(v) * ll2.at(v))),
    \/ := (l1, l2) -> let(ll1:=l1.lambda(), ll2:=l2.lambda(), v:=l1.vars[1].clone(),
    Lambda(v, ll1.at(v) / ll2.at(v))),
    \^ := (l1, l2) -> let(ll1:=l1.lambda(), ll2:=l2.lambda(), v:=l1.vars[1].clone(),
    Lambda(v, ll1.at(v) ^ ll2.at(v))),
    \= := (l1, l2) -> IsLambda(l1) and IsLambda(l2) and l1.lambda().forCmp() = l2.lambda().forCmp(),
    \< := (l1, l2) -> When( IsLambda(l1) and IsLambda(l2),
                              l1.lambda().forCmp() < l2.lambda().forCmp(),
                              BagAddr(l1) < BagAddr(l2) ),
));

#F FDataOfs(<datavar>, <len>, <ofs>)
#F
Class(FDataOfs, Function, rec(
    __call__ := (self, datavar, len, ofs) >> WithBases(self, rec(
    var := datavar,
    operations := PrintOps,
    ofs := toExpArg(ofs),
    len := Checked(IsPosIntSym(len), len)
    )),

# <-Daniele's changes
#    rChildren := self >> [ self.var, self.len, self.ofs],
#    rSetChild := rSetChildFields("var", "len", "ofs"),

    rChildren := self >> [ self.var, self.len, self.ofs, self._domain, self._range],
    rSetChild := rSetChildFields("var", "len", "ofs", "_domain", "_range"),
    from_rChildren := (self, rch) >> ObjId(self)(rch[1], rch[2], rch[3]).setDomain(rch[4]).setRange(rch[5]),
# ->

    domain := self >> self.len,
    print := self >> Print(self.name,"(",self.var,", ",self.len,", ",self.ofs,")"),

    at := (self, n) >> When(IsInt(n) and IsValue(self.ofs) and IsBound(self.var.value),
        self.var.value.v[n + self.ofs.v + 1],
        nth(self.var, n + self.ofs)),

    tolist := self >> List([0..EvalScalar(self.len-1)], i -> nth(self.var, self.ofs+i)),
    lambda := self >> let(x := Ind(self.domain()), Lambda(x, nth(self.var, self.ofs+x))),

    domain := self >> self.len,
    range := self >> When(self._range=false, self.var.t.t, self._range),
    inline := true,
    free := self >> self.ofs.free()
));

#F FData(<datavar>) -- symbolic function i -> datavar[i],
#F
#F domain = datavar.range
#F range = datavar.t
#F
#F
Class(FData, Function, rec(
   __call__ := arg >> let(
       self := arg[1],
       _val := Cond(Length(arg)=2, When(IsList(arg[2]), arg[2], [arg[2]]),
                               Drop(arg, 1)),
       val := When(Length(_val)=1 and IsLoc(_val[1]), _val[1], V(_val)),
       datavar := Cond(IsLoc(val), val,
                   Dat(val.t).setValue(val)),
       WithBases(self, rec(var := datavar, operations := PrintOps))),

   print := self >> Print(self.name, "(", self.var, ")"),
   rChildren := self >> [ self.var ],
   rSetChild := rSetChildFields("var"),

   at := (self, n) >> When(IsInt(n), self.var.value.v[n+1], self.lambda().at(n)),
   tolist := self >> When(IsBound(self.var.value), self.var.value.v, self.lambda().tolist()),
   lambda := self >> let(x := Ind(self.domain()), Lambda(x, nth(self.var, x))),

   domain := self >> self.var.t.size,
   range := self >> self.var.t.t,

   inline := true,
   free := self >> Set([]),
   part := (self, len, ofs) >> FDataOfs(self.var, len, ofs),
));

Class(RCData, Function, rec(
   __call__ := (self, func) >> WithBases(self, rec(
       func := Checked(IsFunction(func) or IsFuncExp(func), func),
       operations := PrintOps)),

   print := self >> Print(self.name, "(", self.func, ")"),
   rChildren := self >> [ self.func ],
   rSetChild := rSetChildFields("func"),

   at := (self, n) >> cond(n mod 2, im(self.func.at(idiv(n,2))), re(self.func.at(idiv(n,2)))),
   tolist := self >> ConcatList(self.func.tolist(), x->[re(x), im(x)]),

   domain := self >> let(d:=self.func.domain(),
       Cond(not IsType(d), d * 2,
            d=TInt,        TInt,
            Error("RCData: Invalid self.func.domain()"))),

   range := self >> self.func.range().realType(),
   # this .lambda corrently handles the case of multirank self.func
   lambda := self >> let(
       x := Ind(self.domain()),
       f := self.func.lambda(),
       jv := DropLast(f.vars, 1),
       Lambda(Concatenation(jv, [x]),
              cond(neq(imod(x, 2),0), im(f.relaxed_at(jv, idiv(x,2))), re(f.relaxed_at(jv, idiv(x,2)))))),

   inline := true,
   free := self >> self.func.free()
));

#F
#F CRData is the opposite of RCData: takes real function which assumed to be interleaved complex data
#F   and packs it into complex function:
#F   [re0, im0, re1, im1, ... ] -> [ (re0, im0), (re1, im1), ... ]
#F

Class(CRData, RCData, rec(

   at     := (self, n) >> cxpack(self.func.at(2*n), self.func.at(2*n+1)),
   tolist := self >> let( l := self.func.tolist(), List([1..Length(l)/2], i->cxpack(l[2*i-1], l[2*i]))),

   lambda := self >> let(
       x  := Ind(self.domain()),
       f  := self.func.lambda(),
       jv := DropLast(f.vars, 1),
       Lambda(Concatenation(jv, [x]), cxpack(f.relaxed_at(jv, 2*x), f.relaxed_at(jv, 2*x+1)))),

   domain := self >> let(d:=self.func.domain(),
       Cond(not IsType(d), div(d, 2),
            d=TInt,        TInt,
            Error("CRData: Invalid self.func.domain()"))),

   range := self >> self.func.range().complexType(),
));

Declare(FConj, FInv, FRConj);

Class(FConj, Function, rec(
   __call__ := (self, func) >> When(ObjId(func)=FConj, func.func, WithBases(self, rec(
       func := Checked(IsFunction(func) or IsFuncExp(func), func),
       operations := PrintOps))),

   print := self >> Print(self.name, "(", self.func, ")"),
   rChildren := self >> [ self.func ],
   rSetChild := rSetChildFields("func"),

   at := (self, n) >> conj(self.func.at(n)),
   tolist := self >> List(self.func.tolist(), conj),

   domain := self >> self.func.domain(),
   range := self >> self.func.range(),
   lambda := self >> let(x := Ind(self.domain()), Lambda(x, self.at(x))),

   inline := true,
   free := self >> self.func.free()
));

Class(FInv, Function, rec(
   __call__ := (self, func) >> When(ObjId(func)=FInv, func.func, WithBases(self, rec(
       func := Checked(IsFunction(func) or IsFuncExp(func), func),
       operations := PrintOps))),

   print := self >> Print(self.name, "(", self.func, ")"),
   rChildren := self >> [ self.func ],
   rSetChild := rSetChildFields("func"),

   at := (self, n) >> 1/(self.func.at(n)),
   tolist := self >> List(self.func.tolist(), x -> 1/x),

   domain := self >> self.func.domain(),
   range := self >> self.func.range(),
   lambda := self >> let(x := Ind(self.domain()), Lambda(x, self.at(x))),

   inline := true,
   free := self >> self.func.free()
));

# conjugation of N complex entries stored as 2*N reals r,i,r,i,...
Class(FRConj, Function, rec(
   __call__ := (self, func) >> When(ObjId(func)=FRConj, func.func, WithBases(self, rec(
       func := Checked(IsFunction(func) or IsFuncExp(func), func),
       operations := PrintOps))),

   print := self >> Print(self.name, "(", self.func, ")"),
   rChildren := self >> [ self.func ],
   rSetChild := rSetChildFields("func"),

   at := (self, n) >> cond(n mod 2, neg(self.func.at(n)), self.func.at(n)),
   tolist := self >> let(ll := self.func.tolist(),
       List([0..Length(ll)-1], i -> When(i mod 2 = 1, neg(ll[1+i]), ll[1+i]))),

   domain := self >> self.func.domain(),
   range := self >> self.func.range(),
   lambda := self >> let(x := Ind(self.domain()), Lambda(x, self.at(x))),

   inline := true,
   free := self >> self.func.free()
));

#F FList(<type>, <list>) -- symbolic function i -> list[i],
#F
#F Example: FList(TInt, [1,2,3])
#F See also: FData
Class(FList, Function, rec(
   __call__ := (self, t, list) >> Checked(IsList(list), WithBases(self, rec(
       list := list,
           t := t, operations := PrintOps))),

   print := self >> Print(self.name, "(", self.t, ", ", self.list, ")"),
   rChildren := self >> [self.t, self.list],
   rSetChild := rSetChildFields("t", "list"),

   domain := self >> Length(self.list),
   range  := self >> self.t,

   at := (self, n) >> When(IsInt(n), self.list[n+1], self.lambda().at(n)),

   _typ := self >> When(IsInt(self.t), TInt, self.t),
   tolist := self >> List(self.list, e->When(not IsSymbolic(e), Value(self._typ(), e), e)),
   tolistval := self >> Value(TArray(self._typ(), Length(self.list)), self.tolist()),
   lambda := self >> let(x := Ind(self.domain()), Lambda(x, nth(self.tolistval(), x))),

   todata := self >> FData(self.list),
   inline := false,
   free := self >> Union(List(self.list, FreeVars))
));

FF := x -> Cond(IsRec(x) and IsBound(x.isFunction), x,
                IsList(x), FList(x),
                Error("<x> must be a function or a list"));

#F FPerm(<perm>) -- function representing permutation
#F
#F Example: FPerm(Perm((1,2,3),4))
Class(FPerm, Function, rec(
   __call__ := (self, perm) >> WithBases(self, rec(
           perm := perm,
           operations := PrintOps)),

   print := self >> Print(self.name, "(", self.perm, ")"),
   rChildren := self >> [self.perm],
   rSetChild := rSetChildFields("perm"),

   domain := self >> self.perm.dims()[1],
   range  := self >> self.perm.dims()[1],

   at := (self, n) >> When(IsInt(n), self.list[n+1], self.lambda().at(n)),

   _typ := self >> TInt,
   tolist := self >> List(self.list, e->When(not IsSymbolic(e), Value(self._typ(), e), e)),
   tolistval := self >> Value(TArray(self._typ(), Length(self.list)), self.tolist()),
   lambda := self >> self.perm.sums().func.lambda(),

   inline := true,
   free := self >> FreeVars
));


text:
--
-- The Java KeyWord Lexer
--
%Options fp=CncKWLexer
%options package=CnCParser
%options template=KeywordTemplateF.gi

%Notice
/.
//
// This file is part of the CNC-C implementation and
// distributed under the Modified BSD License.
// See LICENSE for details.
//
// I AM A GENERATED FILE. PLEASE DO NOT CHANGE ME!!!
//
./
%End

%Include
    KWLexerFoldedCaseMapF.gi
%End

%Export

    T_ENV
    T_UNSIGNED
    T_STRUCT

%End

%Terminals
    a    b    c    d    e    f    g    h    i    j    k    l    m
    n    o    p    q    r    s    t    u    v    w    x    y    z
%End

%Start
    KeyWord
%End

%Rules

    KeyWord ::= e n v /.$setResult($_T_ENV);./
    KeyWord ::= s t r u c t /.$setResult($_T_STRUCT);./
    KeyWord ::= u n s i g n e d /.$setResult($_T_UNSIGNED);./

%End

max_stars_repo_name: pelmers/cnc-ocr

max_issues_repo_name: pelmers/cnc-ocr

max_forks_repo_name: pelmers/cnc-ocr


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

#Don't forget to Import(platforms.intel);

#_svctSizes := (N, p, v) -> Filtered([1..N], i->ForAll(Factors(i), j->j<=p) and IsInt(i/v^2));
#_defaultSizes := (func, default) -> ((arg) -> When(Length(arg) = 1, func(When(IsList(arg[1]), arg[1], [arg[1]])), func(default)));

_cores := 4;
_vlen := 4;
_maxN := 16;
_factors :=  [1, 2, 3, 5];
_parSizes := (N, p, v) -> List(Set(Filtered(Flat(Map(Cartesian(v^2*p^2 * List([1.._maxN], i->2^i), _factors), Product)), i->i<=2^N and i>v^2*p^2)));
_large := List([10..20], i->2^i);

_smpOpts := (intlCmplx, buf, omp) -> rec(use_functions := false, use_openmp := omp, use_buffering := buf, interleavedComplex := intlCmplx,
    simd_opts:=rec(svct:=true, splitL:=false, oddSizes:=false));

benchCore2 := function()
    if LocalConfig.cpuinfo.SIMD().hasSSE3() and LocalConfig.cpuinfo.cores >= 2 then
        return rec(
            OpenMP_2x64f := rec(
                split := rec(
                    medium := _defaultSizes(s->spiral.libgen.doParSimdDft(SSE_2x64f, LocalConfig.cpuinfo.cores, s, _smpOpts(false, false, true)), _parSizes(_maxN, _cores, _vlen)),
                    large := _defaultSizes(s->spiral.libgen.doParSimdDft(SSE_2x64f, LocalConfig.cpuinfo.cores, s, _smpOpts(false, true, true)), _large)
                ),
                interleaved := rec(
                    medium := _defaultSizes(s->spiral.libgen.doParSimdDft(SSE_2x64f, LocalConfig.cpuinfo.cores, s, _smpOpts(true, false, true)), _parSizes(_maxN, _cores, _vlen)),
                    large := _defaultSizes(s->spiral.libgen.doParSimdDft(SSE_2x64f, 1, s, _smpOpts(true, true, true)), _large)
                )
            ),
            OpenMP_4x32f := rec(
                split := rec(
                    medium := _defaultSizes(s->spiral.libgen.doParSimdDft(SSE_4x32f, LocalConfig.cpuinfo.cores, s, _smpOpts(false, false, true)), _parSizes(_maxN, _cores, _vlen)),
                    large := _defaultSizes(s->spiral.libgen.doParSimdDft(SSE_4x32f, LocalConfig.cpuinfo.cores, s, _smpOpts(false, true, true)), _large),

                ),
                interleaved := rec(
                    medium := _defaultSizes(s->spiral.libgen.doParSimdDft(SSE_4x32f, LocalConfig.cpuinfo.cores, s, _smpOpts(true, false, true)), _parSizes(_maxN, _cores, _vlen)),
                    large := _defaultSizes(s->spiral.libgen.doParSimdDft(SSE_4x32f, LocalConfig.cpuinfo.cores, s, _smpOpts(true, true, true)), _large)
                )
            ),
            threads_2x64f := rec(
                split := rec(
                    medium := _defaultSizes(s->spiral.libgen.doParSimdDft(SSE_2x64f, LocalConfig.cpuinfo.cores, s, _smpOpts(false, false, false)), _parSizes(_maxN, _cores, _vlen)),
                    large := _defaultSizes(s->spiral.libgen.doParSimdDft(SSE_2x64f, LocalConfig.cpuinfo.cores, s, _smpOpts(false, true, false)), _large)

                ),
                interleaved := rec(
                    medium := _defaultSizes(s->spiral.libgen.doParSimdDft(SSE_2x64f, LocalConfig.cpuinfo.cores, s, _smpOpts(true, false, false)), _parSizes(_maxN, _cores, _vlen)),
                    large := _defaultSizes(s->spiral.libgen.doParSimdDft(SSE_2x64f, 1, s, _smpOpts(true, true, false)), _large)
                )
            ),
            threads_4x32f := rec(
                split := rec(
                    medium := _defaultSizes(s->spiral.libgen.doParSimdDft(SSE_4x32f, LocalConfig.cpuinfo.cores, s, _smpOpts(false, false, false)), _parSizes(_maxN, _cores, _vlen)),
                    large := _defaultSizes(s->spiral.libgen.doParSimdDft(SSE_4x32f, LocalConfig.cpuinfo.cores, s, _smpOpts(false, true, false)), _large)

                ),
                interleaved := rec(
                    medium := _defaultSizes(s->spiral.libgen.doParSimdDft(SSE_4x32f, LocalConfig.cpuinfo.cores, s, _smpOpts(true, false, false)), _parSizes(_maxN, _cores, _vlen)),
                    large := _defaultSizes(s->spiral.libgen.doParSimdDft(SSE_4x32f, LocalConfig.cpuinfo.cores, s, _smpOpts(true, true, false)), _large)
                )
            )
        );
    else
        return false;
    fi;
end;


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

ImportAll(spiral);

alphai := (n,i,a) ->  Cond(
    i mod 2 = 0, (  a + Int(i/2)) / n,
    i mod 2 = 1, (1-a + Int(i/2)) / n);

#alphai := (n,i,alpha) -> fr(n,i,2*alpha).ev()/2;
cat := Concatenation;
cos := x->CosPi(2*x);
sin := x->SinPi(2*x);

ls := (n,f) -> List([0..n-1], f);

MatRDFT1 := n -> cat(
     [ List([0..n-1], x->1), List([0..n-1], x->(-1)^x) ],
     cat(List([1..n/2-1], i->
                    [ ls(n, j -> cos(j*i/n)),
                      ls(n, j -> sin(j*i/n)) ])));

MatRDFT3 := n ->
     cat(List([0..n/2-1], i->
                    [ ls(n, j -> cos(j*(i+1/2)/n)),
                      ls(n, j -> sin(j*(i+1/2)/n)) ]));

MatRDFT := (N,a) -> let(n:=N/2,
     cat(ls(n, i-> let(aa := alphai(n,i,a), # 2 since we use 2pi*a
                    [ ls(2*n, j -> cos(j*aa)),
                      ls(2*n, j -> sin(j*aa)) ]))));

MatrDFT := (N,a) -> let(n:=N/2,
     cat(ls(n, i-> let(aa := alphai(n,i,a),
                    [ cat(ls(n, j->cos(j*aa)), (-1)^i*ls(n, j -> -sin(j*aa))),
                      cat(ls(n, j->sin(j*aa)), (-1)^i*ls(n, j ->  cos(j*aa))) ]))));

MatBRDFT := (N,a) -> MatSPL(BlockMat(let(n:=N/2,
     ls(n, i->let(phi := [[0,-1],[1,2*cos(alphai(n,i,a))]],
	 ls(n, j->Mat(phi^(2*j))))))));


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

NewRulesFor(Circulant, rec(
    Circulant_PRDFT_FDataNT := rec(
    	forTransposition := false,
    	switch           := true,

    	applicable       := nt -> ObjId(nt.params[2]) = FDataNT and ObjId(nt.params[2].nt) = IPRDFT and
                                  nt.params[2].nt.params = [nt.params[1], 1] and nt.params[1] = -nt.params[3],

    	children      := nt -> let(n := nt.params[1],
    	    [[ spiral.transforms.realdft.IPRDFT(n, 1), spiral.transforms.realdft.PRDFT(n, -1) ]]),

        apply := (nt, C, cnt) -> C[1] * RCDiag(FDataOfs(nt.params[2].var, Cols(nt.params[2].nt), 0)) * C[2]
    )
));

#_isFractionP := (n, p) -> IntDouble(p*n)/p = n;
#
#NewRulesFor(TResample, rec(
#    TResample_TGath := rec(
#    	forTransposition := false,
#    	switch           := true,
#
#    	applicable       := nt -> ForAny(Zip2(nt.params[1], nt.params[2]), l-> l[1]<l[2]),
#
#    	children      := nt -> [[TResample(nt.params[1], nt.params[1], nt.params[3]), TGath(_toBox(nt.params[2], nt.params[1]))]],
#
#        apply := (nt, C, cnt) -> C[1] * C[2]
#    ),
#
#    TResample_TScath := rec(
#    	forTransposition := false,
#    	switch           := true,
#
#    	applicable       := nt -> ForAny(Zip2(nt.params[1], nt.params[2]), l-> l[1]>l[2]),
#
#    	children      := nt -> [[TScat(_toBox(nt.params[1], nt.params[2])), TResample(nt.params[2], nt.params[2], nt.params[3])]],
#
#        apply := (nt, C, cnt) -> C[1] * C[2]
#    ),
#
#    TResample_RC := rec(
#    	forTransposition := false,
#    	switch           := true,
#
#    	applicable       := nt -> Length(nt.params[1]) > 1 and Length(nt.params[2]) > 1 and Length(nt.params[3]) > 1
#                                  and nt.params[1][1] = nt.params[2][1],
#
#    	children      := nt -> [ List(Zip2(Zip2(nt.params[1], nt.params[2]), nt.params[3]), e-> ApplyFunc(TResample, List(Flat(e), i->[i]))) ],
#
#        apply := (nt, C, cnt) -> let(_C:= Reversed(C), n := Length(nt.params), Compose(
#            List([1..n], i->let(lens := List(_C, j->j.dims()[1]), Tensor(I(Product(lens{[1..i-1]})), _C[i], I(Product(lens{[i+1..n]})))))))
#    ),
#
#    TResample_I := rec(
#    	forTransposition := false,
#    	switch           := true,
#
#    	applicable       := nt -> Length(nt.params[1]) = 1 and Length(nt.params[2]) = 1 and Length(nt.params[3]) = 1
#                                  and nt.params[1][1] = nt.params[2][1]  and nt.params[3][1] = 0,
#
#    	children      := nt -> [[]],
#
#        apply := (nt, C, cnt) -> I(nt.dims()[1])
#    ),
#
#    TResample_PRDFT_frac := rec(
#    	forTransposition := false,
#        denominator      := 2,
#    	switch           := true,
#
#    	applicable       := (self,nt) >> Length(nt.params[1]) = 1 and Length(nt.params[2]) = 1 and Length(nt.params[3]) = 1
#                                  and nt.params[1][1] = nt.params[2][1] and nt.params[3][1] <> 0 and _isFractionP(nt.params[3][1], self.denominator),
#
#    	children      := nt -> [[ IPRDFT1(nt.params[1][1], -1), PRDFT1(nt.params[1][1], 1) ]],
#
#        apply := (self, nt, C, cnt) >> let(n := nt.params[1][1], nf := n/2+1, shft := IntDouble(nt.params[3][1] * self.denominator),
#            df := diagMul(fConst(nf, 1/n), fCompose(dOmega(self.denominator*n, shft), dLin(nf, 1, 0, TInt))),
#            dm :=RCDiag(RCData(df)),
#            C[1] *dm * C[2])
#    )
#));
#


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

################################################################################
##
#F DG_S_and_T(<group>,<bool>) Computes Simples, S, T, FPdim(Xi), FPdimC for D(G)
##
InstallGlobalFunction(DG_S_and_T, function(G,display)
  local g, elmsG, R, r, tbl, chi, Simples, x, y, S, s, T, rank, x_idx, y_idx, FPdim, FPdimC,l,pseudounitary_data,sorted_data;
  elmsG := Elements(G);
  if display then
    Print("\nG = ", StructureDescription(G), "              ", IdGroup(G), "\n");
  fi;

  # compute the simples
  R := List(ConjugacyClasses(G), Representative);
  Simples := [];
  for r in R do
    tbl := CharacterTable(Centralizer(G, r));
    for chi in Irr(tbl) do
      Add(Simples, [r, chi, ConjugacyClasses(tbl)]);
    od;
  od;
  rank:=Size(Simples);

  if display then
    Print("Rank is\n\n");
    Display(rank);
    Print("Computing S matrix\n");
  fi;

  # compute the S-matrix
  # since s is symmetric we can half the size of the inner loop
  #S := NullMat(rank,rank);
  S:=[];
#  for x_idx in [1..rank] do
#    x:=Simples[x_idx];
  for x in Simples do
    l:=[];
    #for y_idx in [x_idx..rank] do
#    for y_idx in [1..rank] do
#      y:=Simples[y_idx];
    for y in Simples do
      s := 0;
      for g in elmsG do
        if x[1]*g*y[1]*Inverse(g) = g*y[1]*Inverse(g)*x[1] then
         s := s + ComplexConjugate(x[2][Position(x[3],
                  ConjugacyClass(Centralizer(G, x[1]), g*y[1]*Inverse(g)))])
                * ComplexConjugate(y[2][Position(y[3],
                  ConjugacyClass(Centralizer(G, y[1]), Inverse(g)*x[1]*g))]);
        fi;
      od;
      s := s * (1/(Size(Centralizer(G, x[1])) * Size(Centralizer(G, y[1])))) * Size(G);
      #S[x_idx][y_idx] := s;
      #S[y_idx][x_idx] := s;
      Add(l,s);
     od;
     Add(S,l);
  od;

  # compute the T-matrix
  T := [];
  for x in Simples do
    Add(T, x[2][Position(x[3], ConjugacyClass(Centralizer(G, x[1]), x[1]))]/x[2][1]);
  od;

  if display then
    Print("\n The S-matrix is: \n\n");
    Print(S);
    Print("\n\n The T-matrix is: \n\n");
    Print(T);
    Print("\n\n");
  fi;

  pseudounitary_data := pseudounitary_sort(S,T,Simples);
  Simples := pseudounitary_data[1];
  S := pseudounitary_data[2];
  T := pseudounitary_data[3];
  FPdim := pseudounitary_data[4];
  FPdimC := pseudounitary_data[5];

  # sort by dimensions (if two objects have the same dimension, their
  # relative order is unchanged).
  sorted_data := dimension_sort(S,T,Simples,FPdim);
  Simples:=sorted_data[1];
  S:=sorted_data[2];
  T:=sorted_data[3];
  FPdim:=sorted_data[4];

  return [Simples, S, T, FPdim, FPdimC];
end);

################################################################################
##
#F DG_data(<group>) . . . . . . . . . . . . . . .Computes modular data for D(G).
InstallGlobalFunction(DG_data, function(G)
  local Simples, S, T, FPdim, FPdimC, N, data;
  data:=DG_S_and_T(G,false);
  Simples:=data[1];
  S:=data[2];
  T:=data[3];
  FPdim:=data[4];
  FPdimC:=data[5];
  N:=compute_fusion_rules(S,FPdimC);
  return [Simples,S,T,N,FPdim,FPdimC];
end);

#E DG_data.gi . . . . . . . . . . . . . . . . . . . . . . . . . . . . .ends here

max_stars_repo_name: pnnl/GroupTheoretical

max_issues_repo_name: pnnl/GroupTheoretical

max_forks_repo_name: pnnl/GroupTheoretical


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

# Permutation from standard to gap
# / 1 2 3 4 \
# \ 2 4 3 1 /
# Position(pr,Permuted([2,3,4,2],(2,4,1)))=24;

#......Ross Jordan Block Sizes........Donna & Ross nilpotents........(in GAP notation)
reps:=[[[2,1],[4,1],[3,1],[1,1]],    # F_4
       [[2,1],[4,1],[7,1],[5,1]],    # F_4(a_1) [[2,1],[4,1],[10,1],[1,1]],
       #[[6,1],[10,1],[1,1],[5,1]],   # F_4(a_2) [[9,1],[1,1],[10,1],[4,1]],
       # F_4(a_2): s_11s_9*[4,3,1,24]=[4,5,9,10] also good .. has 2 semisimple elements?
       [[1,1],[4,1],[9,1],[10,1]],   # F_4(a_2) good for all p\neq 2 <C_3A_1
       [[4,1],[3,1],[1,1]],          # C_3      [[1,1],[3,1],[4,1]],
       [[2,1],[4,1],[3,1]],          # B_3      the same
       [[4,1],[6,1],[10,1],[18,1]],  # F_4(a_3) [[4,1],[13,1],[11,1],[12,1]],
       [[4,1],[1,1],[10,1]],         # C_3(a_1) [[1,1],[10,1],[4,1]],
       [[2,1],[3,1],[1,1]],          # sA_2A_1  [[3,1],[1,1],[2,1]],
       [[4,1],[3,1]],                # B_2      the same
       [[2,1],[4,1],[1,1]],          # A_2sA_1  the same
       [[3,1],[1,1]],                # sA_2     the same
       [[2,1],[4,1]],                # A_2      the same
#       [[2,1],[3,1]],                # A_1sA_1  [[2,1],[1,1]],
#       [[3,1]],                      # sA_1     [[1,1]],
#       [[2,1]]                       # A_1      [[2,1]]
       [[2,1],[1,1]],                # A_1sA_1  [[2,1],[1,1]],
       [[1,1]],                      # sA_1     [[1,1]],
       [[2,1]]                       # A_1      [[2,1]]
       ];

I:=[[],     # F_4
    [3],    # F_4(a_1)
    [3,4],  # F_4(a_2)
    [],     # C_3
    [],     # B_3
    [2,3,1],# F_4(a_3)
    [3],    # C_3(a_1)
    [],     # sA_2A_1
    [],     # B_2
    [],     # A_2sA_1
    [],     # sA_2
    [],     # A_2
    [],     # A_1sA_1
    [],     # sA_1
    []      # A_1
    ];

#
# The parabolic where the representative is distinguished
#
II:=[[1,2,3,4],   # F_4
     [1,2,3,4],   # F_4(a_1)
     [1,2,3,4],   # F_4(a_2)
     [4,3,1],     # C_3
     [2,4,3],     # B_3
     [1,2,3,4],   # F_4(a_3)
     [4,3,1],     # C_3(a_1)
     [3,1,2],     # sA_2A_1
     [4,3],       # B_2
     [2,4,1],     # A_2sA_1
     [3,1],       # sA_2
     [2,4],       # A_2
     [2,1],       # A_1sA_1
     [1],         # sA_1
     [2]          # A_1
     ];

labels:=["F_4",
         "F_4(a_1)",
         "F_4(a_2)",
         "C_3",
         "B_3",
         "F_4(a_3)",
         "C_3(a_1)",
         "\\tilde A_2A_1",
         "B_2",
         "A_2\\tilde A_1",
         "\\tilde A_2",
         "A_2",
         "A_1\\tilde A_1",
         "\\tilde A_1",
         "A_1"
         ];

tori:=[
       [2,2,2,2],  # F_4
       [2,2,0,2],  # F_4(a_1)
       [2,0,0,2],  # F_4(a_2)
       [2,-9,2,2], # C_3
       [-6,2,2,2], # B_3
       [0,0,0,2],  # F_4(a_3)
       [2,-5,0,2], # C_3(a_1)
       [2,2,2,-5], # sA_2A_1
       [-3,-4,2,2],# B_2
       [2,2,-3,2], # A_2sA_1
       [2,0,2,-4], # sA_2
       [0,2,-2,2], # A_2
       [2,2,-1,-1],# A_1sA_1
       [2,0,-1,0], # sA_1
       [0,2,0,-1]  # A_1
    ];

#blocks:=[];
blocks:=[[[19, 1], [15, 1], [9, 2]],                      # F_4
         [[9, 5], [7, 1]],                                # F_4(a_1)
         [[9, 3], [7, 1],[6, 2], [3, 2]],                 # F_4(a_2)
         [[9, 2], [8, 2], [6, 2], [3, 1], [1, 3]],        # C_3
         [[9, 2], [7, 4], [3, 1], [1, 3]],                # B_3
         [[7, 1], [6, 2], [5, 3], [3, 6]],                # F_4(a_3)
         [[7, 1], [6, 2], [5, 1], [4, 4], [3, 3], [1, 3]],# C_3(a_1)
         [[3, 16], [2, 2]],                               # sA_2A_1
         [[7, 1], [5, 4], [4, 4], [3, 1], [1, 6]],        # B_2
         [[3, 16], [2, 2]],                               # A_2sA_1
         [[3, 15], [1, 7]],                               # sA_2
         [[3, 15], [1, 7]],                               # A_2
         [[3, 10], [2, 8], [1, 6]],                       # A_1sA_1
         [[3, 7], [2, 8], [1, 15]],                       # sA_1
         [[3, 1], [2, 14], [1, 21]]                       # A_1
         ];

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

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


text:
InstallMethod(Witt,
              "Constructor for Witt Structure",
              [IsPosInt,IsAlgebraicU],
              function(dim,sys)
              local object,
              BPR,xvarnames,yvarnames,xvars,
              p,max_level,generic_poly;

              object:=Objectify(NewType(NewFamily("WittFamily"),
                                        IsAttributeStoringRep and
                                        IsWitt),
                                rec());

              SetDimension(object,dim);
              SetCharacteristic(object,Characteristic(sys));

              SetalgebraicU(object,sys);

              xvarnames:=List([1..dim],i->Concatenation("x_",String(i)));
              yvarnames:=List([1..dim],i->Concatenation("y_",String(i)));
              BPR:=PolynomialRing(ring(sys),Concatenation(xvarnames,yvarnames));
              xvars:=IndeterminatesOfPolynomialRing(BPR){[1..dim]};

              SetchevalleyAdjBPR(object,ChevalleyAdj(sys,BPR));

              SetGenericSum(object,GenericSumOp(object));# this needs to be set after Characteristic

              p:=Characteristic(sys);

              #max_level:=Maximum(List(positiveRoots(simpleAdj(sys)),i->Sum(i)));
              max_level:=20;
              generic_poly:=function(degree)
                  local result,
                  tuples,tuple_degreeP;

                  tuple_degreeP:=function(t)
                      return Sum(List([1..Length(t)],i->t[i]*p^(i-1)));
                  end;

                  #tuples:=Tuples([1..dim+1]-1,dim);
                  tuples:=Tuples([1..max_level]-1,dim);
                  tuples:=Filtered(tuples,i->tuple_degreeP(i)=degree);
                  tuples:=List(tuples,t->List([1..dim],i->xvars[i]^t[i]));
                  tuples:=List(tuples,t->Product(t));

                  return tuples;
              end;
              #SetGenericLevelPolynomials(object,List([1..max_level],i->generic_poly(i)));

              generic_poly:=function(level)
                  local result,
                  i,Unu,tmp;

                  if level=0 then return [List([1..dim],i->0)]; fi;
                  if level<0 then return []; fi;

                  Unu:=DiagonalMat(List([1..dim],i->1));

                  result:=[];
                  for i in [1..dim] do
                      tmp:=generic_poly(level-p^(i-1));
                      if tmp<>[] then
                          result:=Concatenation(result,List(tmp,j->j+Unu[i]));
                      fi;
                  od;

                  tmp:=[];
                  for i in result do
                      if not i in tmp then
                          tmp:=Concatenation(tmp,[i]);
                      fi;
                  od;

                  return tmp;#result;
              end;
              SetGenericLevelPolynomials(object,List([1..max_level],
                                                     i->List(generic_poly(i),
                                                             j->Product([1..dim],
                                                                        k->xvars[k]^j[k]))));

              SetName(object,Concatenation("<witt structure of dimension ",String(dim),
                                           " over ring of characteristic ",String(p),">"));
              return object;
end);

InstallMethod(GenericSumOp,
              "Polynomials for generic sum in the Witt structure",
              [IsWitt],
              function(w)
              local result,
              p,dim,phiX,phiY,psiXY,phiXY,
              localIPR,Xvarnames,Yvarnames,X,Y,ext_reps,
              BPR,xyvars,i,j,k,poly,monom,pos;

              dim:=Dimension(w);
              p:=Characteristic(w);
              BPR:=ring(chevalleyAdjBPR(w));

              Xvarnames:=List([1..dim],i->Concatenation("X",String(i)));
              Yvarnames:=List([1..dim],i->Concatenation("Y",String(i)));
              localIPR:=PolynomialRing(Integers,Concatenation(Xvarnames,Yvarnames));
              X:=IndeterminatesOfPolynomialRing(localIPR){[1..dim]};
              Y:=IndeterminatesOfPolynomialRing(localIPR){dim+[1..dim]};

              phiX:=List([1..Length(X)],i->Sum(List([1..i],j->p^(j-1)*X[j]^(p^(i-j)))) );
              phiY:=List([1..Length(Y)],i->Sum(List([1..i],j->p^(j-1)*Y[j]^(p^(i-j)))) );

              psiXY:=[];
              phiXY:=phiX+phiY;

              for i in [1..Length(phiXY)] do
                  psiXY:=Concatenation(psiXY,[(1/p^(i-1))*(phiXY[i]-Sum(List([1..i-1],j->p^(j-1)*psiXY[j]^(p^(i-j)))))]);
              od;

              ext_reps:=List(Concatenation(X,Y),i->ExtRepPolynomialRatFun(i));
              ext_reps:=List(ext_reps,i->i[1][1]);

              xyvars:=IndeterminatesOfPolynomialRing(BPR);
              psiXY:=List(psiXY,i->List(ExtRepPolynomialRatFun(i)));
              result:=[];
              for i in [1..Length(psiXY)] do
                  poly:=Zero(BPR);
                  for j in [1..Length(psiXY[i])/2] do
                      monom:=One(BPR);
                      for k in [1..Length(psiXY[i][2*j-1])/2] do
                          pos:=Position(ext_reps,psiXY[i][2*j-1][2*k-1]);
                          monom:=monom*xyvars[pos]^psiXY[i][2*j-1][2*k];
                      od;
                      monom:=monom*psiXY[i][2*j];
                      poly:=poly+monom;
                  od;
                  Add(result,poly);
              od;

              return result;

    # adjust the form to the one I need
#    psiXY_:=List([1..dim],i->psiXY[i]-X[i]-Y[i]);
end);


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

Declare(TDCT2, TDCT3, TDCT4);

_DCT_CONST := 2.5;

#######################################################################################
#   tSPL DCT rules

#    PDCT4_CT_SPL := rec(
#    isApplicable     := P -> P > 2 and ForAny(DivisorPairs(2*P), d->IsEvenInt(d[1]) and IsEvenInt(d[2])),
#    forTransposition := false,
#
#    allChildren  := P -> let(N := P,
#        Map2(Filtered(DivisorPairs(2*N), d -> IsEvenInt(d[1]) and IsEvenInt(d[2])),
#        (m,n) -> [ PRDFT3(m,-1).transpose(), PRDFT3(n) ])),
#
#    rule := (P,C) -> let(N := P, m := Rows(C[1]), n := Cols(C[2]),
#        j := Ind(n/2),
#        T := Diag(Twid(2*N, m/2, 1, 1/2, 1/2, j)),
#
#        Prm(Refl(N, 2*N-1, N, L(m*n, m))) *
#        IterDirectSum(j, j.range, Diag(BHN(m)) * C[1] * RC(T)) *
#        RC(L(m*n/4, n/2)) *
#        Tensor(I(m/2), C[2] * Diag(BHN(n))) *
#        Prm(Refl(N, 2*N-1, N, L(m*n, m)))
#        )
#    )
#

#F TDCT4(<size>, tags) - Discrete Cosine Transform, Type IV, non-terminal
#F Definition: (n x n)-matrix  [ cos((k+1/2)*(l+1/2)*pi/n) | k,l = 0...n-1 ]
#F Note:       DCT4 is symmetric
#F Example:    DCT4(8)
Class(TDCT4, TaggedNonTerminal, rec(
    abbrevs := [ N -> Checked(IsInt(N), N >= 1, [N]) ] ,
    dims := self >> [self.params[1], self.params[1]],
    isReal := True,
    terminate := self >> Mat(DCT_IVunscaled(self.params[1])),
    transpose := self >> Copy(self),
    SmallRandom := () -> Random([2, 3, 4, 6, 8, 9, 12, 16, 18, 24, 27, 32]),
    normalizedArithCost := (self) >> let(n := self.params[1], IntDouble(_DCT_CONST * n * d_log(n) / d_log(2)))
));

NewRulesFor(TDCT4, rec(
#    DCT4_CT_tSPL := rec(
#        switch := false,
#
#        applicable := (self, t) >> let(P:=t.params, P[1] > 2 and ForAny(DivisorPairs(2*P[1]), d->IsEvenInt(d[1]) and IsEvenInt(d[2])) and HasTags(t)),
#
#        children := (self, t) >> let(tags := GetTags(t), N := t.params[1],
#                                    Map2(Filtered(DivisorPairs(2*N), d -> IsEvenInt(d[1]) and IsEvenInt(d[2])),
#                                        (m,n) -> [
#                                            SetTag(TCompose([
#                                                TScat(Refl(N, 2*N-1, N, L(m*n, m))),
#                                                TTensorI(Diag(BHN(m)) * PRDFT3(m,-1).transpose(), n/2, APar, APar),
#                                                TRC(TDiag(fPrecompute(
#                                                    fCompose(dOmega(8 * N, 1), diagTensor(dLin(N/m, 2, 1, TInt), dLin(m/2, 2, 1, TInt)))
#                                                ))),
#                                                TRC(TL(m*n/4, n/2, 1, 1)),
#                                                TTensorI(PRDFT3(n) * Diag(BHN(n)), m/2, APar, APar),
#                                                TGath(Refl(N, 2*N-1, N, L(m*n, m))) ]),
#                                                tags)
#                                        ])),
#
#        apply := (self, t, C, Nonterms) >> C[1]
#    ),
#
#    # Variant better suited for vectorization
#    PDCT4_CT_SPL_Vec := rec(
#    isApplicable     := P -> P > 2 and (P mod 2) = 0,
#    forTransposition := false,
#
#    allChildren  := P -> let(N := P,
#        Map2(Filtered(DivisorPairs(2*N), d -> IsEvenInt(d[1]) and IsEvenInt(d[2])),
#        (m,n) -> [ PRDFT3(m,-1).transpose(), PRDFT3(n) ])),
#
#    rule := (P,C) -> let(N := P, m := Rows(C[1]), n := Cols(C[2]),
#            TT := Diag(fPrecompute(fCompose(dOmega(8 * N, 1),
#                        diagTensor(dLin(N/m, 2, 1, TInt), dLin(m/2, 2, 1, TInt))))),
#
#            IJ(N, n/2) *
#            L(N, m) *
#        Tensor(I(n/2),
#                   M(m,m/2) * Diag(BHN(m)) * C[1]) *
#            RC(TT) *
#            L(N, n/2) *
#        Tensor(I(m/2),
#                   L(n, 2) * C[2] * Diag(BHN(n)) * K(n, 2)) *
#            L(N, m/2) *
#            IJ(N, m/2)
#    )
#    )

    DCT4_CT_tSPL := rec(
        switch := false,

        applicable := (self, t) >> let(
            P:=t.params,
            P[1] > 2
            and ForAny(DivisorPairs(2*P[1]), d ->
                IsEvenInt(d[1]) and IsEvenInt(d[2])
            )
            and t.hasTags()
        ),

        children := (self, t) >> let(
            tags := t.getTags(),
            N := t.params[1],
            Map2(
                Filtered(DivisorPairs(2*N), d -> IsEvenInt(d[1]) and IsEvenInt(d[2])),
                (m,n) -> List([
                    TPrm(IJ(N, n/2)),
                    TTensorI(condM(m,m/2) * Diag(BHN(m)) * PRDFT3(m,-1).transpose(), n/2, AVec, AVec),
                    TTensorI(L(n, 2) * PRDFT3(n) * Diag(BHN(n)) * condK(n, 2), m/2, APar, AVec),
                    TPrm(IJ(N, m/2))
                ], i -> i.setTags(tags))
            )
        ),

        apply := (self, t, C, Nonterms) >> let(
            N:=t.params[1],
            n:=2*Rows(Nonterms[1].params[1].params[2]),
            m:=2*Rows(Nonterms[4].params[1].params[2]),
            Grp(
                C[1] * C[2]
                * ConjDiag(RC(Diag(fPrecompute(
                    fCompose(dOmega(8 * N, 1), diagTensor(dLin(N/m, 2, 1, TInt), dLin(m/2, 2, 1, TInt)))
                ))), L(N, m), L(N, n/2))
            )
            * C[3] * C[4]
        )
    )
));

#F TDCT2(<n>) - Discrete Cosine Transform, Type II, non-terminal
#F Definition: (n x n)-matrix [ cos(k*(l+1/2)*pi/n) | k,l = 0...n-1 ]
#F Note:       DCT2 is the transpose of DCT3
#F Example:    DCT2(8)
Class(TDCT2, TaggedNonTerminal, rec(
    abbrevs := [ N -> Checked(IsInt(N), N >= 1, [N]) ] ,
    dims := self >> [self.params[1], self.params[1]],
    isReal := True,
    terminate := self >> Mat(DCT_IIunscaled(self.params[1])),
    transpose := self >> TDCT3(self.params[1]),
    SmallRandom := () -> Random([2,3,4,5,6,8,9,10,12,15,16,18,24,27,30,32]),
    normalizedArithCost := (self) >> let(n := self.params[1], IntDouble(_DCT_CONST * n * d_log(n) / d_log(2)))
));

NewRulesFor(TDCT2, rec(
    DCT2_DCT4_tSPL := rec(
        switch := false,
        applicable := (self, t) >> true,
        children := (self, t) >> [[ TS(t.params[1]).withTags(t.getTags()), TDCT4(t.params[1]).withTags(t.getTags()) ]],
        apply := (self, t, C, Nonterms) >> let(P := t.params[1],
            Diag(Concat([V(2.0)], List([1..P-1], i->V(1.0)))) *
            C[1].transpose() * C[2] *
            Diag(List([0..P - 1], i -> 1/(2 * CosPi((2*i + 1)/(4*P))))))
    ))
);

#F TDCT3(<n>) - Discrete Cosine Transform, Type III, non-terminal (unscaled)
#F Definition: (n x n)-matrix [ cos((k+1/2)*l*pi/n) | k,l = 0...n-1 ]
#F  [scaled]   (n x n)-matrix [ a_l*cos((k+1/2)*l*pi/n) | k,l = 0...n-1 ]
#F                              a_j = 1/sqrt(2) for j = 0 and = 1 else
#F Note:       DCT3 is the transpose of DCT2, scaled NOT supported yet
#F Example:    DCT3(8)
Class(TDCT3, TaggedNonTerminal, rec(
    abbrevs := [ N -> Checked(IsInt(N), N >= 1, [N]) ] ,
    dims := self >> [self.params[1], self.params[1]],
    isReal := True,
    terminate := self >> Mat(DCT_IIIunscaled(self.params[1])),
    transpose := self >> TDCT2(self.params[1]),
    SmallRandom := () -> Random([2,3,4,5,6,8,9,10,12,15,16,18,24,27,30,32]),
    normalizedArithCost := (self) >> let(n := self.params[1], IntDouble(_DCT_CONST * n * d_log(n) / d_log(2)))
));

NewRulesFor(TDCT3, rec(
    DCT3_DCT2_tSPL := rec(
        switch := false,
        applicable := (self, t) >> true,
        children := (self, t) >> [[ TDCT2(t.params[1]).withTags(t.getTags()) ]],
        apply := (self, t, C, Nonterms) >> C[1].transpose()
    ),
    DCT3_DCT4_tSPL := rec(
        switch := false,
        applicable := (self, t) >> true,
        children := (self, t) >> [[ TDCT4(t.params[1]).withTags(t.getTags()), TS(t.params[1]).withTags(t.getTags()) ]],
        apply := (self, t, C, Nonterms) >> let(P := t.params[1],
            Diag(List([0..P - 1], i -> 1/(2 * CosPi((2*i + 1)/(4*P))))) *
            C[1] * C[2] *
            Diag(Concat([V(2.0)], List([1..P-1], i->V(1.0)))))
    )
));


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

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

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

handleUaUbUc:=function()
	local rels, vals,i,v;
	# UaUb=Uc = U(a1+b1,a2+b2,a3+b3+a2*b1)
	rels:=Set(Concatenation(Ua*Ub-Uc));

	vals:=[];
	Append(vals,[[[vars[22]],[vars[2]+vars[12]]]]);
	Append(vals,[[[vars[21]],[vars[1]+vars[11]]]]);
	Append(vals,[[[vars[23]],[vars[3]+vars[13]+vars[2]*vars[11]]]]);

	for i in [1..Length(rels)] do
		for v in vals do
			Print(rels[i],"\n");
			rels[i]:=One(APR)*Value(rels[i],v[1],v[2]);
		od;
	od;
	return vals;
end;

max_stars_repo_name: iuliansimion/groupsSB

max_issues_repo_name: iuliansimion/groupsSB

max_forks_repo_name: iuliansimion/groupsSB