text:
InstallMethod(TotalDegree,
              "TotalDegree of multivariable polynomial",
              [IsPolynomial],
              function(poly)
              local result;

              result:=ExtRepPolynomialRatFun(poly);
              result:=result{2*[1..Length(result)/2]-1};
              result:=Maximum(List(result,j->Sum(j{2*[1..Length(j)/2]})));
              return result;
end);

InstallMethod(ValuePoly,
              "Value of multivariable polynomial",
              [IsPolynomial,IsList,IsList],
              function(poly,vars,vals)
    local result,ext_rep,ext_vars,
    monom,pos,i,j;

    ext_rep:=ExtRepPolynomialRatFun(poly);
    ext_vars:=List(vars,i->ExtRepPolynomialRatFun(i));
    ext_vars:=List(ext_vars,i->i[1][1]);

    result:=Zero(poly);
    for i in 2*[1..Length(ext_rep)/2]-1 do
        monom:=One(poly);
        for j in 2*[1..Length(ext_rep[i])/2]-1 do
            pos:=Position(ext_vars,ext_rep[i][j]);
            monom:=monom*vals[pos]^ext_rep[i][j+1];
        od;
        monom:=monom*((ext_rep[i+1]*One(poly)));
        result:=result+monom;
    od;

    return result;
end);

pRadical:=function(poly,avars)
    local result,
    ext_rep,ext_avars,
    p,deg,deg_max,p_div,monom,i,j;

    #exterior reprezentation of polynomial and variables
    ext_rep:=ExtRepPolynomialRatFun(poly);
    ext_avars:=List(avars,i->ExtRepPolynomialRatFun(i));
    ext_avars:=List(ext_avars,i->i[1][1]);

    # degree of variables in monomials of poly
    deg:=List(ext_rep{2*[1..Length(ext_rep)/2]-1},i->i{2*[1..Length(i)/2]});

    # Adjustment: (poly)^p ---> poly
    p:=Characteristic(poly);
    if Concatenation(deg)=[] then return poly; fi;
    deg_max:=LogInt(Maximum(Concatenation(deg)),p);
    p_div:=Gcd(p^deg_max,Gcd(Concatenation(deg)));
    if p_div>1 then
        result:=Zero(poly);
        for i in 2*[1..Length(ext_rep)/2]-1 do
            monom:=One(poly);
            for j in 2*[1..Length(ext_rep[i])/2]-1 do
                monom:=monom*avars[Position(ext_avars,ext_rep[i][j])]^(ext_rep[i][j+1]/p_div);
            od;
            result:=result+ext_rep[i+1]*monom;
        od;
        return result;
    fi;
    return poly;
end;

Solve:=function(poly,avars)
    local ext_rep,vars,ext_avars,p,
    deg,tot_deg,poz1,pozm,xpoly,xall,x,
    taboo;

    # Adjustment: (poly)^p ---> poly
    poly:=pRadical(poly,avars);

    #exterior reprezentation of polynomial and variables
    ext_rep:=ExtRepPolynomialRatFun(poly);
    ext_avars:=List(avars,i->ExtRepPolynomialRatFun(i));
    ext_avars:=List(ext_avars,i->i[1][1]);

    # variables of poly
    vars:=Concatenation(ext_rep{2*[1..Length(ext_rep)/2]-1});
    vars:=vars{2*[1..Length(vars)/2]-1};
    # degree of variables in monomials of poly
    deg:=List(ext_rep{2*[1..Length(ext_rep)/2]-1},i->i{2*[1..Length(i)/2]});
    tot_deg:=List(deg,i->Sum(i));
    #ignor constants
    tot_deg:=Filtered(tot_deg,i->i<>0);

    # Restrict: x should not be deduced from x+x^p+y
    taboo:=[];
    if Length(Set(vars))>1 then
    for x in Set(vars) do
        xall:=ext_rep{2*[1..Length(ext_rep)/2]-1};
        xall:=Filtered(xall,i->x in i);
        xpoly:=Filtered(xall,i->Length(i)=2);
        if Length(xpoly)>=2 then
        #    if xall=xpoly then
        #        return [Position(ext_avars,x),Zero(poly),avars[Position(ext_avars,x)]];
        #    fi;
            Add(taboo,Positions(vars,x));
        fi;
    od;
    taboo:=Concatenation(taboo);
    fi;

    p:=Characteristic(poly);
#    if Length(vars)=1 then
#        return [Position(ext_avars,vars[1]),RootsOfUPol(poly)[1],true];
#    fi;
#Error("!");
    if 1 in tot_deg then
        # x^p+P(x)^p
        if Length(Set(vars)) = 1 and Length(Set(tot_deg))>=2 then
#        Print(poly,"\n");
            return [Position(ext_avars,vars[1]),Zero(poly),avars[Position(ext_avars,vars[1])]];
#            return fail;
        fi;

##        if Length(Set(vars)) = 1 and not Difference(tot_deg,[1])[1] mod p =0 then
#        if Length(Set(vars)) = 1 and not Set(tot_deg)=[1,p] then
#            Print(poly,"\n");
#            return [Position(ext_avars,vars[1]),Zero(poly),avars[Position(ext_avars,vars[1])]];
##            return fail;
#        fi;

        if Length(Set(vars)) = 1 and Maximum(tot_deg)>1 then return fail; fi; # such relations produce infinite loops

        # position of variables among monomials of poly
        poz1:=Positions(tot_deg,1);
        poz1:=Difference(poz1,taboo);
        # if all variables are taboo return fail
        if Length(poz1)=0 then return fail; fi;
        # else choose one which is not taboo
        poz1:=poz1[1];
        # position of monomial in ext_rep
        pozm:=Position(ext_rep,[vars[poz1],1]);
        return [Position(ext_avars,vars[poz1]),avars[Position(ext_avars,vars[poz1])]-(ext_rep[pozm+1]^(-1))*poly,true];
    fi;
    return fail;
end;

SolveWithTorsion:=function(poly,avars)
    local ext_rep,vars,ext_avars,
    deg,tot_deg,poz1,pozm,xpoly,x,
    taboo;

    # Adjustment: (poly)^p ---> poly
    poly:=pRadical(poly,avars);

    #exterior reprezentation of polynomial and variables
    ext_rep:=ExtRepPolynomialRatFun(poly);
    ext_avars:=List(avars,i->ExtRepPolynomialRatFun(i));
    ext_avars:=List(ext_avars,i->i[1][1]);

    # variables of poly
    vars:=Concatenation(ext_rep{2*[1..Length(ext_rep)/2]-1});
    vars:=vars{2*[1..Length(vars)/2]-1};
    # degree of variables in monomials of poly
    deg:=List(ext_rep{2*[1..Length(ext_rep)/2]-1},i->i{2*[1..Length(i)/2]});
    tot_deg:=List(deg,i->Sum(i));
    #ignor constants
    tot_deg:=Filtered(tot_deg,i->i<>0);

    # Restrict: x should not be deduced from x+x^p+y
    taboo:=[];
    if Length(Set(vars))>1 then
    for x in Set(vars) do
        xpoly:=ext_rep{2*[1..Length(ext_rep)/2]-1};
        xpoly:=Filtered(xpoly,i->x in i);
        xpoly:=Filtered(xpoly,i->Length(i)=2);
        if Length(xpoly)>=2 then
            Add(taboo,Positions(vars,x));
        fi;
    od;
    taboo:=Concatenation(taboo);
    fi;

#    if Length(vars)=1 then
#        return [Position(ext_avars,vars[1]),RootsOfUPol(poly)[1],true];
#    fi;
#Error("!");
    if 1 in tot_deg then
        # x^p+P(x)^p
        if Length(Set(vars)) = 1 and Length(Set(tot_deg))>=2 then
            return [Position(ext_avars,vars[1]),Zero(poly),avars[Position(ext_avars,vars[1])]];
        fi;

        if Length(Set(vars)) = 1 and Maximum(tot_deg)>1 then return fail; fi; # such relations produce infinite loops

        # position of variables among monomials of poly
        poz1:=Positions(tot_deg,1);
        poz1:=Difference(poz1,taboo);
        # if all variables are taboo return fail
        if Length(poz1)=0 then return fail; fi;
        # else choose one which is not taboo
        poz1:=poz1[1];
        # position of monomial in ext_rep
        pozm:=Position(ext_rep,[vars[poz1],1]);
        return [Position(ext_avars,vars[poz1]),avars[Position(ext_avars,vars[poz1])]-(ext_rep[pozm+1]^(-1))*poly,true];
    fi;
    return fail;
end;

ReduceRelations:=function(relations,avars,values,APR)
    local result,#relations,
    relations_len,schimbare,Ring,torsion,GB,
    index,sol,i,modulo;

    result:=ShallowCopy(values);
    relations:=List(relations,i->Value(i,avars,result)*One(APR));
    relations:=Filtered(relations,i->i <> Zero(APR));
    relations:=Set(relations);
#    relations:=List(coefficients(u),i->i[2]);
    relations_len:=Length(relations);

#    i:=1;
#    while i<relations_len do
#        modulo:=Concatenation([1..i-1],[i+1..relations_len]);
#        modulo:=relations{modulo};
#        relations:=List(relations,i->PolynomialReducedRemainder(i,relations,MonomialLexOrdering()));
#        relations:=Set(Filtered(relations,i->i <> Zero(APR)));
#        if Length(relations)<relations_len then i:=1; fi;
#        relations_len:=Length(relations);
#    od;
#    Error("!");
#Print("...relations: ",relations,"\n");
    schimbare:=true;
    while relations<>[] and schimbare do
        #Print("@@@",relations_len," ");
#        Print("..relations: ",relations,"\n");
#        Print("--->result: ",result,"\n");
        #GB:=GroebnerBasis(relations,MonomialLexOrdering()); ## ATENTIE AICI
        result:=List(result,i->PolynomialReducedRemainder(i,relations,MonomialLexOrdering()));
        #result:=List(result,i->PolynomialReducedRemainder(i,GB,MonomialLexOrdering()));
#        Print("fresh result:",result,"\n");
        schimbare:=false;
        #Error("!!");
        #refresh relations
        #relations:=List(relations,i->Value(i,avars,result)*One(APR));
        relations:=List(relations,i->ValuePoly(i,avars,result)*One(APR));
#        Print("fresh rels:",relations,"\n");
        relations:=Filtered(relations,i->i <> Zero(APR));
        schimbare:=Length(relations)<relations_len;
        relations_len:=Length(relations);
    od;
#Error("!");
#Print(".relations: ",relations,"\n");
    relations:=Filtered(relations,i->i <> Zero(APR));

    return [result,relations];
end;

SolveRelations:=function(relations,avars,values,APR)
    local result,#relations,
    relations_len,schimbare,Ring,torsion,
    index,sol,i;

    result:=ShallowCopy(values);
    relations:=List(relations,i->Value(i,avars,result)*One(APR));
    relations:=Filtered(relations,i->i <> Zero(APR));

    relations:=ReduceRelations(relations,avars,result,APR);
#    #Error("!");
    result:=relations[1];
    relations:=relations[2];
#    #relations:=List(coefficients(u),i->i[2]);
    relations_len:=Length(relations);
#        Print("...relations: ",relations,"\n");

    Ring:=CoefficientsRing(APR);
    schimbare:=true;
    torsion:=[];

    #Error("!");
    while relations<>[] and schimbare do
    ##Print("len:=",Length(relations),"\n");
        schimbare:=false;
        index:=1;
        while index<=relations_len do
        #Print("len===",Length(relations),"\n","index===",index,"\n");
            sol:=Solve(relations[index],avars);
            #Print("solved ",relations[index],"\n");
            #Error("!");
            #sol:=fail;
            if sol<>fail then
                #change result
                result[sol[1]]:=sol[2];
                for i in [1..Length(result)] do
                    if not result[i] in Ring and result[i]<>avars[i] then
                       result[i]:=Value(result[i],avars,result)*One(APR);
                    fi;
                od;
                #refresh relations
                relations:=List(relations,i->Value(i,avars,result)*One(APR));
                relations:=Filtered(relations,i->i <> Zero(APR));
                #relations:=ReduceRelations(relations,avars,result,APR);
                #result:=relations[1];
                #relations:=relations[2];
                schimbare:=Length(relations)<relations_len or sol[2]=Zero(APR);
                relations_len:=Length(relations);

                if Length(sol)=3 and sol[3] in avars then Add(torsion,sol[3]); fi;
            fi;
            index:=index+1;
        od;
    od;

    return [result,relations,torsion];
end;


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

#D TL.withTags := (self, tags) >> CopyFields(self, rec(params := [self.params[1], self.params[2], self.params[3], self.params[4], tags]));
#D TL.getTags := self >> self.params[self.tagpos];

DoubleDivisorPairs:=function(n,m)
    local a,b,l;
    l:=[];
    for a in DivisorPairs(n) do
        Append(l,[Concatenation(a,[1,m]),Concatenation(a,[m,1])]);
        for b in DivisorPairs(m) do
            Append(l,[Concatenation(a,b)]);
        od;
    od;
    for b in DivisorPairs(m) do
        Append(l,[Concatenation([n,1],b),Concatenation([1,n],b)]);
    od;
    return l;
end;

VectorDivisible:=function(m,c,k,b,v)
#not yet taken into account : b or c =1
    local leftpart,rightpart;
    leftpart := v in DivisorsInt(c) or m=1 or (b=1 and k=1);
    rightpart:= v in DivisorsInt(b) or c=1 or k=1;
    return leftpart and rightpart;
end;

DoubleDivisorPairsVector:=function(n,m,v)
    return Filtered(DoubleDivisorPairs(n,m) , a -> VectorDivisible(a[1],a[2],a[3],a[4],v));
end;

Blocking:=function(m,c,k,b,v)
    local t, nv, PageSize, CacheLineSize;
    t:=128/v;
    PageSize:=4*1024*8*(2);
    nv:=m*c*k*b*t;
    if (nv>PageSize) then
        return When((PageSize/8<=b*c*t)and(b*c*t<=PageSize),true,false);

        #CacheLineSize:=64*8*(4);
    #else
        #if (nv>CacheLineSize) then
            #return When((CacheLineSize/8<=b*c*t)and(b*c*t<=CacheLineSize),true,false);
        #fi;
    fi;
    return true;
end;

DoubleDivisorPairsVectorBlocking:=function(n,m,v)
    return Filtered(DoubleDivisorPairsVector(n,m,v), a->Blocking(a[1],a[2],a[3],a[4],v));
end;

