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

IsIntSym     := x -> (IsSymbolic(x) and IsOrdT(x.t)) or (IsValue(x) and IsOrdT(x.t)) or IsInt(x);
IsPosIntSym  := x -> (IsSymbolic(x) and IsOrdT(x.t)) or (IsValue(x) and IsOrdT(x.t) and x.v > 0)  or (IsInt(x) and x > 0);
IsPosInt0Sym := x -> (IsSymbolic(x) and IsOrdT(x.t)) or (IsValue(x) and IsOrdT(x.t) and x.v >= 0) or (IsInt(x) and x >= 0);
IsRatSym     := x -> (IsSymbolic(x) and (IsOrdT(x.t) or IsRealT(x.t) or x.t=TUnknown)) or IsRat(x) or (IsValue(x) and (IsOrdT(x.t) or IsRealT(x.t)));
IsBoolSym    := x -> x.t=TBool and (IsSymbolic(x) or IsValue(x));

AnySyms := arg -> ForAny(arg, IsSymbolic);

_ev := x->When(IsRec(x) and IsBound(x.ev), x.ev(), x);

# in functions below
# vec elements could be non-values (i.e. symbolic == expressions)

isValueZero := e ->
    (e.v=0 or IsDouble(e.v) and AbsFloat(e.v)<TReal.cutoff) or
    (IsVecT(e.t) and ForAll(e.v, x->not IsSymbolic(x) and AbsFloat(_ev(x)) < TDouble.cutoff));

isValueOne := e ->
    (e.v=1 or IsDouble(e.v) and AbsFloat(e.v-1)<TReal.cutoff) or
    (IsVecT(e.t) and ForAll(e.v, x->not IsSymbolic(x) and AbsFloat(_ev(x)-1) < TDouble.cutoff));

isValueNegOne := e ->
    (e.v=-1 or IsDouble(e.v) and AbsFloat(e.v+1)<TReal.cutoff) or
    (IsVecT(e.t) and ForAll(e.v, x->not IsSymbolic(x) and AbsFloat(_ev(x)+1) < TReal.cutoff));

_is0none := e -> Cond(IsValue(e), isValueZero(e), ObjId(e)=noneExp or e=0);
_is1 := e -> Cond(IsValue(e), isValueOne(e),  e=1);

# in the patterns below 'e' could be a plan GAP integer or other non-class GAP object
_0    := @(0).cond(e -> Cond(IsValue(e), isValueZero(e), e=0));
_0none:= @(0).cond(_is0none);
_1    := @(0).cond(_is1);
_2    := @(0).cond(x -> x=2);
_neg1 := @(0).cond(e -> Cond(IsValue(e), isValueNegOne(e), e=-1));

_v0    := @(0, Value, isValueZero);
_v0none:= @(0, [Value, noneExp], e -> Cond(ObjId(e)=noneExp, true, isValueZero(e)));
_v1    := @(0, Value, isValueOne);
_vneg1 := @(0, Value, isValueNegOne);

_vtrivial := @(0, Value, e->isValueZero(e) or isValueOne(e) or isValueNegOne(e));

_vtrue    := @(0, Value, e -> e.v=true);
_vfalse    := @(0, Value, e -> e.v=false);

Declare(_divides_rec);

_divides_rec := (d,n) -> Cond(
    ObjId(n) = add and ForAll( n.args, e -> _divides_rec(d, e) ), true,
    ObjId(n) = mul and ForAny( n.args, e -> _divides_rec(d, e) ), true,
    IsSymbolic(d) or IsSymbolic(n), d=n or n=0,
    (EvalScalar(n) mod EvalScalar(d)) = 0);

_divides := (d,n) -> CondPat( d,
    [mul, param, param, ...], # product of distinct params
        ForAll( d.args, p -> IsParam(p) and _divides_rec(p,n) ) and Length(Set(d.args))=Length(d.args),
    [mul, Value, param],
        _divides_rec(d.args[1], n) and _divides_rec(d.args[2], n),
    # else
    _divides_rec(d,n));

_isEven := x -> _divides(2, x);

# We can't know for sure that d|n, because _dividesUnsafe returns true if either d or n is a variable
_dividesUnsafe := (d,n) -> Cond(IsSymbolic(d) or IsSymbolic(n), true,
    (EvalScalar(n) mod EvalScalar(d)) = 0);

_unwrap := n -> Cond(IsValue(n), n.v, n);
Declare(_stripval);
_stripval := n -> Cond( IsValue(n), _stripval(n.v),
                        IsList(n),  List(n, _stripval),
                        n );

_unwrapV := (n) -> let(_n := EvalScalar(n), Cond(IsValue(_n), n.v, _n));


text:
#
# walrus: Computational Methods for Finitely Generated Monoids and Groups
#

## Presentations
InstallGlobalFunction(NewPregroupPresentation,
function(pg, rels)
    local res;

    res := rec( pg := pg );
    res.rels := List([1..Length(rels)], x -> NewPregroupRelator(res, rels[x], x));
    # Orgh. Make sure Locations and Places are all computed for the moment
    Objectify(PregroupPresentationType, res);
    Locations(res);
    Places(res);
    SetVertexTripleCache(res, HashMap());

    return res;
end);

InstallGlobalFunction(PregroupPresentationFromFp,
function(F, rred, rgreen)
    local inv, pg, rgreens, inv_test;

    inv_test := function(x)
        local lr;
        if Length(x) = 2 then
            lr := LetterRepAssocWord(x);
            return lr[1] = lr[2];
        fi;
        return false;
    end;

    # Get the involutions
    inv := Filtered(rgreen, inv_test);
    rgreens := Filtered(rgreen, x -> not inv_test(x));
    pg := PregroupByRedRelators(F, rred, List(inv, x -> LetterRepAssocWord(x)[1]));

    rgreens := List(rgreens, pg!.convert_word);

    return NewPregroupPresentation(pg, rgreens);
end);

InstallMethod(ViewString
             , "for a pregroup presentation"
             , [IsPregroupPresentationRep],
function(pgp)
    # Note that we do not really regard 1 as a generator
    local res;

    return STRINGIFY("<pregroup presentation with "
                    , Size(pgp!.pg)-1, " generators and "
                    , Length(pgp!.rels), " relators>");
end);

InstallMethod(Pregroup
             , "for a pregroup presentation"
             , [IsPregroupPresentation],
             pgp -> pgp!.pg);

InstallMethod(Generators
             , "for a pregroup presentation"
             , [IsPregroupPresentation],
function(pres)
    local gens;

    gens := List(Pregroup(pres), x->x);
    Remove(gens, 1);
    return gens;
end);

InstallMethod(Relators
             , "for a pregroup presentation"
             , [IsPregroupPresentation],
             pgp -> pgp!.rels);

InstallMethod(RelatorsAndInverses
             , "for a pregroup presentation"
             , [IsPregroupPresentation],
function(pgp)
    local r;
    r := Relators(pgp);
    return Concatenation(r, List(r, Inverse));
end);

InstallMethod(Powers
             , "for a pregroup presentation"
             , [IsPregroupPresentation],
             x -> x!.powers);

InstallMethod(RLetters
             , "for a pregroup presentation"
             , [IsPregroupPresentation],
function(pres)
    local rels, gens, x, r, rlett;

    rels := RelatorsAndInverses(pres);
    gens := Generators(pres);
    rlett := Set([]);
    for r in rels do
        for x in gens do
            if x in r then AddSet(rlett, x); fi;
        od;
    od;
    return rlett;
end);

# An R-letter is a letter that occurs in any (interleave) of
# a relation (Definition 7.4)
# XXX: Note that the code below does not do interleaves yet!
InstallGlobalFunction(IsRLetter,
function(pres, x)
    # determine whether x occurs in I(R)
    return x in RLetters(pres);
end);

InstallMethod(Locations, "for a pregroup presentation",
              [IsPregroupPresentation and IsPregroupPresentationRep],
function(pres)
    local rel, locs, w, i, index;

    locs := [];

    for rel in RelatorsAndInverses(pres) do
        Append(locs, Locations(rel));
    od;

    # Hack
    # In particular this ID is ID within presentation
    index := HashMap(x -> HashBasic([__ID(x[1]), __ID(x[2])]));
    for i in [1..Length(locs)] do
        locs[i]![6] := i;

        if IsBound( index[ [ InLetter(locs[i]), OutLetter(locs[i]) ] ] ) then
            Add( index[ [ InLetter(locs[i]), OutLetter(locs[i]) ] ], locs[i] );
        else
            index[ [ InLetter(locs[i]), OutLetter(locs[i]) ] ] := [ locs[i] ];
        fi;
    od;
    SetLocationIndex(pres, index);
    return locs;
end);

# For a location R(i,a,b) find a matching location
# R'(j,b^(-1),c) such that a reduced diagram only
# containing R and R' exists
FindMatchingLocation := function(loc, c)
    local locs, b, binv, loc2;

    b := OutLetter(loc);
    binv := PregroupInverse(b);
    locs := LocationIndex(PregroupPresentationOf(loc));
    locs := locs[ [binv, c] ];

    if locs <> fail then
        for loc2 in locs do
            if CheckReducedDiagram(loc, loc2) then
                return loc2;
            fi;
        od;
    fi;
    return fail;
end;

# Given a pregroup presentation as the input, find all places
# B is true if a place will always occur on the boundary
InstallMethod(Places, "for a pregroup presentation",
              [IsPregroupPresentation and IsPregroupPresentationRep],
function(pres)
    local loc, loc2,
          places, a, b, c,
          rels,
          locs, i;

    # rels := RelatorsAndInverses(pres);
    rels := Relators(pres);
    locs := Locations(pres);
    places := [];

    for loc in locs do
        a := InLetter(loc);
        b := OutLetter(loc);
        for c in Generators(pres) do
            # Colour is "red"
            if IsIntermultPair(PregroupInverse(b), c) then
                Add(places, NewPlace(loc, c, "red"));
            fi;
            # Colour is "green"
            # find location R'(j,b^(-1),c) and check that a diagram that
            # meets at the edge b doesn't collapse.

            if FindMatchingLocation(loc, c) <> fail then
                Add(places, NewPlace(loc, c, "green"));
            fi;
        od;
    od;
    # Hack
    for i in [1..Length(places)] do
        places[i]![5] := i;
    od;
    return places;
end);

InstallMethod(LengthLongestRelator, "for a pregroup presentation",
              [IsPregroupPresentation and IsPregroupPresentationRep],
function(pres)
    return Maximum(List(Relators(pres), Length));
end);

