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

# Cross-compatibility
#

# MAJOR FORWARD DECLARATION HACK
#
# the current module is loaded prior to paradigms.cache, but we rely on the ACache tag
# declared there. To get around this we create a function, which is evaluated on
# execution, to map ACache to a shorter name. This is necessary so that ObjId comparisons
# work at the tag level.
#
# NOTE: this cannot be done with
#
# _ACache := spiral.paradigms.cache.ACach
#
# because assignments are evaluated immediately, and paradigms.cache does not yet exist.
#
# NOTE: we do the same with the CacheSpec and AInplace

_ACache := () -> spiral.paradigms.cache.ACache;

_CacheDesc := () -> spiral.paradigms.cache.CacheDesc;

_AInplace := () -> spiral.paradigms.cache.AInplace;

_AExpRight := () -> spiral.paradigms.cache.AExpRight;

_AIO := () -> spiral.paradigms.common.AIO;

_ABuf := () -> spiral.paradigms.loops.ABuf;
#
## _divisorCacheTags
#
# takes the ACache(N, ...) with the largest N from 't' and
# tests whether the divisor pair 'd' would work well with the cache
#
# We assume a (W x I) (I x W) expansion, so the left side is
# strided. We also assume out-of-place expansion.
#
# At some point, support for an 'Inplace' tag will be added.
#
# the criteria for working well with cache are as follows:
#
# let m = d[1], n = d[2]
#
# LEFT side: kernel size 'm', stride 'n'
# RIGHT side: kernel size 'n', stride 1
#
# cache specs are [B,S,A,R]
# T is the datatype, eg. float, double
#
# let E = B / sizeof(T)
#
# criteria:
#
# --> see the comments in the code
#
# NOTE: additions to support Inplace tag added, but changes incomplete.

_divisorCacheTags := function(d, t)
    local taglist, N, tag, cs, m, n, elems, csize, cond, inplace;

    taglist := t.getTag(_ACache());

    # tagged for inplace?
    inplace := t.getTag(_AInplace());

    # if we don't have any ACache tags, accept this divisor pair.
    if taglist = false then
        return true;
    fi;

    # we need taglist to be a list
    taglist := When(IsList(taglist), taglist, [taglist]);

    # find the max cache level amongst the tags.
    N := Maximum(List(taglist, e -> e.params[1]));

    # use the first tag which matches it.
    tag := Filtered(taglist, e -> e.params[1] = N)[1];

    # extract the cache spec from the tag
    cs := tag.params[2];

    # paranoia
    Constraint(ObjId(cs) = _CacheDesc());

    # d[1] LEFT side expansion
    # d[2] RIGHT side expansion
    m := d[1]; n := d[2];

    # cache sizes
    elems := cs.blksize;
    csize := elems * cs.nsets * cs.assoc;

    # associativity must be at least 2 for outofplace, otherwise read and
    # write overwrite each other
    if inplace = false and cs.assoc < 2 then
        return false;
    fi;

    # if any of these conditions are true -- we allow this split
    cond := [];

    # if the whole thing fits in cache, let it through
    Add(cond, n * m < csize);

    # if right side is larger than cache, left side only uses the assoc so
    # size must be smaller than it.. (also accounting for read/write)
    Add(cond, (n > elems*cs.nsets) and (m <= (cs.assoc / 2)));

    # if right side fits in cache, left side strides partially into the cache,
    # and the size constraint on 'm' is a function of 'n'
    Add(cond, (n <= csize / 2) and (m <= (csize / (2 * n))));

    #PrintLine(d, ": ", cond);

    return ForAny(cond, e -> e);
end;

#
## _dropCacheTags
#
# drop the cache tag if the problem size fits entirely into cache.
#
_dropCacheTags := (t, size, inplace) ->

    (ObjId(t) <> _ACache())           # let other tags through
    or When(inplace = true,  # inplace?
        ((t.params[2].csize) < size), # only let through ACache tags smaller than problem size
        (t.params[2].csize/2 < size)  # (div 2 accounts for rd/wr)
    )
;

#
## _expandRight
#
_expandRight := function(d, t)
    local tag, leftmax;

    tag := t.getTag(_AExpRight());

    # if the tag isn't present, pass through
    if tag = false then
        return true;
    fi;

    # get the limit on left side size, handle multiple tags
    leftmax := When(IsList(tag),
        Minimum(List(tag, e -> e.params[1])),
        tag.params[1]
    );

    # if left expansion exceeds our max, don't allow this breakdown
    return d[1] <= 2^leftmax;
end;

#
## _allowedExpRight
#
#

_allowedExpRight := (nt) -> let(
    t := nt.getTag(_AExpRight()),
    tt := When(IsList(t), t, [t]),
    maxx := Minimum(List(tt, e -> e.params[1])),
    minn := Maximum(List(tt, e -> When(IsBound(e.params[2]) and nt.params[1] > e.params[2], e.params[2], 0))), # optional 2nd param is minimum
#    debugout := Chain(Print(nt.params[1], ": ", minn, "-", maxx, "\n"), 0),
    [minn..maxx]
);

#
## _dropExpandRight
#
# drop all expandright tags at smallest size.
#
# NOTE: this function returns TRUE if we should keep the tag, and FALSE if the
# tag is to be dropped.
#

_dropExpRightTags := (t, size) ->
    ObjId(t) <> _AExpRight()
    or (
        size > 2^Minimum(t.params)
        and not IsBound(t.params[3]) # the presence of a third param says "drop this tag asap"
    )
;

#
#
## DivisorProds
#
# returns lists made up of DivisorsInt(N) which are \leq M and their product = N
#
Declare(_getdivisors);

DivisorProds := function(N, M)
    local l, r;

    # list of valid divisors.
    l := Filtered(DivisorsInt(N), e -> e <= M and e <> 1);

    r := _getdivisors(N, l, [1]);

    # the last element is a 1, we don't want it in the output.
    return List(r, e -> DropLast(e, 1));
end;

_getdivisors := function(N, list, prod)
    local result, e, p;

    p := Product(prod);

    result := List(
        Filtered(list, e -> e * p = N),
        ee -> Concat([ee], prod)
    );

    for e in Filtered(list, i -> i * p < N) do
        Append(result, _getdivisors(N, list, Concat([e], prod)));
    od;

    return result;
end;

_loop := function(p, i)
    local idx, ii, a, c, size;

    # paranoia, 'i' in range.
    Constraint(i >=1 and i <= Length(p));

    idx := [
        XChain([0]),
        GTPar,
        GTVec,
        XChain([1,0,2])
    ];

    #
    a := Product(p{[1..(i-1)]});
    c := Product(p{[i+1..Length(p)]});

    ii := When(a > 1, 1, 0) + When(c > 1, 2, 0) + 1;

    size := [
        [],
        [a],
        [c],
        [a,c]
    ];

    return [idx[ii], size[ii]];
end;

#
## WHT RULES ###############
#