NewRulesFor(TL,rec(
    L_cx_real := rec(
        switch:=false,
        forTransposition := false,
        applicable := t ->  t.isTag(1, AVecRegCx) and let(P := t.params, P[1] = P[2]^2), #FF: NOTE!! can only terminate symmetrical ones...
        children := nt -> let(P:=nt.params, [[TL(P[1], P[2], P[3], P[4]*2).withTags([AVecReg(nt.getTags()[1].isa)])]]),
        apply := (t, C, Nonterms) -> CR(When(t.params[1]=t.params[2]^2, SymSPL(C[1]), C[1]))
    ),

    L_base_vec := rec(
        switch:=false,
        forTransposition := false,
        applicable := t -> t.isTag(1, AVecReg) or t.isTag(1, AVecRegCx),
        apply := (t, C, Nonterms) -> let(
            C1:=When(t.params[3]=1, [], [I(t.params[3])]),
            C2:=When(t.params[4]=1, [], [I(t.params[4])]),
            Tensor(Concat(C1, [L(t.params[1], t.params[2])], C2))
        )
#D        applicable := (self, t) >> FirstTagEq(t, AVecReg) or FirstTagEq(t, AVecRegCx),
#D        apply := (t, C, Nonterms) -> let(C1:=When(t.params[3]=1, [], [I(t.params[3])]),
#D            C2:=When(t.params[4]=1, [], [I(t.params[4])]), Tensor(Concat(C1, [L(t.params[1], t.params[2])], C2)))),
    ),

   #L^mn_m -> (L^mn/v_m x I_v)(I_mn/v2 x L^v2_v)(I_n/v x L^m_m/v x I_v)
    L_mn_m_vec := rec(
        forTransposition := false,
        freedoms := (self, t) >> [],
        applicable := t ->
            t.params[3] = 1
            and t.params[4] = 1
            and (t.isTag(1, AVecReg) or t.isTag(1, AVecRegCx))
            and let(
                v := t.firstTag().v,
                m := t.params[2],
                n := t.params[1] / m,
                IsInt(m/v)
                and IsInt(n/v)
                and not (
                    t.params[1] = v*v
                    and t.params[2] = v
                )
            ),
        child := (nt,freedoms) -> let (v := nt.firstTag().v, [ TL(v*v, v, 1, 1).withTags(nt.getTags()) ] ),

        apply := (nt, C, cnt) -> let(
            v := nt.firstTag().v,
            v2 := v*v,
            m := nt.params[2],
            n := nt.params[1]/m,
            VTensor(L(m*n/v, m), v)
            * SymSPL(BlockVPerm(m*n/v2, v, C[1], L(v^2,v)))
            * VTensor(Tensor(I(n/v), L(m,m/v)), v)
        )

#D        applicable := (self, t) >> t.params[3] = 1 and t.params[4] = 1 and FirstTagEq(t, AVecReg) and
#D                        let (v := GetFirstTag(t).v, m:=t.params[2], n:= t.params[1]/m,
#D                               IsInt(m/v) and IsInt(n/v) and not(t.params[1]=v*v and t.params[2]=v))  ,
#D
#D        child := (nt,freedoms) -> let (v := GetFirstTag(nt).v, [TL(v*v, v, 1, 1, GetTags(nt))]),
#D
#D        apply := (nt,C,cnt) -> let(v:=GetFirstTag(nt).v, v2:=v*v, m:=nt.params[2], n:= nt.params[1]/m,
#D                VTensor(L(m*n/v, m), v) *
#D                SymSPL(BlockVPerm(m*n/v2, v, C[1], L(v^2,v))) *
#D                VTensor(Tensor(I(n/v), L(m,m/v)), v))),
    ),

#SymSPL-BlockVPerm can be replaced by Tensor(I(m*n/v2), C[1])

#L^mn_m could probably also be splitted like that:

#   #L^mn_m -> (L^mn/v_m x I_v)(I_mn/v2 x L^v2_v)(I_n/v x L^m_m/v x I_v)
#   L_mn_m_vec := rec(
#       forTransposition := false,
#       applicable := (self, t) >> t.params[3] = 1 and t.params[4] = 1 and FirstTagEq(t, AVecReg) and
#                        let (v := GetFirstTag(t).v, m:=t.params[2], n:= t.params[1]/m,
#                                  IsInt(m/v) and IsInt(n/v)),

#       freedoms := (self, t) >> [],
#       child := (nt,freedoms) -> let (v := GetFirstTag(nt).v, m:=nt.params[2], n:= nt.params[1]/m,
#         [TL(m*n/v, m, 1, v, GetTags(nt)),
#                TL(v*v, v, m*n/(v*v), 1, GetTags(nt)),
#                   TL(m,m/v,n/v, v, GetTags(nt))]),
#       apply := (nt,C,cnt) -> let(v:=GetFirstTag(nt).v, v2:=v*v, m:=nt.params[2], n:= nt.params[1]/m, C[1] * C[2] * C[3])),

#   #I x L -> BlockVPerm
#   IxLxI_BlockVPerm := rec(
#       forTransposition := false,
#       applicable := (self, t) >> FirstTagEq(t, AVecReg) and (t.params[4]=1)
#                   and let (v := GetFirstTag(t).v, (t.params[1]=v*v) and (t.params[2]=v)),
#       freedoms := (self, t) >> [],
#       child := (nt,freedoms) -> let (v := GetFirstTag(nt).v,[TL(v*v, v, 1, 1, GetTags(nt))]),
#       apply := (t, C, cnt) -> let(v:=GetFirstTag(t).v,
#SymSPL(BlockVPerm(t.params[3], v, C[1], DropTag(cnt[1]))))),

   #A x I_v -> VTensor
   IxLxI_vtensor := rec(
       forTransposition := false,
       applicable := t ->
            (
                t.isTag(1, AVecReg)
                or t.isTag(1, AVecRegCx)
            )
            and IsPosInt(t.params[4] / t.firstTag().v),
       apply := (t, C, Nonterms) -> let(
            v:=t.firstTag().v,
            l1 := When(t.params[3] > 1, [I(t.params[3])], []),
            l2 := When(t.params[4]/v > 1, [I(t.params[4]/v)], []),
            VTensor(Tensor(Concat(l1, [L(t.params[1], t.params[2])], l2)), v)
        )
#D       applicable := (self, t) >> (FirstTagEq(t, AVecReg) or FirstTagEq(t, AVecRegCx)) and
#D                              let(v:=GetFirstTag(t).v, IsPosInt(t.params[4]/v)),
#D       apply := (t, C, Nonterms) -> let(v:=GetFirstTag(t).v, l1:= When(t.params[3] > 1, [I(t.params[3])], []), l2:=When(t.params[4]/v > 1, [I(t.params[4]/v)], []) , VTensor(Tensor(Concat(l1, [L(t.params[1], t.params[2])], l2)), v))),
    ),

#GV1 shouldn't be enabled for stuff that are not strides!
#It is highly disruptive because it doesn't respect the VTensor VPerm nomenclatura
#F L_GV1: TL(kmbc,kb,r,s) = I_r tensor ((L(kbm,bk) tensor I_c) * (I_m tensor TL(bc,b,k,1)) * (I_m tensor L(kc,k) tensor I_b)) tensor I_s
    L_GV1 := rec(
        switch:=false,
        forTransposition := false,
        freedoms := nt -> [
            MapN(
                DoubleDivisorPairsVector(nt.params[1] / nt.params[2], nt.params[2], nt.firstTag().v),
                (m,c,k,b) -> [b,c]
            )
        ],
        child := (nt,freedoms) -> let(b:=freedoms[1][1],c:=freedoms[1][2], [TL(b*c,b,1,1).withTags(nt.getTags())]),
        apply := (nt,C,cnt) -> let(
            n:=nt.params[1],
            b:=cnt[1].params[2],
            c:=cnt[1].params[1]/b,
            k:=nt.params[2]/b,
            m:=nt.params[1]/(k*b*c),
            r:=nt.params[3],
            s:=nt.params[4],
            Tensor(
                I(r),
                Tensor(
                    Tensor(L(k*b*m,b*k),I(c))
                    * Tensor(I(m),Tensor(I(k),C[1]))
                    * Tensor(Tensor(I(m),L(k*c,k)),I(b)),
                    I(s))
            )
        )
#D       freedoms := nt -> [MapN(DoubleDivisorPairsVector(nt.params[1]/nt.params[2],nt.params[2],GetFirstTag(nt).v),(m,c,k,b) -> [b,c])],
#D       child := (nt,freedoms) -> let(b:=freedoms[1][1],c:=freedoms[1][2], [TL(b*c,b,1,1,GetTags(nt))]),
#D       apply := (nt,C,cnt) -> let(
#D       n:=nt.params[1],b:=cnt[1].params[2],c:=cnt[1].params[1]/b,k:=nt.params[2]/b,m:=nt.params[1]/(k*b*c),r:=nt.params[3],s:=nt.params[4],
#D       Tensor(I(r),(Tensor(
#D               Tensor(L(k*b*m,b*k),I(c))*
#D               Tensor(I(m),Tensor(I(k),C[1]))*
#D               Tensor(Tensor(I(m),L(k*c,k)),I(b))
#D               , I(s)))))),
    ),

    L_GV1_vtensor := rec(
        switch:=false,
        forTransposition := false,
        applicable := nt -> let(p := nt.params,
            nt.isTag(1, AVecReg)
            and p[1] <> p[2]
            and p[2] <> 1
            and Length(
                DoubleDivisorPairsVectorBlocking(p[1] / p[2], p[2], nt.firstTag().v)
            ) > 0
        ),
        freedoms := nt -> let(p := nt.params, [
            MapN(
                DoubleDivisorPairsVectorBlocking(p[1] / p[2], p[2], nt.firstTag().v),
                (m,c,k,b) -> [b,c]
            )
        ]),
        child := (nt, freedoms) -> let(
            b:=freedoms[1][1],
            c:=freedoms[1][2],
            [ TL(b*c,b,1,1).withTags(nt.getTags()) ]
        ),
        apply := (nt,C,cnt) -> let(
            n := nt.params[1],
            b := cnt[1].params[2],
            c := cnt[1].params[1]/b,
            k := nt.params[2]/b,
            m := nt.params[1]/(k*b*c),
            r := nt.params[3],
            s := nt.params[4],
            v := nt.firstTag().v,
#            garbage:=fPrint(["n:",n," b:",b," c:",c," k:",k," m:",m," r:",r," s:",s, " v:",v,"\n"]),
            prefactor := When(r = 1, [], [I(r)]),
            firstfactor := When(m = 1 or (b = 1 and k = 1),
                [],
                [When(c=1,L(k*b*m,b*k), VTensor(Tensor(L(k*b*m,b*k),I(c/v)),v))]
            ),
            middlefactor := [ Tensor(I(m), Tensor(I(k), C[1])) ],
#            middlefactor:=When(k>1 and b*c<1024,[Tensor(I(m),I(k/2),BB(Tensor(I(2),C[1])))],[Tensor(I(m),Tensor(I(k),C[1]))]),
            thirdfactor := When(c=1 or k=1,
                [],
                [
                    When(b=1,
                        Tensor(I(m),L(k*c,k)),
                        VTensor(Tensor(I(m),L(k*c,k),I(b/v)),v)
                    )
                ]
            ),
            postfactor := When(r=s, [], [I(s)]),

            Tensor(Concat(
                prefactor,
                [Compose(Concat(
                    firstfactor,
                    middlefactor,
                    thirdfactor
                ))],
                postfactor
            ))
        )
#D       applicable := (self, nt)>> FirstTagEq(nt, AVecReg) and (nt.params[1]<>nt.params[2]) and (nt.params[2]<>1) and
#D           (Length(DoubleDivisorPairsVectorBlocking(nt.params[1]/nt.params[2],nt.params[2],GetFirstTag(nt).v))>0),
#D       freedoms := nt -> [MapN(DoubleDivisorPairsVectorBlocking(nt.params[1]/nt.params[2],nt.params[2],GetFirstTag(nt).v),(m,c,k,b) -> [b,c])],
#D       child := (nt,freedoms) -> let(b:=freedoms[1][1],c:=freedoms[1][2], [TL(b*c,b,1,1,GetTags(nt))]),
#D       apply := (nt,C,cnt) -> let(
#D       n:=nt.params[1],b:=cnt[1].params[2],c:=cnt[1].params[1]/b,k:=nt.params[2]/b,m:=nt.params[1]/(k*b*c),r:=nt.params[3],s:=nt.params[4],v:=GetFirstTag(nt).v,
#D#       garbage:=fPrint(["n:",n," b:",b," c:",c," k:",k," m:",m," r:",r," s:",s, " v:",v,"\n"]),
#D       prefactor  :=When(r=1, [], [I(r)]),
#D       firstfactor:=When(m=1 or (b=1 and k=1), [],[When(c=1,L(k*b*m,b*k), VTensor(Tensor(L(k*b*m,b*k),I(c/v)),v))]),
#D       middlefactor:=[Tensor(I(m),Tensor(I(k),C[1]))],
#D#      middlefactor:=When(k>1 and b*c<1024,[Tensor(I(m),I(k/2),BB(Tensor(I(2),C[1])))],[Tensor(I(m),Tensor(I(k),C[1]))]),
#D       thirdfactor:=When(c=1 or k=1, [], [When(b=1,Tensor(I(m),L(k*c,k)), VTensor(Tensor(I(m),L(k*c,k),I(b/v)),v))]),
#D       postfactor :=When(r=s, [], [I(s)]),
#D       Tensor(Concat(
#D           prefactor,
#D           [Compose(Concat(
#D               firstfactor,
#D               middlefactor,
#D               thirdfactor))
#D           ],
#D           postfactor))))
    )
));

######################################################################
#  OBSOLETE RULES
######################################################################

NewRulesFor(TL,rec(

    # L^nv_n -> (I_n/v x L^v2_v)(L^n_n/v x I_v)
    # superseeded by L_mn_m_vec with n=v
    L_nv_n_vec := rec(
        switch := false,
        forTransposition := false,
        applicable := t ->
            t.params[3] = 1
            and t.params[4] = 1
            and (t.isTag(1, AVecReg) or t.isTag(1, AVecRegCx))
            and let (v := t.firstTag().v,
                t.params[1] = t.params[2] * v
                and t.params[2] <> v
                and IsInt(t.params[2]/v)
            ),
        freedoms := (self, t) >> [],
        child := (nt,freedoms) -> let(
            v := nt.firstTag().v,
            [ TL(v*v, v, 1, 1).withTags(nt.getTags()) ]
        ),
        apply := (nt, C, cnt) -> let (
            v := nt.firstTag().v,
            n := nt.params[2],
            SymSPL(
                BlockVPerm(n/v, v, C[1], L(v^2,v))
            ) * VTensor(L(n, n/v), v)
        )
#D       applicable := (self, t) >> t.params[3] = 1 and t.params[4] = 1 and FirstTagEq(t, AVecReg) and
#D                         let (v := GetFirstTag(t).v, t.params[1] = t.params[2] * v and t.params[2] <> v and IsInt(t.params[2]/v)),
#D       freedoms := (self, t) >> [],
#D       child := (nt,freedoms) -> let (v := GetFirstTag(nt).v, [TL(v*v, v, 1, 1, GetTags(nt))]),
#D       apply := (nt,C,cnt) -> let (v := GetFirstTag(nt).v, n:= nt.params[2], SymSPL(BlockVPerm(n/v, v, C[1], L(v^2,v))) * VTensor(L(n, n/v), v))),
   #SymSPL-BlockVPerm can be replaced by Tensor(I(n/v), C[1])
    ),

   #L^nv_v -> (L^n_v x I_v)(I_n/v x L^v2_v)
   # superseeded by L_mn_m_vec with m=v
   L_nv_v_vec := rec(
       switch := false,
       forTransposition := false,
       applicable := t ->
            t.params[3] = 1
            and t.params[4] = 1
            and (t.isTag(1, AVecReg) or t.isTag(1, AVecRegCx))
            and let(v := t.firstTag().v,
                t.params[2] = v
                and t.params[1] <> v * v
                and IsInt(t.params[1]/(v*v))
            ),
       freedoms := (self, t) >> [],
       child := (nt, freedoms) -> let(v := nt.firstTag().v, [ TL(v*v, v, 1, 1).withTags(nt.getTags()) ]),
       apply := (nt, C, cnt) -> let(
            v := nt.firstTag().v,
            n:= nt.params[1]/v,
            VTensor(L(n, v), v)
            * SymSPL(BlockVPerm(n/v, v, C[1], L(v^2,v)))
        )
#D       applicable := (self, t) >> t.params[3] = 1 and t.params[4] = 1 and FirstTagEq(t, AVecReg) and
#D                        let (v := GetFirstTag(t).v, t.params[2] = v and t.params[1] <> v * v and IsInt(t.params[1]/(v*v))),
#D       freedoms := (self, t) >> [],
#D       child := (nt,freedoms) -> let (v := GetFirstTag(nt).v, [TL(v*v, v, 1, 1, GetTags(nt))]),
#D       apply := (nt,C,cnt) -> let (v := GetFirstTag(nt).v, n:= nt.params[1]/v, VTensor(L(n, v), v) * SymSPL(BlockVPerm(n/v, v, C[1], L(v^2,v))))),
   #SymSPL-BlockVPerm can be replaced by Tensor(I(n/v), C[1])))
    )

));