# Definition 3.3: A diagram is semi-reduced, if no distinct adjacent faces
# are labelled by ww_1 and w_1^{-1}w for a relator ww_1 and have a common
# consolidated edge labelled by w and w^-1
# it is reduced if this also holds for a face incident with itself.
#
# test whether label on Relator(l2) beginning at b^(-1) is equal to inverse
# of label on Relator(l1) ending at b
InstallGlobalFunction(CheckReducedDiagram,
function(l1, l2)
    local i, j, r1, r2, iend, jend;

    r1 := Relator(l1);
    i := Position(l1);
    iend := i + 1;
    if iend > Length(r1) then iend := iend - Length(r1); fi;

    r2 := Relator(l2);
    j := Position(l2) - 1;
    jend := j - 1;
    if jend < 0 then jend := jend + Length(r2); fi;

    # Here inverses should match, or otherwise the passed
    # locations are already incompatible
    if r1[i] <> PregroupInverse(r2[j]) then
        Error("loc1 and loc2 are not compatible");
        return fail;
    fi;

    repeat
        i := i - 1;
        if i = 0 then i := Length(r1); fi;
        j := j + 1;
        if j > Length(r2) then j := 1; fi;

#        # Print("r1[", i, "] = ", r1[i], " r2[", j, "] = ", r2[j], "\n");
        if r1[i] <> PregroupInverse(r2[j]) then
            return true;
        fi;
    until (i = iend) or (j = jend);

    return false;
end);

InstallGlobalFunction(PregroupPresentationToFpGroup,
function(pres)
    local d, F, rels, gens, pg, p, i, j;

    # By Definition 2.2 of the paper
    d := Size(Pregroup(pres)) - 1;
    F := FreeGroup(d);
    pg := Pregroup(pres);

    rels := [];
    # rels := List([1..d], x -> [ F.(x) * F.(pg!.inv(x+1)-1), One(F) ]);
    for i in [2..d+1] do
        for j in [2..d+1] do
            p := pg[i] * pg[j];
            if p <> fail then
                if p = One(pg) then
                    Add(rels, [F.(i-1) * F.(j-1), One(F)]);
                else
                    Add(rels, [F.(i-1) * F.(j-1) * F.(__ID(PregroupInverse(p))-1), One(F)]);
                fi;
            fi;
        od;
    od;
    Append(rels, List(Relators(pres), x -> [Product(List(List(x), y -> F.(__ID(y)-1))), One(F)]));

    return F / rels;
end);

#T Put pregroup relations in (if pregroup is not the pregroup of the free group)
#T Choose sensible generator names and print them (maybe just use x1,x2,... or a,b,c))
if WALRUS_kbmag_available then

InstallGlobalFunction(PregroupPresentationToKBMAG,
    pres -> KBMAGRewritingSystem(PregroupPresentationToFpGroup(pres)));

fi;

# Writes out pregroup as a multiplication table
# and then the relators, intended as a simple exchange
# format
InstallGlobalFunction(PregroupPresentationToSimpleStream,
function(stream, pres)
    local row, col, table, n, i, rel;

    table := MultiplicationTableIDs(Pregroup(pres));
    n := Length(table);
    AppendTo(stream, "pregroup:\n" );

    for row in [1..n] do
        AppendTo(stream, "  ");
        for col in [1..n-1] do
            AppendTo(stream, String(table[row][col]));
            AppendTo(stream, " ");
        od;
        AppendTo(stream, String(table[row][n]));
        AppendTo(stream, "\n");
    od;

    AppendTo(stream, "relators:\n");
    for rel in Relators(pres) do
        AppendTo(stream, "  ");
        for i in [1..Length(rel)-1] do
            AppendTo(stream, String(__ID(rel[i])));
            AppendTo(stream, " ");
        od;
        AppendTo(stream, String(__ID(rel[Length(rel)])));
        AppendTo(stream, "\n");
    od;
end);

# Serialises a pregroup presentation as a GAP record
# that GAP can read back easily
InstallGlobalFunction(PregroupPresentationToStream,
function(stream, pres)
    local res, rel, rels;

    res := rec( table := MultiplicationTableIDs(Pregroup(pres)),
                rels := List(Relators(pres), r -> List(r, x -> __ID(x))));

    PrintTo(stream, res, ";\n");
end);

InstallGlobalFunction(PregroupPresentationFromStream,
function(stream)
    local res, r, pg, rels;
    res := READ_ALL_COMMANDS(stream, false, false, IdFunc);

    if Length(res) = 1
       and res[1][1] = true
       and IsBound(res[1][2])
       and IsRecord(res[1][2]) then
        # FIXME: Store element names, at least don't put arbitarary
        #        restriction o nnumber of pregroup elts here
        r := res[1][2];
        pg := PregroupByTable("1abcdefghijklmnopqrstuvwxyz"{[1..Length(r.table)]}, r.table);
        return NewPregroupPresentation(pg, List(r.rels, x -> pg_word(pg, x)));
    else
        Error("Could not read pregroup presentation from stream");
    fi;
end);

InstallGlobalFunction(PregroupPresentationToFile,
function(filename, pres)
    local outs;

    outs := OutputTextFile(filename, false);
    SetPrintFormattingStatus(outs, false);
    PregroupPresentationToStream(outs, pres);
    CloseStream(outs);
end);

InstallGlobalFunction(PregroupPresentationToSimpleFile,
function(filename, pres)
    local outs;

    outs := OutputTextFile(filename, false);
    SetPrintFormattingStatus(outs, false);
    PregroupPresentationToSimpleStream(outs, pres);
    CloseStream(outs);
end);

ext: gi

max_stars_repo_name: RussWoodroofe/walrus

max_issues_repo_name: RussWoodroofe/walrus

max_forks_repo_name: RussWoodroofe/walrus


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

_orderedUniquify := l-> Flat([l[1]]::List([2..Length(l)], i->When(l[i] in l{[1..i-1]}, [], [l[i]])));

Class(FFTXGenMixin, rec(
    search := (self, t) >> When(IsBound(self.useDP) and self.useDP, DP(t, When(IsBound(self.dpopts), self.dpopts, rec()), self)[1].ruletree, RuleTreeMid(t, self)),
    preProcess := (self, t) >> let(t1 := RulesFFTXPromoteNT(Copy(t)), RulesFFTXPromoteNT_Cleanup(t1)),

    codeSumsCPU := meth(self, ss)
                        local opts2, ss2, c;

                        opts2 := Copy(SpiralDefaults);
                        opts2.useDeref := false;
#                        opts2.codegen := VecRecCodegen;
#                        opts2.unparser := SSEUnparser;

                        ss2 := Copy(ss);

                        ss2 := SubstBottomUp(ss2, @(1, SIMTISum),
                            e -> ISum(@(1).val.var, @(1).val.domain, @(1).val.child(1)));

                        self.debug.ss_cpu := ss2;

                        c := CodeSums(ss2, opts2);
                        c := fixReplicatedData(c, opts2);

                        self.cpu_opts := opts2;
                        c.dimensions := ss.dims();

                        return c;
                    end,

    prettyPrintCPU := meth(self, c)
                       local name;
                       name := "fftx_generated";
                       if (IsBound(c.ruletree) and IsBound(c.ruletree.node.params[2].fname)) then
                           name := c.ruletree.node.params[2].fname;
                       fi;
                       PrintCode(name, c, self.cpu_opts);
                   end,

    fftxGenCPU := meth(self, t)
                   local tt, rt, c, s, r;
				   if Length(t.params) > 1 then
				       r := t.params[2];
					   if IsRec(r) and IsBound(r.fname) then
					      self.cudasubName := r.fname;
					   fi;
				   fi;

                   self.debug := rec();
                   tt := self.preProcess(t);
                   self.debug.tt := tt;
                   rt := self.search(tt);
                   self.debug.rt := rt;
                   s := self.sumsRuleTree(rt);
                   self.debug.ss := s;

                   c := self.codeSumsCPU(s);
                   c := SubstBottomUp(c, @(1,decl), e->let(fr := @(1).val.cmd.free(), decl(Filtered(e.vars, v->v in fr), @(1).val.cmd)));

                   self.debug.c := c;
                   c.ruletree := rt;
                   c.dimensions := s.dims();

                   return c;
               end,

    cmeasureCPU := (self, c) >> CMeasure(c, self.cpu_opts),

    codeSums := meth(self, ss)
                        local c, cc, Xptr, Yptr, plist, plist1, tags;

                        if IsBound(self.Xptr) then
                            Xptr := self.Xptr;
                        else
                            Xptr := X;
                        fi;
                        if IsBound(self.Yptr) then
                            Yptr := self.Yptr;
                        else
                            Yptr := Y;
                        fi;

                        c := CodeSums(ss, self);
                        tags := rec();
                        if IsBound(c.dimensions) then tags.dimensions  := c.dimensions; fi;
                        if IsBound(c.ruletree) then tags.ruletree  := c.ruletree; fi;

                        c := SubstVars(c, rec(X := Xptr, Y := Yptr));
                        cc := fixReplicatedData(c, self);
                        cc := Rewrite(cc, RulesDropSymDecl, self);

                        plist := Collect(cc, @(1,func, e->e.id="transform"))[1].params;
                        plist1 := _orderedUniquify(plist);
                        if Length(plist1) = 2 then
                            plist1 := plist1 :: self.symbol;
                        fi;
                        cc := SubstTopDown(cc, @(1).cond(e->IsList(e) and e = plist), e->plist1);

                        cc := fixReplicatedData(cc, self);
                        cc := SubstTopDown(cc, @(1, Value, e->e.t = TReal and IsDouble(e.v) and let(delta := AbsFloat(e.v - IntDouble(e.v)), delta >0 and delta  < 0.00000000001 )),
                            e->Value(TReal, IntDouble(@(1).val.v))
                        );

                        if IsBound(self.postProcessCode) then
                            cc := self.postProcessCode(cc, self);
                        fi;

                        cc := CopyFields(tags, cc);
                        return cc;
                    end,

    sumsRuleTree := meth(self, rt)
                        local  s, s2, ss, Xptr, Yptr, tags;

                        self.symbol := When(IsBound(rt.node.params[2].params), rt.node.params[2].params, []);

                        if IsBound(rt.node.params[2].XType) then
                            Xptr := var("Xptr", rt.node.params[2].XType);
                        else
                            Xptr := X;
                        fi;
                        if IsBound(rt.node.params[2].YType) then
                            Yptr := var("Yptr", rt.node.params[2].YType);
                        else
                            Yptr := Y;
                        fi;

                        self.Xptr := Xptr;
                        self.Yptr := Yptr;

                        s := SumsRuleTree(rt, self);
                        tags := rec(ruletree := rec());
                        if IsBound(s.ruletree) then tags.ruletree  := s.ruletree; fi;

                        s2 := fixUpSigmaSPL(s, self);

                        ss := SubstVars(s2, rec(X := Xptr, Y := Yptr));
                        ss := SubstTopDown(ss, @(1, [NoDiagPullinLeft, NoDiagPullinRight, Grp]), e->e.child(1));
                        if IsBound(self.postProcessSums) then ss := self.postProcessSums(ss, self); fi;
                        ss.ruletree := tags.ruletree;

                        return ss;
                    end,

    prettyPrint := meth(self, c)
                       local name;
                       name := "fftx_generated";
                       if (IsBound(c.ruletree) and IsBound(c.ruletree.node.params[2].fname)) then
                           name := c.ruletree.node.params[2].fname;
                       fi;
                       PrintCode(name, c, self);
                   end,

    fftxGen := meth(self, t)
                   local tt, rt, c, s, r;
				   if Length(t.params) > 1 then
				       r := t.params[2];
					   if IsRec(r) and IsBound(r.fname) then
					      self.cudasubName := r.fname;
					   fi;
				   fi;
                   tt := self.preProcess(t);
                   rt := self.search(tt);
                   s := self.sumsRuleTree(rt);
                   c := self.codeSums(s);
                   return c;
               end
));

