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

Verifier := rec(
    MaxError := 6,
    IgnoreRoots := false
);

if not IsBound(GenerateErrorReports) then
    GenerateErrorReports := true;
fi;

#F HandleTestError ( <func-name>, <func-args>, <bool> )
#F    Creates a directory and reports the given function call in
#F    'test.g', backtrace in 'backtrace' and current configuration in
#F    'conf.g'
#F
#F    Boolean parameter, if true, tells HandleTestError to all Error()
#F    thus causing an exception, if false, program will continue
#F    normally.
#F
HandleTestError := function (FuncName, FuncArgs, doErr)
    local dir,file,arg,n,k;
    if IsBound(GenerateErrorReports) and GenerateErrorReports = false then
	return;
    fi;
    dir := Concat(Conf("tmp_dir"), Conf("path_sep"), "Error-",
		  String(TimeInSecs()), Conf("path_sep"));
    file := Concat(dir, "test.g");
    SysMkdir(dir);
    # Create a file (truncate if it exists)
    PrintTo(file, "Import(spiral.spl, spiral.nt, spiral.code); \n");
#F    AppendTo(file, "config_update_val(\"remove_temporaries\", SRC_CMDLINE, int_val(0));\n");
#F    AppendTo(file,"config_update_val(\"tmp_dir\", SRC_CMDLINE, str_val(\"./tmp\"));\n");
    AppendTo(file, "GenerateErrorReports := false; \n");
    AppendTo(file, "SPL_DEFAULTS := ", SPL_DEFAULTS, ";\n");
    # write function arguments
    n := 1;
    for arg in FuncArgs do
	AppendTo(file, "arg", String(n), " := ");
	if IsString(arg) then
	    AppendTo(file, "\"", arg, "\"");
	else
	    AppendTo(file, arg);
	fi;
	AppendTo(file, ";\n\n");
	n := n+1;
    od;

    # write function call
    k := 1;
    AppendTo(file, "result := ", FuncName, "( ");
    while k<>n do
	AppendTo(file, "arg", String(k));
	if k <> n-1 then AppendTo(file, ", "); fi;
	k := k + 1;
    od;
    AppendTo(file, " );\n\n");

    # other info
    # NOTE: implement a config dump in sys_conf
#F    AppendTo(Concat(dir, "conf.g"), ConfigProfileList(), ";\n");
    BacktraceTo(Concat(dir, "backtrace.txt"), 100);

    if doErr then
	Error(FuncName, " failed, see ", dir);
    else
	Print(FuncName, " failed, see ", dir, "\n");
    fi;
end;

#F DeriveSPLOptions ( <spl>, <spl-options-record> )
#F    Merges the defaults with <spl-options-record>, derives other
#F    fields, such as dataType, from <spl>, and returns a complete
#F    options record.
#F
DeriveSPLOptions := function (S, R)
    # set options
    R := MergeSPLOptionsRecord(R);

    # check if MPI req'd
#    if IsDMP(S) then
#	R.language := "c.mpi.mpich";
#    fi;

    # if user didn't specify data type determine it from S
    if R.dataType = "no default" then
	if IsRealSPL(S) then R.dataType := "real";
    	else R.dataType := "complex";
	fi;
    else
	;# prevent user from doing nonsense
	#if not IsRealSPL(S) and R.dataType = "real" then
	#    Error("invalid combination: complex <S> and real data type");
	#fi;
    fi;
    return R;
end;

#F DeriveScalarType ( <spl-options-record> )
#F
DeriveScalarType := function(SPLOpts)
    local suffix;
    if IsBound(SPLOpts.customDataType) then return SPLOpts.customDataType;
    elif IsBound(SPLOpts.customReal) and SPLOpts.dataType = "real" then return SPLOpts.customReal;
    elif IsBound(SPLOpts.customComplex) and SPLOpts.dataType = "complex" then return SPLOpts.customComplex;
    else
	if SPLOpts.dataType = "real" then suffix := "";
	elif SPLOpts.dataType = "complex" then suffix := "_cplx";
	else Error("SPLOpts.dataType has invalid value '", SPLOpts.dataType, "'");
	fi;
	if SPLOpts.precision = "single" then return Concat("float",suffix);
	elif SPLOpts.precision = "double" then return Concat("double",suffix);
	elif SPLOpts.precision = "extended" then return Concat("long_double",suffix);
	else Error("SPLOpts.precision has invalid value '", SPLOpts.dataType, "'");
	fi;
    fi;
end;

ProgInputType := rec(
    SPLSource := 0,
    TargetSource := 1,
    ObjFile := 2
);

#F DeriveType ( <spl-options-record> )
DeriveType := (opts)->When(IsBound(opts.vector), opts.vector.isa.ctype,
    DeriveScalarType(opts));

#F ProgSPL ( <spl> , <spl-options-record> )
#F    Convert <spl> to a 'Prog' record used by xxxProg functions.
#F    See gap/src/spiral_spl_prog.c for details.
#F
ProgSPL := function (SPL, Opts)
    local prog;
    Opts := DeriveSPLOptions(SPL, Opts);
    prog := rec();
    prog.profile   := Opts.language;
    prog.type      := ProgInputType.SPLSource;
    prog.data_type := DeriveScalarType(Opts);

    if IsBound(Opts.zeroBits) then prog.zero_bits := Opts.zeroBits;
    else prog.zero_bits := 0; fi;

    prog.dim_rows  := EvalScalar(SPL.dimensions[1]);
    prog.dim_cols  := EvalScalar(SPL.dimensions[2]);
    prog.auto_dim  := 0;

    if IsBound(Opts.compflags) then prog.compiler_flags := Opts.compflags; fi;

    if IsBound(Opts.file) then prog.file := Opts.file;
    else prog.file := SysTmpName(); fi;

    if IsBound(Opts.subName) then prog.sub_name := Opts.subName;
    else prog.sub_name := "sub"; fi;

    prog.spl_file := prog.file;
    return prog;
end;

#F   Valid <compare-type>'s are
#F     "random": compare on random vector (default)
#F     "basis" : compare on standard basis
#F     <int>   : compare on <int> random standard base vectors
#F
VerifierOpts := function (CO)
    local opts;
    opts := Concat("-g -e ", String(Verifier.MaxError), " ");
    if CO = "basis" then return Concat(opts, " -b");
    elif CO = "random" then return Concat(opts, " -r");
    elif IsInt(CO) then return Concat(opts, " -s ", String(CO));
    else Error("<CO> must be an integer, \"random\", or \"basis\"");
    fi;
end;


text:
# Action of a group on a signed point.
BindGlobal("OnSignedPoints", function(dp, signs)
    local p1, p2;
    p1 := Projection(dp, 1);
    p2 := Projection(dp, 2);
    return function(e, g)
        return [OnPoints(e[1], Image(p1, g)),
                Permuted(signs, Image(p2, g))[Position(signs, e[2])]];
    end;
end);

# Action of a product group on vertices of a product graph.
BindGlobal("OnProduct", function(n, dp)
    local pr;
    pr := List([1..n], i -> Projection(dp, i));
    return function(e, g)
        return List([1..n],
            i -> OnPoints(e[i], Image(pr[i], g)));
    end;
end);

# Action of a product group on vertices of a sum graph.
BindGlobal("OnSum", dp -> function(e, g)
        return [e[1], OnPoints(e[2], Image(Projection(dp, e[1]), g))];
    end
);

# Action of a product group on the multiplication table of its factors.
BindGlobal("OnLatinSquare", function(dp)
    local p1, p2;
    p1 := Projection(dp, 1);
    p2 := Projection(dp, 2);
    return function(e, g)
        return [Image(p1, g) * e[1], e[2] * Image(p2, g)];
    end;
end);

# Action of a wreath product on vectors over a ring.
BindGlobal("OnZmodnZVectors", function(d, e)
    local ij;
    ij := Cartesian([1..d], [0..e-1]);
    return function(vec, g)
        local ji, v, w, vw;
        ji := ij{OnTuples([1..d*e], g)};
        vw := List([1..d], i -> ji[Int(vec[i])+(i-1)*e+1]);
        v := List(vw, x -> ZmodnZObj(x[2], e));
        w := List(vw, x -> x[1]);
        return v{List([1..d], i -> Position(w, i))};
    end;
end);

# Action of a wreath product on pairs of vectors with 3 elements.
BindGlobal("OnVectorPairs", function(q)
    local ij;
    ij := Cartesian([1..2], [1..3], [0..q-1]);
    return function(M, g)
        local ji, v, w, vw;
        ji := ij{OnTuples([1..6*q], g)};
        vw := List([1..2],
            i -> List([1..3],
                j -> ji[FFEToInt(M[i][j], q)+(j-1)*q+(i-1)*q*3+1]));
        v := List(vw, r -> List(r, x -> IntToFFE(x[3], q)));
        w := List(vw, r -> List(r, x -> x[2]));
        v := List([1..2], i -> v[i]{List([1..3], j -> Position(w[i], j))});
        w := List([1..2], i -> M[i]{List([1..3], j -> Position(w[i], j))});
        return [v[1] + VectorProduct(w[2], v[2]), v[2]];
    end;
end);

# Action of a wreath product on matrices over a field.
BindGlobal("OnMatrices", function(q, d, e, dp)
    local F, N, p1, p2, pm;
    F := Elements(GF(q));
    N := Position(F, 0*Z(q));
    p1 := Projection(dp, 1);
    p2 := Projection(dp, 2);
    pm := List([0..d-1], i -> List([1..e], j -> Projection(dp, e*i+j+2)));
    return function(M, g)
        return Image(p1, g)*M*Image(p2, g) +
            List(pm, r -> List(r, p -> F[N^Image(p, g)]));
    end;
end);

# Action of a wreath product on Hermitean matrices over a field.
BindGlobal("OnHermiteanMatrices", function(r, d, dp)
    local C, F, K, N, p1, pm;
    F := Elements(GF(r));
    N := Position(F, 0*Z(r));
    K := Elements(GF(r^2));
    C := List(K, x -> Conjugates(GF(r^2), GF(r), x)[2]);
    p1 := Projection(dp, 1);
    pm := List([0..d-1], i -> List([1..d], j -> Projection(dp, d*i+j+1)));
    return function(M, g)
        local H;
        H := Image(p1, g);
        return List(H, r -> List(r, x -> C[Position(K, x)]))*M*TransposedMat(H)
            + ToHermitean(List(pm, r -> List(r, p -> F[N^Image(p, g)])), r);
    end;
end);

# Action of a matrix group on subspaces of a vector space over a finite field.
BindGlobal("OnSubspaces",
    V -> function(S, g)
        return Subspace(V, OnSubspacesByCanonicalBasis(Basis(S), g));
    end);

# Action of a matrix group on a set of subspaces of a vector space over a
# finite field.
BindGlobal("OnSetsSubspaces", function(V)
    local A;
    A := OnSubspaces(V);
    return function(L, g)
        return Set(List(L, S -> A(S, g)));
    end;
end);

# Action on the vertices of doubled Odd graphs.
BindGlobal("OnDoubledOdd", function(n, dp)
    local F, p1, p2;
    F := [x -> x, x -> Difference([1..n], x)];
    p1 := Projection(dp, 1);
    p2 := Projection(dp, 2);
    return function(s, g)
        return F[1^Image(p2, g)](OnSets(s, Image(p1, g)));
    end;
end);

# Action on the vertices of doubled Grassmann graphs.
BindGlobal("OnDoubledGrassmann", function(V, dp)
    local A, F, p1, p2;
    F := [x -> x, OrthogonalSpaceInFullRowSpace];
    p1 := Projection(dp, 1);
    p2 := Projection(dp, 2);
    A := OnSubspaces(V);
    return function(S, g)
        return F[1^Image(p2, g)](A(S, Image(p1, g)));
    end;
end);

# Action on the vertices of Paley graphs.
BindGlobal("OnPaley", function(q, dp)
    local F, N, p1, p2;
    F := Elements(GF(q));
    N := Position(F, 0*Z(q));
    p1 := Projection(dp, 1);
    p2 := Projection(dp, 2);
    return function(x, g)
        return (x + F[N^Image(p1, g)]) * Z(q)^((q-1)^Image(p2, g));
    end;
end);

# Action on the vertices of Preparata graphs.
BindGlobal("OnPreparata", function(q, s, dp)
    local F, N, p1, p2, p3, p4;
    F := Elements(GF(q));
    N := Position(F, 0*Z(q));
    p1 := Projection(dp, 1);
    p2 := Projection(dp, 2);
    p3 := Projection(dp, 3);
    p4 := Projection(dp, 4);
    return function(t, g)
        local w, z;
        z := Z(q)^((q-1)^Image(p1, g));
        w := [z*t[1], t[2] + 2^Image(p2, g),
              t[3]*z^(s+1) + F[N^Image(p3, g)]];
        return List(w, x -> F[Position(F, x)^Image(p4, g)]);
    end;
end);