ext: gi


text:
#!/usr/bin/gap

# time gap -o 8g hierarchy.gap

# Build clifford groups, up to 3 qubits, and
# search for T operators defined by field extensions......

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

# See:
# https://www.mathstat.dal.ca/~selinger/papers/clifford.pdf

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

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

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

for U in Cliff1 do
    found := false;;
    #Udag := Inverse(U);
    #Print("found? ");
    for g in Pauli1 do
        if (U*g*Inverse(U)*Inverse(g) in Pauli1)  then found:=true; break; fi;
    od;
    if not found then Print("Not found\n"); fi;
od;

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

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

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

# Works:
for U in Cliff2 do
    found := false;;
    #Udag := Inverse(U);
    #Print("found? ");
    for g in Pauli2 do
        if (U*g*Inverse(U)*Inverse(g) in Pauli2)  then found:=true; break; fi;
    od;
    if not found then Print("Not found\n"); fi;
od;

Hadamard4 := (1/2)*[
    [0, 1, 0, 1],
    [2, 0, 2, 0],
    [0, 1, 0, -1],
    [2, 0, -2, 0]];
#Print(Hadamard4*Hadamard4, "\n");
#Print(Hadamard4, " in Cliff2 ", Hadamard4 in Cliff2, "\n");

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

# Print(a in Cliff2, "\n"); # true

#hom:=EpimorphismFromFreeGroup(Cliff2:names:=["si", "is", "hi", "ih", "wi", "cz"]);
#Print(PreImagesRepresentative(hom, a), "\n"); # si^3*is^3*hi*si^2*ih*is^2*ih*cz*hi*si
#Print(PreImagesRepresentative(hom, a*a), "\n"); # (si^2*hi)^2
#Print(PreImagesRepresentative(hom, a*a*a), "\n"); # si^2*is^3*hi*(si^2*hi*si)^2*si*ih*is^2*ih*cz*hi*si

cls := ConjugacyClass(Cliff2, a);;
#Print(a,               a in cls, "\n");               # true
#Print(a*a,             a*a in cls, "\n");             # false
#Print(a*a*a,           a*a*a in cls, "\n");           # true
#Print(a*a*a*a,         a*a*a*a in cls, "\n");         # false
#Print(a*a*a*a*a,       a*a*a*a*a in cls, "\n");       # true
#Print(a*a*a*a*a*a,     a*a*a*a*a*a in cls, "\n");     # false
#Print(a*a*a*a*a*a*a,   a*a*a*a*a*a*a in cls, "\n");   # true
#Print(a*a*a*a*a*a*a*a, a*a*a*a*a*a*a*a in cls, "\n"); # false
#Print(cz, cz in cls, "\n"); # false

#b := [
#    [1, 0, 0, 0],
#    [0, 1, 0, 0],
#    [0, 0, 0, 1],
#    [0, 0, -1, 0]];;
#Print(cz, cz in ConjugacyClass(Cliff2, b), "\n"); # false
#for b in ConjugacyClass(Cliff2, cz) do Print(b, "\n"); od;

#a1 := [
#    [0, 0, 1, 0],
#    [0, 0, 0, 1],
#    [0, 1, 0, 0],
#    [-1, 0, 0, 0]];;
#Print(a1 in Pauli2, "\n"); # false
#Print(a1 in Cliff2, "\n"); # true
#Print(a1 in cls, "\n");    # true

#Print(Eigenvalues(Cyclotomics, b), "\n"); # fail

#Print(a*a, "\n");
#Print(a*a*a, "\n");
#Print(a*a*a*a, "\n");

# Print(Order(G2));

Hadamard8 := (1/2)*[
    [ 0, 1, 0,-1, 0, 1, 0,-1],
    [ 1, 0, 1, 0, 1, 0, 1, 0],
    [ 0, 1, 0, 1, 0, 1, 0, 1],
    [-1, 0, 1, 0,-1, 0, 1, 0],
    [ 0, 1, 0,-1, 0,-1, 0, 1],
    [ 1, 0, 1, 0,-1, 0,-1, 0],
    [ 0, 1, 0, 1, 0,-1, 0,-1],
    [-1, 0, 1, 0, 1, 0,-1, 0]];

#Print(Hadamard8, "\n");
#Print(Hadamard8*Hadamard8, "\n"); # = identity

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

in_third_level_2 := function(U2)
    # Is U2 in the third level of the clifford hierarchy ?
    local A;
    for g in Pauli2 do
        A := U2*g*Inverse(U2)*Inverse(g);
        if A in Cliff2 then continue; fi;
        return false; # no
    od;
    return true; # yes
end;;

in_fourth_level_2 := function(U2)
    # Is U2 in the fourth level of the clifford hierarchy ?
    local A;
    for g in Pauli2 do
        A := U2*g*Inverse(U2)*Inverse(g);
        if in_third_level_2(U2) then continue; fi;
        return false; # no
    od;
    return true; # yes
end;;

in_fifth_level_2 := function(U2)
    # Is U2 in the fifth level of the clifford hierarchy ?
    local A;
    for g in Pauli2 do
        A := U2*g*Inverse(U2)*Inverse(g);
        if in_fourth_level_2(U2) then continue; fi;
        return false; # no
    od;
    return true; # yes
end;;

get_level_2 := function(U2)
    if U2 in Phase2 then return "0\c"; fi;
    if U2 in Pauli2 then return "1\c"; fi;
    if U2 in Cliff2 then return "2\c"; fi;
    if in_third_level_2(U2) then return "3\c"; fi;
    if in_fourth_level_2(U2) then return "4\c"; fi;
    if in_fifth_level_2(U2) then return "5\c"; fi;
    return ".\c";
end;;

if true then
    Print("# Where do control 1-qubit Pauli gates live ?\n");
    for g in Pauli1 do
        #Print(g, "\n");
        U2{[3,4]}{[3,4]} := g;
        Print(get_level_2(U2));
    od;
    Print("\n");

    Print("# Where do control 1-qubit clifford gates live ?\n");
    for g in Cliff1 do
        #Print(g, "\n");
        U2{[3,4]}{[3,4]} := g;
        Print(get_level_2(U2));
    od;
    Print("\n");
fi;

# ---------------------------------------------------
#
# 3 qubit gates
#

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

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

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

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

wii := KroneckerProduct(wi, i);;

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

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

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

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

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

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

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

Cliff3 := Group(sii, isi, iis, hii, ihi, iih, wii, icz, czi);; # Order 743178240
Pauli3 := Group(wii, xii, ixi, iix, zii, izi, iiz);; # Order 256

# Print(Order(Pauli3), "\n");;

# ca in Cliff3 = false
# ca in third level of clifford hierarchy = true

Cliff3_12 := Group(sii, isi, hii, ihi, wii, czi);;
Cliff3_13 := Group(sii, iis, hii, iih, wii);; # cz on 1&3 ??
Cliff3_23 := Group(isi, iis, ihi, iih, wii, icz);;

#SmCliff3 := Group(sii, isi, iis, hii, ihi, iih, icz, czi);; # Order 743178240

Print("warming up...\n");
Order(Cliff3);; # need this line otherwise membership test eats all memory...!
Print("Ok\n");

#Print("Hadamard8 in Cliff3 ", Hadamard8 in Cliff3, "\n"); # true

in_first_level := function(U3)
    return U3 in Pauli3;
end;;

in_second_level := function(U3)
    if U3 in Cliff3_12 then return true; fi;
    if U3 in Cliff3_13 then return true; fi;
    if U3 in Cliff3_23 then return true; fi;
    if U3 in Cliff3 then return true; fi;
    return false;
end;;

#in_second_level := function(U3)
#    # Is U3 in the second level of the clifford hierarchy ?
#    local A;
#    for g in Pauli3 do
#        A := U3*g*Inverse(U3)*Inverse(g);
#        if in_first_level(A) then continue; fi;
#        return false; # no
#    od;
#    return true; # yes
#end;;

in_third_level := function(U3)
    # Is U3 in the third level of the clifford hierarchy ?
    local A;
    for g in Pauli3 do
        A := U3*g*Inverse(U3)*Inverse(g);
        if in_second_level(A) then continue; fi;
        return false; # no
    od;
    return true; # yes
end;;

#Print(in_third_level(ca1), "\n"); # true
#quit;

in_fourth_level := function(U3)
    # Is U3 in the fourth level of the clifford hierarchy ?
    local A;
    for g in Pauli3 do
        A := U3*g*Inverse(U3)*Inverse(g);
        if in_third_level(A) then continue; fi;
        return false; # no
    od;
    return true; # yes
end;;

in_fifth_level := function(U3)
    # Is U3 in the fifth level of the clifford hierarchy ?
    local A;
    for g in Pauli3 do
        A := U3*g*Inverse(U3)*Inverse(g);
        if in_fourth_level(A) then continue; fi;
        return false; # no
    od;
    return true; # yes
end;;

get_level := function(U3)
    if U3 in Pauli3 then return "1\c"; fi;
    if in_second_level(U3) then return "2\c"; fi;
    if in_third_level(U3) then return "3\c"; fi;
    if in_fourth_level(U3) then return "4\c"; fi;
    if in_first_level(U3) then return "5\c"; fi;
    return ".\c";
end;;

#Print("in_third_level(d):", in_third_level(cc), "\n"); # true
#Print("in_third_level(ca):", in_third_level(ca), "\n"); # true
#Print("in_third_level(cb):", in_third_level(cb), "\n"); # true
#Print("in_third_level(cc):", in_third_level(cc), "\n"); # true
#Print("in_third_level(Tofolli):", in_third_level(Tofolli), "\n"); # true
#Print("in_third_level(Tofolli*ca):", in_third_level(Tofolli*ca), "\n"); # false

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

Print("# Where do control control 1-qubit Pauli gates live ?\n");
for g in Pauli1 do
    #Print(g, "\n");
    U3{[7,8]}{[7,8]} := g;
    Print(get_level(U3));
od;
Print("\n");

Print("# Where do control control 1-qubit clifford gates live ?\n");
for g in Cliff1 do
    #Print(g, "\n");
    U3{[7,8]}{[7,8]} := g;
    Print(get_level(U3));
od;
Print("\n");

Print("# Where do control 2-qubit Pauli gates live ?\n");
for g in Pauli2 do
    #Print(g, "\n");
    U3{[5,6,7,8]}{[5,6,7,8]} := g;
    Print(get_level(U3));
od;
Print("\n");

Print("# Where do control 2-qubit clifford gates live ?\n");
for g in Cliff2 do
    #Print(g, "\n");
    U3{[5,6,7,8]}{[5,6,7,8]} := g;
    Print(get_level(U3));
od;

Print("\n");

Print("Done.\n");

ext: gap

max_stars_repo_name: punkdit/bruhat

max_issues_repo_name: punkdit/bruhat

max_forks_repo_name: punkdit/bruhat


text:
#
# GreenMachine: The GreenMachine: Hyperbolic Groups in GAP
#
# Implementations
#
InstallGlobalFunction( GreenMachine_Example,
function()
	Print( "This is a placeholder function, replace it with your own code.\n" );
end );


text:
multiply := function(a, b)
    return a*b;
end;

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

# This file contains code relevant to arbitrary bit matrices.
# At some point, Marek and I should make one file of this sort that we can both use. (PM)

##    BitMatrixToInts(P) returns bit matrix converted into integers (0, 1).
##
BitMatrixToInts := meth(x)
   return List([1 .. Length(x)], i-> IntVecFFE(x[i]));
end;

##   PrintBitMatrix(P) prints bit matrix P in a readable form.
##
PrintBitMatrix := meth(x)
   PrintMat(BitMatrixToInts(x));
end;

##   NumberBitRep(x, b) returns a b-bit number, formed as a vector over GF(2),
##       corresponding to decimal number x.
##
NumberBitRep := meth(x, bits)
   return Reversed(List([1..bits], i -> [ Mod(QuoInt(x, 2^(i-1)), 2) * GF(2).one ] ));
end;

##   BitsToNumber(x) returns the decimal number associated with bit vector x.
##
BitsToNumber := meth(x)
   return Sum(List([1..Length(x)], i -> 2^(i-1) * IntVecFFE(Reversed(x)[i])[1]));
end;

##    BasisVecToNum(b) returns a number corresponding to given basis
##       vector b.
##
BasisVecToNum := meth(b)
   return Sum(List([1..Length(b)], i -> Cond(IsList(b[i]), b[i][1], b[i]) * (i-1)));
end;

##   BitsToPermMatrix(P) returns the SPL permutation matrix corresponding
##      to bit matrix P.
##
BitsToPermMatrix := meth(P)
   local x, l;
   x := [];
   l := Dimensions(P)[1];

   x := List([1..2^l], i -> BasisVec(2^l, BitsToNumber(P * NumberBitRep(i-1, l))));

   return Transposed(x);
end;

Declare(TransposeVector);

##   PermMatrixToBits(P) returns the bit matrix associated with
##      permutation matrix P.
##
PermMatrixToBits := function(perm)
   local length, res, i, logLength, pm;

   length := Length(perm);
   logLength := Log2Int(length);

   res := [];
   for i in [0..logLength-1] do
      res[logLength-i] :=
         TransposeVector(
	    NumberBitRep(
	       BasisVecToNum(
	          perm * TransposeVector(BasisVec(length, 2^i))
	       ), logLength
	    )
	 );
   od;

   res := TransposedMat(res);

   pm := BitsToPermMatrix(res);

   if (Length(pm) <> Length(perm) or Maximum(Maximum(pm - perm)) <> 0) then
#       Print("ERROR: Given matrix is not a permutation or is not linear on the bits.\n");
       return -1;
   fi;

   return res;
 end;