_tofAdd := (n, l) -> When(IsList(l), let(mn := Minimum(l), mx := Maximum(l), Checked(ForAll([mn..mx], i->i in l), fAdd(n, mx-mn+1, mn))),
                          When(n <> l, fAdd(n, l, 0), fId(n)));
_toBox := (ns, nps)-> When(IsList(ns), fTensor(List(Zip2(ns, nps), i->ApplyFunc(_tofAdd, i))), _tofAdd(ns, nps));

ZeroEmbedBox := (ns, nps) -> Scat(_toBox(ns, nps));
ExtractBox := (ns, nps) -> Gath(_toBox(ns, nps));
BoxND := (l, stype) -> Cond(IsInt(l), TArray(stype, l),
                             IsList(l) and Length(l) = 1, TArray(stype, l[1]),
                             TArray(BoxND(Drop(l, 1), stype), l[1]));

BoxNDcmaj := (l, stype) -> TColMaj(BoxND(l, stype));

fBox := l -> fTensor(List(l, fId));
fBoxInBox := (l2, l1, l3) ->fTensor(List([1..Length(l1)], _l -> fAdd(l1[_l], l2[_l], let(shft := l3[_l], When(IsInt(shft), shft, IntDouble(shft))))));

#_cplx := (r,i) -> r + E(4)*i;
#_imag := i -> E(4) * i;
_real := r -> r;
M_PI := d_acos(-1);

_cplx := (r,i) -> cxpack(r,i);
_imag := i -> cxpack(0, i);

TMapPar := (krn, dims) -> When(Length(dims) = 0, krn, let(tti := When(IsVar(dims[1]), TTensorInd, TTensorI), ApplyFunc(tti, [TMapPar(krn, Drop(dims, 1)), dims[1], APar, APar])));
TMap := (krn, dims, al, ar) -> TCompose(
    When(al = AVec, [TL(Product(List(dims, i->i.range)) * Rows(krn), Rows(krn), 1, 1)], []) ::
    [ TMapPar(krn, dims)] ::
    When(ar = AVec, [TL(Product(List(dims, i->i.range)) * Cols(krn), Product(List(dims, i->i.range)), 1, 1)], []));

lin_idx := arg -> fTensor(List(arg, fBase)).at(0);

TIterator1D := (A, n, l,  r) -> When(n = 1, A, TTensorI(A, n, l, r));

_gflops := (n, b, t) -> (b*5*n*LogInt(n, 2))/(t*1000000);

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:
##############################################################
##
##           Helps compute the decomposition of
##           the ramification module and the equivariant
##           degree of a multiple of a simple orbit for
##           the modular curves X(p) with automorphism group
##           G = PSL(2,p), p a prime.
##
##           Using these computations and Borne's formula, the
##           G-module structure of the RR spaces of equivariant
##           divisors can be determined explicitly.
##
## Reference: D. Joyner and A. Ksir, "Modular representations
##            on some Riemann-Roch spaces of modular curves
##            $X(N)$, Computational Aspects of Algebraic Curves,
##            (Editor: T. Shaska) Lecture Notes in Computing,
##            WorldScientific, 2005.)
##
##
## 12-30-2004, wdj
###########################################################

## To read this in, use for example:
### Read("/home/wdj/gapfiles/curves/mod_crv_aut_gp3.gap");
## To log your results, use for example:
### LogTo("/home/wdj/gapfiles/mod_crv_aut_gp1.log");

ram_module_X:=function(p)
## p is a prime
## output is [m1,...,mn]
## where n = # conj classes of G=PSL(2,p)
## and mi = mult of pi_i in ram mod of modular curve
## X with AutGp(X) = G.
## Here Irr(G) = [pi_1,...,pi_n] (in that order).
##
local G,i,j,n,n0,H,G1,H_chars,CG,G_chars,w,m,theta,pi_theta_to_the;
  G:=PSL(2,p);
  H:=[];
  H_chars:=[];
  n0:=[ [ Z(p)^0, Z(p) ], [ Z(p)*0, Z(p)^0 ] ];
  CG:=ConjugacyClassesSubgroups(G);
  H[1]:=Representative(CG[2]); # size 2
  H[2]:=Representative(CG[3]); # size 3
  n :=Size(CG);
  for i in [1..n] do
    if Size(Representative(CG[i]))=p then
      H[3]:=Representative(CG[i]);
    fi;
  od;
## H[3]:=Group(n0);       # size p
  for i in [1..Size(H)] do
   H_chars[i]:=Irr(H[i]);
  od;
  G_chars:=Irr(G);
  m:=[];
  m[1]:=[];m[2]:=[];m[3]:=[];
  theta:=List([1..3],i->H_chars[i][2]);
  pi_theta_to_the:=List([1..3],i->Sum([1..(Size(H[i])-1)],
      j->j*InducedClassFunction(theta[i]^(j),G)));
#Print("\n\n pi_theta_to_the = ",Sum(pi_theta_to_the),"\n\n");
  for i in [1..3] do
    m[i]:=List(G_chars, pi->ScalarProduct(pi_theta_to_the[i],pi))/Size(H[i]);
  od;
#Print("\n\n m = ",m,"\n\n");
return Sum(m);
end;

#ram_module_X(5);
#[ 0, 3, 3, 4, 5 ]
#ram_module_X(7);
#[ 0, 4, 3, 6, 7, 8 ]
#ram_module_X(11);
#[ 0, 5, 6, 11, 10, 12, 12, 12 ]
#ram_module_X(13);
#[ 0, 7, 7, 13, 13, 13, 13, 14, 15 ]
#ram_module_X(17);
#[ 0, 9, 9, 18, 17, 17, 17, 18, 18, 19, 19 ]
#ram_module_X(19);
#[ 0, 9, 10, 20, 20, 19, 19, 20, 20, 21, 21, 21 ]
#ram_module_X(23);
#[ 0, 14, 11, 24, 24, 24, 23, 23, 25, 25, 25, 25, 25, 25 ]
#ram_module_X(29);
#[ 0, 16, 16, 30, 31, 31, 30, 30, 30, 30, 31, 32, 32, 32, 31, 31, 31 ]
#ram_module_X(31);
#[ 0, 18, 15, 33, 33, 33, 32, 32, 32, 32, 33, 34, 33, 33, 34, 34, 34, 34 ]
#ram_module_X(37);
#[ 0, 20, 20, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 40, 39, 41, 41, 41, 40, 40, 40 ]

ram_module_X_JK:=function(p)
## p is a prime
## output is (m1,...,mn)
## where n = # conj classes of G=PSL(2,p)
## and mi = "mult of pi_i in ram mod of G" using JK formula
##
local G,i,j,n,n0,H,G1,H_chars,CG,G_chars,A,B,C,D,pi;
  G:=PSL(2,p);
  H:=[];
  H_chars:=[];
  n0:=[ [ Z(p)^0, Z(p) ], [ Z(p)*0, Z(p)^0 ] ];
  CG:=ConjugacyClassesSubgroups(G);
  H[1]:=Representative(CG[2]); # size 2
  H[2]:=Representative(CG[3]); # size 3
  n:=Size(CG);
  for i in [1..n] do
    if Size(Representative(CG[i]))=p then
      H[3]:=Representative(CG[i]);
    fi;
  od;
## H[3]:=Group(n0);       # size p
  for i in [1..Size(H)] do
   H_chars[i]:=Irr(H[i]);
  od;
  G_chars:=Irr(G);
  n:=Length(G_chars);
  A:=[];
  for i in [1..n] do
    pi:=G_chars[i];
    A[i]:=ScalarProduct(H_chars[1][1],RestrictedClassFunction(pi,H[1]));
  od;
  B:=[];
  for i in [1..n] do
    pi:=G_chars[i];
    B[i]:=ScalarProduct(H_chars[2][1],RestrictedClassFunction(pi,H[2]));
  od;
  C:=[];
  for i in [1..n] do
    pi:=G_chars[i];
    C[i]:=ScalarProduct(H_chars[3][1],RestrictedClassFunction(pi,H[3]));
  od;
  D:=[];
  for i in [1..n] do
    pi:=G_chars[i];
    D[i]:=DegreeOfCharacter(pi);
  od;
 return (1/2)*(3*D-A-B-C);
end;

#ram_module_X_JK(5);
#[ 0, 3, 3, 4, 5 ]
#ram_module_X_JK(7);
#[ 0, 7/2, 7/2, 6, 7, 8 ]
#ram_module_X_JK(11);
#[ 0, 11/2, 11/2, 11, 10, 12, 12, 12 ]
#ram_module_X_JK(13);
#[ 0, 7, 7, 13, 13, 13, 13, 15, 14 ]
#ram_module_X_JK(17);
#[ 0, 9, 9, 18, 17, 17, 17, 18, 18, 19, 19 ]
#ram_module_X_JK(19);
#[ 0, 19/2, 19/2, 20, 20, 19, 19, 20, 20, 21, 21, 21 ]
#ram_module_X_JK(23);
#[ 0, 25/2, 25/2, 24, 24, 24, 23, 23, 25, 25, 25, 25, 25, 25 ]
#ram_module_X_JK(29);
#[ 0, 16, 16, 30, 31, 31, 30, 30, 30, 30, 31, 32, 32, 32, 31, 31, 31 ]