NewRulesFor(WHT, rec(

    #F WHT_BinSplit:
    #F
    #F same as WHT_GeneralSplit, but only allows for binary splits
    #F
    #F   WHT_(2^k) =
    #F     WHT_(2^k1)   tensor I_(2^(k-k1)) *
    #F     I_(2^(k-k2)) tensor WHT_(2^k2)
    #F
    #F We use this rule only, when WHT_GeneralSplit is turned off.
    #F
    #F NOTE: propagates tags, handles expand-right tag.
    WHT_BinSplit := rec (
        info             := "WHT_(2^k) -> (WHT_(2^k1) tensor I) (I tensor WHT_(2^k2))",
        forTransposition := false,
        minSize          := 2,
        maxRadix         := 32,
        applicable       := (self, nt) >> nt.params[1] >= self.minSize,
        requiredFirstTag := ANoTag,

        children := (self, nt) >> List([1.. Minimum(self.maxRadix, nt.params[1]-1)],
            i -> [ WHT(i), WHT(nt.params[1]-i) ]
        ),

        apply := (nt, c, cnt) ->
            Tensor(c[1], I(Rows(c[2])))
            * Tensor(I(Rows(c[1])), c[2])
    ),

    # handle the right expanded WHT.
    WHT_ExpandRight := rec (
        info             := "WHT_(2^k) -> (WHT_(2^k1) tensor I) (I tensor WHT_(2^k2))",
        forTransposition := false,
        minSize          := 2,
        applicable       := (self, nt) >> nt.params[1] >= self.minSize and nt.hasTag(_AExpRight()),

        children := nt -> List(Intersection([1..nt.params[1]-1], _allowedExpRight(nt)),
            i -> [
                WHT(i).withTags(
                    Filtered(nt.getTags(), e ->
                        _dropExpRightTags(e, 2^i)
                    )
                ),
                WHT(nt.params[1]-i).withTags(
                    Filtered(nt.getTags(), e ->
                        _dropExpRightTags(e, 2^(nt.params[1]-i))
                    )
                )
            ]
        ),

        apply := (nt, c, cnt) ->
            Tensor(c[1], I(Rows(c[2])))
            * Tensor(I(Rows(c[1])), c[2])
    ),

    # top level bin split rule (requires ATopLevel tag) which inserts some permutations which are
    # linear on the bits into the ruletree.
    # this rule expects ATopLevel param1 to be a size at which to stop
    # propogating the ATopLevel tag and param2 to be a cache specification
    # of the form [E,S,A,R]
    WHT_BinSplitB := rec (
        info             := "WHT_(2^k) -> (WHT_(2^k1) tensor I) (I tensor WHT_(2^k2))",
        forTransposition := false,
        minSize          := 2,
        applicable       := (self, nt) >> nt.hasTag(ATopLevel) and Rows(nt) > self.minSize,

        children := nt -> let(t := nt.getTag(ATopLevel).params[1],
            List([1..(nt.params[1]-1)], i ->
                When(i >= t or (nt.params[1]-i) >= t,
                    [ WHT(i).withTags([ATopLevel(t)]), WHT(nt.params[1]-i).withTags([ATopLevel(t)]) ],
                    [ WHT(i), WHT(nt.params[1]-i) ]
                )
            )
        ),

        apply := (nt, c, cnt) ->
                Tensor(c[1], I(Rows(c[2])))
                * PushL(
                    CL(Rows(nt), Rows(c[2]), nt.getTag(ATopLevel).params[2]))
                * PushR(
                    CL(Rows(nt), Rows(c[2]), nt.getTag(ATopLevel).params[2]).transpose())
                * Tensor(I(Rows(c[1])), c[2]),

        switch := true
    ),

    WHT_BinSplitB_binloops := rec (
        info             := "WHT_(2^k) -> (WHT_(2^k1) tensor I) (I tensor WHT_(2^k2))",
        forTransposition := false,
        minSize          := 2,
        applicable       := (self, nt) >> nt.hasTag(ATopLevel) and Rows(nt) > self.minSize,

        children := nt -> let(t := nt.getTag(ATopLevel).params[1],
            List([1..(nt.params[1]-1)], i ->
                When(i >= t or (nt.params[1]-i) >= t,
                    [ WHT(i).withTags([ATopLevel(t)]), WHT(nt.params[1]-i).withTags([ATopLevel(t)]) ],
                    [ WHT(i), WHT(nt.params[1]-i) ]
                )
            )
        ),

        apply := function(nt, c, cnt)
            local a, b;

            a := WHT_BinSplit_binloops.apply(nt, c, cnt);

            b := CL(Rows(nt), Rows(c[2]), nt.getTag(ATopLevel).params[2]);

            if b = I(Rows(nt)) then
                return a;
            else
                return Grp(a._children[1] * b) * Grp(b.transpose() * a._children[2]);
            fi;
        end,

        switch := true,
    ),
    #######################################################################################################
    #   tSPL WHT rule
    # tSPL WHT_(2^k) -> (WHT_(2^k1) tensor I) (I tensor WHT_(2^k2))
    WHT_tSPL_BinSplit := rec (
        forTransposition := false,
        minSize          := 2,
        applicable       := (self, nt) >>
            nt.hasTags()
            and nt.params[1] >= self.minSize
            and not nt.hasTag(_AInplace()),
        children         := nt -> List([1..nt.params[1] - 1], i -> [ TTensor(WHT(i), WHT(nt.params[1]-i)).withTags(nt.getTags()) ]),
        apply            := (nt, c, cnt) -> c[1],
        switch := true,
    ),

    #F WHT_Base: WHT_1 = F_2
    #F
    WHT_tSPL_Base := rec(
        switch := false,
        applicable := (t) -> Rows(t) = 2 and t.hasTags(),
        children := t -> [[ TTensorI(F(2), 1, AVec, AVec).withTags(t.getTags()) ]],
        apply := (t, C, Nonterms) -> C[1]
    ),

    #######################################################################################################
    #   tSPL Pease WHT rule
    # Pease tSPL WHT_(2^k) -> \Prod((I tensor F_2)L)
    WHT_tSPL_Pease   := rec (
        forTransposition := false,
        minSize          := 2,
        applicable       := (self, nt) >> nt.hasTags()
                                          and nt.params[1] >= self.minSize
                                          and IsBound(nt.firstTag().legal_kernel)
                                          and ForAny(self.radices, i -> nt.firstTag().legal_kernel(2^i)),

        radices          := [1, 2, 3, 4, 5],

        children         := (self, nt) >>
            let(
                k := nt.params[1],
                streamsize := Cond(
                    nt.hasTags() and IsBound(nt.firstTag().bs),
                    nt.firstTag().bs,
                    NULL
                ),
                ap := Filtered(self.radices, n -> IsInt(k/n) and streamsize >= 2^n),
                List(ap, i -> let(
                    r := i,
                    Cond(
                        IsInt(k/r) and nt.hasTags() and IsBound(nt.firstTag().bs)
                            and nt.firstTag().bs >= 2^r,
                        [ TICompose(var("j"), k/r, TTensorI(WHT(i), 2^(k-r), APar, AVec)).withTags(nt.getTags())],
                        []
                    )
                ))
            ),

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

    #######################################################################################################
    #   tSPL Korn-Lambiotte WHT rule
    WHT_tSPL_KornLambiotte   := rec (
        forTransposition := false,
        minSize          := 2,
        applicable       := (self, nt) >> nt.hasTags()
                                and nt.params[1] >= self.minSize
                                and ForAny(self.radices, i->IsInt(nt.params[1]/i)
                                and IsBound(nt.firstTag().legal_kernel)
                                and nt.firstTag().legal_kernel(2^i)
                            ),

        radices          := [1, 2, 3, 4, 5],

        children         := (self, nt) >> let(
                                k := nt.params[1],
                                ap := Filtered(self.radices, n -> IsInt(k/n) and nt.firstTag().legal_kernel(2^n)),
                                List(ap, i ->
                                    [ TICompose(var("j"), k/i,
                                        TTensorI(WHT(i), 2^(k-i), AVec, APar),
                                        nt.getTags()
                                    )]
                                )
                            ),

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

    # WHT GT rule. Direct copy of DFT_GT_CT minus the twiddles
    #
    # WHT_nm -> GT(WHT_n, ...) * GT(WHT_m, ...)
    #
    # supports the ACache tag.
    #
    WHT_GT_CT := rec(
        switch := true,
        maxSize := false,
        minSize := false,
        forTransposition := false,

        applicable := (self, t) >> let(
            n := Rows(t),

            n > 2
            and (self.maxSize=false or n <= self.maxSize)
            and (self.minSize=false or n >= self.minSize)
            and not IsPrime(n))
            and not t.hasTag(_AInplace()),

        children := (self, t) >> Map2(
            # _divisorCacheTags determines if the split is allowed (according to the cache params)
            Filtered(
                DivisorPairs(Rows(t)),
                d -> _divisorCacheTags(d, t) and _expandRight(d, t)
            ),
            (m,n) -> [
                GT(WHT(Log2Int(m)), XChain([0, 1]), XChain([0, 1]), [n]).withTags(
                    Filtered(t.getTags(), e ->
                        _dropCacheTags(e, m, false) and _dropExpRightTags(e, m)
                    )
                ),
                GT(WHT(Log2Int(n)), XChain([1, 0]), XChain([1, 0]), [m]).withTags(
                    Filtered(t.getTags(), e ->
                        _dropCacheTags(e, n, false) and _dropExpRightTags(e, n)
                    )
                )
            ]
        ),

        apply := (self, t, C, Nonterms) >>  C[1] * C[2]
    ),

    #
    ## WHT_TopInplace
    #
    # the input and output strides must be the same, and intermediate
    # arrays are *full length*, which means the code generated is not
    # recursive.
    #
    # each 'stage' is <= kernel size.
    #
    WHT_TopInplace := rec(
        switch := true,
        maxSize := false,
        minSize := false,

        applicable := (self, t) >> let(
            n := Rows(t),

            t.hasTag(_AInplace())
            and n > 2
            and (self.maxSize=false or n <= self.maxSize)
            and (self.minSize=false or n >= self.minSize)
            and not IsPrime(n)),

        children := meth(self, t)
            local iptag, divs, wa, w, p, i, loopdata, newtags;

            iptag := t.getTag(_AInplace());

            divs := DivisorProds(
                Rows(t),
                When(iptag <> false, iptag.params[1], Rows(t))
            );

            wa := [];
            newtags := Filtered(t.getTags(), e -> ObjId(e) <> _AInplace());

            for p in divs do
                w := [];
                for i in [1..Length(p)] do
                    loopdata := _loop(p, i);
                    Add(w, Inplace(
                        GT(
                            WHT(Log2Int(p[i])),
                            loopdata[1],
                            loopdata[1],
                            loopdata[2]
                        ).withTags(newtags)
                    ));
                od;
                Add(wa, w);
            od;

            return wa;
        end,

        apply := (self, t, C, Nonterms) >>  ApplyFunc(Compose, C)
    ),

    #
    ## WHT_GT_Inplace
    #
    # standard Binsplit CT rule except in place, tags ignored and not
    # propagated.

    WHT_GT_Inplace := rec(
        switch := true,
        maxSize := false,
        minSize := false,

        applicable := (self, t) >> let(
            n := Rows(t),

            t.hasTag(_AInplace()) <> false
            and n > 2
            and not IsPrime(n)
        ),

        children := (self, t) >> Map2(
            DivisorPairs(Rows(t)),
            (m,n) -> [
                GT(WHT(Log2Int(m)), XChain([0, 1]), XChain([0, 1]), [n]).withTags(t.getTags()),
                GT(WHT(Log2Int(n)), XChain([1, 0]), XChain([1, 0]), [m]).withTags(t.getTags())
            ]
        ),

        apply := (self, t, C, nt) >> Inplace(ApplyFunc(Compose, C))
    )
));

#
## GT-WHT RULES ########################
#
# rules for GT(WHT, ...)
# these rules are a copy of the rules for the DFT

NewRulesFor(GT, rec(
    # copy and simplification of GT_DFT_Base2
    # matches GT(WHT(n), ...) where n = 1
    GT_WHT_Base := rec(
        forTransposition := false,
        switch           := false,
        applicable := (self, t) >> let(
            rank := Length(t.params[4]),                                # rank is number of outer loops

            rank = 0                                                    # means no outer loops
            and PatternMatch(t, [GT, WHT, @, @, @, @, @], empty_cx())   # t.rChildren(), used by patternmatch appends t.transposed and t.tags, so GT, with its 4 params, actually has 6.
            and Rows(t.params[1])=2                                     # we're looking for a WHT of size 2
        ),

        apply := (t, C, Nonterms) -> F(2)
    ),

    # matches GT(WHT(n), ...) where n > 1
    GT_WHT_Inplace := rec(
        maxSize       := false,
        minSize       := false,
        minRank       := 0,
        maxRank       := 3,
        forTransposition := false,
        switch           := false,

        applicable := (self, t) >> let(
            rank := Length(t.params[4]),
            wht := t.params[1],

            t.getTag(_AInplace()) <> false
            and rank >= self.minRank
            and rank <= self.maxRank
            and rank = 0 # MRT: we use NthLoop for rank > 0
            and (self.maxSize=false or Rows(wht) <= self.maxSize)
            and (self.minSize=false or Rows(wht) >= self.minSize)
            and PatternMatch(t, [GT, WHT, XChain, XChain, @, @, @], empty_cx())
            and WHT_GT_CT.applicable(wht)
        ),

        children := (self, t) >> let(
            wht := t.params[1],
            g := t.params[2],
            s := t.params[3],
            loop_dims := t.params[4],
            nloops := Length(loop_dims),
            tags := t.getTags(),

            Map2(
                DivisorPairs(Rows(wht)),
                (m,n) -> [
                    GT(
                        WHT(Log2Int(m)),
                        s.composeWith(XChain([0, 1])),
                        s.composeWith(XChain([0, 1])),
                        Concatenation([n], loop_dims)
                    ).withTags(Filtered(t.getTags(), e ->
                        ObjId(e) <> _AInplace() or m > e.params[1]
                    )),
                    GT(
                        WHT(Log2Int(n)),
                        g.composeWith(XChain([1, 0])),
                        s.composeWith(XChain([1, 0])),
                        Concatenation([m], loop_dims)
                    ).withTags(Filtered(t.getTags(), e ->
                        ObjId(e) <> _AInplace() or n > e.params[1]
                    ))
                ]
            )
        ),

        apply := (self, t, C, Nonterms) >> Inplace(C[1] * C[2])
    ),

    GT_WHT_CT := rec(
        maxSize       := false,
        minSize       := false,
        minRank       := 0,
        maxRank       := 3,
        codeletSize   := 32,
        forTransposition := false,
        switch           := false,

        applicable := (self, t) >> let(
            rank := Length(t.params[4]),
            wht := t.params[1],

#            not t.hasTag(_ACache())
#            and
            not t.hasTag(_AInplace())
            and not t.hasTag(_ABuf())
            and rank >= self.minRank
            and rank <= self.maxRank
#            and When(rank>0, t.hasTags(), true)
            and rank = 0 # MRT: we use NthLoop for rank > 0
            and (self.maxSize=false or Rows(wht) <= self.maxSize)
            and (self.minSize=false or Rows(wht) >= self.minSize)
            and PatternMatch(t, [GT, WHT, XChain, XChain, @, @, @], empty_cx())
            and WHT_GT_CT.applicable(wht)
        ),

        children := (self, t) >> let(
            wht := t.params[1],
            g := t.params[2],
            s := t.params[3],
            loop_dims := t.params[4],
            nloops := Length(loop_dims),
            tags := t.getTags(),
            inp := t.getTag(_AInplace()) <> false,

            Map2(
                Filtered(
                    DivisorPairs(Rows(wht)),
                    d -> ( d[1] <= self.codeletSize) and _divisorCacheTags(d, t) and _expandRight(d, t)
                ),
                (m,n) -> [
                    GT(
                        WHT(Log2Int(m)),
                        s.composeWith(XChain([0, 1])),
                        s.composeWith(XChain([0, 1])),
                        Concatenation([n], loop_dims)
                    ).withTags(
                        Filtered(t.getTags(), e ->
                            _dropCacheTags(e, m, inp) and _dropExpRightTags(e, m)
                        )
                    ),

                    GT(
                        WHT(Log2Int(n)),
                        g.composeWith(XChain([1, 0])),
                        s.composeWith(XChain([1, 0])),
                        Concatenation([m], loop_dims)
                    ).withTags(
                        Filtered(t.getTags(), e ->
                            _dropCacheTags(e, n, inp) and _dropExpRightTags(e, n)
                        )
                    )
                ]
            )
        ),

        apply := (self, t, C, Nonterms) >> C[1] * C[2]
    ),

    GT_NthLoop_PassTags := rec(
        switch := false,

        applicable := t -> let(
            rank := Length(t.params[4]),

            rank > 0 and t.hasTags()
        ),

        # restrict to innermost first (i.e. loop interchange)
        # to search over loop orders use [1..nloops]

        # Limit tag reduces the number of loop interchanges.
        # it is useful when you don't want the number of potential
        # ruletrees to explode, and you are not overtly concerned
        # with having access to all the potential ones.
        freedoms := t -> let(
            fr := [1..Length(t.params[4])],
            When(t.hasTag(ALimitNthLoop),
                [[1]],
                [fr]
            )
        ),

        child := (t, fr) -> let(
#            codeletSize := 32,
            spl := t.params[1],
            g := t.params[2],
            s := t.params[3],
            loopid := fr[1],
            gt := GT(spl, g.without(loopid),
                s.without(loopid), ListWithout(t.params[4], loopid)
            ),
            [
                gt.withTags(t.getTags()),
                # Filtered(t.getTags(), e ->
                #     ObjId(e) <> _AInplace() or Rows(gt) >= e.params[1]
                # )),
                InfoNt(loopid)
            ]
        ),

        apply := (t, C, Nonterms) -> let(
            loopid := Nonterms[2].params[1],
            dft := Nonterms[1].params[1],
            g := t.params[2],
            s := t.params[3],
            loop_dims := t.params[4],
            i := Ind(loop_dims[loopid]),

            ISum(i, Scat(s.part(loopid, i, Rows(dft), loop_dims))
            * C[1]
            * Gath(g.part(loopid, i, Cols(dft), loop_dims)))
        )
    ),
));

#   (WHTmn x Ik) -> (WHTm x Ikn) (Im x (WHTn x Ik))
NewRulesFor(TTensorI, rec(
    WHTxI_vecrec := rec(
        switch           := false,
        forTransposition := false,
        minKernel := false,
        maxKernel := false,
        applicable := (self, nt) >> nt.hasTags() and IsVecVec(nt.params) and
            (not IsInt(self.minKernel) or nt.params[1].dims()[2] > self.minKernel) and
            (not IsInt(self.maxKernel) or nt.params[1].dims()[2] <= self.maxKernel),
        children := nt -> let(mn := nt.params[1].dims()[2], k := nt.params[2],
            List(Flat(List(DivisorPairs(mn), i -> let(m := i[1], n := i[2],
                [
                    TCompose([
                        TTensorI(WHT(Log2Int(m)), n*k, AVec, AVec),
                        TTensorI(TTensorI(WHT(Log2Int(n)), k, AVec, AVec), m, APar, APar)
                    ]).withTags(nt.getTags()),
                    TCompose([
                        TTensorI(TTensorI(WHT(Log2Int(n)), k, AVec, AVec), m, APar, APar),
                        TTensorI(WHT(Log2Int(m)), n*k, AVec, AVec)
                    ]).withTags(nt.getTags()),
                 ])
            )), j->[j])
        ),
        apply := (nt, c, cnt) -> c[1]
    )
));


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

# Icode objects that are not cuda-specific (those are defined in cuda_icode.gi)

# A simt index is derived from a SIMT dim tag

Class(simtIdx, Exp, rec(
    computeType := (self) >> TInt,

    size := self >> When(Length(self.args)>0, self.args[1], 0),

    count := (self) >> let(rng := self.get_rng(), When(rng<>[], rng[2]-rng[1]+1, 0 ) ),

    set_rng := meth(self, l, h)
        if Length(self.args)>1 then
            self.args[2] := [l, h];
        else
            Add(self.args, [l, h]);
        fi;
        return self;
    end,

    eval := self >> self,

    can_fold := False,

    get_rng := (self) >> Cond(Length(self.args)>1, When(IsValue(self.args[2]), self.args[2].v, self.args[2]), Length(self.args)>0, [0, self.args[1]-1], [])

));

Class(simtBlockIdx, simtIdx);
Class(simtThreadIdx, simtIdx);

Class(simtBlockIdxX, simtBlockIdx);
Class(simtBlockIdxY, simtBlockIdx);
Class(simtBlockIdxZ, simtBlockIdx);

IsSimtBlockIdx := (idx) -> ObjId(idx) in simtBlockIdx._all_descendants();

Class(simtThreadIdxX, simtThreadIdx);
Class(simtThreadIdxY, simtThreadIdx);
Class(simtThreadIdxZ, simtThreadIdx);

IsSimtThreadIdx := (idx) -> ObjId(idx) in simtThreadIdx._all_descendants();

# Used to identify blocks of code that can be mapped to a kernel
Class(simt_block, chain);

Class(simt_loop, loop, rec(
    drop_range1 := false,

    __call__ := meth(self, simt_idx, loopvar, range, cmd)
        local result;
        Constraint(IsCommand(cmd));
        if IsSymbolic(range) then return loopn(loopvar, range, cmd); fi;
        range := toRange(range);
        if self.drop_range1 and range = 1 then
            return SubstBottomUp(Copy(cmd), @(1, var, e->e=loopvar), e->V(0));
        elif range = 0 then
            return skip();
        elif simt_idx = Ignore then
            return loop(loopvar, range, cmd);
        else
            loopvar.setRange(range);
            range := listRange(range);
            result := WithBases(self,
                rec(operations := CmdOps, simt_idx := simt_idx, cmd := cmd, var := loopvar, range := range));
            loopvar.isLoopIndex := true;
            #loopvar.loop := result;
            return result;
        fi;
    end,

    rChildren := self >> [self.simt_idx, self.var, self.range, self.cmd],
    rSetChild := rSetChildFields("simt_idx", "var", "range", "cmd"),

    ));

# A function that allows declaration specifiers (eg, __global__, __device__, etc..)
Class(specifiers_func, func, rec(
    __call__ := (self, decl_specs, ret, id, params, cmd) >> WithBases(self, rec(
            decl_specs := Checked(IsList(decl_specs), decl_specs),
            ret    := Checked(IsType(ret), ret),
            id     := Checked(IsString(id), id),
            params := Checked(IsList(params), params),
            cmd    := Checked(IsCommand(cmd), cmd),
            operations := CmdOps)),

    rChildren := self >> [self.decl_specs, self.ret, self.id, self.params, self.cmd],
    rSetChild := rSetChildFields("decl_specs", "ret", "id", "params", "cmd"),

    print := (self, i, si) >> Print(self.__name__, "(", self.decl_specs, ", ", self.ret, ", \"", self.id, "\", ", self.params, ", \n",
        Blanks(i+si), self.cmd.print(i+si, si), "\n", Blanks(i), ")", self.printA()),
    )
);

# sync primitives. Their relationship is defined by their sync scope (cluster < block < grid).
# A grid is composed of blocks and blocks of clusters (eg, cuda: grid, blocks, warps).

Class(simtSyncOps, CmdOps, rec(
    \< := (v1, v2) -> Cond(
                            ObjId(v1) = ObjId(v2) or ObjId(v1) = simt_syncgrid, false,
                            ObjId(v1) = skip, true,
                            ObjId(v1) = simt_synccluster and ObjId(v2) <> skip, true,
                            ObjId(v1) = simt_syncblock and not ObjId(v2) in [skip, simt_synccluster], true,
                            false
                        )
    ));

Class(simt_sync, call, rec(
    __call__ := (arg) >> let(o := Inherited(), CopyFields(o, rec(operations := simtSyncOps) ) )
    ));

Class(simt_synccluster, simt_sync);

Class(simt_syncblock, simt_sync);

Class(simt_syncgrid, simt_sync);

# unparsed into C-style printf

Class(cprintf, call);

ext: gi

max_stars_repo_name: mikefranusich/spiral-package-simt

max_issues_repo_name: mikefranusich/spiral-package-simt

max_forks_repo_name: mikefranusich/spiral-package-simt


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

########################################################################
#   rules for A x I, I x A, (A x I)L, (I x A)L
NewRulesFor(TTensorI, rec(
#   loop splitting for A x I
    AxI_LS_L := rec(
        info := "(A_rxs x I_n)(L|I) -> Sum(Sum(SAG))",
        forTransposition := false,
        applicable := nt -> let(P := nt.params,
            P[3] = AVec
            and nt.isTag(1, AVecMemL)
            and IsInt(P[2]/nt.firstTag().v)
            and P[2]/nt.firstTag().v>1
        ),
        children := nt -> [[
            When(nt.numTags() = 1,
                nt.params[1],
                nt.params[1].setTags(nt.params[1].withoutFirstTag())
            )
        ]],
        apply := (nt, C, cnt) -> let(
            v := nt.firstTag().v,
            i := Ind(v),
            m_by_v := nt.params[2]/v,
            j := Ind(m_by_v),
            d := C[1].dims(),
            fr := When(nt.params[4] = AVec,
                fTensor(fId(d[2]), fBase(m_by_v, j), fBase(v, i)),
                fTensor(fBase(m_by_v, j), fBase(v, i), fId(d[2]))
            ),

            ISumLS(j, m_by_v,
                Buf(Scat(fTensor(
                    fId(d[1]), fBase(m_by_v, j), fId(v)
                )))
                * ISum(i, v, Scat(fTensor(
                    fId(d[1]), fBase(v, i))) * C[1] * Gath(fr)
                )
            )
        )

#D        isApplicable := P -> P[3].isVec and Length(P[5]) > 0 and P[5][1].isMemL and IsInt(P[2]/P[5][1].v) and P[2]/P[5][1].v>1,
#D        allChildren := P -> let(pv:=P[5], v:=pv[1].v, d:=P[1].dims(), [[When(Length(pv)=1, P[1], P[1].setpv(Drop(pv, 1)))]]),
#D        rule := (P, C) -> let(
#D                v:=P[5][1].v, i:=Ind(v), m_by_v := P[2]/v, j:=Ind(m_by_v), d:=C[1].dims(),
#D                fr:=When(P[4].isVec, fTensor(fId(d[2]), fBase(m_by_v, j), fBase(v, i)), fTensor(fBase(m_by_v, j), fBase(v, i), fId(d[2]))),
#D                ISumLS(j, m_by_v, Buf(Scat(fTensor(fId(d[1]), fBase(m_by_v, j), fId(v)))) *
#D                    ISum(i, v, Scat(fTensor(fId(d[1]), fBase(v, i))) * C[1] * Gath(fr)))
#D            )
    ),
    AxI_LS_R := rec(
        info := "(I|L)(A_rxs x I_n) -> Sum(Sum(SAG))",
        forTransposition := false,
        applicable := nt -> let(P := nt.params,
            P[4] = AVec
            and nt.isTag(1, AVecMemR)
            and IsInt(P[2]/nt.firstTag().v)
            and P[2]/nt.firstTag().v > 1
        ),
        children := nt -> [[
            When(nt.numTags() = 1,
                nt.params[1],
                nt.params[1].setTags(nt.params[1].withoutFirstTag())
            )
        ]],
        apply := (nt, C, cnt) -> let(
            v := nt.firstTag().v,
            i := Ind(v),
            m_by_v := nt.params[2]/v,
            j := Ind(m_by_v),
            d := C[1].dims(),
            fw := When(nt.params[3] = AVec,
                fTensor(fId(d[1]), fBase(m_by_v, j), fBase(v, i)),
                fTensor(fBase(m_by_v, j), fBase(v, i), fId(d[1]))
            ),

            ISumLS(j, m_by_v,
                ISum(i, v, Scat(fw) * C[1] * Gath(fTensor(fId(d[2]), fBase(v, i))))
                * Buf(Gath(fTensor(
                    fId(d[2]), fBase(m_by_v, j), fId(v)
                )))
            )
        )

#D        isApplicable := P -> P[4].isVec and Length(P[5]) > 0 and P[5][1].isMemR and IsInt(P[2]/P[5][1].v) and P[2]/P[5][1].v>1,
#D        allChildren := P -> let(pv:=P[5], v:=pv[1].v, d:=P[1].dims(), [[When(Length(pv)=1, P[1], P[1].setpv(Drop(pv, 1)))]]),
#D        rule := (P, C) -> let(
#D                v:=P[5][1].v, i:=Ind(v), m_by_v := P[2]/v, j:=Ind(m_by_v), d:=C[1].dims(),
#D                fw:=When(P[3].isVec, fTensor(fId(d[1]), fBase(m_by_v, j), fBase(v, i)), fTensor(fBase(m_by_v, j), fBase(v, i), fId(d[1]))),
#D                ISumLS(j, m_by_v,
#D                    ISum(i, v, Scat(fw) * C[1] * Gath(fTensor(fId(d[2]), fBase(v, i)))) *
#D                            Buf(Gath(fTensor(fId(d[2]), fBase(m_by_v, j), fId(v)))))
#D            )
    ),
    LS_drop := rec(
        info := "drop tag",
        forTransposition := false,
        applicable := nt -> let(P := nt.params,
            (nt.isTag(1, AVecMemL) or nt.isTag(1, AVecMemR))
            and (
                (not IsInt(P[2]/nt.firstTag().v))
                or P[2]/nt.firstTag().v = 1
            )
        ),
        children := nt -> let(P := nt.params,
            [[ TTensorI(P[1], P[2], P[3], P[4]) ]]
        ),
        apply := (nt, C, cnt) -> C[1],

#D        isApplicable := P -> Length(P[5]) > 0 and (P[5][1].isMemL or P[5][1].isMemR) and
#D                            ((not IsInt(P[2]/P[5][1].v)) or P[2]/P[5][1].v=1),
#D        allChildren := P -> [[TTensorI(P[1], P[2], P[3], P[4], [])]],
#D        rule := (P, C) -> C[1],
    )
));

########################################################################
#   TCompose rules
NewRulesFor(TCompose, rec(
    TCompose_LS_L := rec(
        info := "TCompose loop splitting left",
        forTransposition := false,
        applicable := nt -> nt.isTag(1, AVecMemL),
        children := nt -> [ Concat(
            [ nt.params[1][1].setTags(nt.getTags()) ],
            Drop(nt.params[1], 1)
        )],
        apply := (nt, C, cnt) -> Compose(C)
#D        isApplicable := P -> Length(P[2]) > 0 and P[2][1].isMemL,
#D        allChildren := P -> [Concat([P[1][1].setpv(P[2])], Drop(P[1], 1))],
#D        rule := (P, C) -> Compose(C)
    ),
    TCompose_LS_R := rec(
        info := "TCompose loop splitting right",
        forTransposition := false,
        applicable := nt -> nt.isTag(1, AVecMemR),
        children := nt -> [ Concat(
            DropLast(nt.params[1], 1),
            [
                nt.params[1][Length(nt.params[1])].setTags(nt.getTags())
            ]
        )],
        apply := (nt, C, cnt) -> Compose(C)

#D        isApplicable := P -> Length(P[2]) > 0 and P[2][1].isMemR,
#D        allChildren := P -> [Concat(DropLast(P[1], 1), [P[1][Length(P[1])].setpv(P[2])])],
#D        rule := (P, C) -> Compose(C)
    )
));


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

# ==========================================================================
# BaseOperation
# ==========================================================================
Class(BaseOperation, ClassSPL, rec(
    #-----------------------------------------------------------------------
    child       := (self, n) >> self._children[n],
    children    := self >> self._children,
    numChildren := self >> Length(self.children()),
    isReal      := self >> ForAll(self.children(), IsRealSPL),
    setChild := meth(self, n, what) self._children[n] := what; end,
    # -------- Transformation rules support ---------------------------------
    rChildren := ~.children,
    rSetChild := ~.setChild,
));


text:
#
# francy: Interactive Discrete Mathematics in GAP
#

#############################################################################
##
#M  PrintObj( <obj> ) . . . . . . . . . . . . .  override for IsFrancyObjects
##
InstallOtherMethod(PrintObj,
  "a francy object",
  true,
  [IsFrancyObject],
  0,
function(object)
  Print(Sanitize(object));
end);

#############################################################################
##
#M  ViewString( <obj> )  . . . . . . . . . . . . override for IsFrancyObjects
##
InstallOtherMethod(ViewString,
  "a francy object",
  true,
  [IsFrancyObject],
  0,
function(object)
  return (Concatenation( "<",
      CategoriesOfObject( object )[1],
      "/", CategoriesOfObject( object )[2], ">"));
end);

#############################################################################
##
#M  JUPYTER_ViewString( <obj> )  . . . . . . . . override for IsFrancyObjects
##
InstallOtherMethod(JUPYTER_ViewString,
  "a francy object",
  true,
  [IsFrancyObject],
  0,
  ViewString
);

#############################################################################
##
#M  Sanitize( <obj> )  . . . . . . . . simple properties clone for FrancyObjects
##
## This method will clone a FrancyObject and return a record, traversing all the
## components and converting when appropriate.
##
## This method removes components of type IsFunction, as they can't be
## converted to JSON string, converts lists of objects into lists of strings
## and everything else that is not a FrancyObject and therefore unknown!
##
InstallMethod(Sanitize,
  "an object",
  true,
  [IsObject],
  0,
function(object)
  return Sanitize(object, rec());
end);

#############################################################################
##
#M  Sanitize( <obj> )  . . . . . . . . simple properties clone for Records
##
## This method will clone a FrancyObject into the given record
##
InstallOtherMethod(Sanitize,
  "an object, a record",
  true,
  [IsObject,
   IsRecord],
  0,
function(object, record)
  local component, copy, tmp;
  copy := StructuralCopy(object);
  for component in NamesOfComponents(copy) do
    if IsRecord(copy!.(component)) or IsFrancyObject(copy!.(component)) then
      record!.(component) := rec();
      Sanitize(copy!.(component), record!.(component));
    elif IsList(copy!.(component)) and not IsString(copy!.(component)) then
      record!.(component) := [];
      Sanitize(copy!.(component), record!.(component));
    elif IsFunction(copy!.(component)) then
      record!.(component) := NameFunction(copy!.(component));
    else
      record!.(component) := copy!.(component);
    fi;
  od;
  return record;
end);

#############################################################################
##
#M  Sanitize( <obj> )  . . . . . . . . simple properties clone for Records
##
## This method will return a sanitized list from another list
##
InstallOtherMethod(Sanitize,
  "a list, another list",
  true,
  [IsList,
   IsList],
  0,
function(list, result)
  local item;
  for item in list do
    # well, everything that is important for the client is in records
    # everything inside arrays we just convert to string...
    # ...if you wonder, these are most likely known arguments that are stored
    # in order to execute callbacks, so the client does nothing with them
    Add(result, String(item));
  od;
  return result;
end);

#############################################################################
##
#O  MergeRecord( <obj>, <obj> )  . . . . . . . . simple properties merge
##
InstallMethod(MergeObjects,
  "an object, another object",
  true,
  [IsFrancyObject, IsFrancyObject],
  0,
function(dst, src)
  local name;
  for name in NamesOfComponents(src) do
    dst!.(name) := src!.(name);
  od;
  return dst;
end);

#############################################################################
##
#O  GenerateID( ) . . . . . . . . . . . Generates sequential ids for objects
##
InstallMethod(GenerateID,
  "",
  true,
  [],
  0,
function()
  FrancyGeneratedID := FrancyGeneratedID + 1;
  return Concatenation("F", String(FrancyGeneratedID));
end);

max_stars_repo_name: LaGuer/francy

max_issues_repo_name: LaGuer/francy

max_forks_repo_name: LaGuer/francy


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

Class(_ExpandInvars, RuleSet);
RewriteRules(_ExpandInvars, rec(
    shift_depends_left_in_add := ARule(add, [@(1).cond(e-> ObjId(e)<>depends), @(2,depends)],
        e->[@(2).val, @(1).val]),

    shift_depends_left_in_mul := ARule(mul, [@(1).cond(e-> ObjId(e)<>depends), @(2,depends)],
        e->[@(2).val, @(1).val]),

    fuse_depends_in_add := ARule(add, [@(1,depends), @(2,depends)],
        e->[depends(add(@(1).val.args[1], @(2).val.args[1]), Union(@(1).val.args[2].v, @(2).val.args[2].v))]),

    fuse_depends_in_mul := ARule(mul, [@(1,depends), @(2,depends)],
        e->[depends(mul(@(1).val.args[1], @(2).val.args[1]), Union(@(1).val.args[2].v, @(2).val.args[2].v))]),

    leq_depends := Rule([leq, @(1, depends), @(2, depends), @(3, depends)],
        e->depends(leq(@(1).val.args[1], @(2).val.args[1], @(3).val.args[1]),
                Union(@(1).val.args[2].v, @(2).val.args[2].v, @(3).val.args[2].v))),

    imod_depends := Rule([imod, @(1, depends), @(2, depends)],
        e-> depends(imod(@(1).val.args[1],@(2).val.args[1]), Union(@(1).val.args[2].v, @(2).val.args[2].v))),

    extract_depends := Rule([@(0, [add,mul]), @(1,depends)],
        e-> @(1).val),

    div_depends := Rule([div, @(1, depends), @(2, depends)],
        e-> depends(div(@(1).val.args[1],@(2).val.args[1]), Union(@(1).val.args[2].v, @(2).val.args[2].v))),

    idiv_depends := Rule([idiv, @(1, depends), @(2, depends)],
        e-> depends(idiv(@(1).val.args[1],@(2).val.args[1]), Union(@(1).val.args[2].v, @(2).val.args[2].v))),

    tcast_depends := Rule([tcast, @(1), @(2, depends)],
        e -> depends(tcast(@(1).val,@(2).val.args[1]), @(2).val.args[2])),

    add_assoc := RulesStrengthReduce.rules.add_assoc,

    mul_assoc := RulesExpensiveStrengthReduce.rules.mul_assoc,

    depends_deref := Rule( [deref, @(1,depends)],   #One needs to introduce depends memory to handle loads and stores differently
        e -> depends_memory(deref(@(1).val.args[1]), @(1).val.args[2])),

    nth_depends := Rule([nth, @(1, depends), @(2, depends)],
        e -> depends_memory(nth(@(1).val.args[1],@(2).val.args[1]),  Union(@(1).val.args[2].v, @(2).val.args[2].v))),

    vdup_depends_memory := Rule([vdup, @(1,depends_memory), @(2)],
        e -> depends_memory(vdup(@(1).val.args[1], @(2).val), @(1).val.args[2])),

    assign_load_depends := Rule([assign, @(1), @(2, depends_memory)], #This is a load, good invariant!
        e -> assign(@(1).val, depends(@(2).val.args[1], @(2).val.args[2]))),

    mul_depends_accu := ARule(mul, [@(1,depends, e->Length(e.args[2].v)=0), @(2,accu)],
        e -> [accu(mul(@(1).val, @(2).val.rChildren()[1]),
            mul(@(1).val, @(2).val.rChildren()[2]),
            mul(@(1).val, @(2).val.rChildren()[3]))]),

    add_depends_accu := ARule(add, [@(1,depends, e->Length(e.args[2].v)=0), @(2,accu)],
        e -> [accu(add(@(1).val, @(2).val.rChildren()[1]),
            @(2).val.rChildren()[2],
            add(@(1).val, @(2).val.rChildren()[3]))]),

    div_depends_accu := Rule([div, @(1, accu), @(2,depends, e->Length(e.args[2].v)=0)],
        e -> accu(div(@(1).val.rChildren()[1], @(2).val),
            div(@(1).val.rChildren()[2], @(2).val),
            div(@(1).val.rChildren()[3], @(2).val))),

    accu_init := Rule( [accu, @(1,depends), @(2), @(3)], #The initialization step doesn't need to be hoisted
        e -> accu(@(1).val.args[1], @(2).val, @(3).val)),        #Since it is gonna be out of the loop

    #These are security checks
    depends_depends := Rule([depends, depends],
        e->Error("Your code is somehow aliased and SimplifyLoop will probably mess it up. Fix it first then come back!")),

    assign_load_depends := Rule([assign, depends_memory, @], #This is a store invariant! Too good to be true...
        e -> Error("A store invariant has been found which is unlikely to be correct")),
));

Class(_ExpandInvarsUnsafeXXX, RuleSet);
RewriteRules(_ExpandInvarsUnsafeXXX, rec(
    unsafe_memory_depends := Rule(@(1,depends_memory), e-> depends(@(1).val.args[1], @(1).val.args[2]))
));

_ExpandInvarsUnsafe := MergedRuleSet(_ExpandInvarsUnsafeXXX, _ExpandInvars);

Declare(_CreateVirtualVars);
_CreateVirtualVars := function(dep, str, csetable, opts)
    local idx, rem, n, a, v, l;
    if Length(dep.args[2].v)=0 then
        v := SReduce(dep.args[1], opts);
        if IsValue(v) or IsVar(v) or ObjId(v)=param then
            return [v, skip()];
        else
            if csetable<>false then
                l := csetable.cseLookup(v);
                if (l<>false) then
                    return [l, skip()];
                else
                    n := var.fresh_t(str, dep.t);
                    csetable.cseAdd(n, v);
                    return [n, assign(n, v)];
                fi;
            else
                n := var.fresh_t(str, dep.t);
                return [n, assign(n, v)];
            fi;
        fi;
    else
        idx := dep.args[2].v[1];
        rem := Drop(dep.args[2].v, 1);
        [n, a] := TransposedMat(List([1.. idx.range],
                i -> _CreateVirtualVars(depends(SubstVars(Copy(dep.args[1]), tab((idx.id) := V(i-1))), rem), str, csetable, opts)));
        return [virtual_var(n, idx), a];
    fi;
end;

_SimplifyLoop := function(c, to_be_unrolled, opts)
    local free_vars, invars, l, v, assign_accs_init, assign_accs_final, x, accu_init, accu_acc, _DisableAccumulatorBound,
          loop_var, accu_bound, virtual_free_vars, csetable, selectedaccus, allaccus, allincs, i, j, t;

    if not(ObjId(c) in [loop,loopn]) then
            Error("SimplifyLoop has to be called on a loop!");
    fi;

    #Make sure the loop is not simplified again later!
    #Note that this is *critical* : we do not do any dependence analysis and thus
    #it is assumed that all free variables are loop independent which is not the case
    #anymore after the loop has been processed
    c.simplified := true;

    #If the loop is not gonna be unrolled, we will convert it to an accumulator loop
    #That is, we will use a while() loop instead of a for loop and this will allow us
    # to manipulate the loop index aka do a loop with multiple indices that start at different
    #values and increments differently
    if not(to_be_unrolled) then
        c := When(ObjId(c)=loopn,
            SubstBottomUp(c, c.var, x->accu(V(0), V(1), V(c.range))),
            SubstBottomUp(c, c.var, x->accu(V(0), V(1), V(Last(c.range)+1))));
    fi;

    #All free variables, params and values are assumed to be pure loop invariants (depends on nothing)
    free_vars := Set(c.free());
    c := SubstBottomUp(c, @(1, [param, Value], e -> not (ObjId(e.t) = TSym)), x->depends(x, Set([])));
    c := SubstBottomUp(c, @(1,var, e-> e in free_vars ), x->depends(x, Set([])));

    #Propagation of these pure loop invariants allow us to propagate inside the loop accumulators and
    #Define them properly.
    #NOTE: we could have add accumulators that actually depend on internal unrolled loops (as explained later) but then
    #it would actually require as many accumulators as the internal loop range which would waste registers. Note that
    #it's just a feeling, it wasn't tested.

    #This is a bad HACK
    #depends_memory exists because of the DGEMM assign_accs, but it
    #also breaks the DFTs. So when there's no assign acc, we turn them off.
    c := When(Length(Collect(c, assign_acc))>0, _ExpandInvars(c), _ExpandInvarsUnsafe(c));

    #Loop indices that will be unrolled later are somehow loop invariant of the external loop without
    #being loop invariants of the internal loop. As we want to hoist these variables anyhow, we introduce them as depends(exp, [idx])
    virtual_free_vars := When(to_be_unrolled,
        Set(List(Collect(c.cmd, loop), l -> l.var)),
        Set(List(Collect(c.cmd, @@(0, loop,(e, cx) -> (IsBound(cx.unroll_cmd)) and (Length(cx.unroll_cmd)>0))), l -> l.var)));
    c := SubstBottomUp(c, @(1,var, e-> e in virtual_free_vars ), x->depends(x, Set([x])));

    #Propagation now extends dependent variables. Note that dependent variables do not mix with accu so this prevents
    #dependent accumulators
    #Again, the hack.
    c := When(Length(Collect(c, assign_acc))>0, _ExpandInvars(c), _ExpandInvarsUnsafe(c));

    #Kick out trivial invariants. e.g.: constant are always invariants, useless ones!
    c := SubstBottomUp(c, [depends, @(1, [var, Value, param, tcast]), @], x->@(1).val);

    #This rule speeds up processing. It is not necessary but it works well... usually
    #Invariants that do not comprise params are probably constants so forget them!
    c := SubstBottomUp(c, @(0, depends, e-> Length(Collect(e, param))=0), x->@(0).val.args[1]);

    #As scalar increments are essentially free, we merge accumulators that have the same
    #increment if the difference of their initial values is a scalar
    allaccus := Set(Collect(c, accu));
    allincs := Set(List(allaccus, x->x.args[2]));
    for i in allincs do
        selectedaccus := Filtered(allaccus, x-> x.args[2] = i);
        for j in DropLast(selectedaccus,1) do
            t := j.args[1]-Last(selectedaccus).args[1];
            if ObjId(t)=Value then
                j.substitute := add(Last(selectedaccus), t);
            fi;
        od;
    od;
    c := SubstBottomUp(c, @(1, accu, e ->IsBound(e.substitute)), e->e.substitute);

    #If we are to transform in a while loop, we only really need one loop counter
    #Therefore, we kick out all loop counters but one
    if not(to_be_unrolled) then
        x := true;
        _DisableAccumulatorBound := function(e)
            if x then
                x:=false;
                e.isLoopBound := true;
                return e;
            else
                return accu(e.rChildren()[1], e.rChildren()[2], V(0));
            fi;
        end;
        SubstBottomUp(c.cmd, accu, e -> _DisableAccumulatorBound(e));
    fi;

    #Create new variables for all invariants
    csetable := CSE.init();
    l := Set(Collect(c, depends));
    [v, invars] := let(a := List(l, x->_CreateVirtualVars(x, "ivr", csetable, opts)),
        When(Length(a)>0, TransposedMat(a), [[],[]]));
    c := SubstBottomUp(c, @(1, depends), e->v[Position(l,e)]);

    #Same for all assign_accs!
    l := Set(List(Collect(c, [assign_acc, depends_memory, @]), x->x.loc));
    [v, assign_accs_init] := let(a := List(l, x->_CreateVirtualVars(x, "acc", false, opts)),
        When(Length(a)>0, TransposedMat(a), [[],[]]));
    assign_accs_final := SubstTopDownNR(Copy(assign_accs_init), assign, x->assign(x.exp, x.loc));
    c := SubstBottomUp(c, @(1, depends_memory), e->v[Position(l,e)]);

    #If the loop is to be unrolled then we're done
    if to_be_unrolled then
        return chain(invars, assign_accs_init, c, assign_accs_final);
    else #if not, then we need to finish the transformation in a while loop

        #drop the loop itself, we will replace it by a doloop later
        c := c.cmd;

        #Replace all accus by new vars
        #retrieve the loop bound and set all bounds to V(0) so they can merge in the Set()
        accu_bound := Copy(Collect(c, @(1, accu, x->IsBound(x.isLoopBound)))[1]);
        [loop_var, accu_bound] := [accu(accu_bound.rChildren()[1], accu_bound.rChildren()[2], V(0)),
            accu_bound.rChildren()[3]];
        c := SubstBottomUp(c, [accu, @(1), @(2), @(3)], x->accu(@(1).val, @(2).val, V(0)));
        l := Set(Collect(c, accu));
        v := List(l, x->var.fresh_t("accu", x.t));
        accu_init := List([1..Length(l)], i->assign(v[i], l[i].rChildren()[1]));
        #Copy here prevents aliasing with previous
        accu_acc := List([1..Length(l)], i->assign_acc(Copy(v[i]), l[i].rChildren()[2]));

        loop_var := let(p:=FirstPosition(l, x->(x=loop_var)), v[p]);
        c := SubstBottomUp(c, @(1, accu), e->v[Position(l,e)]);

        #if we cross the unroll_cmd boundary, we need to flag the new variables as alive
        #so copypropagation doesn't mess with it later
        for x in Collect([assign_accs_init, invars, accu_init], assign) do x.loc.live_out := true ; od;

        c := chain(unroll_cmd(chain(invars, assign_accs_init, accu_init)),
             doloop(loop_var, accu_bound, chain(c, accu_acc)), unroll_cmd(chain(assign_accs_final)));
        return c;
    fi;
end;

SimplifyLoops := function(c, opts)
    c := FlattenCode0(c); #Flatten chains, keep unroll_cmds
    c := SReduce(c, opts); #Simplify index expressions

    c := SubstTopDown(c, @.cond( #Remove all aliases!
            e->IsRec(e) and IsBound(e.rChildren) and IsBound(e.from_rChildren)),
        e -> e.from_rChildren(List(e.rChildren(), Copy)));

    c := SubstTopDown(c, @(1, [loop,loopn], e->not(IsBound(e.simplified))),
        (e, cx) -> _SimplifyLoop(e, (IsBound(cx.unroll_cmd)) and (Length(cx.unroll_cmd)>0), opts));
    return c;
end;


text:
version0000640030014500116100000000001612523536103011462 0ustar  dgentryenggfibertv-40.2
rootfs.squashfs0000640030014500116100000000001212523536621013146 0ustar  dgentryengrootfs.imgvmlinuz0000640030014500116100000000001212523536672011510 0ustar  dgentryengkernel.imgloader.bin0000640030014500116100000000001212523536403012011 0ustar  dgentryengloader.img

max_stars_repo_name: DentonGentry/gfiber-platform

max_issues_repo_name: DentonGentry/gfiber-platform

max_forks_repo_name: DentonGentry/gfiber-platform


text:
# At top level, global variables are declared when they are assigned, so one only writes
global_var := 1;

# In a function, local variables are declared like this
func := function(n)
    local a;
    a := n*n;
    return n + a;
end;

# One can test whether a variable is assigned
IsBound(global_var);
# true;

# And destroy a variable
Unbind(global_var);

# This works with list elements too
u := [11, 12, , 14];
IsBound(u[4]);
# true
IsBound(u[3]);
# false
Unbind(u[4]);

ext: gap

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


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

RewriteRules(RulesStrengthReduce, rec(

################## HALF ####################
    vunpacklo_half := Rule( vunpacklo_half,
        e -> vextract_half(vtrnq_half(e.args[1], e.args[2]), [0])),
    vunpackhi_half := Rule( vunpackhi_half,
        e -> vextract_half(vtrnq_half(e.args[1], e.args[2]), [1])),

################## NEON ####################
    vpacklo_neon := Rule( vpacklo_neon,
        e -> vextract_neon_4x32f(vuzpq_32f(e.args[1], e.args[2]), [0])),
    vpackhi_neon := Rule( vpackhi_neon,
        e -> vextract_neon_4x32f(vuzpq_32f(e.args[1], e.args[2]), [1])),
    vunpacklo_neon := Rule( vunpacklo_neon,
        e -> vextract_neon_4x32f(vzipq_32f(e.args[1], e.args[2]), [0])),
    vunpackhi_neon := Rule( vunpackhi_neon,
        e -> vextract_neon_4x32f(vzipq_32f(e.args[1], e.args[2]), [1])),
    vtransposelo_neon := Rule( vtransposelo_neon,
        e -> vextract_neon_4x32f(vtrnq_32f(e.args[1], e.args[2]), [0])),
    vtransposehi_neon := Rule( vtransposehi_neon,
        e -> vextract_neon_4x32f(vtrnq_32f(e.args[1], e.args[2]), [1])),
    vunpacklolo2_neon := Rule( vunpacklolo2_neon,
        e -> vextract_neon_4x32f(vtrnq_32f(vpacklo_neon(e.args[1], e.args[2]), vpackhi_neon(e.args[1], e.args[2])), [0])),
    vunpackhihi2_neon := Rule( vunpackhihi2_neon,
        e -> vextract_neon_4x32f(vtrnq_32f(vpacklo_neon(e.args[1], e.args[2]), vpackhi_neon(e.args[1], e.args[2])), [1])),

));


text:
# Copyright 2018-2019, Carnegie Mellon University
# See LICENSE for details

Class(TState, VStack, rec(
 dims := self >> [0, FoldL(self._children,(a,b)->When(not IsSymbolic(b.dims()[2]), Max2(a,b.dims()[2]), a), 0) ],
 isCommand := true
));

HCOLSumsGen.TState := (self, o, opts) >> let(
      ch := o.children(),
      ApplyFunc(TState, List([1..Length(ch)], i->OLCompose(
        ScatHUnion(Rows(o), Rows(ch[i]),  Sum(List(ch{[1..i-1]}, child->child.dims()[1])), 1),
	self(ch[i], opts)
      )))
    );

HCOLCodegen.TState := (self, o, y, x, opts) >> List(o.children(), c -> self(c, opts.state, x, opts));

max_stars_repo_name: spiral-software/spiral-package-hcol

max_issues_repo_name: spiral-software/spiral-package-hcol

max_forks_repo_name: spiral-software/spiral-package-hcol


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

# NOTE: too slow
PackCode := function(c)
     local pulled, decls, datas, dims;
     if IsBound(c.dimensions) then dims := c.dimensions; fi;
     pulled := PullBU(Copy(c), decl, e -> e.cmd, e -> e.vars);
     decls := Set(Concatenation(pulled[1]));
     c := pulled[2];

     pulled := Pull(Copy(c), data, e -> e.cmd, e -> [e.var, e.value]);
     datas := pulled[1];
     c := pulled[2];

     c := SubstBottomUp(c, chain, e -> e.flatten());
     c := FoldL(datas, (cod,d) -> data(d[1], d[2], cod), decl(decls, c));
     if IsBound(dims) then c.dimensions := dims; fi;
     return c;
end;

FoldIf := c -> SubstTopDown(c,
    @(1).target(IF).cond(e->IsValue(e.cond)),
    e -> When(e.cond.v=0, e.else_cmd, e.then_cmd));

DeclareHidden := function(code)
    local free, v, dims;
    free := code.free();
    if IsBound(code.dimensions) then dims := code.dimensions; fi;
    for v in free do
        if IsBound(v.value) then code := data(v, v.value, code); fi;
    od;
    if IsBound(dims) then code.dimensions := dims; fi;
    return code;
end;

# This eliminates aliasing between subexpressions of different commands
# (ie pointers can alias, and this is bad)
UntangleChain := c -> SubstTopDownNR(c, chain, e -> chain(List(e.cmds, Copy)));

# reduce complex constants with im(c)=0 to float
_prepConst := e -> Cond(
    e.t <> TComplex, e,
    IsCyc(e.v) and Im(e.v)=0, TDouble.value(e.v),
    IsComplex(e.v) and ImComplex(e.v) = 0, TDouble.value(ReComplex(e.v)),
    e);

Class(_dropTRealInRefs, RuleSet);
RewriteRules(_dropTRealInRefs, rec(
    TReal_left  := ARule(ListClass, [ [@(1, TPtr), @(2, AtomicTyp, x->x=TReal), ...],
        [@(3, TPtr, x->x.alignment=@(1).val.alignment), T_Real, ...]],
        e -> [@(3).val]),
    TReal_right := ARule(ListClass, [ [@(1, TPtr), T_Real, ...],
        [@(2, TPtr), @(3, AtomicTyp, x->x=TReal and @(1).val.alignment=@(2).val.alignment), ...]],
        e -> [@(1).val]),
    TVect_TReal_left  := ARule(ListClass, [ [@(1, TPtr), [TVect, @(2, AtomicTyp, x->x=TReal), @(3)], ...],
        [@(4, TPtr), [TVect, T_Real, @(5).cond(x->x=@(3).val and @(1).val.alignment=@(4).val.alignment)], ...]],
        e -> [@(4).val]),
    TVect_TReal_right := ARule(ListClass, [ [@(1, TPtr), [TVect, T_Real, @(2)], ...],
        [@(3, TPtr), [TVect, @(4, AtomicTyp, x->x=TReal), @(5).cond(x->x=@(2).val and @(1).val.alignment=@(3).val.alignment)], ...]],
        e -> [@(1).val]),
));

#F Compile(<code>)
#F
#F Fully unroll and optimize <code>
#F
Class(Compile, rec(
    datas := 0,
    decls := 0,
    refs := 0,

    status := Ignore,
    timingStatus := Ignore,

    pullDataDeclsRefs := meth(self, code)
        local pulled, c, d, p, id, typecasts, tcast, ref, loopvars, doNotScalarize;
        self.free := code.free();
        pulled := PullBU(code, @(1, [decl, data, nth, deref]),
            e -> When(IsBound(e.cmd), e.cmd, e), e -> e);

        self.decls := Set([]);
        self.datas := tab();
        self.refs  := tab();

        loopvars := Pull(code, @@(1).cond( (x, cx)->IsLoop(x) ),
            e -> e, (cx, e) -> e.var)[1];

        doNotScalarize := (ref) -> ref=false or Collect(ref, @(1, [var, param], x -> (IsLoopIndex(x) and not(x in loopvars)) or x _is param))<>[];

        for p in pulled[1] do
            if   ObjId(p)=decl then UniteSet(self.decls, p.vars);
            elif ObjId(p)=data then self.datas.(p.var.id) := p.value; if IsArrayT(p.var.t) then p.var.value := p.value; fi;
            # nth, deref
            elif IsVar(p.loc) then
                id := p.loc.id;
                if not IsBound(self.refs.(id)) then self.refs.(id) := Set([]); fi;
                AddSet(self.refs.(id), TPtr(p.t));
                if doNotScalarize(p) then # prevent scalarization when there is a dependency on outer loop variable or param
                    AddSet(self.refs.(id), TPtr(TVoid));
                fi;
            fi;
        od;

        # Pull out type casts and figure out granularity of accesses to given variable
        # the complicated condition tries to figure out <..what?..>
        typecasts := Pull(code,
            @@(1,var, (e,cx) -> IsVar(e) and IsArrayT(e.t) and IsBound(cx.tcast) and
                                cx.tcast<>[] and Last(cx.tcast).args[2].t in [e.t, e.t.toPtrType()] ),
            e -> e,
            (cx, e) -> [e, Last(cx.tcast), Cond(IsBound(cx.nth) and Length(cx.nth)>0, Last(cx.nth), false)])[1];

        for p in typecasts do
            [id, tcast, ref] := [p[1].id, p[2], p[3]];
            if not IsBound(self.refs.(id)) then self.refs.(id) := Set([]); fi;
            AddSet(self.refs.(id), tcast.args[1]);
            # below prevents scalarization if pointer arithmetic was used or
            # if there is dependency on outer loop variable or param
            if not IsVar(tcast.args[2]) or doNotScalarize(ref) then
                AddSet(self.refs.(id), TPtr(TVoid));
            fi;
        od;
        # TReal pops up in diffrent places
        for id in UserNSFields(self.refs) do
            self.refs.(id) := _dropTRealInRefs(self.refs.(id));
        od;
        # plugs in data definitions into code, this process could be recursive, since one datavar can depend on another
        return SubstTopDown(pulled[2], @(1,var,e->IsBound(self.datas.(e.id))), e->self.datas.(e.id));
    end,

    scalarizationBarriers := [call, fcall],

    fastScalarize := meth(self,code)
        local d,len,tentry, barriers, typecasts, newlen, newt, newscal;
        barriers := Collect(code, @(1, self.scalarizationBarriers));
        self.doNotScalarize := Union(
            # do not scalarize arrays from barriers (like from function calls)
            Union(List(barriers, c->Collect(c, @(1,var,x->IsArrayT(x.t))))),
            # arrays with vector accesses of different granularity
            Set(Filtered(self.decls, v -> not IsBound(self.refs.(v.id)) or
                    (IsBound(self.refs.(v.id)) and (Length(self.refs.(v.id)) > 1 or ForAny(self.refs.(v.id), IsUnalignedPtrT))) or
                    (IsBound(v.doNotScalarize) and v.doNotScalarize)))
        );

        self.scalarized := Set([]);
        newscal := [];
        for d in Difference(self.decls, self.doNotScalarize) do
            newt := self.refs.(d.id)[1];
            if IsArrayT(d.t) and IsPtrT(newt) then
                newt := newt.t; # base type of ptr
                newlen := d.t.size * When(IsVecT(d.t.t), d.t.t.size, 1) / When(IsVecT(newt), newt.size, 1);
                if IsSymbolic(newlen) or IsValue(newlen) then newlen := newlen.ev(); fi;
                tentry := List([0..newlen-1], i -> var.fresh_t("scal", newt));
                Append(newscal, tentry);
                d.value := Value(TArray(newt, newlen), tentry);
            fi;
        od;
        self.scalarized := Set(newscal);

        return code;
    end,

    # NOTE: make this more general, in particular, if variable does not appear
    # on the left-hand side of the assignment it won't be declared, scalar
    # variables won't be declared if they are only used inside call() (pathological case)
    # Vectors can in fact be only used inside call, so this situtation is handled correctly.
    declareVars := meth(self, code)
        local vars, vects;
        vars  := Filtered(Set(ConcatList(Collect(code, @.cond(e->IsCommand(e) and IsBound(e.op_out))), e->e.op_out())), IsVar);
        vects := Filtered(code.free(), x -> IsArray(x.t));

        IntersectSet(vects, self.decls);
	SubtractSet(vects, self.free);
	SubtractSet(vars,  self.free);

        if Length(vars) > 0 then code := decl(vars, code); fi;
        if Length(vects) > 0 then code := decl(vects, code); fi;
        return code;
    end,

    # Old and slow scalarizer no longer works, but old scalarizer did not make an
    # assumption of "constant geomery code" and the current scalarizer does make that
    # assumption. It is becoming a limitation for certain applications, so we might need
    # to redesign our current fast scalarizer to be more like the old one.
    #
    #D scalarize := (self, c) >> Scalarize(c, self.decls),

    showTimes := meth(self)
        local i;
        Print("times := [\n");
        for i in [1..Length(self.times)] do
	    PrintEval("$1,  $2,\n", i, StringDouble("%.3g", self.times[i]));
	od;
        Print("];\n");
    end,

    __call__ := meth(self, c, opts)
        local dims, root, stage, compileStrategy, t, i;

        if IsBound(c.dimensions) then dims := c.dimensions; fi;
        if IsBound(c.root) then root := c.root; fi;

        if opts.printWebMeasure then Print("web:measure\n"); fi;
        self.curcode := [Copy(c)];

        # self.times keeps compilation time information to be used for profiling
        if not IsBound(self.times) or Length(self.times)<>Length(opts.compileStrategy) then
            self.times := Replicate(Length(opts.compileStrategy), 0.0); fi;

        for i in [1..Length(opts.compileStrategy)] do
            stage := opts.compileStrategy[i];
            if IsMeth(stage) and IsCallableN(stage, 2) then
                [c,t] := UTimedAction(stage(self, c, opts));
            elif IsMeth(stage) then
                [c,t] := UTimedAction(stage(self, c));
            elif IsCallableN(stage, 2) then
                [c,t] := UTimedAction(stage(c, opts));
            else
                [c,t] := UTimedAction(stage(c));
            fi;
            Add(self.curcode, Copy(c));
            self.times[i] := self.times[i] + t;
	    self.timingStatus(i, stage, t); # print timing info

            if opts.printWebMeasure then Print("web:measure\n"); fi;
        od;

        if IsBound(dims) then c.dimensions := dims; fi;
        if IsBound(root) then c.root := root; fi;
        return c;
    end
));


text:
Ackermann := function( m, n )

  if m = 0 then
      return n + 1;
  elif m > 0 and n = 0 then
      return Ackermann( m - 1, 1 );
  fi;

  return Ackermann( m - 1, Ackermann( m, n - 1 ) );

end;

max_stars_repo_name: Mynogs/Algorithm-Implementations

max_issues_repo_name: Mynogs/Algorithm-Implementations

max_forks_repo_name: Mynogs/Algorithm-Implementations


text:
InstallMethod(NilpotentChv,
              "NilpotentChv element in simple Lie type",
              [IsChevalleyAdj,IsList],
              function(sys,coeffs)
              local object,
              B,admat,l;

              if Filtered(coeffs,i->not i[2] in Integers)=[] then
                  coeffs:=List(coeffs,i->[i[1],One(ring(sys))*i[2]]);
              elif Filtered(coeffs,i->not i[2] in ring(sys))<>[] then
                  Error("NilpotentChv: The coefficients ar not in indicated ring.");
              fi;

              object:=Objectify(NewType(NewFamily("NilpotentChvFamily"),
                                        IsAttributeStoringRep and
                                        IsNilpotentChv and
                                        IsAdditiveElementWithInverse and
                                        IsAdditiveElementWithZero),
                                rec());

              SetchevalleyAdj(object,sys);

              # Simplified form of the element (negative root vectors are also allowd "2*Len..")
              coeffs:=Filtered(List([1..2*Length(positiveRoots(sys))],
                                    i->[i,One(ring(sys))*Sum(List(Filtered(coeffs,
                                                                           j->j[1]=i),
                                                                  k->k[2]))]),
                               l->l[2]<>Zero(ring(sys)));

              Setcoefficients(object,coeffs);

              B:=Basis(lieAlgebra(sys));
              if coeffs<>[] then SetLieAlgebraElement(object,Sum(coeffs,i->i[2]*B[i[1]]));
              else SetLieAlgebraElement(object,Zero(lieAlgebra(sys))); fi;
              SetLieAlgebraCoeffs(object,Coefficients(B,LieAlgebraElement(object)));

#              l:=Length(PRoots(R));
#              admat:=AdjointMatrix(B,LieAlgebraElement(e));
#              SetPositiveAde(object,List([1..l],i->admat[i]{[1..l]}));

              SetName(object,Concatenation("<nilpotent element for ",
                                           type(sys),String(rank(sys)),
                                           " in characteristic ",String(Characteristic(sys)),">"));
              return object;
end);

InstallMethod(Ascend,
              "If possible embed the parameters in a polynomial ring with coefficients ring the base ring",
              [IsNilpotentChv,IsPolynomialRing],
              function(e,inel)
              local sys,
              check;

              sys:=chevalleyAdj(e);

              if not IsPolynomialRing(inel) then
                  Error("Ascend NilpotentChv: The new ring is not a polynomial ring to embed in!\n");
              elif CoefficientsRing(inel)<>ring(sys) then
                  Error("Ascend NilpotentChv: CoefficientsRing and ring(sys) do not coincide!\n");
              fi;

              return NilpotentChv(ChevalleyAdj(sys,inel),List(coefficients(e),i->[i[1],One(inel)*i[2]]));
end);

InstallMethod(Ascend,
              "If possible embed the parameters in a polynomial ring with coefficients ring the base ring",
              [IsNilpotentChv,IsChevalleyAdj],
              function(e,sys)
              local syse,
              check;

              syse:=chevalleyAdj(e);

              if not IsPolynomialRing(ring(sys)) then
                  Error("The new ring is not a polynomial ring to embed in.");
              elif CoefficientsRing(ring(sys))<>ring(syse) then
                  Error("CoefficientsRing and ring(syse) do not coincide.");
              fi;

              return NilpotentChv(sys,List(coefficients(e),i->[i[1],One(ring(sys))*i[2]]));
end);

InstallMethod(\*,
              "multiplication by scalar for nilpotent elements s*e",
              [IsRingElement,IsNilpotentChv],
              function(s,e)
              local sys;

              sys:=chevalleyAdj(e);
              if not s in ring(sys) then
                  Error("Scalar*NilpotentChv: Scalar not in ring of nilpotent!\n");
              fi;

              return NilpotentChv(sys,List(coefficients(e),i->[i[1],s*i[2]]));
end);

InstallMethod(\+,
              "addition for nilpotent elements e1+e2",
              [IsNilpotentChv,IsNilpotentChv],
              function(e1,e2)
              local sys1,sys2;

              sys1:=chevalleyAdj(e1);
              sys2:=chevalleyAdj(e2);
              if type(sys1)<>type(sys2) or
                 rank(sys1)<>rank(sys2) or
                 ring(sys1)<>ring(sys2) then
                  Error("NilpotentChv+NilpotentChv: Not in the same family!\n");
              fi;

              return NilpotentChv(sys1,Concatenation(coefficients(e1),coefficients(e2)));
end);

InstallMethod(AdditiveInverseMutable,
              "additive inverse for nilpotent elements -e1",
              [IsNilpotentChv],
              function(e)
              local sys;

              sys:=chevalleyAdj(e);

              return NilpotentChv(sys,List(coefficients(e),i->[i[1],-i[2]]));
end);

InstallMethod(ZeroMutable,
              "Zero vector of the Lie algebra",
              [IsNilpotentChv],
              function(e)
              local sys;

              return NilpotentChv(chevalleyAdj(e),[]);
end);


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

reps:=[[[1,1],[2,1],[3,1],[4,1],[5,1],[6,1]],        # E_6
       [[3,1],[5,1],[1,1],[6,1],[9,1],[10,1],[2,1]], # E_6(a_1)
       [[1,1],[3,1],[4,1],[2,1],[5,1]],              # D_5
       [[13,1],[1,1],[15,1],[6,1],[14,1],[4,1]],     # E_6(a_3)
       [[1,1],[3,1],[8,1],[10,1],[2,1],[5,1]],       # D_5(a_1)
       [[1,1],[3,1],[4,1],[5,1],[6,1]],              # A_5
       [[1,1],[3,1],[4,1],[2,1],[6,1]],              # A_4A_1
       [[3,1],[4,1],[2,1],[5,1]],                    # D_4
       [[1,1],[3,1],[4,1],[2,1]],                    # A_4
       [[3,1],[8,1],[10,1],[2,1],[5,1]],             # D_4(a_1)
       [[1,1],[3,1],[4,1],[6,1]],                    # A_3A_1
       [[1,1],[3,1],[5,1],[6,1],[2,1]],              # A_2^2A_1
       [[1,1],[3,1],[4,1]],                          # A_3
       [[2,1],[4,1],[1,1],[6,1]],                    # A_2A_1^2
       [[1,1],[3,1],[5,1],[6,1]],                    # A_2^2
       [[1,1],[3,1],[2,1]],                          # A_2A_1
       [[1,1],[3,1]],                                # A_2
       [[1,1],[4,1],[6,1]],                          # A_1^3
       [[1,1],[6,1]],                                # A_1^2
       [[1,1]]                                       # A_1
       ];

I:=[[],[4],[],[2,3,5],[4],[],[],[],[],[4],
    [],[],[],[],[],[],[],[],[],[]];

#
# The parabolic where the representative is distinguished
#
II:=[[1,2,3,4,5,6],
     [1,2,3,4,5,6],
     [1,2,3,4,5],
     [1,2,3,4,5,6],
     [1,2,3,4,5],
     [1,3,4,5,6],
     [1,2,3,4,6],
     [2,3,4,5],
     [1,2,3,4],
     [2,3,4,5],
     [1,3,4,6],
     [1,2,3,5,6],
     [1,3,4],
     [1,2,4,6],
     [1,3,5,6],
     [1,2,3],
     [1,3],
     [1,4,6],
     [1,6],
     [1]];

labels:=["E_6",
         "E_6(a_1)",
         "D_5",
         "E_6(a_3)",
         "D_5(a_1)",
         "A_5",
         "A_4A_1",
         "D_4",
         "A_4",
         "D_4(a_1)",
         "A_3A_1",
         "A_2^2A_1",
         "A_3",
         "A_2A_1^2",
         "A_2^2",
         "A_2A_1",
         "A_2",
         "A_1^3",
         "A_1^2",
         "A_1"];

tori:=[
    [2,2,2,2,2,2],
    [2,2,2,0,2,2],
    [2,2,2,2,2,-10],
    [2,0,0,2,0,2],
    [2,2,2,0,2,-7],
    [2,-9,2,2,2,2],
    [2,2,2,2,-7,2],
    [-6,2,2,2,2,-6],
    [2,2,2,2,-6,0],
    [-4,2,2,0,2,-4],
    [2,-3,2,2,-4,2],
    [2,2,2,-5,2,2],
    [2,-3,2,2,-3,0],
    [2,2,-3,2,-3,2],
    [2,0,2,-4,2,2],
    [2,2,2,-3,0,0],
    [2,0,2,-2,0,0],
    [2,-1,-2,2,-2,2],
    [2,0,-1,0,-1,2],
    [2,0,-1,0,0,0]
    ];

#blocks:=[];
blocks:=[[[19, 1], [15, 2], [9, 3], [2, 1]], [[9, 8], [6, 1]], [[9, 7], [7, 2], [1, 1]], [[9, 4], [7, 1], [6, 4], [3, 3], [2, 1]], [[9, 3], [8, 2], [6, 4], [3, 2], [2, 2], [1, 1]], [[9, 3], [8, 2], [6, 4], [3, 2], [2, 1], [1, 3]], [[9, 1], [8, 2], [7, 1], [6, 2], [5, 3], [3, 6], [1, 1]], [[9, 2], [7, 7], [3, 1], [1, 8]], [[9, 1], [7, 5], [5, 3], [3, 5], [1, 4]], [[7, 1], [6, 2], [5, 6], [3, 9], [1, 2]], [[7, 1], [6, 2], [5, 3], [4, 6], [3, 4], [2, 2], [1, 4]], [[3, 24], [2, 3]], [[7, 1], [5, 5], [4, 8], [3, 1], [1, 11]], [[3, 23], [2, 4], [1, 1]], [[3, 23], [2, 1], [1, 7]], [[3, 22], [2, 2], [1, 8]], [[3, 21], [1, 15]], [[3, 13], [2, 14], [1, 11]], [[3, 8], [2, 16], [1, 22]], [[3, 1], [2, 20], [1, 35]]];

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

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


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

Class(RulesSIMTFission, RuleSet);

RewriteRules(RulesSIMTFission, rec(
    FissionSIMTISum := Rule([@(1, SIMTISum), @(2, Compose, e->ForAll(e.children(), c->ObjId(c) = SIMTISum))],
        e-> let(gath := Gath(fTensor(fBase(@(1).val.var), fId(Cols(@(2).val.child(1))))),
                scat := Scat(fTensor(fBase(@(1).val.var), fId(Cols(@(2).val.child(1))))),
                c1 := @(2).val.child(1) * gath,
                cn := scat * Last(@(2).val.children()),
                cms := List(DropLast(Drop(@(2).val.children(), 1), 1), c->scat * c *gath),
                Compose(List([c1]::cms::[cn], c -> SIMTISum(@(1).val.simt_dim, @(1).val.var, @(1).val.domain, c))))),
));

RewriteRules(RulesSums, rec(
    OO_Gath := Rule([@(1, Compose), OO, Gath],
        e -> ApplyFunc(OO, @(1).val.dims())),
    Gath_OO := Rule([@(1, Compose), Scat, OO],
        e -> ApplyFunc(OO, @(1).val.dims())),
    Scat_OO_Gath := Rule([@(1, Compose), Scat, OO, Gath],
        e -> ApplyFunc(OO, @(1).val.dims())),
));

SIMT_NextLoop := (s, simtidx) ->
    SubstBottomUp(s, [@(1, SIMTISum), @(2, Compose, e->ForAll(e.children(), c->ObjId(c) in [ISum, SUM]))],
        e->let(sx1c := @(2).val.children(),
            doms := List(sx1c, c->When(ObjId(c) = ISum, c.domain, Length(c.children()))),
            mdom := Maximum(doms),
            ranges := List(List(sx1c, c->[0, When(ObjId(c) = ISum, c.domain, Length(c.children()))-1])),
            newc := List([1..Length(sx1c)],
                i-> When(ObjId(sx1c[i]) = ISum,
                        SIMTISum(simtidx(mdom, ranges[i]), sx1c[i].var, sx1c[i].domain, sx1c[i].child(1)),
                        SIMTSUM(simtidx(mdom, ranges[i]), sx1c[i].children()))),
            SIMTISum(@(1).val.simt_dim, @(1).val.var, @(1).val.domain, Compose(newc))
        )
    );

_fBaseVar := function(s, i, v)
    local cands;

    cands := Collect(s, @(0, fBase, e->e.params[1]=i.range and IsInt(e.params[2]) and e.params[2] = v));
    if cands = [] then return s; fi;

    s := SubstTopDown(s, @(0, fBase, e->e.params[1]=i.range and IsInt(e.params[2]) and e.params[2] = v),
        e->fBase(i));
    return s;
end;

FixUpCUDASigmaSPL := function(s, opts)
    local srt, kernels, i, _s, newv;

    if IsBound(s.ruletree) then srt := s.ruletree; fi;

    s := SIMT_NextLoop(s, ASIMTBlockDimY);
    s := SIMT_NextLoop(s, ASIMTBlockDimX);

    s:= SubstTopDown(s, @(1, SIMTSUM),
        e->let(
            its := @(1).val.simt_dim.params[1],
            ii := Ind(its),
            ch := @(1).val.children(),
            nch := Length(ch),
            SIMTISum(@(1).val.simt_dim, ii, nch, SUM(List([0..nch-1],
               i->COND(eq(ii, V(i)), _fBaseVar(ch[i+1], ii, i), ApplyFunc(OO, ch[i+1].dims()))))))
    );

    if Length(Collect(s,  @(1, SIMTISum, e->ObjId(e.simt_dim) = ASIMTBlockDimX))) > 0 then
        s := SubstTopDown(s, [@@(1, SIMTISum, (e,cx)->ObjId(e.simt_dim) = ASIMTBlockDimY), @(2, [BB, SUM])],
            (e,cx)->let(xdimcds := Collect(Filtered(cx.SIMTISum, s->ObjId(s.simt_dim) = ASIMTKernelFlag), @(0, SIMTISum, (e)->ObjId(e.simt_dim) = ASIMTBlockDimX)),
                xdim := When(xdimcds <> [], xdimcds[1].domain, @@(1).val.domain),
                SIMTISum(@@(1).val.simt_dim, @@(1).val.var, @@(1).val.domain,
                    SIMTISum(ASIMTBlockDimX(xdim,[0,0]), Ind(1), 1, @@(1).val.child(1))))
        );

        s := ApplyStrategy(s, [RulesSIMTFission, RulesSums], BUA, opts);
#        s := RulesSums(RulesSIMTFission(s));

        s := SubstTopDown(s, @@(1, SIMTSUM, (e,cx)->ObjId(e.simt_dim) = ASIMTBlockDimY and ForAll(e.children(), c->not ObjId(c) in [SIMTISum, SIMTSUM])),
            (e,cx)->let(xdimcds := Collect(Filtered(cx.SIMTISum, s->ObjId(s.simt_dim) = ASIMTBlockDimZ), @(0, SIMTISum, (e)->ObjId(e.simt_dim) = ASIMTBlockDimX)),
                xdim := When(xdimcds <> [], xdimcds[1].domain, When(IsBound(@@(1).val.domain), @@(1).val.domain, 1)),
                SIMTSUM(@@(1).val.simt_dim,
                    List(@@(1).val.children(), c->SIMTISum(ASIMTBlockDimX(xdim, [0,0]), Ind(1), 1, c))))
        );

        s := SubstTopDown(s, [@(1, SIMTISum, e->ObjId(e.simt_dim) = ASIMTBlockDimX and Set(Flat(e.simt_dim.simt_all_ranges())) = [0]), @(2, ISum)],
            e -> SIMTISum(ASIMTBlockDimX(@(2).val.domain), @(2).val.var, @(2).val.domain, @(2).val.child(1))
        );

        if ObjId(s) = Compose then
            kernels := s.children();

            for i in [1..Length(kernels)] do
                _s := kernels[i];
                if Length(Collect(_s,  @(1, SIMTISum, e->ObjId(e.simt_dim) = ASIMTBlockDimX))) > 0 then
                    newv := Maximum(List(List(Collect(_s, @(1, [SIMTSUM, SIMTISum], e->ObjId(e.simt_dim) = ASIMTBlockDimX)), e-> e.simt_dim), i->i.params[1]));
                    _s := SubstTopDown(_s, @(1, SIMTISum, e->ObjId(e.simt_dim) = ASIMTBlockDimX),
                        e->SIMTISum(ApplyFunc(ASIMTBlockDimX, [newv]::Drop(@(1).val.simt_dim.params, 1)), @(1).val.var, @(1).val.domain, @(1).val.child(1))
                    );

                    _s := SubstTopDown(_s, @(1, SIMTSUM, e->ObjId(e.simt_dim) = ASIMTBlockDimX),
                        e->SIMTSUM(ApplyFunc(ASIMTBlockDimX, [newv]::Drop(@(1).val.simt_dim.params, 1)), @(1).val.children())
                    );
                fi;
                kernels[i] := _s;
            od;
            s := Compose(kernels);
        fi;

        s := SubstTopDown(s, [@(1, SIMTSUM, e->ObjId(e.simt_dim) = ASIMTBlockDimX), @(2,BB), @(3, ISum)],
            e->SIMTSUM(ApplyFunc(ObjId(@(1).val.simt_dim),  [@(1).val.simt_dim.params[1], [0, @(3).val.domain-1], [@(3).val.domain,@(3).val.domain]]),
                SIMTISum(ApplyFunc(ObjId(@(1).val.simt_dim),  [@(1).val.simt_dim.params[1], [0, @(3).val.domain-1]]), @(3).val.var, @(3).val.domain, @(3).val.child(1)),
                @(2).val)
        );
    fi;
    # somehow SIMTSUM by itself does not create the necessary loop dimension (?)
    s := SubstBottomUp(s, [@(1,SIMTISum, e->ObjId(e.simt_dim) = ASIMTBlockDimZ), [BB, @(2,SUM)], ...],
        e->let(ch := @(2).val.children(), nch := Length(ch),
            ii := Ind(nch),
            SIMTISum(@(1).val.simt_dim, @(1).val.var, @(1).val.domain,
                SIMTISum(ASIMTBlockDimY(nch), ii, nch, SUM(List([0..nch-1], i->COND(eq(ii, V(i)), BB(_fBaseVar(ch[i+1], ii, i)), ApplyFunc(OO, ch[i+1].dims())))))))
    );

    s := SubstTopDown(s, @(1, Grp), e->e.child(1));

    s := SubstBottomUp(s, @(1, SIMTISum, e->ObjId(e.simt_dim) = ASIMTBlockDimZ),
        e->let(SIMTISum(ASIMTGridDimZ(@(1).val.simt_dim.params[1]), @(1).val.var, @(1).val.domain, @(1).val.child(1)))
    );

    if IsBound(opts.fixUpZ) and opts.fixUpZ then
        s := SubstBottomUp(s, @(1, SIMTISum, e->ObjId(e.simt_dim) = ASIMTGridDimZ),
            e -> let(
                orig_zdim := @(1).val.domain,
                xdim := Maximum([1]::List(Collect(@(1).val.child(1), @(2, [SIMTSUM, SIMTISum], e->ObjId(e.simt_dim) = ASIMTBlockDimX)), i->i.simt_dim.params[1])),
                ydim := Maximum([1]::List(Collect(@(1).val.child(1), @(2, [SIMTSUM, SIMTISum], e->ObjId(e.simt_dim) = ASIMTBlockDimY)), i->i.simt_dim.params[1])),
                zdim_inner := Maximum([1]::Filtered(DivisorsInt(orig_zdim), i -> i * xdim * ydim <= (opts.max_threads))),
                zdim_outer := orig_zdim / zdim_inner,
                io := Ind(zdim_outer),
                sdo := ASIMTGridDimZ(zdim_outer),
                ii := Ind(zdim_inner),
                sdi := ASIMTBlockDimZ(zdim_inner),
                krn := SubstVars(Copy(@(1).val.child(1)), rec((@(1).val.var.id) := (io * zdim_inner + ii))),
                SIMTISum(sdo, io, zdim_outer,
                    SIMTISum(sdi, ii, zdim_inner, krn))
            )
        );
    fi;

    if IsBound(s.ruletree) then s.ruletree := srt; fi;

    return s;
end;

FixUpCUDASigmaSPL_3Stage := function(ss, opts)
    local kernels, _s, newv;

    # drop grp
    ss := SubstTopDown(ss, @(1, Grp), e->e.child(1));

    # parallelize and flatten loop
    ss:= let(simtidx := ASIMTBlockDimX,
        SubstBottomUp(ss, [@(1, SIMTISum), @(2, Compose, e->ForAll(e.children(), c->ObjId(c) = ISum))],
            e->let(sx1c := @(2).val.children(),
                    doms := List(sx1c, c->c.domain),
                    mdom := Maximum(doms),
                    ranges := List(List(sx1c, c->[0, c.domain-1])),
                    newc := List([1..Length(sx1c)],
                        i-> SIMTISum(@(1).val.simt_dim, @(1).val.var, @(1).val.domain, SIMTISum(simtidx(mdom, ranges[i]), sx1c[i].var, sx1c[i].domain, sx1c[i].child(1)))),
                ApplyFunc(Compose, newc)
            ))
    );

    #ll := Collect(ss, [@(1, SIMTISum, e->ObjId(e.simt_dim) = ASIMTBlockDimX), [@(2, SIMTISum, e->ObjId(e.simt_dim) = ASIMTBlockDimX), BB]]);
    #
    #s1 := ll[1];
    #i1 := s1.var;
    #i2 := s1.child(1).var;
    #
    #rng := i1.range * i2.range;
    #ii := Ind(rng);
    #sr := rec(
    #    (i1.id) := idiv(ii, i2.range),
    #    (i2.id) := imod(ii, i2.range)
    #);
    #
    #sdim := ASIMTBlockDimX(rng);
    #
    #ss2 := SIMTISum(sdim, ii, ii.range, SubstVars(s1.child(1).child(1), sr));

    # normalize loop
    ss := SubstTopDown(ss,
        [@(1, SIMTISum, e->ObjId(e.simt_dim) = ASIMTBlockDimX), [@(2, SIMTISum, e->ObjId(e.simt_dim) = ASIMTBlockDimX), BB]],
        e->let(s1 := @(1).val,
            i1 := s1.var,
            i2 := s1.child(1).var,
            rng := i1.range * i2.range,
            ii := Ind(rng),
            sr := rec(
                (i1.id) := idiv(ii, i2.range),
                (i2.id) := imod(ii, i2.range)
            ),
            sdim := ASIMTBlockDimX(rng),
            SIMTISum(sdim, ii, ii.range, SubstVars(s1.child(1).child(1), sr))
        )
    );

    # fix loop iterations
    if ObjId(ss) = Compose then
        kernels := ss.children();

        for i in [1..Length(kernels)] do
            _s := kernels[i];
            if Length(Collect(_s,  @(1, SIMTISum, e->ObjId(e.simt_dim) = ASIMTBlockDimX))) > 0 then
                newv := Maximum(List(List(Collect(_s, @(1, [SIMTSUM, SIMTISum], e->ObjId(e.simt_dim) = ASIMTBlockDimX)), e-> e.simt_dim), i->i.params[1]));
                _s := SubstTopDown(_s, @(1, SIMTISum, e->ObjId(e.simt_dim) = ASIMTBlockDimX),
                    e->SIMTISum(ApplyFunc(ASIMTBlockDimX, [newv]::Drop(@(1).val.simt_dim.params, 1)), @(1).val.var, @(1).val.domain, @(1).val.child(1))
                );

                _s := SubstTopDown(_s, @(1, SIMTSUM, e->ObjId(e.simt_dim) = ASIMTBlockDimX),
                    e->SIMTSUM(ApplyFunc(ASIMTBlockDimX, [newv]::Drop(@(1).val.simt_dim.params, 1)), @(1).val.children())
                );
            fi;
            kernels[i] := _s;
        od;
        ss := Compose(kernels);
    fi;

    return ss;
end;

max_stars_repo_name: broderickpt/spiral-package-fftx

max_issues_repo_name: broderickpt/spiral-package-fftx

max_forks_repo_name: broderickpt/spiral-package-fftx


text:
~% gap

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

     Information at:  http://www.gap-system.org
     Try '?help' for help. See also  '?copyright' and  '?authors'

   Loading the library. Please be patient, this may take a while.
GAP4, Version: 4.4.12 of 17-Dec-2008, x86_64-unknown-linux-gnu-gcc
Components:  small 2.1, small2 2.0, small3 2.0, small4 1.0, small5 1.0, small6 1.0, small7 1.0, small8 1.0,
             small9 1.0, small10 0.2, id2 3.0, id3 2.1, id4 1.0, id5 1.0, id6 1.0, id9 1.0, id10 0.1, trans 1.0,
             prim 2.1  loaded.
Packages:    AClib 1.1, Polycyclic 2.6, Alnuth 2.2.5, AutPGrp 1.4, CrystCat 1.1.3, Cryst 4.1.6, CRISP 1.3.2,
             CTblLib 1.1.3, TomLib 1.1.4, FactInt 1.5.2, GAPDoc 1.2, FGA 1.1.0.1, IRREDSOL 1.1.2, LAGUNA 3.5.0,
             Sophus 1.23, Polenta 1.2.7, ResClasses 2.5.3  loaded.
gap> join := function(a, b, sep)
>   return Concatenation(a, sep, sep, b);
> end;
function( a, b, sep ) ... end
gap>
gap> join("Rosetta", "Code", ":");
"Rosetta::Code"
gap>