# Action on the vertices of Kasami graphs.
BindGlobal("OnKasami", function(q, s, dp)
    local Fq, Fs, Nq, Ns, p1, p2, p3;
    Fq := Elements(GF(q));
    Nq := Position(Fq, 0*Z(q));
    Fs := Elements(GF(s));
    Ns := Position(Fs, 0*Z(s));
    p1 := Projection(dp, 1);
    p2 := Projection(dp, 2);
    p3 := Projection(dp, 3);
    return function(v, g)
        local w;
        w := [v[1] + Fq[Nq^Image(p1, g)],
              v[2] + Fs[Ns^Image(p2, g)]];
        return List(w, x -> Fs[Position(Fs, x)^Image(p3, g)]);
    end;
end);

# Action on vertices of additive symplectic covers of complete graphs.
BindGlobal("OnAdditiveSymplecticCover", function(q, m, B, dp)
    local F, N, p1, p2, pr;
    F := Elements(GF(q));
    N := Position(F, 0*Z(q));
    p1 := Projection(dp, 1);
    p2 := Projection(dp, 2);
    pr := List([1..m], i -> Projection(dp, i+2));
    return function(p, g)
        local z;
        z := List([1..m], i -> F[N^Image(pr[i], g)]);
        return [p[1] + F[N^Image(p2, g)] + p[2]*B*z,
                p[2]*Image(p1, g) + z];
    end;
end);

# Action on vertices of multiplicative symplectic covers of complete graphs.
BindGlobal("OnMultiplicativeSymplecticCover", function(q, dp)
    local F, p1, p2, p3;
    F := Elements(GF(q));
    p1 := Projection(dp, 1);
    p2 := Projection(dp, 2);
    p3 := Projection(dp, 3);
    return function(s, g)
        local g2, g3;
        g2 := Image(p2, g);
        g3 := Image(p3, g);
        return OnSets(Set(List(s, p -> List(Permuted([1,2], g2),
                                            i -> F[Position(F, p[i])^g3]))),
                      Image(p1, g));
    end;
end);

# Action on the vertices of affine polar graphs.
BindGlobal("OnAffine", function(q, d, dp)
    local F, N, p1, pr;
    F := Elements(GF(q));
    N := Position(F, 0*Z(q));
    p1 := Projection(dp, 1);
    pr := List([1..d], i -> Projection(dp, i+1));
    return function(v, g)
        return v^Image(p1, g) + List([1..d], i -> F[N^Image(pr[i], g)]);;
    end;
end);

# Action on the roots of E_8
BindGlobal("OnRoots", function(v, g)
    return (-1)^(10^g) * v{OnTuples([1..Length(v)], g)};
end);

# Action of a matrix group on normalized vectors over a semifield.
BindGlobal("OnSemifieldVectors", function(div)
    local norm;
    norm := NormalizeSemifieldVector(div);
    return function(s, g)
        return norm(s*g);
    end;
end);

# Action of a matrix group on sets of normalized vectors over a semifield.
BindGlobal("OnSetsSemifieldVectors", function(div)
    local norm;
    norm := NormalizeSemifieldVector(div);
    return function(S, g)
        return Set(List(S*g, norm));
    end;
end);

# Action on points and lines of Desarguesian projective planes.
BindGlobal("OnProjectivePlane", function(V, dp)
    local F, P, p1, p2;
    F := OnSubspaces(V);
    P := [x -> x, OrthogonalSpaceInFullRowSpace];
    p1 := Projection(dp, 1);
    p2 := Projection(dp, 2);
    return function(S, g)
        return P[2^Image(p2, g)](F(S, Image(p1, g)));
    end;
end);

# Action on points or lines of a point-line geometry.
BindGlobal("OnPointsOrLines", function(act, line)
    return function(x, g)
        if line(x) then
            return Set(List(x, p -> act(p, g)));
        else
            return act(x, g);
        fi;
    end;
end);

# Action on points and lines of Hall planes.
BindGlobal("OnHallPlane", function(q, dp)
    local p5, p6, p7, p8, pr, A, F, N;
    pr := List([0,1], i -> List([1,2], j -> Projection(dp, 2*i+j)));
    p5 := Projection(dp, 5);
    p6 := Projection(dp, 6);
    p7 := Projection(dp, 7);
    F := Elements(GF(q));
    N := Position(F, 0*Z(q));
    A := function(p, g)
        local s, M;
        M := Image(p7, g);
        p := List(p, z -> [z[1] + F[N^Image(p5, g)]*z[2],
                            z[2]*(Z(q)^((q-1)^Image(p6, g)))]);
        if p = [] then
            if IsZero(M[2][1]) then
                return [];
            else
                return [[M[2][2]/M[1][2], 0*Z(q)]];
            fi;
        elif Length(p) = 1 then
            if not IsZero(p[1][2]) then
                return p;
            fi;
            s := M[1][1] + p[1][1]*M[1][2];
            if IsZero(s) then
                return [];
            else
                return [[(M[2][1] + p[1][1]*M[2][2])/s, 0*Z(q)]];
            fi;
        else
            return M*p + List(pr, r -> List(r, e -> F[N^Image(e, g)]));
        fi;
    end;
    return OnPointsOrLines(A, x -> Length(x) > 2);
end);

# Action on points and lines of Hughes planes.
BindGlobal("OnHughesPlane", function(q, div, dp)
    local p1, p2, act, F;
    p1 := Projection(dp, 1);
    p2 := Projection(dp, 2);
    act := OnSemifieldVectors(div);
    F := [Inverse, TransposedMat];
    return function(p, g)
        return [p[1]^Image(p2, g), act(p[2], F[p[1]](Image(p1, g)))];
    end;
end);

ext: gap

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-2021, Carnegie Mellon University
# See LICENSE for details

# Symbolic Scalars
# ================

#F Scalars
#F =======
#F

#F The class "Scalar" represents a symbolic scalar, that is not
#F evaluated. The main motivation is to
#F deal with trigonometric expressions (e.g., cosines) efficiently.
#F The usual numbers provided by gap will be called numbers.
#F
#F A Scalar is a GAP record with the following fields mandatory
#F for all types
#F
#F   isScalar    (=true)
#F   operations  (=ScalarOps)
#F   domain      (=Scalars)
#F   val
#F
#F where val is the actual symbolic expression represented as the delayed
#F evaluation object constructed with Delay() function.
#F

#F CanBeScalar( <obj> )
#F    This functions returns true if an object can become a Scalar record.
#F    Integers and rationals can't become a Scalar record, although they
#F    are considered to be part of Scalars field. This implies that this
#F    function is *misleading*. (NOTE: function def + description)
CanBeScalar :=  obj ->
    IsFloat(obj) or IsComplex(obj) or obj in FieldElements;

Declare(ScalarsOps, IsScalar);

_abs := x->Cond(IsComplex(x), AbsComplex(x),
                IsCyc(x), AbsComplex(Complex(x)),
        AbsFloat(x));

inf_norm := mat ->
    Maximum(List(mat, r -> Sum(r, x -> Cond(IsExp(x),   _abs(x.ev()),
                                    IsValue(x), _abs(x.v), _abs(x)))));
evmat := mat ->
    List(mat, r -> List(r, x -> Cond(IsExp(x), x.ev(), IsValue(x), x.v, x)));
absmat := mat ->
    List(mat, r -> List(r, x -> Cond(IsExp(x), _abs(x.ev()), IsValue(x), _abs(x.v), _abs(x))));
remat := mat -> MapMat(mat, x->ReComplex(Complex(x)));

col := (mat, c) -> List(mat, row->[row[c]]);

InfinityNormMat := inf_norm;

RealMatComplexMat := function(mat)
    local res, row, row1, row2, r, i, e, ee;
    res := []; # list of rows
    for row in mat do
        row1 := []; row2 := [];
    for e in row do
       ee := When(IsExp(e), e.eval(), e); r := re(ee).eval(); i := im(ee).eval();
       Append(row1, [r, -i]);
       Append(row2, [Copy(i), Copy(r)]);
    od;
    Add(res, row1); Add(res, row2);
    od;
    return res;
end;

RealVMatComplexVMat := function(vmat)
    local res, row, row1, row2, r, i, e, ee, el;
    res := []; # list of rows
    for row in vmat do
        row1 := []; row2 := [];
        for el in row do
            if IsRec(el) and IsBound(el.t) and IsArrayT(el.t) then el := el.v; fi;
            if IsList(el) then
                r := [];
                i := [];
                for e in el do
                    ee := When(IsExp(e), e.eval(), e);
                    Append(r, [re(ee).eval()]);
                    Append(i, [im(ee).eval()]);
                od;
            else
                Constraint((IsExp(el) or IsValue(el)) and IsVecT(el.t) and el.t.t<>TComplex);
                r := el;
                i := el.t.value(0);
            fi;
            Append(row1, [r, -i]);
            Append(row2, [Copy(i), Copy(r)]);
        od;
        Add(res, row1); Add(res, row2);
    od;
    return res;
end;

RCMat := RealMatComplexMat;

RCMatCyc := function ( mat )
    local  res, row, row1, row2, r, i, e, ee;
    res := [];
    for row  in mat  do
        row1 := []; row2 := [];
        for e  in row  do
            ee := When(IsExp(e), e.ev(), e);
            r := Re(ee); i := Im(ee);
            Append(row1, [ r, -i ]);
            Append(row2, [ i, r ]);
        od;
        Add(res, row1); Add(res, row2);
    od;
    return res;
end;

getReIm := function(c)
    local r, i;
    if IsCyc(c) then c:=ComplexCyc(c); fi;
    if IsInt(c) or IsFloat (c) then r:=c; i:=0;
    else if IsComplex(c) then r:=ReComplex(c); i:=ImComplex(c);
    else Error("bad matrix entry type");
    fi; fi;
    return [r,i];
end;

#############################################################################
#V  Scalars  . . . . . . . . . . . . . . . . . . . . . . field of all scalars
#V  ScalarsOps  . . . . . . . . . . . . . .  operations record for the domain
##
Class(Scalars, rec(
    isDomain   := true,
    isField    := true,
    #generators := [ ]
    zero       := V(0),
    one        := V(1),
    size       := "infinity",
    char       := 0,
    field      := 0,
    isFinite   := false,
));

Scalars.operations := Class(ScalarsOps, OpsOps, rec(
    \in := ( x, Scalars ) -> IsScalar(x) or IsRat(x) or IsCyc(x) or CanBeScalar(x),
    Field  := elms -> Scalars,
    DefaultField := elms ->Scalars
));

#F IsScalar( <obj> )
#F   tests whether <obj> is a scalar and returns true in the affirmative
#F   case. Else, false is returned.
#F
IsScalar := s -> IsRec(s) and (IsExp(s) or IsValue(s) or IsLoc(s));

#F IsScalarOrNum( <obj> )
#F
IsScalarOrNum := x -> IsScalar(x) or CanBeScalar(x);

#F Scalar( <d> )
#F   adds fields common to all scalars to a record to form a scalar
#F   from the delayed evaluation object <d>
#F
Scalar := d -> When(not IsDelay(d) and IsList(d),
    List(d, Scalar),
    toExpArg(d));

# Construction macro (allows to skip D())
SS := UnevalArgs( x -> Scalar(x) );

# we use lower case function names for optical convenience in SPLs
normPi := function(r)
    # move r into the interval [0,2)
    r := 2 * (r/2 - Int(r/2));
    if r < 0 then
    r := r + 2;
    fi;
    return r;
end;

ExportBagSPL := x->x;

TranslateGapFuncToSPL := rec(
    CosPi := cc -> Print("cos(", ExportBagSPL(cc[2]), "*pi)"),
    SinPi := cc -> Print("sin(", ExportBagSPL(cc[2]), "*pi)"),
    Sqrt  := cc -> Print("sqrt(", ExportBagSPL(cc[2]), ")")
);

ExportBagSPL := function(o)
    local t, cc;
    t := BagType(o);
    cc := Children(o);
    cc := cc{[2..Length(cc)]};

    Cond(t = T_PROD, Print("(", ExportBagSPL(cc[1]), " * ", ExportBagSPL(cc[2]), ")"),
     t = T_QUO,  Print("(", ExportBagSPL(cc[1]), " / ", ExportBagSPL(cc[2]), ")"),
     t = T_SUM,  Print("(", ExportBagSPL(cc[1]), " + ", ExportBagSPL(cc[2]), ")"),
     t = T_DIFF, Print("(", ExportBagSPL(cc[1]), " - ", ExportBagSPL(cc[2]), ")"),
     t = T_FUNCCALL, When(BagType(cc[1]) = T_VAR,
                       let(f  := NameOf(cc[1]),
			   tr := TranslateGapFuncToSPL,
                          When(IsBound(tr.(f)), tr.(f)(cc),
                   Error("Don't know how to translate '", f, "' to SPL")) ),
               Error("Can't translate unnamed function to SPL")),

     t = T_CYC,  let(c := ComplexCyc(o),
                     Print("(", String(ReComplex(c)), ",", String(ImComplex(c)), ")")),
     Print(o));
    return "";