pieces_ram_module_X:=function(p)
## p is a prime
## output is (m1,...,mn)
## where n = # conj classes of G=PSL(2,p)
## and mi = mult of pi_i in ram mod of G
##
local G,i,j,n,n0,H,G1,H_chars,CG,G_chars,w,m,theta,pi_theta_to_the;
  G:=PSL(2,p);
  H:=[];
  H_chars:=[];
  n0:=[ [ Z(p)^0, Z(p) ], [ Z(p)*0, Z(p)^0 ] ];
  CG:=ConjugacyClassesSubgroups(G);
  H[1]:=Representative(CG[2]); # size 2
  H[2]:=Representative(CG[3]); # size 3
  n :=Size(CG);
  for i in [1..n] do
    if Size(Representative(CG[i]))=p then
      H[3]:=Representative(CG[i]);
    fi;
  od;
## H[3]:=Group(n0);       # size p
  for i in [1..Size(H)] do
   H_chars[i]:=Irr(H[i]);
  od;
  G_chars:=Irr(G);
  m:=[];
  m[1]:=[];m[2]:=[];m[3]:=[];
  theta:=List([1..3],i->H_chars[i][2]);
  pi_theta_to_the:=List([1..3],i->Sum([1..(Size(H[i])-1)],
      j->j*InducedClassFunction(theta[i]^(j),G)));
#Print("\n\n pi_theta_to_the = ",Sum(pi_theta_to_the),"\n\n");
  for i in [1..3] do
    m[i]:=List(G_chars, pi->ScalarProduct(pi_theta_to_the[i],pi));
  od;
#Print("\n\n m = ",m,"\n\n");
return m;
end;

equiv_deg_module_X:=function(p,ii,r)
## p is a prime
## ii  = 1 for H[1], size 2
## ii  = 2 for H[2], size 3
## ii  = 3 for H[3], size p
## output is (m1,...,mn)
## where n = # conj classes of G=PSL(2,p)
## and mi = mult of pi_i in deg_equiv module of G
##
local G,i,j,n,n0,H,G1,H_chars,CG,G_chars,w,m,theta,pi_theta,pi_theta_to_the;
G:=PSL(2,p);
H:=[]; # 3 decomp gps
H_chars:=[];
n0:=[ [ Z(p)^0, Z(p) ], [ Z(p)*0, Z(p)^0 ] ];
CG:=ConjugacyClassesSubgroups(G);
H[1]:=Representative(CG[2]); # size 2
H[2]:=Representative(CG[3]); # size 3
n:=Size(CG);
for i in [1..n] do
  if Size(Representative(CG[i]))=p then
    H[3]:=Representative(CG[i]);
  fi;
od;
## H[3]:=Group(n0);       # size p
for i in [1..3] do
 H_chars[i]:=Irr(H[i]);
od;
G_chars:=Irr(G);
m:=[];
m[1]:=[];m[2]:=[];m[3]:=[];
theta:=List([1..3],i->H_chars[i][2]);
pi_theta_to_the:=List([1..3],i->List([1..r],j->InducedClassFunction(theta[i]^(-j),G)));
for i in [1..3] do
 for j in [1..r] do
  m[i][j]:=List(G_chars, pi->ScalarProduct(pi_theta_to_the[i][j],pi));
 od;
od;
return Sum(m[ii]);
end;

# equiv_deg_module_X(7,2,1);
#[ 0, 1, 1, 2, 2, 3 ]
# equiv_deg_module_X(7,1,1);
#[ 0, 2, 2, 2, 4, 4 ]
# equiv_deg_module_X(7,3,1);
#[ 0, 1, 0, 1, 1, 1 ]
# equiv_deg_module_X(7,3,7);
#[ 1, 3, 3, 6, 7, 8 ]
# equiv_deg_module_X(7,1,2);
#[ 1, 3, 3, 6, 7, 8 ]

ext: gap

max_stars_repo_name: bopopescu/sage

max_issues_repo_name: Ivo-Maffei/sage

max_forks_repo_name: dimpase/sage


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

# Recursive data structure for SPL programs
# =========================================
# needs AREP, approx.g, aux.g, scalar.g, optrec.g

#F SPLs
#F ====
#F
#F *** WARNING ***
#F *** INFORMATION BELOW IS OBSOLETE ***
#F
#F This file contains a recursive data structure to
#F represent SPL programs in GAP. The data structure is
#F also named SPL. An spl is a GAP record with certain
#F fields as explained below.
#F There are two main applications of spl:
#F 1. They are used by the rule-based formula generator to
#F    generate SPL programs within GAP.
#F 2. As an interface between amats (cf. amat.g in AREP) and
#F    SPL programs.
#F
#F We define a SPL object (SPL) recursively in BNF as
#F the disjoint union of the following cases
#F
#F <SPL> ::=
#F ; atomic cases
#F   | <mat>                                 ; "mat"
#F   | <sparse>                              ; "sparse"
#F   | <mon>                                 ; "mon"  (invertible)
#F   | <diag>                                ; "diag"
#F   | <perm>                                ; "perm" (invertible)
#F
#F ; parametrized symbols
#F   | Symbol(<name>, <params>)              ; "symbol"
#F
#F ; non-terminal symbols
#F   | NonTerminal(<name>, <params>)         ; "nonTerminal"
#F
#F ; composed cases
#F   | <SPL> * .. * <SPL>                    ; "compose"
#F   | DirectSum(<SPL>, .., <SPL>)           ; "directSum"
#F   | TensorProduct(<SPL>, .., <SPL>)       ; "tensor"
#F   | <scalar> * <SPL>                      ; "scalarMultiple"
#F   | <SPL> ^ <permSPL>                     ; "conjugate"
#F   | <SPL> ^ <positive-int>                ; "power"
#F
#F The design of SPL is similar to that of AMat with the
#F exception of the types "symbol" and "nonTerminal".
#F The type "symbol" has been introduced to represent parametrized
#F matrices in an efficient way. E.g., the Fourier transform F_n of
#F size n is presented as
#F   SPLSymbol( "F", [ 3 ] )  and in SPL as  (F 3)
#F without any matrix being stored. The symbol should also be part of SPL
#F (the language). The symbols are contained in the Symboltable (below)
#F together with symbol-specific functions. The symboltable is meant to be
#F extended.
#F
#F Note: you can use scalars to build spls (see scalar.g) but currently
#F only for "scalarMultiple" and "diag".
#F
#F An SPL is a GAP-Rec T with the following mandatory
#F fields common to all types of SPL
#F
#F   isSPL        := true, identifies spls
#F   operations   := SPLOps, GAP operations record for T
#F   type         : a string identifying the type of T
#F   dimensions   : size of the matrix represented (= [rows, columns])
#F
#F The characteristic of the base field is always = 0.
#F
#F The following fields are optional for any type of spl S:
#F
#F   splOptions: a list of strings supplying options for the spl compiler
#F     "unrolled": specifies that S is to be unrolled
#F     when exported to spl. In this case this option is disregarded for
#F     all subtrees (if an spl is to be unrolled, then all its subtrees have
#F     to be unrolled).
#F     "pcl": a pcl label should be set in the spl program
#F   pclLabel: a label for pcl profiling
#F   root: contains the spl that has been expanded to obtain S.
#F
#F The following fields are mandatory to the type of T (T.type):
#F
#F type = "mat"
#F   T.element: matrix represented
#F type = "sparse"
#F   T.element: list of triples [i, j, a_(i,j)]
#F type = "mon"
#F   T.element: mon represented
#F type = "diag"
#F   T.element: list containing the diagonal
#F type = "perm"
#F   T.element: perm represented
#F type = "symbol"
#F   T.symbol: string representing the symbol
#F   T.params: list of parameters for the symbol
#F type = "nonTerminal"
#F   T.symbol: string representing the symbol
#F   T.isRealTransform: true if the transform is real, false else;
#F     a transform is real iff the matrix is real for all choices of parameters
#F   T.params: list of parameters for the symbol
#F   T.index: the index of the non-terminal in NonTerminalTableSPL
#F   T.ruleinds: the indices of applicable rules from RuleTable
#F type = "compose"
#F   T.children(): list of spls, the factors
#F type = "directSum"
#F   T.children(): list of spls, the summands
#F type = "tensor"
#F   T.children(): list of spls, the tensor factors
#F type = "scalarMultiple"
#F   T.scalar: the scalar
#F   T.element: the SPL multiplied
#F type = "conjugate"
#F   T.element: the SPL to be conjugated
#F   T.conjugation: the conjugating SPL
#F type = "power"
#F   T.element: the SPL raised to the exponent
#F   T.exponent: the exponent, an integer > 0
#F
#F *** WARNING ***
#F *** INFORMATION ABOVE IS OBSOLETE ***

# Internal Functions for SPLs
# ----------------------------
#

# the operations record
SPLOps := CantCopy(OperationsRecord("SPLOps"));

# SPLOps.IsDegree( <perm>, <degree> )
#   tests whether <degree> is a degree consistent with <perm>

SPLOps.IsDegree := function ( p, n )
  if not IsPerm(p) then
    Error("usage: SPLOps.IsDegree( <perm>, <degree> )");
  fi;
  return
    IsInt(n) and
    n >= 0 and
    ( p = ( ) or n >= LargestMovedPointPerm(p) );
end;

# SPLOps.IsDimension( <dims/deg> )
#   tests whether <dims/deg> is a valid dimension for a
#   matrix, i.e. a positive integer or a pair of them and
#   returns a pair representing the dimension.

SPLOps.IsDimension := function ( D )
  if IsList(D) then
    if not (
      Length(D) = 2 and
      ForAll(D, n -> IsInt(n) and n > 0)
    ) then
      Error("<D> must be a list of two positive ints");
    fi;
    return D;
  else
    if not ( IsInt(D) and D > 0 ) then
      Error("<D> must be a positive int");
    fi;
    return [D, D];
  fi;
end;

#F Fundamental Constructors and Tests for SPLs
#F --------------------------------------------
#F

#F IsSPL( <obj> )
#F   tests whether <obj> is an SPL. The particular case of SPL
#F   is provided by SPL.type and the additional fields are available
#F   directly from the record.
#F
IsSPL := S -> IsRec(S) and IsBound(S.isSPL) and S.isSPL;