# ===============================================
# A permutation represented as a linear mapping on the bits
# LinearBits(matrix over GF(2), <SPL Permutation object>)
# ===============================================
Class(LinearBits, BaseMat, SumsBase, rec(
	new := (self, m, p) >>
	   SPL(WithBases(self,
		   rec(
		       dimensions := [2^Length(m), 2^Length(m)],
		       _children  := [m, p]
		   )
	)),

	rSetChild := meth(self, n, what)
	   self._children[n] := what;
	end,

	rChildren  := self >> self._children,
	child      := (self, n) >> self._children[n],
	children   := self >> self._children,

	createCode := self >> self,
        print      := (self,i,is) >> Print(self.name, "(", self._children[1], ")"),

	toAMat     := self >> AMatMat(BitsToPermMatrix(self.child(1))),
	dims       := self >> [2^(Length(self._children[1])), 2^(Length(self._children[1]))],

	sums       := self >> self
));

TL.permBits :=
    meth(self)
        local stride, logStride, size, logSize, res, row, i, j, a, b;
	stride    := self.params[2];
	logStride := Log2Int(stride);
	size      := self.params[1];
	logSize   := Log2Int(size);
	a         := Log2Int(self.params[3]);
	b         := Log2Int(self.params[4]);
	res       := [];

	for i in [1 .. a] do
	   row := Concatenation([ List([1 .. i-1],     i -> 0*GF(2).one),
		                  [ GF(2).one ],
				  List([i+1 .. a],     i -> 0*GF(2).one),
				  List([1 .. logSize], i -> 0*GF(2).one),
				  List([1 .. b],       i -> 0*GF(2).one) ]);

	   Append(res, [row]);
	od;

	for i in [logSize-logStride+1 .. logSize] do
	   row := Concatenation([ List([1 .. a],         i -> 0*GF(2).one),
		                  List([1 .. i-1],       i -> 0*GF(2).one),
				  [ GF(2).one ],
				  List([i+1 .. logSize], i -> 0*GF(2).one),
				  List([1 .. b],         i -> 0*GF(2).one)]);

	   Append(res, [row]);
	od;

	for i in [1 .. logSize-logStride] do
	   row := Concatenation([ List([1 .. a],   i -> 0*GF(2).one),
		                  List([1 .. i-1], i -> 0*GF(2).one),
				  [ GF(2).one ],
				  List([i+1 .. logSize], i -> 0*GF(2).one),
				  List([1 .. b],         i -> 0*GF(2).one)]);

	   Append(res, [row]);
	od;

	for i in [1 .. b] do
	   row := Concatenation([ List([1 .. a],       i -> 0*GF(2).one),
		                  List([1 .. logSize], i -> 0*GF(2).one),
				  List([1 .. i-1],     i -> 0*GF(2).one),
				  [ GF(2).one ],
				  List([i+1 .. b],     i -> 0*GF(2).one)]);

	   Append(res, [row]);
	od;

	return res;
end;

DR.permBits := meth(self)
       local k,  n, res;
       k := LogInt(self.params[2], 2);
       n := LogInt(self.domain(), 2);
       return MatSPL(Tensor(J(n/k), I(k))) * GF(2).one;
end;


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

# dummy to restart doc string
_dummy := 0;

#F ExhaustiveSearch(spl, opts)
#F

ExhaustiveSearch := function(spl, opts)
	local idxList, resList, bestRec;

	if not IsSPL(spl) then
		Error("invalid SPL");
	fi;

	idxList := [1 .. NofRuleTrees(spl, opts)];

	resList := TimeRuleTrees(spl, opts, idxList);

	bestRec := BestTimedRuleTree(resList);
	bestRec.ruletree := RuleTreeN(spl, bestRec.index, opts);

	return bestRec;
end;


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

DimLength := (l) -> When(IsList(l), Length(l), 1);

Class(CrossBase, BaseOperation, rec(
  abbrevs := [ arg -> [Flat(arg)] ],

   new := meth(self, L)
        if Length(L)=1 then return L[1]; fi;
        return SPL(WithBases(self, rec(_children:=L, dimensions := TransposedMat(List(L, e -> e.dims())))));
   end,

   area := self >> Sum(self.children(), x->x.area()),

   normalizedArithCost := self >> Sum(self.children(), x -> x.normalizedArithCost()),

));

Declare(Cross);

Class(Cross, CrossBase, rec(
   codeletName:="C",

   isPermutation := self >> false,

   dmn:=self>>let(li:=[],Flat(List(self._children, x-> x.dmn()))),

   rng:=self>>let(li:=[],Flat(List(self._children, x-> x.rng()))),

   sums:= self>>ObjId(self)(List(self.rChildren(),x->x.sums())),

   transpose := self >> self,

   isIdentity := self >> ForAll(self._children, IsIdentitySPL),

   isBlockTransitive := true,
));

Declare(fCross);
Class(fCross, FuncClassOper, rec(
   _perm := true,

   isFunction := true,
   skipOneChild := true,

   range := (self) >> StripList(List(self.rChildren(), (x) -> x.range())),

   advrange := (self) >> List(self.rChildren(), (x) -> x.advrange()[1]),

   domain := (self) >> StripList(List(self.rChildren(), (x) -> x.domain())),

   advdomain := (self) >> List(self.rChildren(), (x) -> x.advdomain()[1]),

   toSpl := (self, allinds, cols) >>
       fCross(List([1..Length(self.rChildren())],
               i -> self.rChildren()[i].toSpl(allinds,cols[i]))),

   downRank := (self, loopid, ind) >>
       fCross(List(self.rChildren(),
           c -> c.downRank(loopid, ind))),

   downRankFull := (self, allinds) >>
       fCross(List(self.rChildren(),
           c -> c.downRankFull(allinds))),

   isIdentity := self >> ForAll(self._children, IsIdentity),

   transpose  := self >> ObjId(self)(List(self._children, e -> e.transpose())),
));

Class(ExplicitGath, Gath);
Class(ExplicitCross, Cross);

Class(ExplicitOLRuleset, RuleSet);
RewriteRules(ExplicitOLRuleset, rec(
   Gath_f2DTrExpR := Rule([Gath, [@(2,fCompose), @(3,f2DTrExplicitR), ...]],
      e-> Gath(fCompose(Drop(@(2).val.rChildren(),1)))*ExplicitGath(ApplyFunc(f2DTr,@(3).val.rChildren()))),

   Gath_f2DTrExpLR := Rule([Gath, [@(2,fCompose), @(3,f2DTrExplicitL), @(4,f2DTrExplicitR), ...]],
      e-> Gath(fCompose([ApplyFunc(f2DTr,@(4).val.rChildren())]::Drop(@(2).val.rChildren(),2)))
      *ExplicitGath(ApplyFunc(f2DTr,@(3).val.rChildren()))),

   Cross_Compose_ExplicitGath_1 := Rule([Cross, [@(1,Compose), ..., @(2,ExplicitGath)], @(3) ],
       e ->
       Compose(Cross(Compose(DropLast(@(1).val.rChildren(),1)), @(3).val),
       ExplicitCross(ApplyFunc(Gath, @(2).val.rChildren()), ApplyFunc(2DI, @(3).val.advdims()[2][1])))),

   Cross_Compose_ExplicitGath_2 := Rule([Cross, @(3), [@(1,Compose), ..., @(2,ExplicitGath)] ],
       e ->
       Compose(Cross(@(3).val, Compose(DropLast(@(1).val.rChildren(),1))),
       ExplicitCross(ApplyFunc(2DI, @(3).val.advdims()[2][1]), ApplyFunc(Gath, @(2).val.rChildren())))),
));

Class(FinalOLRuleset, RuleSet);
RewriteRules(FinalOLRuleset, rec(
   ExplicitCross := Rule(@(1,ExplicitCross),
       e -> ApplyFunc(Cross, @(1).val.rChildren()))
));


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

Class(SMP_Unparser,      SMP_UnparseMixin, CUnparserProg);
Class(SMP_MacroUnparser, SMP_UnparseMixin, CMacroUnparserProg);

Class(OpenMP_Unparser,      OpenMP_UnparseMixin, CUnparserProg);
Class(OpenMP_MacroUnparser, OpenMP_UnparseMixin, CMacroUnparserProg);

# suggested values: bufIters=64 (16 for older machines), maxRank=1 (larger value increases search space)
# Example: opts := InitGTLibgen(64, 1)
#
InitGTLibgen := function(bufIters, maxRank, useComplex)
    local opts;
    LibgenHardcodeStrides();

    opts := CopyFields(InitLibgen(When(useComplex, CplxLibgenDefaults, LibgenDefaults)),
        rec(
            useDeref := true,
            breakdownRules := rec(
                GT  := [ CopyFields(GT_Base, rec(maxSize := 32)),
                         CopyFields(GT_BufReshape, rec(bufIters := bufIters)),
                         CopyFields(GT_DFT_CT, rec(minSize := 33, maxRank := maxRank)),
                         GT_NthLoop, GT_Par ],
                DFT := [ CopyFields(DFT_CT, rec(maxSize:=32)),
                         CopyFields(DFT_GT_CT, rec(minSize:=32)),
                         DFT_Base ],
                InfoNt := [Info_Base])
        ));
    opts.formulaStrategies.preRC := [ HfuncSumsRules ];
    return opts;
end;

InitSMPGTLibgen := function(bufIters, maxRank, useComplex, useOpenMP)
    local opts;
    opts := CopyFields(InitGTLibgen(bufIters, maxRank, useComplex), rec(
            unparser := Cond(
                useComplex and useOpenMP,         OpenMP_MacroUnparser,
                useComplex and not useOpenMP,        SMP_MacroUnparser,
                not useComplex and useOpenMP,     OpenMP_Unparser,
                not useComplex and not useOpenMP,    SMP_Unparser)));

    opts.formulaStrategies.sigmaSpl := [ MergedRuleSet(StandardSumsRules,RulesSMP), HfuncSumsRules ];
    opts.formulaStrategies.rc := opts.formulaStrategies.sigmaSpl;

    if not useOpenMP then
        opts.subParams := [var("num_threads", TInt), var("tid", TInt)];
        opts.profile := When(LocalConfig.osinfo.isWindows(),
            LocalConfig.cpuinfo.profile.threads(),
            profiler.default_profiles.linux_x86_threads
        );
    fi;
    return opts;
end;


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

Class(APar, ASingletonTag);
Class(AVec, ASingletonTag);

IsParPar := P -> P[3] = APar and P[4] = APar;
IsVecVec := P -> P[3] = AVec and P[4] = AVec;
IsVecPar := P -> P[3] = AVec and P[4] = APar;
IsParVec := P -> P[3] = APar and P[4] = AVec;
IsParXXX := P -> P[3] = APar;
IsVecXXX := P -> P[3] = AVec;
IsXXXPar := P -> P[4] = APar;
IsXXXVec := P -> P[4] = AVec;

BaseOperation.doNotMeasure := true;
PermClass.doNotMeasure := true;

###################################################################
#F tSPL - base class for all tSPL constructs
#F *OBSOLETE*. Use Tagged_tSPL
#F

#D Class(tSPL, NonTerminal, rec(
#D    isTSPL := true
#D ));

###################################################################
#F Tagged_tSPL - base class for all tSPL constructs with tag support
#F
#F *!!* Should be used for all new constructs instead of tSPL
#F
Class(Tagged_tSPL, TaggedNonTerminal, rec(
    isTSPL := true
));

###################################################################
#F Tagged_tSPL_Container - base class for a generic tSPL container
#F
Class(Tagged_tSPL_Container, Tagged_tSPL, rec(
    abbrevs :=  [ nt -> Checked(IsSPL(nt), [nt]) ],
    dims := self >> self.params[1].dims(),
    advdims := self >> self.params[1].advdims(),
    terminate := self >> self.params[1].terminate(),
    transpose := self >> ObjId(self)(self.params[1].transpose())
                             .withTags(self.getTags()),
    isReal := self >> self.params[1].isReal(),
    isInplace := self >> self.params[1].isInplace(),
    normalizedArithCost := self >> self.params[1].normalizedArithCost(),
    doNotMeasure := true,
));

IsTSPL := o -> IsRec(o) and IsBound(o.isTSPL) and o.isTSPL;

Declare(TGrp, TTensor, TTensorI, TTensorInd, TDirectSum, TInplace);

#F TGrp(<s>)
#F  tSPL grouping container (i.e. parenthesis)
#F
Class(TGrp, Tagged_tSPL_Container, rec(
    abbrevs :=  [ s -> Checked(IsSPL(s), [s]) ],
    HashId := self >> let(
	p := self.params[1],
	h := When(IsBound(p.HashId), p.HashId(), p),
        [ TGrp, h ] :: When(IsBound(self.tags), self.tags, [])
    )
));

#F TTensor(<A>, <B>)
#F    Nonterminal for Tensor(A, B)
#F
Class(TTensor, Tagged_tSPL_Container, rec(
    abbrevs :=  [ (A, B) -> Checked(IsNonTerminal(A), IsNonTerminal(B), [A,B]) ],

    dims := self >> let(
	a := self.params[1].dims(), b := self.params[2].dims(),
	[ a[1]*b[1], a[2]*b[2] ] ),

    terminate := self >> Tensor(self.params[1], self.params[2]),

    transpose := self >>
        TTensor(self.params[1].transpose(), self.params[2].transpose())
           .withTags(self.getTags()),

    isReal := self >> self.params[1].isReal() and self.params[2].isReal(),

    normalizedArithCost := self >>
        self.params[1].normalizedArithCost() * Rows(self.params[2]) +
	self.params[2].normalizedArithCost() * Cols(self.params[1]),

    HashId := self >> let(
	h := List(self.params, i -> When(IsBound(i.HashId), i.HashId(), i)),
        When(IsBound(self.tags), Concatenation(h, self.tags), h))
));

#F TDirectSum(<A>, <B>)
#F    Nonterminal for DirectSum(A, B)
#F
Class(TDirectSum, Tagged_tSPL_Container, rec(
    abbrevs :=  [ (A, B) -> Checked(IsSPL(A), IsSPL(B), [ A, B ]) ],
    dims := self >> let(
	a:=self.params[1].dims(), b:=self.params[2].dims(),
	[ a[1]+b[1], a[2]+b[2] ]),

    terminate := self >> DirectSum(self.params[1], self.params[2]),

    transpose := self >>
        TDirectSum(self.params[1].transpose(), self.params[2].transpose())
	    .withTags(self.getTags()),

    isReal := self >> self.params[1].isReal() and self.params[2].isReal(),

    doNotMeasure := true,

    HashId := self >> let(
	h := List(self.params, i->When(IsBound(i.HashId), i.HashId(), i)),
        When(IsBound(self.tags), Concatenation(h, self.tags), h))
));