end;

# PrintScalarToSPLNC ( <scalar> )
#   prints <scalar> in spl syntax. Used for exporting (spl.g).
#
PrintScalarToSPLNC := s -> s.printSPL();

Value.printSPL := self >> spiral.spl.PrintNumberToSPLNC(self.v);
Exp.printSPL := self >> Print(self.name, "(",
    DoForAllButLast(self.args, a->Chain(a.printSPL(), Print(", "))),
    Last(self.args).printSPL(), ")");

Loc.printSPL := self >> self.print();

fdiv.printSPL := self >> Print("(", self.args[1].printSPL(), "/", self.args[2].printSPL(), ")");
mul.printSPL := self >> Print("(", self.args[1].printSPL(), "*", self.args[2].printSPL(), ")");
add.printSPL := self >> Print("(", self.args[1].printSPL(), "+", self.args[2].printSPL(), ")");
sub.printSPL := self >> Print("(", self.args[1].printSPL(), "-", self.args[2].printSPL(), ")");
neg.printSPL := self >> Print("(-", self.args[1].printSPL(), ")");

cospi.printSPL := self >> Print("cos(", self.args[1].printSPL(), "*pi)");
sinpi.printSPL := self >> Print("sin(", self.args[1].printSPL(), "*pi)");
ScalarIsCos := x -> ObjId(x) = cospi;
ScalarIsSin := x -> ObjId(x) = sinpi;
ScalarCosArg := x -> x.args[1];

#F Functions for Scalars
#F ---------------------
#F

#F EvalScalar ( <scalar>/<number> )
#F   returns the number represented by <scalar>. If a <number>
#F   is given, then <number> is returned.
#F
EvalScalar := s -> _unwrap(
    Cond(s _is funcExp,
	    s.eval(),
	 IsScalar(s) and IsBound(s.ev),
	    s.ev(),
	 Eval(s)));

#F <s1> = <s2>
#F   returns true if both <s1> and <s2> are scalars, not of
#F   type "float" or "complex", and represent the same number.
#F
#Ops.\= := function ( s1, s2 )
#    return IsScalar(s1) and IsScalar(s2) and
#           s1.ev() = s2.ev();
#end;

#Exp.domain := Scalars;
#Value.domain := Scalars;
#Loc.domain := Scalars;

ext: gi


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

# mn=16, m=4, n=4, h=3
ImportAll(spiral);

# even k, P permutation on pairs
 P := (km,k) -> (i -> let(m := km/k, h := Int((m+2)/2), mc:=Int((m+1)/2),
	 Cond(
	     i < h, k*i,
	     i >= km/2-mc+1, k/2 + k*(i-(km/2-mc+1)),
	     ((i-h) mod m) mod 2 = 0, k*((i-h) mod m)/2 + 1 + Int((i-h)/m),
		              k*((i-h) mod m)/2 + 1 + Int(k/2-1-(i-h)/m)
	 )));

# even k, P permutation on single el-s
 rP := (km,k) -> (i -> let(m := km/k,
	 Cond(
	     i = 0, 0,
	     1 <= i and i < m, -1 + 2*k*Int((i+1)/2) + ((i-1) mod 2),
	     i >= km-m, k + 2*k*Int((i-(km-m))/2) - 1 + ((i-m) mod 2),
	     (Int((i-m)/2) mod m) mod 2 = 0, k*(Int((i-m)/2) mod m) + 1 + 2*Int((i-m)/(2*m)) + (i-m) mod 2,
	                                     k*(Int((i-m)/2) mod m) + 1 + 2*Int((km-m-i)/(2*m)) + (i-m) mod 2
	 )));

# odd k, \hat{P} permutation on pairs
 Phat := (km,k) -> (i -> let(m := km/k, h := Int((m+2)/2),
	 Cond(
	     i < h, k*i,
	     ((i-h) mod m) mod 2 = 0, Int(k*((i-h) mod m)/2) + 1 + Int((i-h)/m),
		                      Int(k*((i-h) mod m)/2) +  Int((k+1)/2-(i-h)/m)
	 )));

# odd k
 rPhat := (km,k) -> (i -> let(m := km/k,
	 Cond(
	     i = 0, 0,
	     1 <= i and i < m, -1 + 2*k*Int((i+1)/2) + ((i-1) mod 2),
	     (Int((i-m)/2) mod m) mod 2 = 0, k*(Int((i-m)/2) mod m) + 1 + 2*Int((i-m)/(2*m)) + (i-m) mod 2,
		                             k*(Int((i-m)/2) mod m) +     2*Int((km-i)/(2*m)) + (i-m) mod 2
	 )));

# odd k
 rQhat := (km,k) -> (i -> let(m := km/k,
	 Cond(
	     i >= km-m, k-1 + 2*k*Int((i-(km-m))/2) + (i-(km-m)) mod 2,
	     (Int(i/2) mod m) mod 2 = 0, k*(Int(i/2) mod m) +     2*Int(i/(2*m)) + i mod 2,
	                                 k*(Int(i/2) mod m) + 1 + 2*Int((km-m-i)/(2*m)) + i mod 2
	 )));

# Pshow := (km,k)->let(N := Int(km/2)+1, pm(Perm(PermFunc(P(km,k),N),N)));
# rPshow := (km,k)->let(N := km, pm(Perm(PermFunc(rP(km,k),N),N)));
# Phatshow := (km,k)->let(N := Int((km+2)/2), pm(Perm(PermFunc(Phat(km,k),N),N)));

PermP := (m, k) -> Checked(IsEvenInt(k), PermFunc(rP(k*m,k),k*m)^-1);
PermPhat := (m, k) -> Checked(IsOddInt(k), PermFunc(rPhat(k*m,k),k*m)^-1);
PermQhat := (m, k) -> Checked(IsOddInt(k), PermFunc(rQhat(k*m,k),k*m)^-1);

p_4_4 := PermP(4,4);
p_4_8 := PermP(4,8);
p_8_4 := PermP(8,4);
p_8_8 := PermP(8,8);

p_3_4 := PermP(3,4);
p_5_4 := PermP(5,4);
p_5_8 := PermP(5,8);

phat_4_3 := PermPhat(4,3);
phat_4_5 := PermPhat(4,5);
phat_8_5 := PermPhat(8,5);
phat_8_7 := PermPhat(8,7);
phat_7_7 := PermPhat(7,7);

qhat_4_3 := PermQhat(4,3);
qhat_4_5 := PermQhat(4,5);
qhat_8_5 := PermQhat(8,5);
qhat_8_7 := PermQhat(8,7);
qhat_7_7 := PermQhat(7,7);


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

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

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

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

fftCUDADeviceOpts := function(arg) # specific to FFT size 100...
    local opts;

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

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

    opts.globalUnrolling := 27;
    opts.codegen.GathPtr := fftx.codegen.MultiPtrCodegenMixin.GathPtr;
    opts.codegen.ScatPtr := fftx.codegen.MultiPtrCodegenMixin.ScatPtr;
    opts.codegen.OO := (self, o, y, x, opts) >> skip();

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

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

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

    return opts;
end;

confFFTCUDADevice := rec(
    defaultName := "confFFTCUDADevice",
    defaultOpts := (arg) >> rec(useFFTCUDADevice := true),
    confHandler := fftCUDADeviceOpts
);

fftx.FFTXGlobals.registerConf(confFFTCUDADevice);

max_stars_repo_name: mikefranusich/spiral-package-fftx

max_issues_repo_name: mikefranusich/spiral-package-fftx

max_forks_repo_name: mikefranusich/spiral-package-fftx


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

ImportAll(paradigms.vector);

Declare(ASPF);

Class(flipBot, PermClass, rec(
    def := k -> rec(size := 2*k),
    lambda := self >> let(k := self.params[1],
	i := Ind(2*k),
	Lambda(i, cond(leq(i, k-1), i, i + 1 - 2*imod(i,2))))
));