#F IsPermutationSPL( <spl> )
#F   returns true, if <spl> is of type "perm" or a symbol that
#F   represents a permutation.
#F
IsPermutationSPL := function ( S )
  if not IsSPL(S) then
      Error("<S> must be an spl");
  elif IsBound(S.isPermutation) then
      return S.isPermutation();
  else
      return false;
  fi;
end;

#F Comparison of SPLs
#F ------------------
#F

#F SPLOps.\=( <spl1>, <spl2> )
#F   returns true if <spl1> and <spl2> are structurally equal.
#F
SPLOps.\= := (S1, S2) -> Cond(IsRec(S1) and IsBound(S1.equals), S1.equals(S2), false);
# NOTE: this function will not yield a good order is Set(..) is used
SPLOps.\< := (S1, S2) -> Cond(IsRec(S1) and IsBound(S1.lessThan), S1.lessThan(S2), BagAddr(S1) < BagAddr(S2));

#F IsIdenticalSPL ( <spl1>, <spl2> )
#F    returns true if <spl1> and <spl2> are spls with identical
#F    syntax trees.
#F
IsIdenticalSPL := (S1, S2) ->
   Checked(IsSPL(S1), IsSPL(S2),
       Cond(S1.dimensions <> S2.dimensions, false,
	    ObjId(S1) <> ObjId(S2), false,
	    S1.equals(S2)));

IsEqualObj := (o1,o2) -> Cond(
    Same(o1,o2), true,
    not IsRec(o1) or not IsBound(o1.equals), o1=o2,
    o1.equals(o2)
);

#F IsHashIdenticalSPL ( <spl1>, <spl2> )
#F    returns true if <spl1> and <spl2> are identical for the
#F    hashing purposes (see spiral/search).
#F

IsHashIdenticalSPL := function( S1, S2 )
    if not ( IsSPL(S1) and IsSPL(S2) ) then
	Error("usage: IsHashIdenticalSPL ( <spl1>, <spl2> )");
    fi;

    if IsBound(S1.hashId) and IsBound(S2.hashId) then
       return Same(ObjId(S1), ObjId(S2)) and
              S1.dimensions = S2.dimensions and
              S1.transposed = S2.transposed and
              S1.hashId = S2.hashId;
    elif IsBound(S1.HashId) and IsBound(S2.HashId) then #NOTE: IS THAT RIGHT???
       return Same(ObjId(S1), ObjId(S2)) and
              When(ForAll(Concat(S1.dimensions, S2.dimensions), IsInt), S1.dimensions = S2.dimensions, true) and
              S1.transposed = S2.transposed and
              S1.HashId() = S2.HashId();
    else
       return IsIdenticalSPL( S1, S2 );
    fi;
end;

#F Converting SPLs
#F ---------------

#F MatSparseSPL ( <"sparse"-spl> )
#F   returns the matrix represented by the spl of type "sparse".
#F
MatSparseSPL := function ( S )
  local M, t;
  #Constraint(ObjId(S) = Sparse);
  M := NullMat(S.dimensions[1], S.dimensions[2]);
  for t in S.element do
    M[t[1]][t[2]] := t[3];
  od;

  return M;
end;

#F AMatSPL( <SPL> )
#F   converts <SPL> into an amat.
#F
AMatSPL := S -> S.toAMat();

#F MatSPL( <spl> )
#F   returns the matrix represented by <spl>.
#F
MatSPL := S -> MatAMat(AMatSPL(S));

#F PermSPL( <spl> )
#F   returns the permutation represented by <spl>,
#F   or false if <spl> is no permutation.
#F
PermSPL := S -> PermAMat(AMatSPL(S));

#F Other Functions
#F ---------------
#F

#F TerminateSPL( <spl> )
#F   returns an spl equivalent to <spl> by replacing
#F   all non-terminals with terminal spls.
#F
TerminateSPL := S -> S.terminate();

#F TransposedSPL ( <SPL> )
#F   returns an spl representing the transpose of <SPL>.
#F
TransposedSPL := S -> S.transpose();

# ArithmeticCostScalarComplex ( <scalar> )
#   returns the arithmetic cost of multiplying <scalar>
#   to a complex number in the form
#     [adds, mults, sign-changes]
#   <scalar> can be a scalar (see scalar.g) or a number.
#
ArithmeticCostScalarComplex := function ( s )
  local n;
  # convert scalar into number (scalar.g)
  s := EvalScalar(s);

  if s = 1 then     return [0, 0, 0];
  elif s = -1 then  return [0, 0, 2];
  # s is real
  elif GaloisCyc(s, -1) = s then return [0, 2, 0];
  else
    n := NofCyc(s);

    # determine exponent of s as root of unity
    if n mod 2 = 1 and s ^ n = -1 then
      n := 2 * n;
    fi;

    # s is no root of unity
    if s^n <> 1 then
      return [2, 4, 0];
    fi;

    # special roots of unity
    if n = 4 then  return [0, 0, 1];
    elif n = 8 then return [2, 2, 0];
    else return [2, 4, 0];
    fi;
  fi;
end;

IsRealNumber := function (num)
   if IsValue(num) then num := num.v; fi;

   if IsSymbolic(num) then
       return When(IsVecT(num.t) or IsArrayT(num.t), not IsComplexT(num.t.t), not IsComplexT(num.t));
   elif IsComplex(num) then return ImComplex(num)=0;
   elif IsDouble(num) then return true;
   elif IsCyc(num) then return GaloisCyc(num, -1) = num;
   elif IsList(num) then return ForAll(num, IsRealNumber);
   else return true;
   fi;
end;

#F IsRealSPL( <spl> )
#F   returns true if <spl> contains no complex subexpressions
#F   and false else. Note that this function does not determine whether
#F   <terminal spl> represents a real matrix! E.g.,
#F     E(4) * Mat([[-E(4)]])
#F   is a real matrix, but contains complex expressions. The result is
#F   stored in .isReal
#F   If the field nonTerminal is set in <spl> then true is returned iff
#F   the nonTerminal is a real transform. This assumes that a real transform
#F   is never expanded using complex expressions!
#F   The function is used to determine whether the compiler option -R
#F   should be used on exporting to spl.
#F
IsRealSPL     := S -> S.isReal();

IsIdentity    := S -> IsBound(S.isIdentity) and S.isIdentity();
IsIdentitySPL := S -> IsSPL(S) and IsIdentity(S);
IsIdentityObj := S -> IsBound(S.isIdentity) and S.isIdentity = True;


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

# The Wire class implements OL connection system
# usage: Wire(inputs, interconnection matrix)
# ex: Wire([1,2,3,4],MatSPL(L(4,2)));
 Class(Wire, BaseMat, rec(
   abbrevs := [ (inp, mat) ->  Checked(IsList(inp), IsMat(mat), [[inp, mat]])],
   new := meth(self, p)
       local res;
       res := SPL(WithBases(self, rec(params := p, TType:=Replicate(Length(p[1]),TUnknown))));
       return res;
   end,
   isPermutation := self >> false,
   #now just make it work, assuming all inputs have the same type
   dmn:=self >> List(self.params[1],x->TArray(self.TType[1],x)),
   rng:=self >> List(self.params[2]*self.params[1],x->TArray(self.TType[1],x)),
   sums:= self>>self,
   numops:=self>>0,
   rChildren := self >> [],
   print := meth(self,i,is)
     Print(self.name, "(", PrintCS(self.params), ")");
     return;
   end
 ));

 Class(LinkIO, BaseMat, rec(
   abbrevs := [ (sums, mat) ->  Checked(IsMat(mat), [[sums, mat]])],
   new := meth(self, p)
       local res;
       res := SPL(WithBases(self, rec(params := p, TType:=Replicate(Length(p[1]),TUnknown))));
       return res;
   end,
   isPermutation := self >> false,
   dmn:=self >> self.params[1].dmn(),
   rng:=self >> self.params[1].rng(),
   sums:= self>>self,
   numops:=self>>0,
   rChildren := self >> self.params,
   rSetChild := meth ( self, n, newChild )
     if n = 1  then
         self.params[1] := newChild;
     elif n = 2 then
         self.params[2] := newChild;
    else
        Error("<n> not in range");
     fi;
   end,
   print := meth(self,i,is)
     Print(self.name, "(", PrintCS(self.params), ")");
     return;
   end
 ));


text:
#############################################################################
####
##
#W  general.gi                 ACE Package                   Alexander Hulpke
#W                                                                Greg Gamble
##
##  This file installs mainly non-interactive ACE  variables  and  functions.
##  Though Alexander will barely recognise it,  some  of his ideas are  still
##  present.
##
#Y  Copyright (C) 2000  Centre for Discrete Mathematics and Computing
#Y                      Department of Information Technology & Electrical Eng.
#Y                      University of Queensland, Australia.
##

#############################################################################
####
##
#V  ACETCENUM . . . . . . . .  The ACE version of the coset enumerator TCENUM
##  . . . .  CosetTableFromGensAndRels is set to ACECosetTableFromGensAndRels
##
InstallValue(ACETCENUM, rec(
  name := "ACE (Advanced Coset Enumerator)",
  CosetTableFromGensAndRels := ACECosetTableFromGensAndRels
));

#############################################################################
####
##
#F  InfoACELevel . . . . . . . . . . . . . . .  Get the InfoLevel for InfoACE
##
##
InstallGlobalFunction(InfoACELevel, function()
  return InfoLevel(InfoACE);
end);

#############################################################################
####
##
#F  SetInfoACELevel . . . . . . . . . . . . . . Set the InfoLevel for InfoACE
##
##
InstallGlobalFunction(SetInfoACELevel, function(arg)
  if IsEmpty(arg) then
    SetInfoLevel(InfoACE, 1);     # Set to default level
  else
    SetInfoLevel(InfoACE, arg[1]);
  fi;
end);

#############################################################################
####
##
#F  ACEPackageVersion()
##
##  returns the version number of the current ACE package.
##
InstallGlobalFunction(ACEPackageVersion, function()

  return GAPInfo.PackagesInfo.ace[1].Version;
end);