#F TTensorI(<nt>, <s>, <l>, <r>)
#F
#F tSPL nonterminal for A x I, I x A, (A x I)L, (I x A)L
#F
#F Example: TTensorI(F(2), 32, APar, APar)
#F
Class(TTensorI, Tagged_tSPL_Container, rec(
    abbrevs :=  [ (nt, s, l, r) -> Checked(
        IsSPL(nt), IsPosInt(s), l in [APar, AVec], r in [APar, AVec],
	[nt, s, l, r])
    ],

    dims := self >> self.params[1].dims()*self.params[2],

    SPLtSPL := (self, nt, P) >> let(
	A := nt[1], n := P[2], l := P[3], r := P[4],
        Cond(IsParPar(P), Tensor(I(n), A),
             IsVecVec(P), Tensor(A, I(n)),
             IsParVec(P), Tensor(I(n), A) * L(A.dims()[2]*n, n),
             IsVecPar(P), let(m:=A.dims()[2], Tensor(A, I(n)) * L(m*n, m))
    )),

    terminate := self >> let(
	P := self.params, A:=P[1].terminate(), n:= P[2], l:=P[3],  r:=P[4],
        Cond(IsParPar(P), Tensor(I(n), A),
             IsVecVec(P), Tensor(A, I(n)),
             IsParVec(P), Tensor(I(n), A) * L(A.dims()[2]*n, n),
             IsVecPar(P), let(m:=A.dims()[2], Tensor(A, I(n)) * L(m*n, m)))),

    transpose := self >> let(p := self.params,
        TTensorI(p[1].transpose(), p[2], p[4], p[3]).withTags(self.getTags())),

    isReal := self >> self.params[1].isReal(),

    normalizedArithCost := self >>
        self.params[1].normalizedArithCost() * self.params[2],

    doNotMeasure := true,

    HashId := self >> let(
	p := self.params,
	h := When(IsBound(p[1].HashId), p[1].HashId(), p[1]),
        [h, p[2], p[3], p[4]] :: When(IsBound(self.tags), self.tags, [])
    ),
));

Class(TTensorInd, Tagged_tSPL_Container, rec(
    abbrevs :=  [ (nt, idx, l, r) -> Checked(
        IsSPL(nt), IsVar(idx), l in [APar, AVec], r in [APar, AVec],
	[nt, idx, l, r])
    ],

    dims := self >> self.params[1].dims() * self.params[2].range,

    SPLtSPL := (self, nt, P) >> let(
	A := nt[1], idx := P[2], m := A.dims()[2], n := idx.range,
        dsA := IDirSum(idx, A),
        Cond(IsParPar(P), dsA,
             IsVecVec(P), L(m*n, m) * dsA * L(m*n, n),
             IsParVec(P), dsA * L(m*n, n),
             IsVecPar(P), L(m*n, m) * dsA)
    ),

    terminate := self >> let(
	P := self.params, A := P[1], idx:= P[2],
        lstA := List([0..idx.range-1],
	    i -> SubstVars(Copy(A), rec((idx.id) := V(i))).terminate()),
        d := A.dims(),
	n := idx.range,
        dsA := DirectSum(lstA),
        Cond(IsParPar(P), dsA,
             IsVecVec(P), L(d[1]*n, d[1]) * dsA * L(d[2]*n, n),
             IsParVec(P), dsA * L(m*n, n),
             IsVecPar(P), L(m*n, m) * dsA)
    ),

    transpose := self >> TTensorInd(
       self.params[1].transpose(), self.params[2], self.params[4], self.params[3])
          .withTags(self.getTags()),

    isReal := self >> self.params[1].isReal(),

    normalizedArithCost := self >>
        self.params[1].normalizedArithCost() * self.params[2].range,

#    doNotMeasure := true,

    HashId := self >> let(
	p := self.params,
	h := When(IsBound(p[1].HashId), p[1].HashId(), p[1]),
        [h, p[2].range, p[3], p[4]] :: When(IsBound(self.tags), self.tags, [])
    )
));

Declare(TL);

#F  tSPL I x L x I
#F  TL(m, n, k, j)
#F  I(k), L(m, n), I(j)
Class(TL, Tagged_tSPL_Container, rec(
    abbrevs :=  [ (size, stride) -> Checked(ForAll([size, stride], IsPosIntSym), [size, stride, 1, 1]),
                  (size, stride, left, right) -> Checked(ForAll([size, stride, left, right], IsPosIntSym), [size, stride, left, right]) ],
    __call__ := arg >> let(self := arg[1], args := Drop(arg, 1),
        Cond(args=[1,1,1,1], I(1),
             ApplyFunc(Inherited,args))),

    dims := self >> Replicate(2, self.params[1]*self.params[3]*self.params[4]),
    terminate := self >> Tensor(I(self.params[3]), L(self.params[1], self.params[2]), I(self.params[4])),
    transpose := self >> TL(self.params[1], self.params[1]/self.params[2], self.params[3], self.params[4]).withTags(self.getTags()),
    isReal := self >> true,
    noCodelet := true,
    # need that for DPBench!
    doNotMeasure := true,
    transposeSymmetric := False,
    isSymmetric := self >> self.params[2]^2 = self.params[1],
    hashAs := self >> ObjId(self)(self.params[1], self.params[2], self.params[3], self.params[4]).withTags(self.getTags()),
));

#F  tSPL Bit Reversal
Class(TBR, Tagged_tSPL_Container, rec(
    abbrevs :=  [ n -> Checked(IsInt(n), [n])],
    dims := self >> [self.params[1], self.params[1]],
    terminate := self >> BR(self.params[1]),
    transpose := self >> self,
    isReal := True,
    doNotMeasure := true,
));

Declare(TRDiag);

#F TRDiag(<nrows>, <ncols>, <diagfunc>)
#F
#F TRDiag is a rectangular tall matrix, which is a diagonal
#F zero-padded on the bottom.
#F
#F To obtain wide matrix (diagonal zero-padded on the right),
#F use TRDiag(...).transpose()
#F
#F Example:
#F  spiral> PrintMat(MatSPL(  TRDiag(6, 3, FList(TInt, [1,2,3])) ));
#F  [ [ 1,  ,   ],
#F    [  , 2,   ],
#F    [  ,  , 3 ],
#F    [  ,  ,   ],
#F    [  ,  ,   ],
#F    [  ,  ,   ] ]
#F
Class(TRDiag, Tagged_tSPL, rec(
    abbrevs :=  [ (m,n,diag) -> Checked(IsPosIntSym(m), IsPosIntSym(n),
	    AnySyms(m, n) or (m>=n),
	    IsFunction(diag) or IsFuncExp(diag),
	    [m, n, diag])],

    dims :=  self >> let( d := [ self.params[1], self.params[2] ], Cond(self.transposed, Reversed(d), d)),

    # tells autolib's ConstrainPSums that params[2] must be = params[3].domain()
    # NOTE: get rid of params[2]
    paramConstraints := self >> [[self.params[2], self.params[3].domain()]],

    terminate := self >> let(
        res := RI(self.params[1], self.params[2]) * Diag(self.params[3]),
	When(self.transposed, res.transpose(), res)),

    conjTranspose := self >> CopyFields(self,
	rec(params := [self.params[1], self.params[2], FConj(self.params[3])]))
            .transpose(),

    isReal := self >> not IsComplexT(self.params[3].range()),

    doNotMeasure := true,

    hashAs := self >> let(
	t := ObjId(self)(self.params[1], self.params[2],
	                 fUnk(self.params[3].range(), self.params[3].domain())),
	tt := t.withTags(self.getTags()),
	When(self.transposed, tt.transpose(), tt))
));

#F  tSPL Digit Reversal
Class(TDR, Tagged_tSPL_Container, rec(
    abbrevs :=  [ (n, r) -> Checked(IsInt(n) and  IsInt(r) and n=r^LogInt(n,r), [n, r])],
    dims := self >> [self.params[1], self.params[1]],
    terminate := self >> DR(self.params[1], self.params[2]),
    transpose := self >> self,
    isReal := True,
    doNotMeasure := true,
));

Declare(TReflL);
#F  tSPL Refl(L(mn,m))
Class(TReflL, Tagged_tSPL_Container, rec(
    abbrevs :=  [ (mn, m) -> Checked(IsInt(mn) and  IsInt(m), [mn, m]) ],
    dims := self >> self.terminate().dims(),
    terminate := self >> let(mn := self.params[1], m:=self.params[2], N := mn/2,
                             Refl(N, 2*N-1, N, L(mn, m))),
    transpose := self >> TReflL(self.params[1], self.params[1]/self.params[2]).withTags(self.getTags()),
    isReal := True,
    doNotMeasure := true,
));

Declare(TRC, TCompose, TDiag, TICompose);

#F   tSPL RC(.)
Class(TRC, Tagged_tSPL_Container, rec(
    _short_print := true,
    abbrevs :=  [ (A) -> Checked(IsNonTerminal(A) or IsSPL(A), [A]) ],
    dims := self >> 2*self.params[1].dims(),
    terminate := self >> Mat(MatSPL(RC(self.params[1]))),

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

    conjTranspose := self >> self.transpose(),

    isReal := self >> true,

    # Do not use doNotMeasure, this will prevent TRC_By_Def from ever being found!
    doNotMeasure := false,
    normalizedArithCost := self >> self.params[1].normalizedArithCost(),

    HashId := self >> let(
	h := [ When(IsBound(self.params[1].HashId), self.params[1].HashId(),
		    self.params[1]) ],
        When(IsBound(self.tags), Concatenation(h, self.tags), h))
));

Class(TMat, Tagged_tSPL, rec(
    abbrevs :=  [ m -> Checked(IsMat(m), [m]) ],
    dims := self >> DimensionsMat(self.params[1]),
    terminate := self >> Mat(self.params[1]),
    transpose := self >> ObjId(self)(TransposedMat(self.params[1])).withTags(self.getTags()),
    isReal := self >> IsRealMat(self.params[1]),
    isInplace := self >> false,
    normalizedArithCost := self >> 2 * Rows(self.params[1]) * Cols(self.params[1]),
    doNotMeasure := false,
));

#F   tSPL Diag(d)
Class(TDiag, Tagged_tSPL_Container, rec(
    abbrevs :=  [ (D) -> [D] ],
    dims := self >> let(n:=self.params[1].domain(), [n,n]),
    terminate := self >> Diag(self.params[1]),
    transpose := self >> Copy(self),
    conjTranspose := self >> ObjId(self)(FConj(self.params[1])).withTags(self.getTags()),

    isReal := self >> Diag(self.params[1]).isReal(),

    doNotMeasure := true,
    noCodelet := true,
    hashAs := self >> ObjId(self)(FUnk(self.params[1].domain())).withTags(self.getTags()),
    #NOTE: depends on wether the whole transform is complex...
    normalizedArithCost := self >> When(self.params[1].isReal(), 2*Rows(self), 6*Rows(self))
));

#F   tSPL RCDiag(d)
Class(TRCDiag, Tagged_tSPL_Container, rec(
    abbrevs :=  [ (D) -> [D] ],
    dims := self >> let(n:=self.params[1].domain(), [n,n]),
    terminate := self >> RCDiag(self.params[1]),
    # NOTE: transpose on real == conjugate; transpose on complex == id. What to use?
#    transpose := self >> let(d := self.params[1].var.value, t := d.t.t, size := d.t.size,
#        newd := List(d.v, i-> t.value(Complex(ReComplex(i.v), -ImComplex(i.v)))),
#        FData(newd)),
    transpose := self >> CopyFields(self, rec(transposed := not self.transposed)),
    isReal := True,
    doNotMeasure := true,
    noCodelet := true,
    hashAs := self >> ObjId(self)(FUnk(self.params[1].domain())).withTags(self.getTags()),
    normalizedArithCost := self >> 3*Rows(self)
));

#F   tSPL Id(f)
Class(TId, TDiag);

Declare(TCond);
#F   tSPL Cond
Class(TCond, Tagged_tSPL_Container, rec(
    abbrevs   := [ (f, a, b) -> Checked(IsBound(f.ev), [f,a,b]) ],
    dims      := self >> self.params[2].dims(),
    terminate := self >> COND(self.params[1], self.params[2].terminate(), self.params[3].terminate()),
    transpose := self >> TCond(self.params[1], self.params[2].transpose(), self.params[3].transpose()), # maybe
));

Declare(TRaderMid);
#F   tSPL RaderMid(N, k, root, transp)
Class(TRaderMid, Tagged_tSPL_Container, rec(
    # -- Rader infrastructure --
    raderDiag := (N, k, root) -> SubList(ComplexFFT(List([0..N-2], x->ComplexW(N, k*root^x mod N))), 2) / (N-1),
    raderMid := (self, N, k, root) >>
                 DirectSum(Mat([[1, 1], [1, -1/(N-1)]]),
                   Diag(FData(self.raderDiag(N, k, root)))),
    # -- tSPL infrastructure --
    abbrevs :=  [ (N, k, root) -> [N, k, root, false],
                  (N, k, root, transp) -> [N, k, root, transp] ],
    dims := self >> let(N:=self.params[1], [N,N]),
    terminate := self >> let(spl := self.raderMid(self.params[1], self.params[2], self.params[3]),
        When(self.params[4], spl.transpose(), spl)),
    transpose := self >> TRaderMid(self.params[1], self.params[2], self.params[3], not self.params[4]).withTags(self.getTags()),
    isReal := False,
    doNotMeasure := true,
    noCodelet := true,
));

#F   tSPL Compose([A,B,...])
Class(TCompose, Tagged_tSPL_Container, rec(
    abbrevs :=  [ (l) -> Checked(IsList(l), [l]) ],
    dims := self >> [self.params[1][1].dims()[1], self.params[1][Length(self.params[1])].dims()[2]],
    terminate := self >> Compose(List(self.params[1], i->i.terminate())),
    transpose := self >> TCompose(Reversed(List(self.params[1], i->i.transpose()))).withTags(self.getTags()),
    isReal := self >> ForAll(self.params[1], i->i.isReal()),
    doNotMeasure := true,
    normalizedArithCost := self >> Sum(List(self.params[1], i->i.normalizedArithCost())),

    HashId := self >> let(h := List(self.params[1], i->When(IsBound(i.HashId), i.HashId(), i)),
        When(IsBound(self.tags), Concatenation(h, self.tags), h))
));

#F   tSPL ICompose(<var>, <domain>, <nonterm>, <attrib>)
Class(TICompose, Tagged_tSPL_Container, rec(
    abbrevs :=  [ (var,domain,l) -> [var, domain, l] ],
    dims := self >> self.params[3].dims(),
    terminate := self >> ICompose(self.params[1], self.params[2], self.params[3]).unroll(),

    unroll := self >> let(
        var := self.params[1],
        domain := self.params[2],
        expr := self.params[3],
        tag := self.getTags(),
        TCompose(
            List([ 0 .. domain - 1 ],
                index_value -> SubstBottomUp(
                    Copy(expr), var,
                    e -> V(index_value))))),

    isReal := self >> ICompose(self.params[1], self.params[2], self.params[3]).unroll().isReal(),

    transpose := self >> let(
        var := self.params[1],
        dom := self.params[2],
        ch := self.params[3],
        newvar := var.clone(),
        tags := self.getTags(),
        TICompose(newvar, dom,
            SubstVars(Copy(ch.transpose()), tab((var.id) := dom-1-newvar))).withTags(tags)
    )
));