Class(ASPF, TaggedNonTerminal, rec(
   _short_print := true,
   abbrevs := [ (alg, tbasis, fbasis) -> Checked(
           IsASPAlgebra(alg), IsASPTimeBasis(tbasis), IsASPFreqBasis(fbasis),
           [alg, tbasis, fbasis]) ],

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

   dims := self >> let(n:=self.params[1].n, [n, n]),

   print := (self, i, is) >> let(base_print := NonTerminal.print,
       Print(base_print(self, i, is),
           When(self.tags<>[], Print(".withTags(", self.tags, ")")))),

   hashAs := self >> let(p := self.params,
       t := ObjId(self)(p[1].hashAs(), p[2], p[3].hashAs()).withTags(self.getTags()),
       When(self.transposed, t.transpose(), t)),

   HashId := self >> [ # NOTE: why is this function necessary?
       [ObjId(self.params[1]), self.params[1].n], self.params[2], ObjId(self.params[3]),
        self.getTags(), self.transposed ],

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

   _scale12 := (self, s1, s2) >> let(n := self.params[1].n, Cond(
       n<=2,         s1,
       IsEvenInt(n), Diag(diagDirsum(fConst(TReal, 2, s1), fConst(TReal, n-2, s2))),
       <# else #>    Diag(diagDirsum(fConst(TReal, 1, s1), fConst(TReal, n-1, s2))))),

   _scale34 := (self, s1, s2) >> let(n := self.params[1].n, Cond(IsEvenInt(n),
       s2,
       Diag(diagDirsum(fConst(TReal, n-1, s2), fConst(TReal, 1, s1))))),

   terminate := self >> let(
       A := self.params[1], T := self.params[2], F := self.params[3], s1 := F.scale1d,
       n := EvalScalar(A.n),
       res := CondPat(self,
           [ASPF, XN_min_1,  Time_TX, Freq_1,  ...], s1 * DFT (n, A.rot),
           [ASPF, XN_min_1,  Time_TX, Freq_1H, ...], s1 * DFT2(n, A.rot),
           [ASPF, XN_plus_1, Time_TX, Freq_1,  ...], s1 * DFT3(n, A.rot),
           [ASPF, XN_plus_1, Time_TX, Freq_1H, ...], s1 * DFT4(n, A.rot),

           [ASPF, XN_skew,  @,  @(1, [Freq_1,Freq_1H]), ...], let(
               aa := EvalScalar(A.a),
               p := When(aa > 1/2, J(n), I(n)) * LIJ(n).transpose(),
               p * s1 * Cond(ObjId(F)=Freq_1, BSkewDFT3(n, aa, A.rot), BSkewDFT4(n, aa, A.rot)) * T.toX(A)),

           [ASPF, XN_min_1, ...],
               DirectSum(List(A.rspectrum(), a -> F.from1X(a))) * BRDFT1(n, A.rot) * T.toX(A),
           [ASPF, XN_min_1U, ...],
               DirectSum(List(A.rspectrum(), a -> F.from1X(a))) * UBRDFT1(n, A.rot) * T.toX(A),
           [ASPF, XN_plus_1, ...],
               DirectSum(List(A.rspectrum(), a -> F.from1X(a))) * BRDFT3(n, 1/4, A.rot) * T.toX(A),
           [ASPF, XN_skew, ...], let(aa := EvalScalar(A.a),
               p := When(aa > 1/2, RC(J(n/2)), RC(I(n/2))) * RC(LIJ(n/2).transpose()),
               DirectSum(List(A.rspectrum(), s -> F.from1X(s))) * p * BRDFT3(n, aa, A.rot) * T.toX(A)),

           Error("not implemented")),
       Cond(self.transposed, TerminateSPL(res).transpose(), TerminateSPL(res))),

   transpose := self >> Cond(
       ObjId(self.params[1]) = XN_min_1 and self.params[2] = Time_TX and self.params[3]=Freq_1(1),
           self, # DFT
       CopyFields(self, rec(
               transposed := not self.transposed,
               dimensions := Reversed(self.dimensions) ))),

   conjTranspose := self >> Cond(
       self.isReal(),
           self.transpose(),
       CopyFields(self, rec(
	       params := [self.params[1].conj(), self.params[2], self.params[3]],
               transposed := not self.transposed,
               dimensions := Reversed(self.dimensions) ))),

   normalizedArithCost := self >> let(n := self.params[1].n,
       nlogn := n * log(n) / log(2),
       When(self.params[3].dim() = 1, floor(5*nlogn), floor(2.5*nlogn)))
));

Class(ASP, rec(
  rot := 1,
  __call__ := (self, rot) >> WithBases(self, rec(rot := rot, operations := PrintOps)),
  print := self >> Print(self.__name__, "(", self.rot, ")"),

  DFT    := (self, n) >> ASPF(XN_min_1(n,self.rot),  Time_TX, Freq_1(1)),
  RDFT   := (self, n) >> ASPF(XN_min_1(n,self.rot),  Time_TX, Freq_E(1,1)),
  URDFT  := (self, n) >> ASPF(XN_min_1U(n,self.rot), Time_TX, Freq_E(1,1)),
  DHT    := (self, n) >> ASPF(XN_min_1(n,self.rot),  Time_TX, Freq_H(1,1)),
  UDHT   := (self, n) >> ASPF(XN_min_1U(n,self.rot), Time_TX, Freq_H(1,1)),

  BRDFT  := (self, n) >> ASPF(XN_min_1(n,self.rot),  Time_TX, Freq_T(1,1)),
  MBRDFT := (self, n) >> ASPF(XN_min_1(n,self.rot),  Time_SX, Freq_S(1,1)), # Murakami BRDFT variant

  IRDFT  := (self, n) >> ASPF(XN_min_1(n,self.rot),  Time_TX, Freq_E(1,2)).transpose(),
  IURDFT := (self, n) >> ASPF(XN_min_1U(n,self.rot), Time_TX, Freq_E(1,2)).transpose(),

  DFT2   := (self, n) >> ASPF(XN_min_1(n,self.rot), Time_TX,  Freq_1H(1)),
  RDFT2  := (self, n) >> ASPF(XN_min_1(n,self.rot), Time_TX,  Freq_EH(1,1)),
  DHT2   := (self, n) >> ASPF(XN_min_1(n,self.rot), Time_TX,  Freq_HH(1,1)),
  BRDFT2 := (self, n) >> ASPF(XN_min_1(n,self.rot), Time_TX,  Freq_TH(1,1)),
  SBRDFT2:= (self, n) >> ASPF(XN_min_1(n,self.rot), Time_TX,  Freq_THU(1,1)),

  DFT3   := (self, n) >> ASPF(XN_plus_1(n,self.rot), Time_TX, Freq_1(1)),
  RDFT3  := (self, n) >> ASPF(XN_plus_1(n,self.rot), Time_TX, Freq_E(1,1)),
  DHT3   := (self, n) >> ASPF(XN_plus_1(n,self.rot), Time_TX, Freq_H(1,1)),
  BRDFT3 := (self, n) >> ASPF(XN_plus_1(n,self.rot), Time_TX,  Freq_T(1,1)),
  MBRDFT3:= (self, n) >> ASPF(XN_plus_1(n,self.rot), Time_SX, Freq_S(1,1)), # Murakami BRDFT variant

  BDFT   := (self, n,a) >> ASPF(XN_skew(n, a,self.rot), Time_EX, Freq_1(1)),
  rDFT   := (self, n,a) >> ASPF(XN_skew(n, a,self.rot), Time_EX, Freq_E(1,1)),
  rDFTII := (self, n,a) >> ASPF(XN_skew(n, a,self.rot), Time_EX, Freq_EH(1,1)),
  rDHT   := (self, n,a) >> ASPF(XN_skew(n, a,self.rot), Time_HX, Freq_H(1,1)),
  rDHTII := (self, n,a) >> ASPF(XN_skew(n, a,self.rot), Time_HX, Freq_HH(1,1)),

  brDFT   := (self, n,a) >> ASPF(XN_skew(n, a,self.rot), Time_SX, Freq_E(1,1)),
  brDFTII := (self, n,a) >> ASPF(XN_skew(n, a,self.rot), Time_SX, Freq_EH(1,1)),
  brDHT   := (self, n,a) >> ASPF(XN_skew(n, a,self.rot), Time_SX, Freq_H(1,1)),
  brDHTII := (self, n,a) >> ASPF(XN_skew(n, a,self.rot), Time_SX, Freq_HH(1,1)),

  bRDFT3 := (self, n,a) >> ASPF(XN_skew(n,a,self.rot),  Time_SX, Freq_S(1,1)),

  skewSS := (self, n,a) >> ASPF(XN_skew(n,a,self.rot), Time_SX, Freq_S(1,1)),
  skewTT := (self, n,a) >> ASPF(XN_skew(n,a,self.rot), Time_TX, Freq_T(1,1)),

  DFT4   := (self, n) >> ASPF(XN_plus_1(n,self.rot), Time_TX, Freq_1H(1)),
  RDFT4  := (self, n) >> ASPF(XN_plus_1(n,self.rot), Time_TX, Freq_EH(1,1)),
  DHT4   := (self, n) >> ASPF(XN_plus_1(n,self.rot), Time_TX, Freq_HH(1,1)),

  BRDFT4 := (self, n,a) >> ASPF(XN_skew(n,a, self.rot), Time_TX, Freq_TH(1,1)),
  SBRDFT4:= (self, n,a) >> ASPF(XN_skew(n,a, self.rot), Time_TX, Freq_THU(1,1)),
));

Class(TCodelet, Tagged_tSPL_Container);

_mid := F -> When(ObjId(F)=Freq_1, Freq_E, ObjId(F));

Freq_H.fixbot := k -> Prm(flipBot(k));
Freq_E.fixbot := k -> Diag(BHD(k, 1.0, -1.0));
Freq_1.fixbot := k -> I(2*k);
Freq_S.fixbot := k -> I(2*k);
Freq_T.fixbot := k -> I(2*k);
# NOTE: what is it for others?
# NOTE: explain fixbot

nsFiltered := (lst, func) -> Cond(
    IsValue(lst), Filtered(lst.v, func),
    IsSymbolic(lst), lst,
    Filtered(lst, func)
);

ASPF.isType1 := self >> ObjId(self.params[1]) in [XN_min_1, XN_min_1U];
ASPF.isType3 := self >> ObjId(self.params[1]) in [XN_plus_1];
ASPF.isSkew  := self >> ObjId(self.params[1]) in [XN_skew];
ASPF.size := self >> self.params[1].n;

ASPF_Breakdown_Rule := rec(
    a := rec(
        extraLeftTags := [],
        extraRightTags := [],
        inplaceTag := false,
	maxRadix := -1,
    ),

    apply    := (t, C, N) -> C[1],
    inplace := (self, x) >> When(self.getA("inplaceTag"), Inplace(x), x),
    leftTags := (self, t) >> Concatenation(t.tags, self.getA("extraLeftTags", [])),
    forTransposition := true,
);

Class(even, AutoFoldExp, rec(
    ev := self >> let(a := self.args[1].ev(),
        Cond(not IsInt(a), Error("even(<n>) works only with integer <n>"),
             IsEvenInt(a), 1,
             0)),
    computeType := self >> TInt
));

Class(odd, AutoFoldExp, rec(
    ev := self >> let(a := self.args[1].ev(),
        Cond(not IsInt(a), Error("odd(<n>) works only with integer <n>"),
             IsOddInt(a), 1,
             0)),
    computeType := self >> TInt
));

Class(firstEltRft, FuncClass, rec(
    def := (n,f) -> Checked(IsIntSym(n), IsFunction(f), let(d:=f.domain(), r:=f.range(),
                   rec(n := n, N := 2*r-odd(n)))),
    domain := self >> self.params[1],
    range  := self >> 2*self.params[2].range() - odd(self.params[1]),

    lambda := self >> let(
        i := Ind(self.domain()), n := self.params[1], f := self.params[2],
        fst := 1+even(n), o := odd(n),
        Lambda(i, cond(lt(i,fst), i,  imod(i+o, 2) + 2*f.at(idiv(i+o,2)) - o))),

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

Class(ASPF_CT1_DFT_Mat, TaggedNonTerminal, rec(
    isAuxNonTerminal := true,
    abbrevs := [ k -> Checked(IsPosIntSym(k), [k]) ],
    dims := self >> [2*self.params[1], 2*self.params[1]],
    isReal := self >> true,
    terminate := self >> let(k := self.params[1],
        mat := Mat([[1,1],[-E(4), E(4)]]),
        res := L(2*k, 2) * RC(Tensor(I(k/2), mat) * MM(2,k/2)) * Diag(BHD(k,1.0,-1.0)),
        Cond(self.transposed, res.transpose(), res)
    )
));

NewRulesFor(ASPF_CT1_DFT_Mat, rec(
    ASPF_CT1_DFT_Mat_terminate := rec(
        applicable := t -> true,
        freedoms := t -> [],
        child := (t, fr) -> [],
        apply := (self, t, C, Nonterms) >> t.terminate()
    )
));

NewRulesFor(ASPF, rec(
    # Type 1
    #
    ASPF_CT1_URFT := CopyFields(ASPF_Breakdown_Rule, rec(
        applicable := t -> let(n := t.size(),
            logic_and(t.isType1(), logic_and(eq(n mod 2, 0), logic_neg(isPrime(idiv(n,2)))))),

        freedoms := (self, t) >> let(maxR := self.getA("maxRadix", -1),
            [ nsFiltered(divisorsIntNonTriv(t.size()/2), x -> (maxR < 0) or (x <= maxR/2)) ]),

        child := (self, t, fr) >> let(
            ltags := self.leftTags(t), rtags := self.getA("extraRightTags", []),
            A := t.params[1],  F    := t.params[3],  T    := t.params[2],
            k := fr[1],        Fmid := _mid(F),      Tmid := Fmid.timeBasis(),
            m := A.n / (2*k),  j    := Ind(m-1),     aj   := fdiv(j+1, 2*m),

            C1 := ASPF(ObjId(A) (2*k,     A.rot),  Time_TX, F),
            C2 := ASPF(XN_skew  (2*k, aj, A.rot),  Tmid,    F),
            C3 := ASPF(XN_min_1U(2*m,     A.rot),  T,       Fmid(1,1)),

	    [ GT(C1.withTags(ltags), fId(Cols(C1)), fId(Rows(C1)), []),
	      C2.withTags(ltags),
	      GT(C3.withTags(rtags), GTVec, GTVec, [k]).withTags(t.tags),
	      InfoNt(j) ]
        ),

	apply := (self, t, C, Nonterms) >> let(
            A := t.params[1],  F := t.params[3], j := Nonterms[4].params[1],
	    k := Cols(Nonterms[1])/2, m := Cols(Nonterms[3].params[1])/2,
            When(F.dim()=1, Scat(Refl0_u(m, 2*k)),      # if transform is complex
                            RC(Scat(Refl0_u(m, k)))) *  # if it is real
            DirectSum(
                C[1],
		IDirSum(j, F.fixbot(k) * C[2])) *
            self.inplace(C[3])
        )
    )),

    ASPF_CT1Prm_URFT := CopyFields(~.ASPF_CT1_URFT, rec(
        child := (self, t, fr) >> let(
            ltags := self.leftTags(t), rtags := self.getA("extraRightTags", []),
            A := t.params[1],    T := t.params[2],    F := t.params[3],
            k := fr[1],          Fmid := _mid(F),     Tmid := Fmid.timeBasis(),
            m := div(A.n, 2*k),  j := Ind(m-1),       aj := fdiv(j+1, 2*m),

            C1 := ASPF(ObjId(A) (2*k,     A.rot),  Time_TX, F),
            C2 := ASPF(XN_skew  (2*k, aj, A.rot),  Tmid,    F),
            C3 := ASPF(XN_min_1U(2*m,     A.rot),  T,       Fmid(1,1)),

	    [ GT(C1.withTags(ltags), fId(Cols(C1)), fId(Rows(C1)), []),
	      C2.withTags(ltags),
	      GT(C3.withTags(rtags), GTVec, GTPar, [k]).withTags(t.tags),
	      InfoNt(j) ]
	),
        apply := (self, t, C, Nonterms) >> let(
            A := t.params[1],  F := t.params[3], j := Nonterms[4].params[1],
            N := t.params[1].n,          k := div(Cols(C[1]), 2),   tr := Tr(k, 2),
            j := Nonterms[4].params[1],  m := div(N, 2*k),
            fixbot := t.params[3].fixbot(k),

            When(F.dim()=1, Scat(Refl0_u(m, 2*k)),      # if transform is complex
                            RC(Scat(Refl0_u(m, k)))) *  # if it is real
	    DirectSum(
		C[1] * tr,
		IDirSum(j, F.fixbot(k) * C[2] * tr)) *

            RC(Tr(k, m)) * C[3]
        )
    )),

    ASPF_CT1_DFT := CopyFields(ASPF_Breakdown_Rule, rec(
        applicable := t -> let(n := t.size(),
            logic_and(t.isType1(), logic_and(eq(n mod 4, 0), logic_neg(isPrime(idiv(n,4)))))),

        freedoms := (self, t) >> let(maxR := self.getA("maxRadix", -1),
            [ nsFiltered(divisorsIntNonTriv(t.size()/4), x -> (maxR < 0) or (x <= maxR/2)) ]),

        child := (self, t, fr) >> [let(
            ltags := self.leftTags(t),
            A := t.params[1],  F    := t.params[3],  T    := t.params[2],
            k := 2*fr[1],      Fmid := _mid(F),      Tmid := Fmid.timeBasis(),
            m := A.n / (2*k),  j    := Ind(m-1),     aj   := fdiv(j+1, 2*m),

            mat := Mat([[1,1],[-E(4), E(4)]]),
            C1 := ASPF(ObjId(A) (2*k,     A.rot),  Time_TX, F),
            C2 := ASPF(XN_skew  (2*k, aj, A.rot),  Tmid,    CopyFields(F, rec(scale2d:=1/2*F.scale2d))),
            C3 := UDFT(2*m, A.rot),

            P := When(F.dim()=1, Scat(Refl0_u(m, 2*k)),      # if transform is complex
                                 RC(Scat(Refl0_u(m, k)))),   # if it is real
            Pt := When(F.dim()=1, Prm(Refl0_u(m, 2*k)),      # if transform is complex
                                 RC(Prm(Refl0_u(m, k)))),   # if it is real
            self.inplace(
                P *
                DirectSum(
                    C1.withTags(ltags) * RC(L(k, 2)),
                    IDirSum(j, F.fixbot(k) * C2.withTags(ltags) * PushL(ASPF_CT1_DFT_Mat(k)))
                ) *
                Pt
            ) *
            RC(Tensor(I(k/2), C3) *
               Tr(2*m,k/2))
        )]
    )),

    ASPF_CT1Odd_RFT := CopyFields(~.ASPF_CT1_URFT, rec(
        applicable := t -> let(n := t.size(),
            logic_and(t.isType1(), logic_and(logic_neg(isPrime(n)), hasOddDivisors(n)))),

        freedoms := (self, t) >> let(maxR := self.getA("maxRadix", -1),
            [ nsFiltered(oddDivisorsIntNonTriv(t.size()), x -> (maxR < 0) or (n/x <= maxR/2)) ]),

        # m is odd here
        child := (self, t, fr) >> [let(
            ltags := self.leftTags(t),
            A := t.params[1],  T    := t.params[2],  F    := t.params[3],       N  := A.n,
            m := fr[1],        Fmid := _mid(F),      Tmid := Fmid.timeBasis(),  mf := (m-1)/2,
            k := N / m,        j    := ind(mf,1),    aj   := fdiv(j+1, m),      Nc := idiv(N+1, 2),
            kc:= idiv(k+1,2),
            C1 := ASPF(ObjId(A) (  k,     A.rot), Time_TX, F).withTags(ltags),
            C2 := ASPF(XN_skew  (2*k, aj, A.rot), Tmid,    F).withTags(ltags),
            C3 := ASPF(XN_min_1 (  m,     A.rot), T,       Fmid(1,1)),
            fst := 1 + even(k),
            cpx := F.dim()=1,
#           rp := fDirsum(fId(1+even(k)), fTensor(Refl0_odd(mf, k, 1), fId(2))),
#           When(F.dim()=1, Scat(Refl0_odd(mf, 2*k, 0)), Scat(rp)) *
            SUM(
                Scat(Cond(cpx, HH(N,k,0,[m]), firstEltRft(k, HH(Nc,kc,0,[m])))) * C1 * Gath(HH(N, k, 0, [1])),

                GT(F.fixbot(k) * C2,   # reflect = N-2  in BH compensates for fAdd shift
                   HH(N, 2*k, k, [1, 2*k]),
                   Cond(cpx, BHH(N, 2*k, 1, [m,1], 2*N),
                             fCompose(HH(N, N-fst, fst, [1]), fTensor(BHH(Nc-1,k,0,[m,1], N-2), fId(2)))),
                   [mf])) *
            self.inplace(
                GT(C3, GTVec, GTVec, [k]).withTags(t.tags))
        )]
    )),

#     ASPF_CT1Inp_URFT := CopyFields(~.ASPF_CT1_URFT, rec(
#         child := (self, t, fr) >> let(
#             A := t.params[1],    T := t.params[2],    F := t.params[3],         r := A.rot,
#             k := fr[1],          Fmid := _mid(F),     Tmid := Fmid.timeBasis(),
#             m := div(A.n, 2*k),  j := Ind(m-1),       aj := fdiv(j+1, 2*m),

#             [ ASPF(ObjId(A) (2*k, r),     Time_TX, F).withTags([ANoRecurse()]),
#               ASPF(XN_skew  (2*k, aj, r), Tmid,    F).withTags([ANoRecurse()]),
#               ASPF(XN_min_1U(2*m, r),     T,       Fmid(1,1)),
#               InfoNt(j) ]),

#         apply := (self, t, C, Nonterms) >> let(
#             N := t.params[1].n,         k := div(Cols(C[1]), 2),
#             j := Nonterms[4].params[1], m := div(N, 2*k),

#             RC(Scat(Refl0_u(m, k))) *
#             DirectSum(BB(C[1]*Tr(k,2)),
#                       IDirSum(j, BB(Diag(BHD(k,1,-1))*C[2]*Tr(k,2)*RC(MM(2,k/2))))) *
#             RC(Gath(Refl0_u(m, k))) *
#             RC(condIOS(k*m, m))*
#             Tensor(I(k), C[3]) * Tr(2*m, k)
#         )
#     )),

    # Type 3
    #
    ASPF_CT3_RFT := CopyFields(ASPF_Breakdown_Rule, rec(
        aj := (j, a, m) -> fdiv(j+1/2, 2*m),

        libApplicable := t -> logic_and(eq(imod(t.params[1].n, 2),0), logic_neg(isPrime(div(t.params[1].n,2)))),
        applicable := t -> let(A := ObjId(t.params[1]), n := t.params[1].n,
            A = XN_plus_1 and (IsSymbolic(n) or (n > 4 and IsEvenInt(n)))),

        freedoms := (self, t) >> let(n := t.params[1].n, maxR := self.getA("maxRadix", -1),
            [ When(IsSymbolic(n), divisorsIntNonTriv(div(n,2)),
                                  Filtered(divisorsIntNonTriv(div(n,2)).ev(), x -> (maxR < 0) or (x <= maxR/2))) ]),

        child := (self, t, fr) >> let(
            A := t.params[1],     T := t.params[2],  F := t.params[3],         r := A.rot,
            k := fr[1],           Fmid := _mid(F),   Tmid := Fmid.timeBasis(), a := A.a,
            m := div(A.n, 2*k),   j := Ind(m),       aj := self.aj(j, a, m),

            C1 := ASPF(XN_skew(2*k, aj, r),                 Tmid, F),
            C2 := ASPF(CopyFields(A, rec(n:=_unwrap(2*m))), T,    Fmid(1,1)),

            [ C1.withTags(t.tags), GT(C2, GTVec, GTVec, [k]).withTags(t.tags), InfoNt(j) ]),

        apply := (self, t, C, Nonterms) >> let(
            N := t.params[1].n,         k := div(Cols(C[1]), 2),
            j := Nonterms[3].params[1], m := div(N, 2*k),
            When(t.params[3].dim()=1,
                 Scat(Refl1(m, 2*k)) *   IDirSum(j,                   C[1]) * C[2],
                 RC(Scat(Refl1(m, k))) * IDirSum(j, Diag(BHD(k,1.0,-1.0))*C[1]) * C[2]))
    )),

    # Skew
    #
    ASPF_CTSkew_RFT := CopyFields(~.ASPF_CT3_RFT, rec(
        aj := (j, a, m) -> fdiv(j+a, m),

        applicable := t -> let(A := ObjId(t.params[1]), n := t.params[1].n,
            A = XN_skew and (IsSymbolic(n) or (n > 4 and IsEvenInt(n)))),

        apply := (self, t, C, Nonterms) >> let(
            N := t.params[1].n,         k := div(Cols(C[1]),2),
            j := Nonterms[3].params[1], m := div(N, 2*k),
            When(t.params[3].dim()=1,
                 Tr(m, 2*k)   * IDirSum(j, C[1]) * C[2],
                 RC(Tr(m, k)) * IDirSum(j, C[1]) * C[2]))
    ))
));

NewRulesFor(ASPF, rec(
    ASPF_CT1_URFT_LftInplace := CopyFields(ASPF_Breakdown_Rule, rec(
        switch := false,

        applicable := t -> let(A := ObjId(t.params[1]), n := t.params[1].n,
            A in [XN_min_1, XN_min_1U] and
            (IsSymbolic(n) or (n > 2 and IsEvenInt(n))) and
            ObjId(t.params[3]) = Freq_1),

        freedoms := (self, t) >> let(n := t.params[1].n, maxR := self.getA("maxRadix", -1),
            [ Filtered(DivisorsIntNonTrivSym(div(n,2)), x -> (maxR < 0) or (x <= maxR/2)) ]),

        child := ASPF_CT1_URFT.child,

        apply := (self, t, C, Nonterms) >> let(
            N := t.params[1].n, k := Cols(C[1])/2, j := Nonterms[4].params[1],
	    m := div(N, 2*k),
            cplx := ObjId(t.params[3])=Freq_1,
	    jj := Ind( Int(j.range/2) ),
	    pp := When(IsOddInt(m), DirectSum(I(1), condM(m-1, (m-1)/2)),
		                    condMp(m, m/2)),
            # NOTE:  add a GT_IJ' transform (ie with IJ' on either side), and an inplace rule for it
            #        also  GT_IJ transforms (ie with IJ on either side), and an inplace rule for it
	    When(not cplx, Error("not implemented"), #  RC(condKp(div(n,2), k))
                Inplace(
		    DirectSum(
		        C[1],
		        When(IsOddInt(m), [], Data(j, V( (j.range-1)/2), C[2])),
		        When(jj.range = 0, [],
			    IDirSum(jj,
			        condM(4*k, 2)*L(4*k,2*k) *
			        DirectSum(Data(j, jj, C[2]), Data(j, j.range-1-jj, C[2])))))
		    ^ (L(N, m) * Tensor(I(2*k), pp))) *
		Grp(L(N, 2*k) * C[3])))
    ))
));

#
# Base cases
#
NewRulesFor(ASPF, rec(
    ASPF_Base2 := rec(
        requiredFirstTag := ANoTag,
        forTransposition := true,
        applicable := t -> t.params[1].n = 2,
        apply := (t, C, Nonterms) -> t.terminate()
    ),

    # This rules provides the base case for any complex transform (ie for any of the algebras)
    # via the real transform.
    #
    # The trouble here, is that in different cases there are some 1d spectral components
    # in the real transform already, and that makes permutations different.
    #
    ASPF_SmallCpx := rec(
        requiredFirstTag := ANoTag,
        forTransposition := true,
        a := rec(maxSize := 512),
        applicable := (self, t) >> ObjId(t.params[3]) = Freq_1 and t.params[1].n <= self.getA("maxRadix", -1),
        children := t -> let(F := t.params[3], s := F.scale1d,
            [[ ASPF(t.params[1], t.params[2], Freq_E(s, s)) ]]),

        apply := (t, C, Nonterms) -> let(
            A := t.params[1], n := EvalScalar(A.n),
            bcK := nn -> K(nn, 2) * Tensor(I(nn/2), Mat([[1, E(4)], [1, -E(4)]])),
            bcL := nn -> L(nn, 2) * Tensor(I(nn/2), Mat([[1, E(4)], [1, -E(4)]])),
                 # x^n-1
            Cond(ObjId(A) = XN_min_1 and IsEvenInt(n),
		   When(n=2, I(2), DirectSum(I(1), Z(n/2, 1), I(n/2-1)) * DirectSum(I(2), bcK(n-2))) * C[1],
                 ObjId(A) = XN_min_1 and IsOddInt(n),
                   DirectSum(I(1), bcK(n-1)) * C[1],
                 # x^n+1
                 ObjId(A) = XN_plus_1 and IsEvenInt(n),
                   bcK(n) * C[1],
                 ObjId(A) = XN_plus_1 and IsOddInt(n),
                   DirectSum(I((n-1)/2), Z((n+1)/2, -1)) * DirectSum(bcK(n-1), I(1)) * C[1],
                 # skew
                 ObjId(A) = XN_skew,  bcL(n) * C[1]
            ))
    ),

    ASPF_NonSkew_Base_VecN := rec(
       requiredFirstTag := [AVecReg, AVecRegCx],
       forTransposition := false,
       applicable := t -> let(v:=t.firstTag().v, alg := t.params[1],
	   ObjId(alg) in [XN_min_1, XN_min_1U, XN_plus_1] and
	   2 <= alg.n and alg.n <= 2*4*v
       ),
       apply := (t, C, Nonterms) -> VectorizedMatSPL(t.firstTag().isa, t)
    ),
    # need a separate rule because VectorizedMatSPL returns a non-transposeable result,
    ASPF_NonSkew_Base_VecN_tr := rec(
       requiredFirstTag := [AVecReg, AVecRegCx],
       forTransposition := false,
       transposed := true,
       applicable := t -> let(v:=t.firstTag().v, alg := t.params[1],
	   ObjId(alg) in [XN_min_1, XN_min_1U, XN_plus_1] and
	   2 <= alg.n and alg.n <= 2*4*v
       ),
       apply := (t, C, Nonterms) -> VectorizedMatSPL(t.firstTag().isa, t)
    ),

    ASPF_BRDFT3_Base4 := rec(
        requiredFirstTag := ANoTag,
        forTransposition := true,
        applicable := t -> t.hashAs() = ASP.bRDFT3(4,1/16),  # skew versionm but rule is probably invalid up to a permutation
        apply := (t, C, Nonterms) -> let(
	   rot := t.params[1].rot, a := t.params[1].a, D := Dat1d(TReal, 2),
           s := t.params[3].scale2d, dd := When(s=1, I(2), Diag(s, 1)),
	   Data(D, fPrecompute(FList(TReal, [2*s*cospi(rot*a), s*(4*cospi(rot*a)^2 - 1)])),
           L(4,2) * VStack(
	       F(2) * dd * Mat([[1,   0,   -1,        0    ],
		                [0,   0,    0,    -nth(D,0)]]),
	       F(2)   *    Mat([[0,   s,    0,     nth(D,1)],
		                [0,   0,  nth(D,0),   0    ]]))))
    ),

    ASPF_URDFT_Base4 := rec(
        requiredFirstTag := ANoTag,
        forTransposition := true,
        applicable := t -> t.hashAs() = ASP.URDFT(4),
        apply := (t, C, Nonterms) -> let(rot := EvalScalar(t.params[1].rot) mod 4,
            #NoPull
            BB(
                Cond(rot=1, I(4), Diag(1,1,1,-1)) *
                Tensor(t.params[3].scale2d * F(2), I(2))))
    ),

    ASPF_URDFT_Base4_Vec2 := rec(
        requiredFirstTag := [AVecReg, AVecRegCx],
        forTransposition := true,
        applicable := t -> t.hashAs() = ASP.URDFT(4).withTags(t.tags) and t.firstTag().v = 2,
        apply := (t, C, Nonterms) ->
            When(EvalScalar(t.params[1].rot)=1,
                VTensor(t.params[3].scale2d * F(2), 2),
                DirectSum(VBase(I(2), 2), VDiag(FList(TReal, [1,-1]), 2)) *
                VTensor(t.params[3].scale2d * F(2), 2))
    ),

    ASPF_RDFT1_Base4 := rec(
        requiredFirstTag := ANoTag,
        forTransposition := true,
        applicable := t -> t.hashAs() = ASP.RDFT(4),
        apply := (t, C, Nonterms) -> let(rot := t.params[1].rot,
            DirectSum(F(2), When(rot=1, I(2), Diag(1,-1))) *
            Tensor(Diag(t.params[3].scale1d, t.params[3].scale2d) * F(2), I(2)))
    ),

    ASPF_RDFT_toPRDFT := rec(
        requiredFirstTag := ANoTag,
        forTransposition := true,
        a := rec(maxSize := false),

        applicable := (self, t) >> let(n:=t.params[1].n,
            not Is2Power(n) and
            (self.getA("maxSize")=false or n <= self.getA("maxSize")) and
            When(IsEvenInt(n), t.hashAs() in [ASP.RDFT(n), ASP.URDFT(n)],
                               t.hashAs() = ASP.RDFT(n))),

        freedoms := (self, t) >> [],
        child := (self, t, fr) >> [ PRDFT(EvalScalar(t.params[1].n), EvalScalar(t.params[1].rot)) ],

        apply := (t, C, Nonterms) -> let(n := EvalScalar(t.params[1].n),
            rdft := Perm_CCS(n) * C[1],
            When(IsEvenInt(n) and t.hashAs() = ASP.URDFT(n),
                 DirectSum((1/2)*F(2), I(n-2)) * rdft,
                 rdft))
    ),

    ASPF_RDFT1_Base4_Vec2 := rec(
        requiredFirstTag := [AVecReg, AVecRegCx],
        forTransposition := true,
        applicable := t -> t.hashAs() = ASP.RDFT(4).withTags(t.tags) and t.firstTag().v = 2,
        apply := (t, C, Nonterms) -> let(rot := t.params[1].rot, vt := TVect(TReal, 2),
            DirectSum(VBlk( [[ vt.value([1, -1]), vt.value([1,1]) ]], 2) * _VVStack([VTensor(I(1), 2), VIxJ2(2)], 2),
                      When(rot=1, VBase(I(2), 2), VDiag(FList(TReal, [1,-1]), 2))) *
            VTensor(Diag(t.params[3].scale1d, t.params[3].scale2d) * F(2), 2))
    ),

    ASPF_RDFT1_toTRDFT := rec(
        requiredFirstTag := [AVecReg, AVecRegCx],
        forTransposition := true,
	a := rec(maxSize := 512),
        applicable := (self, t) >> let(n:=t.params[1].n,
	    t.hashAs() = ASP.RDFT(n).withTags(t.tags) and n <= self.getA("maxSize", 512)),
	freedoms := (self, t) >> [],
	child := (self, t, fr) >> [ TRDFT(t.params[1].n, t.params[1].rot).withTags(t.tags) ],
        apply := (t, C, Nonterms) -> let(
	    n := t.params[1].n,
            s  := Double(t.params[3].scale1d),
            s2 := Double(t.params[3].scale2d),
	    Cond(s=1 and s2=1, C[1],
		 IsEvenInt(n), Diag(diagDirsum(fConst(TReal, 2, s), fConst(TReal, n-2, s2))) * C[1],
		 IsOddInt(n),  Diag(diagDirsum(fConst(TReal, 1, s), fConst(TReal, n-1, s2))) * C[1]
	    )
	)
    ),

    ASPF_rDFT_Base4 := rec(
       requiredFirstTag := ANoTag,
       forTransposition := true,
       applicable := t -> t.hashAs() = ASP.rDFT(4,1/16),
       apply := (t, C, Nonterms) -> let(
	   a := t.params[1].a,
           rot := t.params[1].rot,
           s := t.params[3].scale2d,
           NoPull(           #Diag(1,1,1,-1) *
               Tensor(F(2), I(2)) *
               DirectSum(s*I(2), RCDiag(fPrecompute(FList(TReal, [s*cospi(rot*a), s*sinpi(rot*a)])))) *
               L(4,2)))
    ),

    ASPF_rDFT_toSkewDFT := rec(
        requiredFirstTag := [ANoTag, ATwidOnline],
        forTransposition := true,
        a := rec(maxSize := 512),
        applicable := (self, t) >> let(n:=t.params[1].n,
            IsEvenInt(n) and n >= 4 and n <= self.getA("maxSize", 512) and t.hashAs().setTags([]) = ASP.rDFT(n,1/16)),

        freedoms := (self, t) >> [],
        child := (self, t, fr) >> let(alg := t.params[1],
            [ SkewDFT(alg.n/2, alg.a, alg.rot).withTags(t.getTags()) ]),

        apply := (t, C, Nonterms) -> let(n:=t.params[1].n, s := t.params[3].scale2d,
            RC(s * C[1]) * L(n, n/2))
    ),

    ASPF_rDFT_BaseN := rec(
       requiredFirstTag := ANoTag,
       forTransposition := true,
       a := rec(maxSize := 512),
       applicable := (self, t) >> let(n:=t.params[1].n,
           IsEvenInt(n) and n >= 4 and n <= self.getA("maxSize", 512) and t.hashAs() = ASP.rDFT(n,1/16)),

       freedoms := (self, t) >> [],
       child := (self, t, fr) >> [ DFT(EvalScalar(t.params[1].n/2), EvalScalar(t.params[1].rot)) ],
       apply := (t, C, Nonterms) -> let(
           n := t.params[1].n,
	   a := t.params[1].a,
           rot := t.params[1].rot,
           s := Double(t.params[3].scale2d), # NOTE: if Double() is not used, vector code uses _mm_set1_epi32(..) intead of _ps(..)
           j := Ind(n-2),
           exp := rot*a*(fdiv(4,n)),
           twid := Lambda(j, cond(neq(imod(j, 2),0), s*sinpi(exp*idiv(j+2,2)), s*cospi(exp*idiv(j+2,2)))),
           #NoPull
	   (
               RC(C[1]) *
               # XXX NOTE: Buf below delays sucking in, terrible hack,
               # due to VJamData performance problems, fix it right now
               # NOTE: fPrecompute inside diagDirsum causes twiddles to be vpacked on the fly in the code
               Buf(RCDiag(diagDirsum(fConst(TReal, 1, s), fConst(TReal, 1, 0.0), #FList(TReal, [s,0]),
                                      fPrecompute(twid)))) *
               L(n,n/2)))
    ),

    ASPF_rDHT_BaseN := rec(
 #      requiredFirstTag := ANoTag,
       forTransposition := true,
       a := rec(maxSize := 512),
       applicable := (self, t) >> let(n:=t.params[1].n,
           IsEvenInt(n) and n >= 4 and n <= self.getA("maxSize", 512) and t.setTags([]).hashAs() = ASP.rDHT(n,1/16)),

       freedoms := (self, t) >> [],
       child := (self, t, fr) >> [ DFT(EvalScalar(t.params[1].n/2), EvalScalar(t.params[1].rot)) ],
       apply := (t, C, Nonterms) -> let(
           n := t.params[1].n,
	   a := t.params[1].a,
           rot := t.params[1].rot,
           s := t.params[3].scale2d,
           j := Ind(n-2),
           exp := rot*a*(fdiv(4,n)),
           twid := Lambda(j, cond(neq(imod(j, 2),0), s*sinpi(exp*idiv(j+2,2)), s*cospi(exp*idiv(j+2,2)))),

           NoPull(
               Tensor(I(n/2), Diag(1,-1)) *
               RC(C[1]) *
               DirectSum(s*I(2), RCDiag(fPrecompute(twid))) *
               L(n,n/2) * Diag(diagDirsum(fConst(TReal, n/2, 1), fConst(TReal, n/2, -1)))))
    ),

    ASPF_DHT1_Base4 := rec(
#        requiredFirstTag := ANoTag,
        forTransposition := true,
        applicable := t -> t.setTags([]).hashAs() = ASP.DHT(4),
        apply := (t, C, Nonterms) -> let(rot := t.params[1].rot,
            DirectSum(F(2), F(2)*When(rot=1, I(2), Diag(1,-1))) *
            Tensor(Diag(t.params[3].scale1d, t.params[3].scale2d) * F(2), I(2)))
    ),

    ASPF_UDHT1_Base4 := rec(
#        requiredFirstTag := ANoTag,
        forTransposition := true,
        applicable := t -> t.setTags([]).hashAs() = ASP.UDHT(4),
        apply := (t, C, Nonterms) -> let(rot := t.params[1].rot,
            NoPull(
                DirectSum(I(2), F(2)*When(rot=1, I(2), Diag(1,-1))) *
                Tensor(Diag(t.params[3].scale1d, t.params[3].scale2d) * F(2), I(2))))
    ),

    ASPF_rDHT_Base4 := rec(
       requiredFirstTag := ANoTag,
       forTransposition := true,
       applicable := t -> t.hashAs() = ASP.rDHT(4,1/16),
       apply := (t, C, Nonterms) -> let(
	   a := t.params[1].a,
           rot := t.params[1].rot,
           s := t.params[3].scale2d,
           NoPull(Diag(1,-1,1,-1) *
               Tensor(F(2), I(2)) *
               DirectSum(s*I(2), RCDiag(fPrecompute(FList(TReal, [s*cospi(rot*a), s*sinpi(rot*a)])))) *
               L(4,2)*Diag(1,1,-1,-1)))),

    ASPF_rDFT_Base4_Vec2 := rec(
       requiredFirstTag := [AVecReg, AVecRegCx],
       forTransposition := true,
       applicable := t -> t.hashAs() = ASP.rDFT(4,1/16).withTags(t.tags) and t.firstTag().v = 2,
       freedoms := t -> [],
       child := (t, fr) -> [TL(4,2,1,1).withTags(t.getTags())],
       apply := (t, C, Nonterms) -> let(
	   a := t.params[1].a,
           rot := t.params[1].rot,
           s := t.params[3].scale2d,
           #VDiag(FList(TReal, [1,1,1,-1]), 2) *
           VTensor(F(2), 2) *
           C[1] *
           VRCDiag(fPrecompute(VData(FList(TReal, [s, s*cospi(rot*a), 0, s*sinpi(rot*a)]), 2)), 2))
    ),

    ASPF_rDFT_Base_VecN_Drop := rec(
	requiredFirstTag := [AVecReg, AVecRegCx],
	forTransposition := true,
	applicable := t -> IsEvenInt(t.params[1].n) and t.hashAs() = ASP.rDFT(t.params[1].n, 1/16).withTags(t.tags),

	freedoms := t -> [],
	child := (t, fr) -> [ t.withoutFirstTag() ],
	apply := (t, C, Nonterms) -> C[1]
    ),

    ASPF_rDFT_Base_VecN := rec(
	requiredFirstTag := [AVecReg, AVecRegCx],
	forTransposition := false,
	applicable := t -> IsEvenInt(t.params[1].n) and t.hashAs() = ASP.rDFT(t.params[1].n, 1/16).withTags(t.tags),
	freedoms := t -> [],
	child := (t, fr) -> let(
            n := t.params[1].n,
	    a := t.params[1].a,
            rot := t.params[1].rot,
            s := Double(t.params[3].scale2d), # NOTE: if Double() is not used, vector code uses _mm_set1_epi32(..) intead of _ps(..)
            j := Ind(n/2-1),
            exp := rot*a*(fdiv(4,n)),
            twid := Lambda(j, s*omegapi(exp*(j+1))),
	    [ TConj(
                  TRC(DFT(n/2, rot) *
		      TDiag(fPrecompute(diagDirsum(fConst(TReal, 1, s), twid)))),  # NOTE: fPrecompute inside diagDirsum causes twiddles to be vpacked on the fly in the code
	          fId(n), L(n, n/2)
	      ).withTags(t.tags) ]
        ),
	apply := (t, C, Nonterms) -> C[1]
    ),

    ASPF_rDFT_Base_VecN_tr := rec(
	requiredFirstTag := [AVecReg, AVecRegCx],
	forTransposition := false,
	transposed := true,
	applicable := t -> IsEvenInt(t.params[1].n) and t.hashAs() = ASP.rDFT(t.params[1].n, 1/16).withTags(t.tags).transpose(),
	freedoms := t -> [],
	child := (t, fr) -> let(
            n := t.params[1].n,
	    a := t.params[1].a,
            rot := t.params[1].rot,
            s := Double(t.params[3].scale2d), # NOTE: if Double() is not used, vector code uses _mm_set1_epi32(..) intead of _ps(..)
            j := Ind(n/2-1),
            exp := rot*a*(fdiv(4,n)),
            twid := Lambda(j, s*omegapi(exp*(j+1))),
	    [ TConj(
                  TRC( TDiag(fPrecompute(diagDirsum(fConst(TReal, 1, s), twid))).conjTranspose() *
                       DFT(n/2, rot).conjTranspose()
                       ),  # NOTE: fPrecompute inside diagDirsum causes twiddles to be vpacked on the fly in the code
	          L(n, 2), fId(n)
	      ).withTags(t.tags) ]
        ),
	apply := (t, C, Nonterms) -> C[1]
    ),

    ASPF_Cpx_rDFT_Base_VecN := rec(
	requiredFirstTag := [AVecReg, AVecRegCx],
	forTransposition := true,
	applicable := t -> IsEvenInt(t.params[1].n) and
	    ObjId(t.params[1])=XN_skew and t.params[2]=Time_EX and
	    ObjId(t.params[3])=Freq_1,

	freedoms := t -> [],
	child := (t, fr) -> let(
            n := t.params[1].n,
	    a := t.params[1].a,
            rot := t.params[1].rot,
            s := Double(t.params[3].scale1d), # NOTE: if Double() is not used, vector code uses _mm_set1_epi32(..) intead of _ps(..)
            j := Ind(n/2-1),
            exp := rot*a*(fdiv(4,n)),
            twid := Lambda(j, s*omegapi(exp*(j+1))),
	    mat := Mat([[1, E(4)], [1, -E(4)]]),
	    [ GT(mat, GTVec, GTVec, [n/2]).withTags(t.tags) *
	      TConj(
                  TRC(DFT(n/2, rot) *
		      TDiag(fPrecompute(diagDirsum(fConst(TReal, 1, s), twid)))),  # NOTE: fPrecompute inside diagDirsum causes twiddles to be vpacked on the fly in the code
	          L(n, 2), L(n, n/2)
	      ).withTags(t.tags) ]
        ),
	apply := (t, C, Nonterms) -> C[1]
    ),

#     ASPF_rDFT_Base4_Vec2b := rec(
#        requiredFirstTag := [AVecReg, AVecRegCx],
#        forTransposition := true,
#        applicable := t -> t.hashAs() = ASP.rDFT(4,1/16).withTags(t.tags) and t.firstTag().v = 2,
#        apply := (t, C, Nonterms) -> let(
# 	   a := t.params[1].a,
#            rot := t.params[1].rot,
#            s := t.params[3].scale2d,
#            VTensor(F(2), 2) *
#            DirectSum(VBase(I(2), 2), VBase(J(2), 2)*VDiag(FList(TReal, [1,-1]), 2)) *
# 	   VTensor(Mat([[1,0],[0,1]]),2) *
#            VRCDiag(fPrecompute(VData(FList(TReal, [s, s*sinpi(rot*a), 0, -s*cospi(rot*a)]), 2)), 2))
#     ),

    ASPF_RDFT3_Base4 := rec(
       requiredFirstTag := ANoTag,
       forTransposition := true,
       applicable := t -> t.hashAs() = ASP.RDFT3(4),
       apply := (t, C, Nonterms) -> let(
	   D := Dat1d(TReal, 2), d := Diag(1,-1), s := t.params[3].scale2d,
           rot := t.params[1].rot, pm := s*sinpi(rot*1/2), # pm = +/- s
           NoPull(
	   Data(D, fPrecompute(FList(TReal, [s*cospi(rot*1/4), s*sinpi(rot*1/4) ])),
               L(4,2) *
               DirectSum(F(2)*Diag(1,nth(D,0)), d*F(2)*Diag(1, nth(D,1))) *
               Mat([[s,  0,  0,   0 ],
                    [0,  1,  0,  -1 ],
                    [0,  0,  pm,  0 ],
                    [0,  1,  0,   1 ]]))))
    )
));


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

NewRulesFor(WHT, rec(
	WHT_tSPL_CT_Mem := rec(
		switch := true,
		applicable := (self, nt) >> nt.params[1] > 2
			and not IsPrime(nt.params[1])
			and nt.hasTag(AMem) 						# make sure we're tagged
			and not IsList(nt.getTag(AMem)) 			# make sure we have only ONE memory tag (single level of cache)
			and nt.params[1] > nt.getTag(AMem).size,	# problem size must be larger than the cache

		# atm, just handle a single AMem
		children := nt -> let(
			memtag := nt.getTag(AMem),
			Map2(
				Filtered(DivisorPairs(nt.params[1]), (m,n) -> n mod memtag.block = 0, m <= memblock.assoc),
				(m,n) -> [
					TCompose([
						TTensorI(WHT(m), n, AVec, AVec),
						TTensor(WHT(n), m, APar, APar)
					])
				]
			)
		),

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


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

BACValue:=[
[IntHexString("5601"), 1, 1, 1],
[IntHexString("3401"), 2, 6, 0],
[IntHexString("1801"), 3, 9, 0],
[IntHexString("0AC1"), 4, 12, 0],
[IntHexString("0521"), 5, 29, 0],
[IntHexString("0221"), 38, 33, 0],
[IntHexString("5601"), 7, 6, 1],
[IntHexString("5401"), 8, 14, 0],
[IntHexString("4801"), 9, 14, 0],
[IntHexString("3801"), 10, 14, 0],
[IntHexString("3001"), 11, 17, 0],
[IntHexString("2401"), 12, 18, 0],
[IntHexString("1C01"), 13, 20, 0],
[IntHexString("1601"), 29, 21, 0],
[IntHexString("5601"), 15, 14, 1],
[IntHexString("5401"), 16, 14, 0],
[IntHexString("5101"), 17, 15, 0],
[IntHexString("4801"), 18, 16, 0],
[IntHexString("3801"), 19, 17, 0],
[IntHexString("3401"), 20, 18, 0],
[IntHexString("3001"), 21, 19, 0],
[IntHexString("2801"), 22, 19, 0],
[IntHexString("2401"), 23, 20, 0],
[IntHexString("2201"), 24, 21, 0],
[IntHexString("1C01"), 25, 22, 0],
[IntHexString("1801"), 26, 23, 0],
[IntHexString("1601"), 27, 24, 0],
[IntHexString("1401"), 28, 25, 0],
[IntHexString("1201"), 29, 26, 0],
[IntHexString("1101"), 30, 27, 0],
[IntHexString("0AC1"), 31, 28, 0],
[IntHexString("09C1"), 32, 29, 0],
[IntHexString("08A1"), 33, 30, 0],
[IntHexString("0521"), 34, 31, 0],
[IntHexString("0441"), 35, 32, 0],
[IntHexString("02A1"), 36, 33, 0],
[IntHexString("0221"), 37, 34, 0],
[IntHexString("0141"), 38, 35, 0],
[IntHexString("0111"), 39, 36, 0],
[IntHexString("0085"), 40, 37, 0],
[IntHexString("0049"), 41, 38, 0],
[IntHexString("0025"), 42, 39, 0],
[IntHexString("0015"), 43, 40, 0],
[IntHexString("0009"), 44, 41, 0],
[IntHexString("0005"), 45, 42, 0],
[IntHexString("0001"), 45, 43, 0],
[IntHexString("5601"), 46, 46, 0]];

#index -> Qe, NMPS, NLPS, SWITCH
BACTable := FList(TArray(TInt, 4), [
[IntHexString("5601"), 1, 1, 1],
[IntHexString("3401"), 2, 6, 0],
[IntHexString("1801"), 3, 9, 0],
[IntHexString("0AC1"), 4, 12, 0],
[IntHexString("0521"), 5, 29, 0],
[IntHexString("0221"), 38, 33, 0],
[IntHexString("5601"), 7, 6, 1],
[IntHexString("5401"), 8, 14, 0],
[IntHexString("4801"), 9, 14, 0],
[IntHexString("3801"), 10, 14, 0],
[IntHexString("3001"), 11, 17, 0],
[IntHexString("2401"), 12, 18, 0],
[IntHexString("1C01"), 13, 20, 0],
[IntHexString("1601"), 29, 21, 0],
[IntHexString("5601"), 15, 14, 1],
[IntHexString("5401"), 16, 14, 0],
[IntHexString("5101"), 17, 15, 0],
[IntHexString("4801"), 18, 16, 0],
[IntHexString("3801"), 19, 17, 0],
[IntHexString("3401"), 20, 18, 0],
[IntHexString("3001"), 21, 19, 0],
[IntHexString("2801"), 22, 19, 0],
[IntHexString("2401"), 23, 20, 0],
[IntHexString("2201"), 24, 21, 0],
[IntHexString("1C01"), 25, 22, 0],
[IntHexString("1801"), 26, 23, 0],
[IntHexString("1601"), 27, 24, 0],
[IntHexString("1401"), 28, 25, 0],
[IntHexString("1201"), 29, 26, 0],
[IntHexString("1101"), 30, 27, 0],
[IntHexString("0AC1"), 31, 28, 0],
[IntHexString("09C1"), 32, 29, 0],
[IntHexString("08A1"), 33, 30, 0],
[IntHexString("0521"), 34, 31, 0],
[IntHexString("0441"), 35, 32, 0],
[IntHexString("02A1"), 36, 33, 0],
[IntHexString("0221"), 37, 34, 0],
[IntHexString("0141"), 38, 35, 0],
[IntHexString("0111"), 39, 36, 0],
[IntHexString("0085"), 40, 37, 0],
[IntHexString("0049"), 41, 38, 0],
[IntHexString("0025"), 42, 39, 0],
[IntHexString("0015"), 43, 40, 0],
[IntHexString("0009"), 44, 41, 0],
[IntHexString("0005"), 45, 42, 0],
[IntHexString("0001"), 45, 43, 0],
[IntHexString("5601"), 46, 46, 0]]);


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

Class(OLBase, SumsBase, rec(
    isOL := true,
    visitAs := "OLBase",
));

## Multiplication operator
## OLMultiplication(1, n) is I(n)
## OLMultiplication(2, n) is a point-wise multiplication of two vectors of size n
Class(OLMultiplication, RewritableObject, BaseMat, OLBase, rec(
    dims := self >> [self.rChildren()[2], StripList(Replicate(self.rChildren()[1], self.rChildren()[2]))],
));

## OLConjMultiplication(n) is a point-wise multiplication of two vectors of size n where second vector is complex conjugated
##
Class(OLConjMultiplication, RewritableObject, BaseMat, OLBase, rec(
    dims := self >> [self.rChildren()[2], StripList(Replicate(self.rChildren()[1], self.rChildren()[2]))],
));

#F RCOLMultiplication is RC(OLMultiplication(..))
#F
Class(RCOLMultiplication, RewritableObject, BaseMat, OLBase, rec(
    dims := self >> [2*self.rChildren()[2], StripList(Replicate(self.rChildren()[1], 2*self.rChildren()[2]))],
));

## RCOLConjMultiplication(n) is RC(OLConjMultiplication(n))
##
Class(RCOLConjMultiplication, RewritableObject, BaseMat, OLBase, rec(
    dims := self >> [2*self.rChildren()[2], StripList(Replicate(self.rChildren()[1], 2*self.rChildren()[2]))],
));

Class( 2DI, SumsBase, RewritableObject, AttrMixin, rec(
   isSPL := true,
   isIdentity := True,
   advdims := self >> let(a:=self.rChildren(), [[a], [a]]),
   dims := self >> let(a:=Product(self.rChildren()), [a, a]),
   rng := self >> [ TArray(TUnknown, Product(self.rChildren()))],
   dmn := self >> [ TArray(TUnknown, Product(self.rChildren()))],
   children := self >> [],
   arity := ClassSPL.arity,
   updateParams := meth(self)
       self.func := fId(Product(self.params));
       Inherited();
   end,
));

Class(LeftOver, RewritableObject, ClassSPL, rec(
    isSums := true,
    isSPL := true,

    __call__ := (self, cond, spl) >> Cond(
	cond=true, spl,
	cond=false, let(d := spl.dims(), VirtualPad(d[1], d[2], I(0))),
	Inherited(cond, spl)),

    rng := self >> self.params[2].rng(),
    dmn := self >> self.params[2].dmn(),
    dims := self >> [ StripList(List(self.rng(), (l) -> l.size)), StripList(List(self.dmn(), (l) -> l.size)) ],
    advdims := self >> self.params[2].advdims(),

    isInplace := self >> self.params[2].isInplace(),
    transpose := self >> ObjId(self)(self.params[1], self.params[2].transpose()),
    conjTranspose := self >> ObjId(self)(self.params[1], self.params[2].conjTranspose()),

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

    a := rec(),

    print:= arg >> let(
	self := arg[1],
	i := Cond(Length(arg)>=3, arg[2], 0),
	is := Cond(Length(arg)>=3, arg[3], 4),
	Print(self.__name__, "(", self.params[1], "\n",
	    Blanks(i+is), self.params[2].print(i+is, is), "\n", Blanks(i), ")"))
));

Class( OLDup, RewritableObject, BaseMat, OLBase, rec(
    dims := self >> [Replicate(self.params[1], self.params[2]), self.params[2]],
));

# ParSeqWrap is a helper class used dduring codegen stage to ease code generation
#

Class(ParSeqWrap, BaseContainer, rec(
    __call__ := (self, p, ci, y, x) >>
        WithBases(self, rec(p := p, ci := ci, y := y, x := x, dimensions := [p.dimsCompL(p.child(ci)), p.dimsCompR(p.child(ci))])),
    dims := self >> self.dimensions,
    children := self >> [self.p.child(self.ci)],
));

#  ParSeq(<N>, <spl>, <spl>, ...) operator implements simultaneous sequentional (first <N> inputs/outputs) and
#    parallel (all other inputs/outputs) data flow:
#
#        +----+-<- +
#        |    |
#  * <- A <- B <- *
#        |    |
#  + <---+----+
#
# For example: ParSeq(1, Addition(2,1), Addition(2,1), Addition(2,1));
#              Y[0] := X1[0] + X2[0] + X2[0] + X2[0];

Class(ParSeq, SumsBase, BaseOperation, rec(
    area := self >> Sum(self.children(), x->x.area()),
    abbrevs := [ arg -> Checked( Length(arg)>1 and IsPosInt0(arg[1]),
                                 [ arg[1], Flat(Drop(arg, 1)) ] )
               ],

    # filter list leaving elements with positions of Compose inputs/outputs
    filtCompL := (self, lst) >> lst{[1..self.fb_cnt]},
    filtCompR := (self, lst) >> lst{[1..self.fb_cnt]},
    filtSUML  := (self, lst) >> lst{[self.fb_cnt+1..Length(lst)]},
    filtSUMR  := (self, lst) >> lst{[self.fb_cnt+1..Length(lst)]},

    dimsCompL := (self, child) >> self.filtCompL(Flat([self.dims()[1]])),
    dimsCompR := (self, child) >> self.filtCompR(Flat([self.dims()[2]])),
    dimsSUML  := (self, child) >> self.filtSUML(Flat([self.dims()[1]])),
    dimsSUMR  := (self, child) >> self.filtSUMR(Flat([self.dims()[2]])),

    checkDimsCompose := (self) >> let(chdims := List(self._children, c -> [self.dimsCompL(c), self.dimsCompR(c)]),
        DoForAll([1..Length(chdims)-1], i ->
            DoForAll(Zip2(chdims[i][2], chdims[i+1][1]), x->
                When( not(IsSymbolic(x[1]) or IsSymbolic(x[2])) and (x[1] <> x[2]),
                      Error("Dimensions of children ",i," and ",i+1," do not match (",x[1]," <> ",x[2],") in ", self._children), 0)))),

    checkDimsSUM := (self) >> let(
        chdims := TransposedMat(List(self._children, c -> [self.dimsSUML(c), self.dimsSUMR(c)])),
        dims := List(TransposedMat(chdims[1]) :: TransposedMat(chdims[2]), d -> Set(Filtered(d, e -> not IsSymbolic(e)))),
        When(not ForAll(dims, d -> Length(d) in [0,1]),
            Error("Dimensions of summands do not match"), 0)),

    #-----------------------------------------------------------------------
    new := meth(self, fb_cnt, C)
        local dims, obj, a;
        Constraint(Length(C) >= 1); Constraint(ForAll(C, IsSPL));
        a := C[1].arity();
        Constraint(ForAll(C, e -> e.arity()=a));
        Constraint(IsInt(fb_cnt) and fb_cnt>=0 and fb_cnt <= a[1] and fb_cnt <= a[2]);

        if Length(C) = 1 then return C[1]; fi;
        if fb_cnt=a[1] and a[1]=a[2] then return Compose(C); fi;
        if fb_cnt=0 then return SUM(C); fi;

        # check 'Compose' dims
        obj := SPL(WithBases(self, rec( _children := C, fb_cnt := fb_cnt)));
        obj.checkDimsCompose();
        obj.checkDimsSUM();
        return obj.setDims();
    end,

    rng:= self >> self._children[1].rng(),
    dmn:= self >> Last(self._children).dmn(),

    isPermutation := self >> false,

    from_rChildren := (self, rch) >> ApplyFunc(ObjId(self), [self.fb_cnt] :: rch).appendAobj(self),

    print := (self, i, is) >> self._print([self.fb_cnt] :: self.rChildren(), i, is),

));


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

@_RecursBase := NewMakeRule(
    rec(applicable       := nt -> true,
	forTransposition := false,
	children         := nt -> [[ ]],
	apply := (nt,ch,nonterms) -> RecursStep(nt)),
        "@_RecursBase", @);

#F CreateRecursBaseHash(<smallHash>)
#F
#F Returns new hash table with T->@_RecursBase(T) ruletrees for each T in <smallHash>
#F
CreateRecursBaseHash := function(smallHash)
    local e, f, t, basehash;
    basehash := HashTableDP();
    for e in smallHash.entries do
	if e<>[] then
            for f in e do
              if f.data<>[] then
                  HashAdd(basehash, f.key, When(IsBound(f.data[1].measuredGP),[ rec(ruletree := @_RecursBase(f.key),measured := f.data[1].measuredGP) ],[ rec(ruletree := @_RecursBase(f.key))]));
              fi;
            od;
	fi;
    od;
    return basehash;
end;

Class(Codelet, BaseContainer);

_lookupFormula := function(smallHash, transform, opts)
    local lkup;
    lkup := HashLookup(smallHash, transform);
    if lkup = false
	then Error("Transform ", transform, " missing in opts.libgen.baseBench");
    else
	if not IsBound(lkup[1].sums) then
	    VPrintLine(opts.verbosity, "Sigma-splizing ", lkup[1].ruletree.node, "...");
	    # NOTE: generateComplexCode
	    lkup[1].sums := SumsRuleTree(lkup[1].ruletree, CopyFields(opts, rec(generateComplexCode := true)));
	fi;
	return Copy(lkup[1].sums);
    fi;
end;

Class(RecursStepCall, BaseOperation, rec(
    new := (self, dimensions, func, bindings) >> SPL(WithBases(self,
        rec(dimensions := dimensions,  func := func, bindings := bindings))),

    dims := self >> self.dimensions,
    rChildren := self >> [self.dimensions, self.func, self.bindings],
    rSetChild := rSetChildFields("dimensions", "func", "bindings"),
    numops := self >> 10^300,

    print := (self, i, is) >> Print(
        self.name, "(", self.dimensions, ", \"", self.func, "\", ", self.bindings, ")")
));

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

    RecursStepToCall := Rule(RecursStep,
        meth(e, cx)
            local rs, opts, bindings, func, shape, lkup, params, clrec, smallhash, preimpl;
            opts := cx.opts;
            e := e.child(1);
	    e := SubstTopDown(e, BB, x->x.child(1));
            shape := CodeletShape(e);
            func := CodeletName(shape);
            rs := MkCodelet(e);
            params := CodeletParams(rs);
            bindings := Zip2(params, CodeletParams(e));

            lkup := HashLookup(opts.libgen.codeletTab, shape);
            smallhash := opts.libgen.baseBench.exp.bases.hashTable;
#should it be NoPull(_lookupFormula(...))??
            if lkup = false then
                preimpl := Copy(rs);
                rs := SubstBottomUpRules(rs, [
                    [RTWrap, r -> _lookupFormula(smallhash, r.rt.node, opts) ],
                    [@.cond(IsNonTerminal), r -> _lookupFormula(smallhash, r, opts) ]]);
                rs := ApplyStrategy(rs, opts.libgen.terminateStrategy, UntilDone, opts);
                clrec := rec(name := func, shape := shape, sums := rs, preimpl := preimpl, params := params);
                HashAdd(opts.libgen.codeletTab, shape, clrec);
                HashAdd(opts.libgen.codeletTab, func, clrec);
            fi;
            return RecursStepCall(rs.dimensions, func, bindings);
        end)
));

CreateCodeletHashTable := () ->
    HashTable((val,size)->InternalHash(CodeletName(val)));

CompileCodelets := function(c, opts)
    local calls, call, clet_code, lkup, res, funcs,outputs,inputs, u, t;
    calls := Collect(c, RecursStepCall);
    res := []; funcs := Set([]);
    for call in calls do
        if not call.func in funcs then
	    Add(funcs, call.func);
	    lkup := HashLookup(opts.libgen.codeletTab, call.func);
	    if lkup = false then Error(call.func, " is missing from codelet hash. ",
                                                  "Run SumsRuleTree again?"); fi;
	    if not IsBound(lkup.code) then
		VPrintLine(opts.verbosity, "Compiling ", call.func, "...");
                t := TPtr(Cond(opts.generateComplexCode, TComplex, TReal));
                outputs:=StripList(List([1..Length(Flat([lkup.sums.dims()[1]]))],x->var(Concat("Y",When(x>1,String(x),"")), t)));
                inputs:=StripList(List([1..Length(Flat([lkup.sums.dims()[2]]))],x->var(Concat("X",When(x>1,String(x),"")), t)));
                u:=opts.libgen.basesUnrolling;
                if IsBound(lkup.unrolling) then
                    u:=lkup.unrolling;
                fi;
		lkup.code := opts.codegen(Codelet(lkup.sums), outputs, inputs,
		    CopyFields(opts, rec(globalUnrolling := u)));
	    fi;
            #if Collect(lkup.code, errExp) <> [] then
            #    lkup.badcode := lkup.code; Unbind(lkup.code); Error("errExp"); fi;
	    Add(res, lkup);
	fi;
    od;
    return res;
end;


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

;; This buffer is for notes you don't want to save, and for Lisp evaluation.
;; If you want to create a file, visit that file with C-x C-f,
;; then enter the text in that file's own buffer.

n := 20;

opts := CopyFields(PRDFTDefaults, rec(
    globalUnrolling := 1000,
    compileStrategy := NoCSE,
    compflags := "-O1"
));

r  := RandomRuleTree(BRDFT3(4,1/16), opts);
rr := RandomRuleTree(RC(SkewDFT(2,1/16)), opts);
s  := SumsSPL(Tensor(I(n), SumsRuleTree(r, opts)), opts);
ss := SumsSPL(Tensor(I(n), SumsRuleTree(rr, opts)), opts);
c  := CodeSums(1000, s, opts);;
cc := CodeSums(1000, ss, opts);;
CMeasure(c, opts);
CMeasure(cc, opts);

r  := RandomRuleTree(BRDFT1(4), opts);
rr := RandomRuleTree(SRDFT1(4), opts);
s  := SumsSPL(Tensor(I(n), SumsRuleTree(r, opts)), opts);
ss := SumsSPL(Tensor(I(n), SumsRuleTree(rr, opts)), opts);
c  := CodeSums(1000, s, opts);;
cc := CodeSums(1000, ss, opts);;
CMeasure(c, opts);
CMeasure(cc, opts);

copts := CopyFields(opts, rec(
   generateComplexCode := true,
   XType := TComplex,
   YType := TComplex,
   dataType := "complex",
   unparser := CMacroUnparser,
   includes := ["<include/complex_gcc_sse3.h>"]
));

opts := copts;

r  := RandomRuleTree(Tensor(I(n), BRDFT3(4,1/16)), opts);
rr := RandomRuleTree(Tensor(I(n), RC(SkewDFT(2,1/16))), opts);

c  := CodeRuleTree(r, opts);;
cc := CodeRuleTree(rr, opts);;


text:
ITER_POLY_WARN:=false;
Read(Filename(home_dir,"lib/io.gi"));
Read(Filename(home_dir,"handle.gi"));

sys:=ChevalleyAdj("E",8,GF(5));
alg:=AlgebraicU(sys);
orbs:=UnipotentClasses(alg,"Mix");

labels:=List(AllClasses(orbs),x->Label(x));