#############################################################################
####
##
#F  CALL_ACE . . . . . . . . . Called by ACECosetTable, ACEStats and ACEStart
##
##
InstallGlobalFunction(CALL_ACE, function(ACEfname, fgens, rels, sgens)
local optnames, echo, errmsg, onbreakmsg, infile, datarec, ToACE, gens,
      standard, ignored;

  if ValueOption("aceexampleoptions") = true and
     IsBound(ACEData.aceexampleoptions) then
    SANITISE_ACE_OPTIONS(OptionsStack[ Length(OptionsStack) ],
                         ACEData.aceexampleoptions);
    PushOptions(ACEData.aceexampleoptions);
    Unbind(ACEData.aceexampleoptions);
    ACEData.options := OptionsStack[ Length(OptionsStack) ];
    PopOptions();
    OptionsStack[ Length(OptionsStack) ] := ACEData.options;
    Unbind(ACEData.options);
  fi;
  optnames := ACE_OPT_NAMES();
  # We have hijacked ACE's echo option ... we don't actually pass it to ACE
  echo := ACE_VALUE_ECHO(optnames);

  ECHO_ACE_ARGS( echo, ACEfname, rec(fgens := fgens,
                                     rels  := rels,
                                     sgens := sgens) );
  # Check arguments are valid
  while IsEmpty(fgens) do
    errmsg :=
        ["fgens (arg[1]) must be a non-empty list of group generators ..."];
    onbreakmsg :=
        ["Type: 'quit;' to quit to outer loop, or",
         "type: 'fgens := <val>; return;' to assign <val> to fgens to continue."
        ];
    Error(ACE_ERROR(errmsg, onbreakmsg), "\n");
  od;
  fgens := ACE_FGENS_ARG_CHK(fgens);
  rels  := ACE_WORDS_ARG_CHK(fgens, rels, "relators");
  sgens := ACE_WORDS_ARG_CHK(fgens, sgens, "subgp gen'rs");

  infile  := VALUE_ACE_OPTION(optnames, fail, "aceinfile");
  if ACEfname = "ACECosetTableFromGensAndRels" and infile <> fail then
    datarec := rec(
        infile  := infile,
        outfile := VALUE_ACE_OPTION(optnames, ACEData.outfile, "aceoutfile"),
        stream  := OutputTextFile(infile, false) );
    ToACE := function(list) WRITE_LIST_TO_ACE_STREAM(datarec.stream, list); end;
  else
    datarec := rec(
        stream := InputOutputLocalProcess(ACEData.tmpdir, ACEData.binary, []) );
    if datarec.stream = fail then
      Error("sorry! Run out of pseudo-ttys. Can't initiate stream\n");
    fi;
    if ACEfname <> ACEStart then
      datarec.procId := 0;
      ACEData.ni := datarec;
    fi;
    FLUSH_ACE_STREAM_UNTIL(datarec.stream, 3, 3, ACE_READ_NEXT_LINE,
                           line -> IsMatchingSublist(line, "name", 3));
    ToACE := function(list)
                 INTERACT_TO_ACE_WITH_ERRCHK(datarec, list);
             end;
  fi;
  datarec.args    := rec(fgens := fgens, rels := rels, sgens := sgens);
  datarec.options := ACE_OPTIONS();
  standard := ACE_COSET_TABLE_STANDARD( datarec.options );

  # Define the group generators ACE will use
  gens := TO_ACE_GENS(fgens);
  ToACE([ "Group Generators: ", gens.toace, ";"]);
  datarec.acegens := gens.acegens;

  # Define the group relators ACE will use
  datarec.enforceAsis :=
      (ACEfname <> "ACEStats") and (standard = "lenlex") and
      not IsACEGeneratorsInPreferredOrder(fgens, rels, "noargchk");
  ToACE([ "Group Relators: ",
          ACE_RELS(rels, fgens, datarec.acegens, datarec.enforceAsis), ";" ]);

  # Define the subgroup generators ACE will use
  ToACE([ "Subgroup Generators: ",
          ACE_WORDS(sgens, fgens, datarec.acegens), ";" ]);

  if ACEfname  = "ACECosetTableFromGensAndRels" then
    ignored := [ ];
  else
    ignored := [ "aceinfile" ];
  fi;
  if ACEfname  = "ACEStart" then
    Add(ignored, "aceoutfile");
  fi;
  if datarec.enforceAsis then
    Add(ignored, "asis");
    ToACE([ "Asis: 1;" ]);
  fi;

  PROCESS_ACE_OPTIONS(
      ACEfname, optnames, optnames, echo, datarec,
      rec(group      := ACE_ERRORS.argnotopt, # disallowed options
          generators := ACE_ERRORS.argnotopt,
          relators   := ACE_ERRORS.argnotopt),
      ignored
      );

  if not IsInputOutputStream(datarec.stream) then
    if VALUE_ACE_OPTION(optnames, fail, ["start", "aep", "rep"]) = fail then
      # if the user hasn't issued there own enumeration initiation directive
      # ... initiate the enumeration
      ToACE([ "Start;" ]);
    fi;
    if ACEfname = "ACECosetTableFromGensAndRels" then
      if standard = "lenlex" then
        ToACE([ "Standard;" ]);
      fi;
      ToACE([ "Print Table;" ]);
    fi;
    CloseStream(datarec.stream);
  elif ACEfname <> "ACEStart" then
    if VALUE_ACE_OPTION(optnames, fail, ["start", "aep", "rep"]) = fail then
      ACE_MODE( "Start", datarec );
    fi;
    if ACEfname = "ACECosetTableFromGensAndRels" and standard = "lenlex" then
      ToACE([ "Standard;" ]);
    fi;
  fi;

  if ACEfname = "ACEStart" then
    datarec.procId := Length(ACEData.io) + 1;
    Add(ACEData.io, datarec);
    return Length(ACEData.io);
  elif ACEfname = "ACECosetTableFromGensAndRels" then
    datarec.silent := VALUE_ACE_OPTION(optnames, false, "silent");
  fi;
  return datarec;
end);

#############################################################################
####
##
#F  ACECosetTableFromGensAndRels . . . . . . .  Non-interactive ACECosetTable
##
##
InstallGlobalFunction(ACECosetTableFromGensAndRels, function(fgens, rels, sgens)
  # Use ACECosetTable non-interactively
  return ACECosetTable(fgens, rels, sgens);
end);

#############################################################################
####
##
#F  IsACEStandardCosetTable . . . . . . Returns true if table is standardised
##  . . . . . . . . . . . . . . . . . . according to GAP's default scheme  or
##  . . . . . . . . . . . . . . . . . . with the lenlex option, according  to
##  . . . . . . . . . . . . . . . . . . . . the lenlex standardisation scheme
##
InstallGlobalFunction(IsACEStandardCosetTable, function(table)
local standard, geninvIndices, index, next, j, i;

  standard := ACE_COSET_TABLE_STANDARD( ACE_OPTIONS() );
  if standard in ["lenlex", "GAPlenlex"] then
    geninvIndices := [1 .. Length(table)];
  elif standard in ["semilenlex", "GAPsemilenlex"] then
    geninvIndices := [1, 3 .. Length(table) - 1];
  else
    return IsStandardized(table); # Should only get here with GAP 4.3+
                                  # ... by which time `IsStandardized'
                                  # will hopefully have been generalised
                                  # to cope with any other standardisation
                                  # schemes
  fi;

  index := Length( table[1] );
  next := 2;
  for j in [1 .. index - 1] do
    for i in geninvIndices do
      if table[i][j] >= next then
        if table[i][j] > next then
          return false;
        fi;
        next := next + 1;
      fi;
    od;
  od;
  return true;
end);

#############################################################################
####
##
#F  IsACEGeneratorsInPreferredOrder . . . . . Returns true if the  generators
##  . . . . . . . . . . . . . . . . . . . . . gens are already in  the  order
##  . . . . . . . . . . . . . . . . . . . . . . . . . . . .  preferred by ACE
##
##  For a presentation with more than one generator, the first  generator  of
##  which is an involution and the second is not, ACE prefers to  switch  the
##  first two generators. IsACEGeneratorsInPreferredOrder returns true if the
##  order of the generators gens would not  be  changed  by  ACE  and  false,
##  otherwise. When necessary, the argument rels (the  relators)  is  scanned
##  for relators that determine  whether  or  not  gens[1]  and  gens[2]  are
##  involutions.
##
##  If IsACEGeneratorsInPreferredOrder would return false, it is possible  to
##  enforce a user's order of the generators within ACE, by  enforcing  ACE's
##  `asis' option and passing the relator,  that  determines  gens[1]  is  an
##  involution,  explicitly  to  ACE   as:   gens[1]*gens[1]   (rather   than
##  gens[1]^2).
##
InstallGlobalFunction(IsACEGeneratorsInPreferredOrder, function(arg)
local ioIndex, gens, rels;

  if Length(arg) < 2 then
    ioIndex := CallFuncList(ACEProcessIndex, arg);
    gens := ACEGroupGenerators(ioIndex);
    rels := ACERelators(ioIndex);
  elif Length(arg) = 2 then
    gens := ACE_FGENS_ARG_CHK(arg[1]);
    rels := ACE_WORDS_ARG_CHK(gens, arg[2], "relators");
  elif Length(arg) = 3 and arg[3] = "noargchk" then
    # This scenario only intended for use internally
    gens := arg[1];
    rels := arg[2];
  else
    Error("expected at most 2 arguments, not ", Length(arg), " arguments.\n");
  fi;

  if Length(gens) = 1 or not ForAny(rels, rel -> rel = gens[1]^2) then
    return true;
  else
    return ForAny(rels, rel -> rel = gens[2]^2);
  fi;
end);

#############################################################################
####
##
#F  ACE_READ_AS_FUNC  . . . . . . . . . . . . . . . Variant of ReadAsFunction
##  . . . . . . . . . . . .  that allows the passing of the function argument
##  . . . . . . . . . . . . . . . . . . . . ACEfunc to ReadAsFunction(file)()
##
##
InstallGlobalFunction(ACE_READ_AS_FUNC, function(file, ACEfunc)
local line, instream, rest;

  instream := InputTextFile(file);
  # We don't want the user to see this ... so we flush at InfoACE level 10.
  line := FLUSH_ACE_STREAM_UNTIL( instream, 10, 10, ReadLine,
                                  line -> IsMatchingSublist(line, "local") );
  rest := ReadAll(instream);
  CloseStream(instream);
  return ReadAsFunction(
             InputTextString(
                 Concatenation([ ReplacedString(line, ";", ", ACEfunc;"),
                                 "ACEfunc := ", NameFunction(ACEfunc), ";",
                                 rest ]) ) )();
end);