Declare(TGath, TScat, TPrm, TConj, TS);
#F   tSPL Gath(d)
Class(TGath, Tagged_tSPL_Container, rec(
    abbrevs :=  [ (tspl) -> [tspl] ],
    dims := self >> [self.params[1].domain(), self.params[1].range()],
    terminate := self >> Gath(self.params[1]),
    transpose := self >> TScat(self.params[1]).withTags(self.getTags()),
    isReal := True,
    doNotMeasure := true,
    noCodelet := true,
    normalizedArithCost := self >> 0
));

Class(TScat, Tagged_tSPL_Container, rec(
    abbrevs :=  [ (tspl) -> [tspl] ],
    dims := self >> [self.params[1].range(),self.params[1].domain()],
    terminate := self >> Scat(self.params[1]),
    transpose := self >> TGath(self.params[1]).withTags(self.getTags()),
    isReal := True,
    doNotMeasure := true,
    noCodelet := true,
    normalizedArithCost := self >> 0
));

Class(TPrm, Tagged_tSPL_Container, rec(
    abbrevs :=  [ (tspl) -> [tspl] ],
    dims := self >> Cond(IsBound(self.params[1].domain), [self.params[1].domain(), self.params[1].range()], self.params[1].dims()),
    terminate := self >> Prm(self.params[1]),
    transpose := self >> TPrm(self.params[1].transpose()).withTags(self.getTags()),
    isReal := True,
    doNotMeasure := true,
    noCodelet := true,
    normalizedArithCost := self >> 0
));

Class(TConj, Tagged_tSPL_Container, rec(
    abbrevs :=  [ (tspl, prm) -> [tspl,prm.transpose(),prm],
                  (tspl, prml, prmr) -> [tspl,prml,prmr] ],
    dims := self >> [self.params[2].domain(), self.params[3].range()],
    terminate := self >> ConjLR(self.params[1].terminate(), self.params[2], self.params[3]),
    transpose := self >> TConj(self.params[1].transpose(), self.params[3].transpose(), self.params[2].transpose()).withTags(self.getTags()),
    isReal := self >> self.params[1].isReal(),
    doNotMeasure := false,
    noCodelet := true,
    normalizedArithCost := self >> self.params[1].normalizedArithCost()
));

Class(TS, Tagged_tSPL_Container, rec(
    abbrevs :=  [ (N) -> [N, false],
                  (N, t) -> [N, t] ],
    dims := self >> [self.params[1], self.params[1]],
    terminate := self >> When(self.params[2], S(self.params[1]).transpose(), S(self.params[1])),
    transpose := self >> TS(self.params[1], not self.params[2]).withTags(self.getTags()),
    isReal := True,
    doNotMeasure := true,
    noCodelet := true,
    normalizedArithCost := self >> self.params[1] - 1
));

Class(TInplace, Tagged_tSPL_Container, rec(
    abbrevs :=  [ (tspl) -> [tspl] ],
    dims := self >> self.params[1].dims(),
    terminate := self >> self.params[1].terminate(),
    transpose := self >> TInplace(self.params[1].transpose()).withTags(self.getTags()),
    isReal := self >> self.params[1].isReal(),
    isInplace := True,
    doNotMeasure := true,
    noCodelet := true,
    normalizedArithCost := self >> self.params[1].normalizedArithCost()
));

InterleavedComplexT := t -> TRC(t);
SplitComplexT := t -> let(r:= Rows(t), c:=Cols(t), TConj(TRC(t), L(2*r, 2), L(2*c, c)));
ComplexT := (t,opts) -> When((not IsBound(opts.interleavedComplex)) or (IsBound(opts.interleavedComplex) and opts.interleavedComplex), InterleavedComplexT(t), SplitComplexT(t));

InterleavedSplitComplexT := t -> let(r:= Rows(t), c:=Cols(t), TConj(TRC(t), fId(2*r), L(2*c, c)));


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

# !!!! what's missing: CRT of 2 primes perm pullout

RewriteRules(RulesFuncSimp, rec(
 # ===================================================================
 # OddStride
 # ===================================================================
 OS_Id := Rule([OS, @N, 1], e->fId(@N.val)),

 OS_OS := ARule(fCompose, [[OS, @N, @r], [OS, @, @s]],
     e -> [ OS(@N.val, @r.val * @s.val mod @N.val) ]),

 OS_RR_PhaseShift := ARule( fCompose, [[OS, @N, @r], [RR, @, @phi, @g]],
     e -> [ RR(@N.val, @phi.val * @r.val, @g.val) ]),

 OS_RM_PhaseShift := ARule( fCompose, [[OS, @N, @r], [RM, @, @phi, @g]],
     e -> [ RM(@N.val, @phi.val * @r.val, @g.val) ]),

 # ===================================================================
 # Cyclic shift (Z) propagation
 # Z o (fBase X f)  ->  (Z o f) X (Z o fBase)
 Z_Propagate := ARule(fCompose, [[Z, @, @z], [fTensor, @F, [fBase, @m, @j]]],
     (e, cx) -> [ fTensor(
       fCompose(Z(range(@F.val), memo(cx, "g", idiv(@j.val+@z.val, @m.val))), @F.val),
       fCompose(Z(@m.val, @z.val), fBase(@m.val, @j.val))) ]),

 # ===================================================================
 # CRT / gammaTensor
 # ===================================================================
#   NOTE: This rule triggers problems when fired too early, as gammatensor cannot be transposed -- what to do??
# CRT_toGammaTensor := Rule(CRT, e->e.toGammaTensor()),

 GammaTensorAssoc := ARule(gammaTensor, [@(1,gammaTensor)], e->@(1).val.children()),

 # gammaTensor o fTensor -> gammaTensor( ... o ..., ... o ... , ...)
 GammaTensorMerge := ARule(fCompose,
      [ @(1,gammaTensor), @(2,fTensor,e->compat_tensor_chains(
                         @(1).val.children(), e.children(), domain, range)) ],
 e -> [ gammaTensor(merge_tensor_chains(
       @(1).val.children(), @(2).val.children(), fCompose, x->x, x->x, domain, range))]),

 ### fTensor(..., Y, ...) o X -> fTensor(..., Y o X, ...)
 GammaTensorComposeMergeRight := ARule(fCompose,
      [ @(1,gammaTensor), @(2).cond(e->compat_tensor_chains(
                         @(1).val.children(), [e], domain, range)) ],
 e -> [ gammaTensor(merge_tensor_chains(
       @(1).val.children(), [@(2).val], fCompose, x->x, x->x, domain, range)) ] ),

 ### X o fTensor(..., Y, ...) -> fTensor(..., X o Y, ...)
 GammaTensorComposeMergeLeft := ARule(fCompose,
      [ @(1), @(2,gammaTensor,e->compat_tensor_chains(
                         [@(1).val], e.children(), domain, range)) ],
 e -> [ gammaTensor(merge_tensor_chains(
       [@(1).val], @(2).val.children(), fCompose, x->x, x->x, domain, range)) ] ),

 # ===================================================================
 # Rader permutations
 # ===================================================================
 # NB: RR = fStack(1,RM), RR plugged into gammaTensor, eg. F_15 -> F_3, F_5

 # RR o fAdd(N,1,0) -> (0)_N
 RR_1stPoint_fAdd := ARule(fCompose, [[RR, @N, @phi, @g], [fAdd, @, 1, 0]],
     e -> [ H(@N.val, 1, 0, 1) ]),
 RR_1stPoint_H := ARule(fCompose, [[RR, @N, @phi, @g], [H, @, 1, 0, @]],
     e -> [ H(@N.val, 1, 0, 1) ]),

 # RR o fAdd(N,N-1,1) -> RM
 RR_toRM_fAdd := ARule(fCompose,  [[RR, @N, @phi, @g], [fAdd, @, @(1).cond(e->e=@N.val-1), 1]],
     e -> [ RM(@N.val, @N.val-1, @phi.val, @g.val) ]),
 RR_toRM_H := ARule(fCompose,  [[RR, @N, @phi, @g], [H, @, @(1).cond(e->e=@N.val-1), 1, 1]],
     e -> [ RM(@N.val, @N.val-1, @phi.val, @g.val) ])
));

Class(RulesRaderTensor, RuleSet);

RewriteRules(RulesRaderTensor, rec(
 # RM o fTensor -> RM
 RM_Tensor1 := ARule(fCompose, [ [RM,@N,@n,@phi,@g], [fTensor, @(1,fId), [fBase,@m,@j]] ],
     (e, cx) -> [ RM(@N.val, domain(@1.val),
           memo(cx, "g", powmod(@phi.val,@g.val,@j.val,@N.val)),
           powmod(1,@g.val,@m.val,@N.val).eval()) ] ),

 RM_Tensor2 := ARule(fCompose, [ [RM,@N,@n,@phi,@g], [fTensor, [fBase,@m,@j], @(1,fId)] ],
     (e, cx) -> [ RM(@N.val, domain(@1.val),
           memo(cx, "g", powmod(@phi.val, powmod(1,@g.val,@m.val,@N.val).eval(), @j.val, @N.val)),
           @g.val) ])
));


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

UnparseQASM := function (spl)
    PrintTo("./qspiralout", spl);
    Exec( "./namespaces/packages/quantum/unparser/bin/unparser ./qspiralout");
end;

max_stars_repo_name: spiral-software/spiral-package-quantum

max_issues_repo_name: spiral-software/spiral-package-quantum

max_forks_repo_name: spiral-software/spiral-package-quantum


text:
RequirePackage("grape");

V := GF(3)^7;
P := Filtered(List(Elements(Subspaces(V, 1)), v -> Elements(v)[2]), u -> u*u = Z(3));;

G := Graph(Group(()), P, function(x,y) return x; end, function(x,y) return x*y = 0*Z(3); end, true);;
GlobalParameters(G);

max_stars_repo_name: jaanos/gap-graphs

max_issues_repo_name: jaanos/gap-graphs

max_forks_repo_name: jaanos/gap-graphs


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

# Radix n kernel for FFTE, targeting ARM SVE

buildSVEKernel_a := function(n, kk, conf, opts)
    local vlen, name, suffix, filesuffix, t, j, k, l, m, tmp1, tmp2, tmp3,
        gi, g, gc, si, ss, s, sc, dft, rts, rt, opcnts, dfts, dftc, i, twt, twf, twl, tws, twc,
        cl, cc, c, clp, useFMA, _FMA, lp;

    name := DFTnameStr(n, kk);
    suffix := "a";
    filesuffix := "sve_";
#    suffix := "_j";
    Print("\n\n// == Generating \"", name, "\" ==================================================\n\n");

    vlen := TInt;
    useFMA := When(IsBound(conf.useFMA), conf.useFMA, true);
    _FMA := When(useFMA, FMA, p -> p);

    # variables
    t := var.fresh_t("t", TInt);
    j := var.fresh_t("j", TInt);
    k := V(0);
    l := var.fresh_t("l", TInt);
    m := V(1);
    lp := var.fresh_t("lp", TPtr(TInt));
#    mp := var.fresh_t("mp", TPtr(TInt));

    # temp arrays
    tmp1 := var.fresh_t("R", TArray(TVectSVE, 2*n));
    tmp2 := var.fresh_t("S", TArray(TVectSVE, 2*n));
    tmp3 := var.fresh_t("T", TArray(TVectSVE, 2*n));

    # gather
    gi := _i -> VGath_L(1, vlen, 2*(k + j*m + _i*l*m));
    g := RulesSums(SumsSPL(VStack(List([0..n-1], gi)), opts));
    gc := Compile(opts.codegen(g, tmp1, X, opts), opts);

    # scatter
    si := _i -> ScatH(1, 1, 2*n*j+_i, 2*n);
    ss := SumsSPL(HStack(List([0..2*n-1], si)), opts);
    s := RulesSums(SubstTopDown(ss, @(1, ScatAcc), e->Scat(@(1).val.func)));
    sc := Compile(opts.codegen(s, Y, tmp3, opts), opts);

    # generate DFT kernel
    dft := RC(DFT(n, kk));
    rts := AllRuleTrees(dft, opts);
    opcnts := List(rts, r -> [ Length(Collect(CodeRuleTree(r, opts), @(1, [add, sub, mul], e->e.t=TReal))), r]);
    rt := Minimum(opcnts)[2];
    dfts := SumsRuleTree(rt, opts);
    dftc := BlockUnroll(opts.codegen(dfts, tmp2, tmp1, opts), opts);

    # twiddles as lookup table
    i := var.fresh_t("i", 2*n);
    twt := var.fresh_t("TW", TPtr(TReal));
    twf := Lambda(i, cond(
        eq(i, 0), 1,
        eq(i, 1), 0,
        sve_gath(nth(twt, 2*((n-1)*j) + 2 * idiv(i, 2) + (imod(i, 2) - 2)).toPtr(TReal), TInt, opts.pg, 2*(n-1))
    ));

    # debug twiddles
    twl := Map(twf.tolist(), RulesStrengthReduce);
    DoForAll([0..Length(twl)-1], _i->Unparse(assign(nth(Y, _i), twl[_i+1]), CUnparser, 0, 1));
    tws := RCDiag(twf);
    twc :=BlockUnroll(unroll_cmd(opts.codegen(tws, tmp3, tmp2, opts)), opts);

    # stitch code fragments together
    cl := func(TVoid, "transform", [Y, X, twt, j, k, l, m, opts.pg],
            decl([tmp1, tmp2, tmp3], chain(gc, dftc, twc, sc)));
    cc := Compile(cl, opts);

    # fixup to push decls inside (?)
    c := func(TVoid, "transform", [Y, X, twt, j, l, opts.pg],
            decl(cc.vars, _FMA(cc.cmd.cmds[1].cmd)));

    # print Radix N FFTE kernel as function
    #PrintCode(name, c, opts);
    #PrintTo(name::".c", PrintCode(name, c, opts));

    ########################################################
    # version with j loop vectorized as m == 1, k == 0
    clp := func(TVoid, "transform", [Y, X, twt, lp],
        decl(c.cmd.vars::c.free()::[opts.pg, l],
            chain(
                assign(l, deref(lp)),
                sve_loopn(j, l,
                    c.cmd.cmd
                )
            )
        )
    );

    # print Radix N FFTE kernel as function
    PrintCode(opts.FortranIze(name::suffix), clp, opts);
    PrintTo(name::filesuffix::suffix::".c", PrintCode(opts.FortranIze(name::suffix), clp, opts));
end;

max_stars_repo_name: spiral-software/spiral-package-ffte

max_issues_repo_name: spiral-software/spiral-package-ffte

max_forks_repo_name: spiral-software/spiral-package-ffte


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

################################################
#   propagate VTensor into products and sums
_VConsRightBlk := [BlockVPerm, VPerm, VDiag_x_I, VDiag, VRCDiag, RCVDiagSplit, RCVDiag, VGath, VGath_u, VGath_sv, VGath_pc, IxVGath_pc, VGath_dup, VPrm_x_I, Conj, ConjL, ConjR, ConjLR, RCVGath_sv, FormatPrm];
_VConsLeftBlk  := [BlockVPerm, VPerm, VDiag_x_I, VDiag, VRCDiag, RCVDiagSplit, RCVDiag, VScat, VScatAcc, VScatAcc_u, VScat_u, VScat_sv, VScat_svAcc, VScat_pc, IxVScat_pc, VScat_pcAcc, VPrm_x_I, Conj, ConjL, ConjR, ConjLR, RCVScat_sv, FormatPrm];

# !!!! do not put VPerm in here - VPerms CANNOT be pulled into ISums!!!!
_VConsDiag := [ VDiag_x_I, VDiag, VRCDiag, RCVDiag, RCVDiagSplit ];

_VConsRightNoDiag := [VGath, VGath_u, VGath_sv, VPrm_x_I];
_VConsLeftNoDiag  := [VScat, VScat_u, VScat_sv, VPrm_x_I, VScatAcc, VScatAcc_u, VScat_svAcc, VScat_pcAcc, VDiag_x_I];

_VConsRight := _VConsRightNoDiag :: _VConsDiag :: [IxVGath_pc ];
_VConsLeft  := _VConsLeftNoDiag :: _VConsDiag :: [IxVScat_pc ];

# things to not pull into Vcontainer:
_VContDontPullIn := [VContainer, Cross, Cvt, TCvt, RC, BB, ISum, SUM, RecursStep];

_PullInLeft  := [RecursStep, Grp, BB, SUM, Buf, ISum, Data, COND, NeedInterleavedComplex];
_PullInRight := _PullInLeft :: [SUMAcc, ISumAcc];

Class(RulesPropagate, RuleSet);
RewriteRules(RulesPropagate, rec(

 VecPullInLeft := ARule( Compose, [ @(1, _VConsLeft), @(2, _PullInLeft :: [NoDiagPullinRight]) ],
  e -> [ CopyFields(@(2).val, rec(
             _children :=  List(@(2).val._children, c -> @(1).val * c),
             dimensions := [Rows(@(1).val), Cols(@(2).val)] )) ]),

 VecPullInRight := ARule( Compose, [ @(1, _PullInRight :: [NoDiagPullinLeft]), @(2, _VConsRight) ],
     e -> [ CopyFields(@(1).val, rec(
                _children := List(@(1).val._children, c -> c * @(2).val),
                dimensions := [Rows(@(1).val), Cols(@(2).val)] )) ]),

 VecPullInLeftNoDiag := ARule( Compose, [ @(1, _VConsLeftNoDiag), @(2, [NoDiagPullinLeft, NoDiagPullin]) ],
  e -> [ CopyFields(@(2).val, rec(
             _children :=  List(@(2).val._children, c -> @(1).val * c),
             dimensions := [Rows(@(1).val), Cols(@(2).val)] )) ]),

 VecPullInRightNoDiag := ARule( Compose, [ @(1, [NoDiagPullinRight, NoDiagPullin]), @(2, _VConsRightNoDiag) ],
     e -> [ CopyFields(@(1).val, rec(
                _children := List(@(1).val._children, c -> c * @(2).val),
                dimensions := [Rows(@(1).val), Cols(@(2).val)] )) ]),

#######################################
    VTensorXXX := Rule([@(1, VTensor), @(2, [SUM, SUMAcc, Compose, BB, Buf, Inplace, RecursStep, ISum, ISumAcc, ISumLS, Data, NeedInterleavedComplex, SMPBarrier, SMPSum])],
    e -> let(s := @(2).val, CopyFields(s,
         rec(_children := List(s.children(), c->VTensor(c, @(1).val.vlen)),
             dimensions := @(1).val.dimensions)))),

    VTensorIndXXX := Rule([@(1, VTensorInd), @(2, [SUM, SUMAcc, Compose, BB, Buf, Inplace, RecursStep, ISum, ISumAcc, ISumLS, Data, NeedInterleavedComplex, SMPBarrier, SMPSum]), @(3)],
    e -> let(s := @(2).val, CopyFields(s,
         rec(_children := List(s.children(), c->VTensorInd(c, @(3).val)),
             dimensions := @(1).val.dimensions)))),

    # Free variables in breakdown rules seems to be broken...
    VTensorInd_VTensor := Rule([@(1, VTensorInd), @(2).cond(e->ObjId(e)<>RTWrap), @(3).cond(e->not(e in @(2).val.free()))],
        e -> VTensor(@(2).val, @(3).val.range)),

#(s) -> SubstTopDownNR(RulesRC(s), @@(1, RCDiag, (x, cx) -> IsBound(cx.VRCLR) and cx.VRCLR <> [  ] or x.element.range() = TComplex), (x) -> x.toloop().sums())

    # NOTE (handle xofs and yofs correctly)
    MergeRS := ARule(Compose, [@(1,RecursStep), @(2,RecursStep)],
    e -> [ RecursStep(0,0,@(1).val.child(1)*@(2).val.child(1)) ]),

    # NOTE (handle xofs and yofs correctly)
    # V * OP(a) -> OP(V*a), including the case V=BlockVPerm
    XXX_VConsBlk := ARule(Compose, [@(1,[RecursStep, Buf, Inplace, BB, NeedInterleavedComplex]), @(2,_VConsRightBlk)],
    e -> [ ObjId(@(1).val)(@(1).val.child(1)*@(2).val) ]),
    XXX_VConsBlk_VCont := ARule(Compose, [@(1,[RecursStep, Buf, Inplace, BB, NeedInterleavedComplex]), [@(2, VContainer), @(3,_VConsRightBlk)]],
    e -> [ ObjId(@(1).val)(@(1).val.child(1)*@(2).val) ]),

    VConsBlk_XXX := ARule(Compose, [@(1,_VConsLeftBlk), @(2, [RecursStep, Buf, Inplace, BB, NeedInterleavedComplex])],
    e -> [ ObjId(@(2).val)(@(1).val*@(2).val.child(1)) ]),
    VConsBlk_VCont_XXX := ARule(Compose, [[@(1, VContainer), @(3,_VConsLeftBlk)], @(2, [RecursStep, Buf, Inplace, BB, NeedInterleavedComplex])],
    e -> [ ObjId(@(2).val)(@(1).val*@(2).val.child(1)) ]),

    VConsBlk_XXX := ARule(Compose, [@(1,_VConsLeftBlk), @(2, [RecursStep, Buf, Inplace, BB, NeedInterleavedComplex])],
    e -> [ ObjId(@(2).val)(@(1).val*@(2).val.child(1)) ]),

    VContainer_XXX := ARule(Compose, [@(1,VContainer), @(2).cond(x -> not (ObjId(x) in _VContDontPullIn))],
        e -> [ ObjId(@(1).val)(@(1).val.child(1)*@(2).val, @(1).val.isa) ]),

    XXX_VContainer := ARule(Compose, [@(1).cond(x -> not (ObjId(x) in _VContDontPullIn)), @(2, VContainer)],
        e -> [ ObjId(@(2).val)(@(1).val*@(2).val.child(1), @(2).val.isa) ]),

    Compose_VContainers := ARule(Compose, [ @(1, VContainer), @(2, VContainer, x -> x.isa=@(1).val.isa and x.isa.isCplx()=@(1).val.isa.isCplx())],
        e -> [VContainer(@(1).val.child(1) * @(2).val.child(1), @(1).val.isa)]),

    ISum_VContainer := Rule( [ISum, [@(1, VContainer), @(2)]],
        e -> VContainer(ISum(e.var, e.domain, @(2).val), @(1).val.isa)),
    IParSeq_VContainer := Rule( [IParSeq, [@(1, VContainer), @(2)]],
        e -> VContainer(IParSeq(e.var, e.domain, e.fb_cnt, @(2).val), @(1).val.isa)),
    SUM_VContainer := Rule( [@(1, SUM), @(2, VContainer, x -> ForAll(@(1).val.children(), a -> ObjId(a)=VContainer)), ...],
        e -> VContainer( ApplyFunc(SUM, List(e.children(), l -> l.child(1))), @(2).val.isa)),
    VContainer_PullFrom := Rule( [@(1, [NoDiagPullin, NoDiagPullinLeft, NoDiagPullinRight, CR, SymSPL]), @(2, VContainer)],
        e -> VContainer( ObjId(@(1).val)(@(2).val.child(1)), @(2).val.isa )),
    BlockVPerm_VContainer := Rule( [BlockVPerm, @(1, VContainer)],
        e -> VContainer( BlockVPerm(e.n, e.vlen, @(2).val.child(1), e.perm), @(1).val.isa )),
    # WrappedVCons_XXX and XXX_WrappedVCons rules duplicate VCons_XXX and XXX_VCons to suck
    # Vconstructs wrapped into VContainer.
    # V * SUM(a,b,...) -> SUM(V*a, V*b, ...), do not suck in BlockVPerm
    WrappedVCons_XXX := ARule(Compose, [ [@(1, VContainer), @(0, _VConsLeft)], @(2, [SUM, SMPSum, SMPBarrier, ISum, Data]) ], e -> let(s:=@(2).val,
    [ CopyFields(s, rec(_children := List(s.children(), c -> @(1).val * c ),
                    dimensions := [@(1).val.dimensions[1], @(2).val.dimensions[2]])) ])),

    XXX_WrappedVCons := ARule(Compose, [ @(1, [SUM, SUMAcc, SMPSum, SMPBarrier, ISum, Data]), [@(2, VContainer), @(0, _VConsRight)]], e -> let(s:=@(1).val,
    [ CopyFields(s, rec(_children := List(s.children(), c -> c * @(2).val ),
                    dimensions := [@(1).val.dimensions[1], @(2).val.dimensions[2]])) ])),

    # V * SUM(a,b,...) -> SUM(V*a, V*b, ...), do not suck in BlockVPerm
    VCons_XXX := ARule(Compose, [ @(1, _VConsLeft), @(2, [SUM, SMPSum, SMPBarrier, ISum, Data]) ], e -> let(s:=@(2).val,
    [ CopyFields(s, rec(_children := List(s.children(), c -> @(1).val * c ),
                    dimensions := [@(1).val.dimensions[1], @(2).val.dimensions[2]])) ])),

    XXX_VCons := ARule(Compose, [ @(1, [SUM, SUMAcc, SMPSum, SMPBarrier, ISum, Data]), @(2, _VConsRight)], e -> let(s:=@(1).val,
    [ CopyFields(s, rec(_children := List(s.children(), c -> c * @(2).val ),
                    dimensions := [@(1).val.dimensions[1], @(2).val.dimensions[2]])) ])),

    COND_VCons := ARule(Compose, [@(1, COND), @(2, _VConsRight)],
        e -> [COND(@(1).val.cond, @(1).val.child(1)*@(2).val, @(1).val.child(2)*@(2).val)]),

    VCons_COND := ARule(Compose, [@(1, _VConsLeft), @(2, COND)],
        e -> [COND(@(2).val.cond, @(1).val*@(2).val.child(1), @(1).val*@(2).val.child(2))]),

    # BlockVPerm * ISumLS
    BlockVPermISum := ARule(Compose,  [ @(1, BlockVPerm), @(2, ISumLS) ],
         e -> [ ISum(@(2).val.var, @(2).val.domain, @(1).val * @(2).val.child(1)) ]),

    ISumBlockVPerm  := ARule(Compose, [ @(1, ISumLS), @(2, BlockVPerm) ],
        e -> [ ISum(@(1).val.var, @(1).val.domain, @(1).val.child(1) * @(2).val) ]),

    # Replace by vector constructs
    VectGath := Rule([@(1, VTensor), @(2, Gath), ...], e->VGath(@(2).val.func, @(1).val.vlen)),
    VectScat := Rule([@(1, VTensor), @(2, Scat), ...], e->VScat(@(2).val.func, @(1).val.vlen)),
    VectScatAcc := Rule([@(1, VTensor), @(2, ScatAcc), ...], e->VScatAcc(@(2).val.func, @(1).val.vlen)),

    # ----------------------
    # Combine Gath/Scat
    #
    ComposeVGathVGath := ARule(Compose, [ @(1, VGath), @(2, [VGath, VGath_u, VGath_dup], e->@(1).val.v=e.v) ], # o 1-> 2->
        e -> [ ObjId(@(2).val)(fCompose(@(2).val.func, @(1).val.func), @(1).val.v) ]),
    # there is a copy of the rule below in autolib, why?
    ComposeVGath_dup_Gath := ARule(Compose, [@(1, VGath_dup), @(2, [Gath, Prm])],
    e -> [ VGath_dup(fCompose(Cond(ObjId(@2.val)=Prm,@(2).val.func,@(2).val.func), @(1).val.func), @(1).val.v) ]),

    ComposeVScatVScat := ARule(Compose, [ @(1, [VScat_u, VScat, VScatAcc, VScatAcc_u]), @(2, VScat, e->@(1).val.v=e.v) ], # <-1 <-2 o
        e -> [ ObjId(@(1).val)(fCompose(@(1).val.func, @(2).val.func), @(1).val.v) ]),

    ComposeVScatVScatAcc := ARule(Compose, [ @(1, [VScat_u, VScat, VScatAcc, VScatAcc_u]), @(2, VScatAcc, e->@(1).val.v=e.v) ], # <-1 <-2 o
        e -> let(scat := Cond(@(1).val _is VScat,   VScatAcc,
		              @(1).val _is VScat_u, VScatAcc_u,
			      @(1).val),
	    [ scat(fCompose(@(1).val.func, @(2).val.func), @(1).val.v) ])),

    ComposeVScatPrm  := ARule(Compose, [@(1, [VScat, VScat_u, VScatAcc, VScatAcc_u]),   @(2, VPrm_x_I)], # 1-> <-2 o
        e -> [ ObjId(@(1).val)(fCompose(@(1).val.func, @(2).val.func.transpose()), @(1).val.v) ]),
    ComposeVGathPrm  := ARule(Compose, [ @(1, VGath), @(2, VPrm_x_I) ], # o 1-> 2->
        e -> [ VGath(fCompose(@(2).val.func, @(1).val.func), @(1).val.v) ]),
    ComposePrmVScat  := ARule(Compose, [ @(1, VPrm_x_I), @(2, [VScat, VScatAcc]) ], # 1-> <-2 o
        e -> [ ObjId(@(2).val)(fCompose(@(1).val.func.transpose(), @(2).val.func), @(2).val.v) ]),
    ComposePrmVGath  := ARule(Compose, [ @(1, VPrm_x_I), @(2, [VGath, VGath_u])], # o 1-> 2->
        e -> [ ObjId(@(2).val)(fCompose(@(2).val.func, @(1).val.func), @(2).val.v) ]),

    ComposeVTensorPrm  := ARule(Compose, [@(1, VTensor),   @(2, Prm)],
        e -> [ @(1).val, VGath_sv(@(2).val.func, @(1).val.vlen, 1) ]),
    ComposePrmVTensor  := ARule(Compose, [ @(1, Prm), @(2, VTensor) ],
        e -> [  VScat_sv(@(1).val.func.transpose(), @(2).val.vlen, 1), @(2).val ]),

## NOTE: These rules seem broken. There seems to be an assumption on the size of the Scatter. look at Compose_IxVGath_pc__Gath which was fixed for BG/Q 3D FFT
    ComposeIxVGath_pc__Prm  := ARule(Compose, [ @(1, IxVGath_pc), @(2, [Prm, DelayedPrm])], # o 1-> 2->
        e -> let(g:=@(1).val, [ VStretchGath(@(2).val.func, g.k, g.v) ])),
    Compose_IxVGath_pc__Gath  := ARule(Compose, [ @(1, IxVGath_pc), @(2, Gath)], # o 1-> 2->
        e -> let(g:=@(1).val, [ VStretchGath(fCompose(@(2).val.func, fTensor(fId(g.k), fAdd(g.N, g.n, g.ofs))), g.k, g.v) ])),
    Compose_IxVGath_pc__VGath_sv  := ARule(Compose, [ @(1, IxVGath_pc), @(2, VGath_sv)], # o 1-> 2->
        e -> let(g:=@(1).val, [ VStretchGath(fTensor(@(2).val.func, fId(@(2).val.sv)), g.k, g.v) ])),
    Compose_IxVGath_pc__VGath  := ARule(Compose, [ @(1, IxVGath_pc), @(2, VGath)], # o 1-> 2->
        e -> let(g:=@(1).val, [ VStretchGath(fTensor(@(2).val.func, fId(g.v)), g.k, g.v) ])),
    Compose_IxRCVGath_pc__RCGath  := ARule(Compose, [ @(1, IxRCVGath_pc), @(2, VGath_sv, e->e.sv=2)], # o 1-> 2->
        e -> let(g:=@(1).val, [ RCVStretchGath(@(2).val.func, g.k, g.v) ])),

    Compose_VGath__IxVGath_pc := ARule(Compose, [ @(1, VGath, e->e.func.domain()<=2), @(2, IxVGath_pc, e->e.n <= e.v)],
        e -> let(g := @(1).val, gpc := @(2).val,
            Cond(g.func.domain()=1,
                [ VGath_pc(gpc.k*gpc.N, gpc.n, fCompose(fTensor(fId(gpc.k), fAdd(gpc.N, gpc.n, gpc.ofs)), fTensor(g.func, fId(gpc.n))).at(0), gpc.v) ],
                [ SUM(
                    VScat(fBase(V(2),V(0)), g.v) * VGath_pc(gpc.k*gpc.N, gpc.n, fCompose(fTensor(fId(gpc.k), fAdd(gpc.N, gpc.n, gpc.ofs)), fTensor(fCompose(g.func, fBase(V(2),V(0))), fId(gpc.n))).at(0), gpc.v),
                    VScat(fBase(V(2),V(1)), g.v) * VGath_pc(gpc.k*gpc.N, gpc.n, fCompose(fTensor(fId(gpc.k), fAdd(gpc.N, gpc.n, gpc.ofs)), fTensor(fCompose(g.func, fBase(V(2),V(1))), fId(gpc.n))).at(0), gpc.v))
                ]
        ))),

## NOTE: These rules seem broken. There seems to be an assumption on the size of the Scatter. look at Compose_IxVGath_pc__Gath which was fixed for BG/Q 3D FFT
    ComposePrmIxVScat_pc  := ARule(Compose, [ @(1, [Prm, DelayedPrm]), @(2, IxVScat_pc) ], # 1-> <-2 o
        e -> let(s:=@(2).val, [ VStretchScat(@(1).val.func.transpose(), s.k, s.v) ])),
    Compose_Scat__IxVScat_pc  := ARule(Compose, [ @(1, Scat), @(2, IxVScat_pc) ], # 1-> <-2 o
        e -> let(s:=@(2).val, [ VStretchScat(fCompose(@(1).val.func, fTensor(fId(s.k), fAdd(s.N, s.n, s.ofs))), s.k, s.v) ])),
    Compose_VScat_sv__IxVScat_pc  := ARule(Compose, [ @(1, VScat_sv), @(2, IxVScat_pc) ], # 1-> <-2 o
        e -> let(s:=@(2).val, [ VStretchScat(fTensor(@(1).val.func, fId(@(1).val.sv)), s.k, s.v) ])),
    Compose_VScat__IxVScat_pc  := ARule(Compose, [ @(1, VScat), @(2, IxVScat_pc) ], # 1-> <-2 o
        e -> let(s:=@(2).val, [ VStretchScat(fTensor(@(1).val.func, fId(s.v)), s.k, s.v) ])),
    Compose_RCScat__IxRCVScat_pc  := ARule(Compose, [ @(1, VScat_sv, e->e.sv=2), @(2, IxRCVScat_pc) ], # 1-> <-2 o
        e -> let(s:=@(2).val, [ RCVStretchScat(@(1).val.func, s.k, s.v) ])),

    Compose_IxVScat_pc__VScat := ARule(Compose, [ @(1, IxVScat_pc, e->e.n <= e.v), @(2, VScat, e->e.func.domain()<=2)],
        e -> let(s := @(2).val, spc := @(1).val,
            Cond(s.func.domain()=1,
                [ VScat_pc(spc.k*spc.N, spc.n, fCompose(fTensor(fId(spc.k), fAdd(spc.N, spc.n, spc.ofs)), fTensor(s.func, fId(spc.n))).at(0), spc.v) ],
                [ SUM(
                    VScat_pc(spc.k*spc.N, spc.n, fCompose(fTensor(fId(spc.k), fAdd(spc.N, spc.n, spc.ofs)), fTensor(fCompose(s.func, fBase(V(2),V(0))), fId(spc.n))).at(0), spc.v) * VGath(fBase(V(2),V(0)), s.v),
                    VScat_pc(spc.k*spc.N, spc.n, fCompose(fTensor(fId(spc.k), fAdd(spc.N, spc.n, spc.ofs)), fTensor(fCompose(s.func, fBase(V(2),V(1))), fId(spc.n))).at(0), spc.v) * VGath(fBase(V(2),V(1)), s.v))
                ]
        ))),

    ComposeVStretchGath__Prm  := ARule(Compose, [ @(1, VStretchGath), @(2, [Prm,DelayedPrm])], # o 1-> 2->
        e -> let(g:=@(1).val, [ VStretchGath(fCompose(@(2).val.func, g.func), g.part, g.v) ])),
    Compose_VStretchGath__Gath  := ARule(Compose, [ @(1, VStretchGath), @(2, Gath)], # o 1-> 2->
        e -> let(g:=@(1).val, [ VStretchGath(fCompose(@(2).val.func, g.func), g.part, g.v) ])),
    Compose_VStretchGath__VGath  := ARule(Compose, [ @(1, VStretchGath), @(2, VGath)], # o 1-> 2->
        e -> let(g:=@(1).val, [ VStretchGath(fCompose(fTensor(@(2).val.func, fId(g.v)), g.func), g.part, g.v) ])),

    ComposePrmVStretchScat  := ARule(Compose, [ @(1, [Prm,DelayedPrm]), @(2, VStretchScat) ], # 1-> <-2 o
        e -> let(s:=@(2).val, [ VStretchScat(fCompose(@(1).val.func.transpose(), s.func), s.part, s.v) ])),
    Compose_Scat__VStretchScat  := ARule(Compose, [ @(1, Scat), @(2, VStretchScat) ], # 1-> <-2 o
        e -> let(s:=@(2).val, [ VStretchScat(fCompose(@(1).val.func, s.func), s.part, s.v) ])),
    Compose_VScat__VStretchScat  := ARule(Compose, [ @(1, VScat), @(2, VStretchScat) ], # 1-> <-2 o
        e -> let(s:=@(2).val, [ VStretchScat(fCompose(fTensor(@(1).val.func, fId(s.v)), s.func), s.part, s.v) ])),

    Drop_GathScat := ARule(Compose, [ @(1), @(2, VGath),  @(3, VScat, x -> x.transpose()=@(2).val), @(4)],
        e -> [@(1).val, @(4).val]),
# uncombinable gather/scatter CANNOT be pulled in!!!
# NOTE: only pull till BB(), then stop -> new rewriting stage _after_ MarkBB
#    XXX_ISum := ARule(Compose,  [ @(1, [IxVScat_pc, IxRCVScat_pc, VStretchScat, RCVStretchScat]), @(2, ISum) ],
#         e -> [ ObjId(@(2).val)(@(2).val.var, @(2).val.domain, @(1).val * @(2).val.child(1)).attrs(@(2).val) ]),
#
#    ISum_XXX  := ARule(Compose, [ @(1, [ISum, ISumAcc]), @(2, [IxVGath_pc, IxRCVGath_pc, VStretchGath, RCVStretchGath]) ],
#        e -> [ ObjId(@(1).val)(@(1).val.var, @(1).val.domain, @(1).val.child(1) * @(2).val).attrs(@(1).val) ]),

# uncombinable gather/scatter CANNOT be pulled in!!!
# NOTE: only pull till BB(), then stop -> new rewriting stage _after_ MarkBB
#    VScatSUM := ARule(Compose, [ @(1, [IxVScat_pc, IxRCVScat_pc, VStretchScat, RCVStretchScat]), @(2, SUM) ],
#     e -> [ ApplyFunc(ObjId(@(2).val),
#                  List(@(2).val.children(), c -> @(1).val * c)) ]),
#
#    VGathSUM  := ARule(Compose, [ @(1, SUM), @(2, [IxVGath_pc, IxRCVGath_pc, VStretchGath, RCVStretchGath])],
#     e -> [ ApplyFunc(ObjId(@(1).val),
#                  List(@(1).val.children(), c -> c * @(2).val)) ]),

################################################
# NOTE!!!!!!!!!!!!!!!!!!
# GUARD MISSING: gather/scatter must be f x fId(v^2)
#    ComposeBlockVPermVScat := ARule(Compose,[ @(1, BlockVPerm), @(2, VScat)],
#        e->[@2.val, BlockVPerm(@2.val.dimensions[2]/@1.val.child(1).dims()[1], @1.val.vlen, @1.val.child(1), @1.val.perm)]
#    ),
#    ComposeVGathBlockVPerm := ARule(Compose,[ @(1, VGath), @(2, BlockVPerm)],
#        e->[BlockVPerm(@1.val.dimensions[1]/@2.val.child(1).dims()[2], @2.val.vlen, @2.val.child(1), @2.val.perm), @1.val ]
#    ),

    ScatH_VScat_pc := ARule(Compose, [ [ @(1, [Scat,ScatAcc]), fId ], @(2, [VScat_pc, VScat_pcAcc]) ],
	e -> let(s   := @(2).val,
	         oid := Cond((@(1).val _is Scat) and (@(2).val _is VScat_pc), VScat_pc, VScat_pcAcc),
	         [ oid(s.N, s.n, s.ofs, s.v) ])),

    ScatId_VScat_pc := ARule(Compose, [ [ @(1, [Scat,ScatAcc]), [ @(2,H), @, @, @, 1 ] ], @(3, [VScat_pc, VScat_pcAcc]) ],
	e -> let(h   := @(2).val,
	         s   := @(3).val,
	         oid := Cond((@(1).val _is Scat) and (@(3).val _is VScat_pc), VScat_pc, VScat_pcAcc),
	         [ oid(h.params[1], s.n, s.ofs + h.params[3], s.v) ])),

    VGath_pc__toVGath := Rule(@(1, VGath_pc, e->let(v := e.v, (e.N mod v = 0) and (e.n mod v = 0) and (e.ofs mod v = 0))),
        e->let(g := @(1).val, v := g.v, VGath(fAdd(g.N/v, g.n/v, g.ofs/v), v))),

    VScat_pc__toVScat := Rule(@(1, VScat_pc, e->let(v := e.v, (e.N mod v = 0) and (e.n mod v = 0) and (e.ofs mod v = 0))),
        e->let(s := @(1).val, v := s.v, VScat(fAdd(s.N/v, s.n/v, s.ofs/v), v))),

    VScat_pcAcc__toVScatAcc := Rule(@(1, VScat_pcAcc, e->let(v := e.v, (e.N mod v = 0) and (e.n mod v = 0) and (e.ofs mod v = 0))),
        e->let(s := @(1).val, v := s.v, VScatAcc(fAdd(s.N/v, s.n/v, s.ofs/v), v)))

));