#############################################################################
####
##
#F  ACEExample( )
#F  ACEExample( <file>[, <ACEfunc>] )
##
##  With no arguments, or with single argument "index", or a string  that  is
##  not a filename  in  the  `examples'  directory,  an  index  of  available
##  examples is displayed.
##
##  With argument <file> that is a filename in the `examples' directory other
##  than "index" the example is displayed as it would  be  when  called  with
##  <ACEfunc> (or `ACEStats', if  the  2nd  argument  is  omitted)  and  then
##  executed via a call to `ReadAsFunction' and a little internal  ``magic''.
##  <ACEfunc>    must    be    one    of    `ACEStats'     (the     default),
##  `ACECosetTableFromGensAndRels'  (or  equivalently   `ACECosetTable',   or
##  `ACEStart'.
##
InstallGlobalFunction(ACEExample, function(arg)
local name, file, instream, line, ACEfunc,
      EnquoteIfString, optnames, lastoptname, optname;

  if IsEmpty(arg) then
    name := "index";
  else
    name := arg[1];
    if Length(arg) > 1 then
      ACEfunc := arg[2];
    else
      ACEfunc := ACEStats;
    fi;
    if not IsEmpty(OptionsStack) then
      ACEData.aceexampleoptions := OptionsStack[ Length(OptionsStack) ];
      PopOptions();
      PushOptions( rec(aceexampleoptions := true) );
    fi;
  fi;
  file := Filename( DirectoriesPackageLibrary( "ace", "examples"), name );
  if file = fail then
    Info(InfoACE + InfoWarning, 1,
         "Sorry! There is no ACE example file with name `", name, "'");
    name := "index";
    file := Filename( DirectoriesPackageLibrary("ace", "examples"), name );
  fi;
  # Display file ... after a few minor modifications
  instream := InputTextFile(file);
  if name <> "index" then
    line := FLUSH_ACE_STREAM_UNTIL( instream, 1, 10, ReadLine,
                                    line -> IsMatchingSublist(line, "local") );
    Info(InfoACE, 1,
         "#", line{[Position(line, ' ')..Position(line, ';') - 1]},
         " are local to ACEExample");
    line := FLUSH_ACE_STREAM_UNTIL( instream, 1, 10, ReadLine,
                                    line -> IsMatchingSublist(line, "return") );
    Info(InfoACE, 1,
         Chomp(ReplacedString(line, "return ACEfunc", NameFunction(ACEfunc)))
         );
    if IsBound(ACEData.aceexampleoptions) then
      line := FLUSH_ACE_STREAM_UNTIL(
                  instream, 1, 10, ReadLine,
                  line -> PositionSublist(line, ");") <> fail );
      Info(InfoACE, 1, Chomp(ReplacedString(line, ");", ", ")));
      Info(InfoACE, 1, "    # User Options");
      optnames := ShallowCopy( RecNames(ACEData.aceexampleoptions) );
      lastoptname := optnames[ Length(optnames) ];
      Unbind(optnames[ Length(optnames) ]);

      EnquoteIfString := function(optval)
      # Puts quotes around optval if it's a string
        if IsString(optval) then
          return Concatenation(["\"", optval, "\""]);
        else
          return optval;
        fi;
      end;

      for optname in optnames do
        Info(InfoACE, 1, "      ", optname, " := ",
                         EnquoteIfString(
                             ACEData.aceexampleoptions.(optname) ), ",");
      od;
      Info(InfoACE, 1, "      ", lastoptname, " := ",
                       EnquoteIfString(
                           ACEData.aceexampleoptions.(lastoptname) ), ");");
    fi;
  fi;
  FLUSH_ACE_STREAM_UNTIL( instream, 1, 10, ReadLine, line -> line = fail );
  CloseStream(instream);
  if name <> "index" then
    return ACE_READ_AS_FUNC(file, ACEfunc);
  fi;
end);

#############################################################################
####
##
#F  ACEReadResearchExample  . . . . . . . .  Read  an  ACE  research  example
##  . . . . . . . . . . . . . . . . . . . .  from the res-examples directory.
##
##  Currently,  all  examples  in  the  res-examples  directory   depend   on
##  pgrelfind.g, which with Info text doubles as an index. This  function  is
##  essentially equivalent to doing a Read of its argument or  "pgrelfind.g",
##  if there is no argument.
##
InstallGlobalFunction(ACEReadResearchExample, function(arg)
local name, file;

  if IsEmpty(arg) then
    name := "pgrelfind.g"; # If there is ever more than one key research
                           # example, we should replace this with an index
  else
    name := arg[1];
  fi;
  file := Filename( DirectoriesPackageLibrary( "ace", "res-examples"), name );
  if file = fail then
    Error("ACEReadResearchExample: Sorry! There is no ACE research example\n",
          "file with name \"", name, "\"\n");
  else
    Read(file);
  fi;
end);

#############################################################################
####
##
#F  ACEPrintResearchExample . . . . . . . . Print  an  ACE  research  example
##  . . . . . . . . . . . . . . . . . . . . from the  res-examples  directory
##  . . . . . . . . . . . . . . . . . . . . to the terminal  or  to  a  file,
##  . . . . . . . . . . . . . . . . . . . . . . . minus header and Info text.
##
##  ACEPrintResearchExample(examplefile)
##      prints examplefile in res-examples directory to the terminal
##
##  ACEPrintResearchExample(examplefile, outfile)
##      prints examplefile in res-examples directory to outfile
##
InstallGlobalFunction(ACEPrintResearchExample, function(arg)
local outstream, print, file, instream, line;

  if IsEmpty(arg) then
    Error("expected 1 or 2 arguments\n");
  fi;

  file := Filename( DirectoriesPackageLibrary( "ace", "res-examples"), arg[1] );
  if file = fail then
    Error("ACEPrintResearchExample: Sorry! There is no ACE research example ",
          "file with name `", arg[1], "'\n");
  fi;

  if Length(arg) > 1 then
    outstream := OutputTextFile(arg[2], false);
    print := function(line) WriteAll(outstream, line); end;
  else
    print := Print;
  fi;

  instream := InputTextFile(file);
  repeat
    line := ReadLine(instream);
  until IsMatchingSublist(line, "## Begin");
  line := ReadLine(instream);
  while not IsMatchingSublist(line, "## End") do
    print(line);
    line := ReadLine(instream);
  od;
  CloseStream(instream);

  if print <> Print then
    CloseStream(outstream);
  fi;
end);

#############################################################################
####
##
#F  ACEDirectoryTemporary( <dir> )
##
##  calls the UNIX command `mkdir' to create <dir>, which must be  a  string,
##  and if successful a directory  object  for  <dir>  is  both  assigned  to
##  `ACEData.tmpdir'  and   returned.   The   fields   `ACEData.infile'   and
##  `ACEData.outfile' are also set to be files in  `ACEData.tmpdir',  and  on
##  exit from {\GAP} <dir> is removed.
##
InstallGlobalFunction(ACEDirectoryTemporary, function(dir)
local created;

  # check arguments
  if not IsString(dir) then
    Error("usage: ACEDirectoryTemporary( <dir> ) : <dir> must be a string.\n");
  fi;

  # create temporary directory
  created := Process(DirectoryCurrent(),
                     Filename( DirectoriesSystemPrograms(), "sh" ),
                     InputTextUser(),
                     OutputTextUser(),
                     [ "-c", Concatenation("mkdir ", dir) ]);
  if created = fail then
    return fail;
  fi;

  Add( GAPInfo.DirectoriesTemporary, dir );
  ACEData.tmpdir := Directory(dir);
  ACEData.infile := Filename(ACEData.tmpdir, "in");
  ACEData.outfile := Filename(ACEData.tmpdir, "out");
  return ACEData.tmpdir;
end);

#############################################################################
####
##
#F  ACE_ERROR(<errmsg>, <onbreakmsg>)
##
##  sets `OnBreakMessage' to  print  <onbreakmsg>  in  order  to  generate  a
##  one-off user-friendly message of how the user may recover from the error,
##  and returns an error message formed from <errmsg> to be used by Error.
##
##  <errmsg> and <onbreakmsg> should be lists of strings;  <onbreakmsg>  must
##  be non-empty and its first member must not be a null string.
##
InstallGlobalFunction(ACE_ERROR, function(errmsg, onbreakmsg)
local NormalOnBreak, NormalOnBreakMessage;

  errmsg := JoinStringsWithSeparator(errmsg, "\n ");
  NormalOnBreakMessage := OnBreakMessage;
  onbreakmsg[1]{[1]} := LowercaseString( onbreakmsg[1]{[1]} );
  OnBreakMessage := function()
    local s;

    for s in onbreakmsg do
      Print(" ", s, "\n");
    od;
    OnBreakMessage := NormalOnBreakMessage;
  end;

  return errmsg;
end);

#############################################################################
####
##
#F  CallACE . . . . . . . . . . . . . . . . . . . . . . . . . . .  deprecated
##
InstallGlobalFunction(CallACE, function(arg)

  Error("CallACE is deprecated: Use `ACECosetTableFromGensAndRels' or\n",
        "`ACECosetTable'.\n");
end);

#E  general.gi  . . . . . . . . . . . . . . . . . . . . . . . . . . ends here

max_stars_repo_name: wilfwilson/ace

max_issues_repo_name: wilfwilson/ace

max_forks_repo_name: wilfwilson/ace


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

#IsDigit := c -> c in "0123456789";
#var.prefix := meth( self )
#    local p;
#    p := FirstPosition(self.id,IsDigit);
#    p := SplitAt(self.id,When(p>0,p-1,0))[1];
#    Constraint(p <> "");
#    return p;
#end;

Class(VarGenBase, rec(
    __call__ := self >> WithBases(self, rec(counter:=rec(), pfx := "")),

    suffix := n -> StringInt(n),

    reset := meth(self) self.counter := rec(); end,

    nextu := (self, t) >> When(self.pfx="", self.next(t, "u"), self.next(t, "")),

    next := meth(self, t, pfx)
        local p;
	Constraint(IsString(pfx)); Constraint(IsType(t));
	pfx := Concat(self.pfx, pfx);
	if not IsBound(self.counter.(pfx)) then
            # first will be plain "pfx", then "pfx2"
	    p := self.var(t, pfx);
	    self.counter.(pfx) := 2;
	else
	    p := self.var(t, Concat(pfx,  self.suffix(self.counter.(pfx))));
	    self.counter.(pfx) := self.counter.(pfx) + 1;
	fi;
	return p;
    end,

    nextName := (self, pfx) >> self.next(TInt, pfx).id,

    var := (t, pfx) -> Error("must be implemented in subclasses"),

    # creates a cloned generator with aliased .counter, and an extra prefix
    # that will be added to all params created with .next()
    withPrefix := (self, pfx) >> When(
        pfx="i",
        Error("i prefix is RESERVED for loop variables!"),
        WithBases(self, rec(pfx := Concat(self.pfx, pfx))))
));