Class(RulesTerm, RuleSet);
RewriteRules(RulesTerm, rec(
    ISumLS_Term := Rule(@(1,ISumLS), e->ISum(@(1).val.var, @(1).val.domain, @(1).val.child(1)).attrs(@(1).val)),
    BlockVPerm_Term := Rule(@(1,[BlockVPerm, BlockVPerm2]), e-> Tensor(I(@1.val.n), @1.val._children[1]).sums()),
#   BROKEN RULES; to be redone for loop code
#    IxVGath_pc_term := Rule(@(1,IxVGath_pc, e->e.k>1),
#        e->let(i := Ind(), k:= @1.val.k, v:=@1.val.v,
#            ISum(i, k,
#                VScat(fTensor(fBase(k,i), fId(@1.val.nv/v)), v) *
#                IxVGath_pc(1, k*@1.val.N, @1.val.n, add(mul(i,k), @1.val.ofs), v)
#            ))),
#    IxVScat_pc_term := Rule(@(1,IxVScat_pc, e->e.k>1),
#        e->let(i := Ind(), k:= @1.val.k, v:=@1.val.v,
#            ISum(i, k,
#                IxVScat_pc(1, k*@1.val.N, @1.val.n, add(mul(i,k), @1.val.ofs), v) *
#                VGath(fTensor(fBase(k,i), fId(@1.val.nv/v)), v)
#            ))),
#
    term_PushLR := Rule(@(1, PushLR), e->e.child(1).sums()),

));

Class(RulesTermGrp, RuleSet);
RewriteRules(RulesTermGrp, rec(
    Grp_Term := Rule(@(1,Grp), e-> @1.val._children[1])
));

Class(RulesKickout, RuleSet);
RewriteRules(RulesKickout, rec(
    IxVGath_pc__IxVScat_pc_kickout := ARule(Compose, [@(1, IxVGath_pc), @(2, IxVScat_pc, g->let(s:=@(1).val, s.k=g.k and s.n=g.n and s.N=g.N and s.ofs=g.ofs and s.v=g.v))],
            e->let(g:=@(1).val, [VGath(fId(Rows(g)/g.v), g.v), VScat(fId(Rows(g)/g.v), g.v)])),

    IxVScat_VScat_kickout := Rule(@(1,IxVScat_pc, s-> IsInt(s.n/s.v) and IsInt(s.N/s.v) and IsInt(s.ofs/s.v)),
            e->let(s:=@(1).val, v:=s.v, VScat(fTensor(fId(s.k), fAdd(s.N/v, s.n/v, s.ofs/v)), v))),

    IxVGath_VGath_kickout := Rule(@(1,IxVGath_pc, s-> IsInt(s.n/s.v) and IsInt(s.N/s.v) and IsInt(s.ofs/s.v)),
            e->let(s:=@(1).val, v:=s.v, VGath(fTensor(fId(s.k), fAdd(s.N/v, s.n/v, s.ofs/v)), v))),
));

Class(TerminateSymSPL, RuleSet);
RewriteRules(TerminateSymSPL, rec(
    terminateSymSPL := Rule(@(1, SymSPL), e->@(1).val.child(1))
));

Class(TerminateDPrm, RuleSet);
RewriteRules(TerminateSymSPL, rec(
    terminateDPrm := Rule(@(1, DelayedPrm), e->Prm(@(1).val.func))
));