Class(VarGenNumeric, VarGenBase, rec(
    var := (t, name) -> var(name, t)
));

Class(VarGenSymbolic, VarGenBase, rec(
    var := (t, name) -> var(name, t),
    suffix := VarNameInt
));

Class(ParamGenNumeric, VarGenBase, rec(
    var := (t, name) -> param(t, name)
));

Class(ParamGenSymbolic, VarGenBase, rec(
    var := (t, name) -> param(t, name),
    suffix := VarNameInt
));

Class(NewInt,rec(
	__call__ := meth(self)
	    self.n := self.n + 1;
	    return self.n;
	end,

	n := 0,
));

Class(VarMapper, rec(
    __call__ := (self, vargen) >> WithBases(self, rec(
	    sn := NewInt(),
	    bindings := tab(),
	    vargen := vargen)),

    reset := meth(self)
	self.vargen.reset();
	self.bindings := tab();
	self.sn := NewInt();
    end,

    ignore := (self, var) >> false,

    alreadyMapped := (self,var) >> IsBound(var.sn) and var.sn = self.sn,

    map := meth(self, var)
        local newvar;
	if self.ignore(var) or (IsBound(var.NoReMap) and var.NoReMap) or self.alreadyMapped(var) then
	    return var;
	fi;
        if IsBound(self.bindings.(var.id)) then
	    return self.bindings.(var.id);
	else
	    newvar := self.vargen.next(var.t, var.id{[1]});
	    newvar.sn := self.sn;
	    self.bindings.(var.id):=newvar;
	    return newvar;
	fi;
    end
));

ProperName := function(o)
    o.properName := true;
    return o;
end;

_RemapVars := function(c, ignore_list, varGen )
    local vmapper;
    vmapper := CopyFields(VarMapper(varGen),
        rec(ignore := (self, var) >> IsBound(var.properName) and var.properName=true
                                     or var.id in ignore_list));
    return SubstTopDownRules(c, [
	    [var,v->vmapper.map(v)],
	    [decl,d->decl(List(d.vars,v->vmapper.map(v)),d.cmd)],
	    [data,d->data(vmapper.map(d.var),d.value,d.cmd)] ]);
end;

RemapVars := c -> _RemapVars(c, ["X", "Y"], VarGenNumeric);

RemapVarsIgnore := (c, ignore_list) -> _RemapVars(c, Concatenation(["X", "Y"], ignore_list));

RemapVarsSafe := function(c, ignore_list)
    local free, args, init;
    free := c.free();
    args := Set(ConcatList( Collect(c, @(1, [func])), x->x.params));
    init := Collect(c, @(1, [var,param], x->IsBound(x.value) or IsBound(x.init))); # .value and .init: this is a hack!
    return  _RemapVars(c, List(free::args, x->x.id) :: init :: ignore_list, ParamGenSymbolic());
end;

Class(RemoveUnusedVars,rec(__call__:=function(code)
    local usedvars,declvars,unusedvars;

    usedvars := Set(Flat([
		  	Collect(code,var)	# all used vars in operands of assign-cmds
		]));

    declvars := Set(Flat([
##  			 X, Y, # implicitly they are arguments of DFT kernel
			 List(Collect(code,decl),d->d.vars)
		]));

    unusedvars := declvars;
    SubtractSet(unusedvars,usedvars);

    # remove unused vars
	code := SubstTopDown(code,decl,d->decl(Difference(d.vars, unusedvars),d.cmd));

    return code;
end));

##  # Maybe, it would be usefull to add someting else to remove unused loops, data...
##  Class(RemoveUnused...,rec(__call__:=function(code)
##      local usedvars,declvars,unusedvars,v;
##      declvars := Set(Flat([
##  			 List(Collect(code,data),d->d.var),
##  			 List(Collect(code,loop),l->l.var)
##  		]));
##      unusedvars := declvars;
##      SubtractSet(unusedvars,usedvars);
##
##			# NOTE: prefixes for data and loop var-names?
##      for v in unusedvars do
##        if v.id[1] = 'D' then
##          code := SubstTopDown(code,data,d->Cond(d.var=v,d.cmd,d));
##        elif v.id[1] = 'i' then
##          code := SubstTopDown(code,loop,d->Cond(d.var=v,d.cmd,d));
##        else # decls
##          code := SubstTopDown(code,decl,d->decl(RemoveList(d.vars,v),d.cmd));
##        fi;
##      od;


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

Local.qmat := Mat(
    [[1,1,1],
     [1,-1,0],
     [1,0,-1]]
);
Local.wmat := Mat(
    [[1,0,1],
     [0,1,0],
     [0,0,1]]);

## C3 = DCT-3
## C4 = DCT-4
## S3 = DST-3
## S4 = DST-4

RulesFor(DTT, rec(
    #
    # DCT-3
    #
    DTT_C3_Base2 := SimpleRule(
	[ DTT_C3, 2, @, @ ],
	(P, C) -> F(2) * Diag(1, cospi(P[3]))
    ),
    DTT_IC3_Base2 := SimpleRule(
	[ DTT_IC3, 2, @, @ ],
	(P, C) -> Diag(1, 1/(2*cospi(P[3]))) * F(2)
    ),
    DTT_C3_Base3 := SimpleRule(
	[ DTT_C3, 3, @, @ ],
	(P, C) -> let(r:=P[3],
	    qmat * DirectSum(I(1), Mat([[ cospi((1+r)/3),  cospi((1-2*r)/3) ],
			             [ cospi((1-r)/3),  cospi((1+2*r)/3) ]])))
    ),
    DTT_C3_Base3 := SimpleRule(
	[ DTT_C3, 3, @, @ ],
	(P, C) -> let(r:=P[3],
	    DirectSum(I(1),
		      1/sinpi(r) *
		      Mat([[ cospi(( 1-2*r)/6),  sinpi((-1+2*r)/3) ],
			   [ cospi((-1+2*r)/6),  cospi(( 1-2*r)/3) ]]) *
		      Mat([[0, sqrt(3)/2],
			   [1, -1/2]]))
	    * qmat)
    ),

    #
    # DST-3
    #
    DTT_S3_Base2 := SimpleRule(
	[ DTT_S3, 2, @, false ],
	(P, C) -> let(r:=P[3],  F(2) * Diag(sinpi(r/2), sinpi(r)))
    ),
    DTT_S3Poly_Base2 := SimpleRule(
	[ DTT_S3, 2, @, true ],
	(P, C) -> let(r:=P[3],  F(2) * Diag(1, 2*cospi(r/2)))
    ),
    DTT_S3Poly_Base3 := SimpleRule(
	[ DTT_S3, 3, @, true ],
	(P, C) -> let(r:=P[3],
	    qmat * DirectSum(I(1), 2* Mat([[ cospi((1+r)/3),  cospi((1-2*r)/3) ],
	       	                           [ cospi((1-r)/3),  cospi((1+2*r)/3) ]]) * wmat))
    ),
    DTT_S3_Base2 := SimpleRule(
	[ DTT_S3, 2, @, false ],
	(P, C) -> let(r:=P[3],  F(2) * Diag(sinpi(r/2), sinpi(r)))
    ),

    #
    # DCT-4
    #
    DTT_C4_Base2 := SimpleRule(
	[ DTT_C4, 2, @, false ],
	(P, C) -> let(r:=P[3],  Mat([[cospi(r/4), cospi(3*r/4)],
		                     [sinpi(r/4), -sinpi(3*r/4)]]))
    ),
    DTT_C4_Base2_LS := SimpleRule(
	[ DTT_C4, 2, @, false ],
	(P, C) -> let(r:=P[3],  Diag(cospi(r/4), sinpi(r/4)) * F(2) * Mat([[1,-1],[0,2*cospi(r/2)]]))
    ),
    DTT_C4Poly_Base2 := SimpleRule(
	[ DTT_C4, 2, @, true ],
	(P, C) -> let(r:=P[3],  F(2) * Mat([[1,-1],[0,2*cospi(r/2)]]))
    ),

    #
    # DST-4
    #
    DTT_S4_Base2 := SimpleRule(
	[ DTT_S4, 2, @, false ],
	(P, C) -> let(r:=P[3],  Mat([[sinpi(r/4), sinpi(3*r/4)],
		                     [cospi(r/4), -cospi(3*r/4)]]))
    ),
    DTT_S4_Base2_LS := SimpleRule(
	[ DTT_S4, 2, @, false ],
	(P, C) -> let(r:=P[3],  Diag(sinpi(r/4), cospi(r/4)) * F(2) * Mat([[1,1],[0,2*cospi(r/2)]]))
    ),
    DTT_S4Poly_Base2 := SimpleRule(
	[ DTT_S4, 2, @, true ],
	(P, C) -> let(r:=P[3],  F(2) * Mat([[1,1],[0,2*cospi(r/2)]]))
    ),

    #
    # Recursive rules
    #
    DTT_DIF_T := rec(
	isApplicable     := P -> not IsPrime(P[2]) and P[2] > 2,
	forTransposition := false,

	baseChange := (self, P, m) >> let(t := P[1], mn := P[2],
	    Cond(Same(t, DTT_C3), B_DCT3_U(mn, m),
		 Same(t, DTT_C4), B_DCT4_U(mn, m),
		 Same(t, DTT_S3), B_DST3_U(mn, m),
		 Same(t, DTT_S4), B_DST4_U(mn, m))),

	recursType := (self, t) >>
	    Cond(Same(t, DTT_C3), DTT_S3,
		 Same(t, DTT_C4), DTT_S3,
		 Same(t, DTT_S3), DTT_S3,
		 Same(t, DTT_S4), DTT_S3),

	allChildren      := (self,P) >> let(
	    type := P[1], MN := P[2], r := P[3],  poly := P[4],
	    divs := DivisorPairs(MN),
	    List(divs, d -> let(j := Ind(d[2]),
		[ DTT(type,                     d[1], fr(d[2], j, r), poly),
		  DTT(self.recursType(type),    d[2], r,              true) ]))),

	rule := (self, P, C, Nonterms) >> let(MN := P[2], N := Rows(C[1]), M := Rows(C[2]),
	                 j := Nonterms[1].params[3].args[2],
	    K(MN, N) *
	    IterDirectSum(j, j.range, C[1]) *
	    Tensor(C[2], I(N)) *
	    self.baseChange(P, M)
	)
    )
));