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 := (ss, opts) -> SubstTopDown(ss, @(1, Grp), e->e.child(1));
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( #Error(),
                    simtloop := @(1).val,
                    sx1c := @(2).val.children(),
                    doms := List(sx1c, c->c.domain),
                    i := @(1).val.var,
                    mdom := Maximum(doms),
                    ranges := List(List(sx1c, c->[0, c.domain-1])),
                    nch := [sx1c[1] * Gath(fTensor(fBase(i), fId(Cols(sx1c[1]))))] ::
                        List(sx1c{[2..Length(sx1c)-1]}, c -> Scat(fTensor(fBase(i), fId(Rows(c)))) * c * Gath(fTensor(fBase(i), fId(Cols(c))))) ::
                        [ Scat(fTensor(fBase(i), fId(Rows(Last(sx1c))))) * Last(sx1c)],
                    newc := List([1..Length(nch)],
                        j-> SIMTISum(simtloop.simt_dim, simtloop.var, simtloop.domain,
                                SIMTISum(simtidx(mdom, ranges[j]), sx1c[j].var, sx1c[j].domain, ApplyStrategy(nch[j], opts.formulaStrategies.sigmaSpl, BUA, opts).child(1)))),
                ApplyFunc(Compose, newc)
            ))
    );
    ss := ApplyStrategy(ss, opts.formulaStrategies.sigmaSpl, BUA, opts);

#    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),
#                    i := @(1).val.var,
#                    mdom := Maximum(doms),
#                    ranges := List(List(sx1c, c->[0, c.domain-1])),
#                    newc := List([1..Length(sx1c)],
#                        j-> SIMTISum(@(1).val.simt_dim, @(1).val.var, @(1).val.domain, SIMTISum(simtidx(mdom, ranges[j]), sx1c[j].var, sx1c[j].domain, sx1c[j].child(1)))),
#                ApplyFunc(Compose, newc)
#            ))
#    );

    # loop distribution Y(X*X)
    ss := SubstBottomUp(ss, [@(1, SIMTISum, e->ObjId(e.simt_dim) = ASIMTBlockDimY), Compose],
        e -> let(ch := @(1).val.child(1).children(), i := @(1).val.var,
            nch := [ch[1] * Gath(fTensor(fBase(i), fId(Cols(ch[1]))))] ::
                List(ch{[2..Length(ch)-1]}, c -> Scat(fTensor(fBase(i), fId(Rows(c)))) * c * Gath(fTensor(fBase(i), fId(Cols(c))))) ::
                [ Scat(fTensor(fBase(i), fId(Rows(Last(ch))))) * Last(ch)],
            Compose(List(nch, c -> SIMTISum(@(1).val.simt_dim, @(1).val.var, @(1).val.var.range, c)))
            ));
    ss := ApplyStrategy(ss, opts.formulaStrategies.sigmaSpl, BUA, opts);

    # privatize temporaries -- that should be a no-op now
    ss := SubstTopDown(ss, [@(1, SIMTISum), @(2, Compose, e->ForAll(e.children(), c->ObjId(c) = ISum))],
        e-> let(ch := @(1).val.child(1).children(),
            i := @(1).val.var,
            nch := [ch[1] * Gath(fTensor(fBase(i), fId(Cols(ch[1]))))] ::
                List(ch{[2..Length(ch)-1]}, c -> Scat(fTensor(fBase(i), fId(Rows(c)))) * c * Gath(fTensor(fBase(i), fId(Cols(c))))) ::
                [ Scat(fTensor(fBase(i), fId(Rows(Last(ch))))) * Last(ch)],
            SIMTISum(@(1).val.simt_dim, @(1).val.var, @(1).val.domain, Compose(List(nch, c->ApplyStrategy(c, opts.formulaStrategies.sigmaSpl, BUA, opts))))
        )
    );

    # flatten X/X -> X loops
    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))
        )
    );

    # flatten Y/X -> X loops
    ss := SubstTopDown(ss,
        [@(1, SIMTISum, e->ObjId(e.simt_dim) = ASIMTBlockDimY), [@(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();
    else
        kernels := [ss];
    fi;

    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;

    if ObjId(ss) = Compose then
        ss := Compose(kernels);
    else
        ss := kernels[1];
    fi;

#    # prepare for PingPong: privatize first ping so it can be mapped to Y without conflict
#    ss := SubstBottomUp(ss, [@(1, SIMTISum), [@(2, Compose), @(3, SIMTISum), @(4, SIMTISum), @(5, SIMTISum)]],
#        e->let(
#            g := Gath(fTensor(fBase(@(1).val.var), fId(Cols(@(4).val)))),
#            s := Scat(fTensor(fBase(@(1).val.var), fId(Rows(@(5).val)))),
#            krn := @(3).val * Grp(@(4).val * g) * Grp(s * @(5).val),
#            CopyFields(@(1).val, rec(_children := [krn])))
#    );
#    ss := ApplyStrategy(ss, opts.formulaStrategies.sigmaSpl, BUA, opts);
#    ss := SubstTopDown(ss, @(1, Grp), e->e.child(1));

    return ss;
end;

PingPong_3Stages := function(c, opts)
    local cands, candvars, outvars, invars, in_loops, substvars, link_loops, linkvars, link_func, subst_list, substrec, shift, shiftvar, v;

    # need to guard that better...
    cands := Collect(c, @(1, specifiers_func, e -> (Collect(e, @(2, decl, f -> Length(f.vars) <= 2 and ForAll(f.vars, k -> ObjId(k.t) = TArray))) <> [])) );

    candvars := chain(cands).free()::[X, Y];

    if Length(cands) > 0 then
#    Error();
        outvars := Filtered([Y]::Flat(List(Collect(c, @@(1, decl,
            (e, cx) -> (not IsBound(cx.specifiers_func) or cx.specifiers_func = []) and
                       (not IsBound(cx.func) or cx.func = []))), x->x.vars)), p -> p in candvars);
#        outvars := List(Collect(c, @(1, cu_call)), e->e.args[2]);

        invars := Filtered([X]::Flat(List(Collect(c, @@(1, decl,
            (e, cx) -> (not IsBound(cx.specifiers_func) or cx.specifiers_func = []) and
                       (not IsBound(cx.func) or cx.func = []))), x->x.vars)), p -> p in candvars);
#        invars := List(Collect(c, @(1, cu_call)), e->e.args[1]);

        in_loops := Collect(cands, @(1, simt_loop, e -> Collect(e.cmd, @(2,simt_loop)) = [] and
            Collect(e.cmd, [assign, @(3), @(4, nth, f->f.loc in invars), ...]) <> []));

        substvars := Set(List(Collect(in_loops, [assign, @(1, nth, e -> ObjId(e.loc.t) = TArray), ...]), f -> f.loc.loc));

        link_loops := Collect(cands, @(1, simt_loop, e -> Collect(e.cmd, @(2,simt_loop)) = [] and
            Collect(e.cmd, [assign, @(3), @(4, nth, f->f.loc in substvars), ...]) <> []));

        linkvars := Set(List(Collect(link_loops, [assign, @(1, nth, e -> ObjId(e.loc.t) = TArray), ...]), f -> f.loc.loc));

        link_func := v -> Set(Collect(Filtered(Collect(c, @(1, specifiers_func)),
            e -> Collect(e.cmd, [assign, @(3), @(4, nth, f->f.loc = v), ...]) <> []), @(1, var, e-> e in outvars)))[1];

        for v in substvars do
            shiftvar := Collect(Collect(c, @(1, decl, e-> v in e.vars)), @(2, simt_loop, f->ObjId(f.simt_idx) = simtBlockIdxX))[1].var;
            shift := v.t.size;
            SubstBottomUp(c, @(1, nth, e->e.loc = v),
                e-> nth(@(1).val.loc, shiftvar * shift + @(1).val.idx));
        od;

        subst_list := List(substvars, vv -> [vv, link_func(vv)]);

        substrec := FoldR(subst_list, (a,b) -> CopyFields(a, rec((b[1].id) := b[2])), rec());

        c := SubstBottomUp(c, @(1, decl), e-> decl(Filtered(@(1).val.vars, v -> not v in substvars), @(1).val.cmd));
        c := SubstVars(c, substrec);
    fi;

    return c;
end;

#FixUpCUDASigmaSPL_3Stage := (ss, opts) -> _FixUpCUDASigmaSPL_3Stage(_FixUpCUDASigmaSPL_3Stage(ss, opts), opts);

FixUpTeslaV_Code := function (c, opts)
    local kernels, kernel_inits, globals, var_decls, var_dels, cx, v, heap_init, dptr;

    if IsBound(opts.fixUpTeslaV_Code) and opts.fixUpTeslaV_Code then
        kernels := List(Collect(c, specifiers_func), k->k.id);

        dptr := var.fresh_t("hp", TPtr(TReal));

#        kernel_inits := List(kernels, k-> call(rec(id := "INIT_KERNEL"), k));
#        #define INIT_KERNEL(k) cudaFuncSetCacheConfig(k, cudaFuncCachePreferEqual)
        kernel_inits := List(kernels, k-> call(rec(id := "cudaFuncSetCacheConfig"), k, "cudaFuncCachePreferL1"));

#       cudaDeviceSetLimit(cudaLimitMallocHeapSize, (DEV_MIN_HEAP_SIZE));
        heap_init := call(rec(id := "cudaDeviceSetLimit"), "cudaLimitMallocHeapSize", opts.max_heap);

        globals := Flat(List(Collect(c, @@(1, decl, (e, cx) -> (not IsBound(cx.specifiers_func) or cx.specifiers_func = []) and
                               (not IsBound(cx.func) or cx.func = []))), x->x.vars));

#        var_decls := List(globals, v -> call(rec(id := "DECLARE_DEVICE_ARRAY"), v.id, sizeof(v.t.t) * v.t.size));
#        double *hp;\
#        cudaMalloc((void**)&hp,(sz));\
#        cudaMemcpyToSymbol(p,&hp,sizeof(double*));\

        var_decls := chain(Flat(List(globals, v -> [
                call(rec(id := "cudaMalloc"), tcast(TPtr(TPtr(TVoid)), addrof(dptr)), sizeof(v.t.t) * v.t.size),
                call(rec(id := "cudaMemcpyToSymbol"), v, addrof(dptr), sizeof(dptr.t))
            ])));

#        var_dels := List(globals, v -> call(rec(id := "DELETE_DEVICE_ARRAY"), v.id));
#        double *hp;\
#        cudaMemcpyFromSymbol(&hp,p,sizeof(double*));\
#        cudaFree((void**)&hp);\

        var_dels := decl(dptr, chain(Flat(List(globals, v -> [
                call(rec(id := "cudaMemcpyFromSymbol"), addrof(dptr), v, sizeof(dptr.t)),
                call(rec(id := "cudaFree"), dptr)
            ]))));

        cx := chain([heap_init] :: kernel_inits :: [var_decls]);
        for v in globals do
            v.t := TPtr(v.t.t, ["__device__"]);
        od;

        c := SubstBottomUp(c, @(1, func, f -> f.id = "init"),
            e -> CopyFields(@(1).val, rec(cmd := decl(dptr, chain(cx, @(1).val.cmd))))
        );
        c := SubstBottomUp(c, @(1, func, f -> f.id = "destroy"),
            e -> CopyFields(@(1).val, rec(cmd := chain(var_dels, @(1).val.cmd)))
        );
        c := SubstBottomUp(c, @(1, chain, e -> ForAny(e.cmds, e -> ObjId(e) = func and e.id = "init")),
            e -> chain(Filtered(@(1).val.cmds, e-> ObjId(e) <> func or e.id <> "init") :: Filtered(@(1).val.cmds, e -> ObjId(e) = func and e.id = "init"))
        );

    fi;

    return c;
end;

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

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

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

#    ss := SIMT_NextLoop(ss, ASIMTBlockDimY);
    ss := SIMT_NextLoop(ss, ASIMTBlockDimX);

    ss:= SubstTopDown(ss, @(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()))))))
    );

    # loop distribution X(X*X)
    ss := SubstBottomUp(ss, [@(1, SIMTISum, e->ObjId(e.simt_dim) = ASIMTBlockDimX), Compose],
        e -> let(ch := @(1).val.child(1).children(), i := @(1).val.var,
            nch := [ch[1] * Gath(fTensor(fBase(i), fId(Cols(ch[1]))))] ::
                List(ch{[2..Length(ch)-1]}, c -> Scat(fTensor(fBase(i), fId(Rows(c)))) * c * Gath(fTensor(fBase(i), fId(Cols(c))))) ::
                [ Scat(fTensor(fBase(i), fId(Rows(Last(ch))))) * Last(ch)],
            Compose(List(nch, c -> SIMTISum(@(1).val.simt_dim, @(1).val.var, @(1).val.var.range, c)))
            ));
    ss := ApplyStrategy(ss, opts.formulaStrategies.sigmaSpl, BUA, opts);

    # flatten X/X -> X loops
    ss := SubstTopDown(ss,
        [@(1, SIMTISum, e->ObjId(e.simt_dim) = ASIMTBlockDimX), [@(2, SIMTISum, e->ObjId(e.simt_dim) = ASIMTBlockDimX), @(3, [BB, SUM]) ]],
        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))
        )
    );

    # loop distribution Y(X*X)
    ss := SubstBottomUp(ss, [@(1, SIMTISum, e->ObjId(e.simt_dim) = ASIMTBlockDimY), Compose],
        e -> let(ch := @(1).val.child(1).children(), i := @(1).val.var,
            nch := [ch[1] * Gath(fTensor(fBase(i), fId(Cols(ch[1]))))] ::
                List(ch{[2..Length(ch)-1]}, c -> Scat(fTensor(fBase(i), fId(Rows(c)))) * c * Gath(fTensor(fBase(i), fId(Cols(c))))) ::
                [ Scat(fTensor(fBase(i), fId(Rows(Last(ch))))) * Last(ch)],
            Compose(List(nch, c -> SIMTISum(@(1).val.simt_dim, @(1).val.var, @(1).val.var.range, c)))
            ));
    ss := ApplyStrategy(ss, opts.formulaStrategies.sigmaSpl, BUA, opts);

    # flatten Y/X -> X loops
    ss := SubstTopDown(ss,
        [@(1, SIMTISum, e->ObjId(e.simt_dim) = ASIMTBlockDimY), [@(2, SIMTISum, e->ObjId(e.simt_dim) = ASIMTBlockDimX), @(3, [BB, SUM, Data, ISum])]],
        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))
        )
    );

    # flatten GridDim X/X -> X loops
    ss := SubstTopDown(ss,
        [@(1, SIMTISum, e->ObjId(e.simt_dim) = ASIMTBlockDimX), [@(2, SIMTISum, e->ObjId(e.simt_dim) = ASIMTBlockDimX), @(3, [BB, SUM]) ]],
        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();
    else
        kernels := [ss];
    fi;

    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;

    if ObjId(ss) = Compose then
        ss := Compose(kernels);
    else
        ss := kernels[1];
    fi;

    return ss;
end;

max_stars_repo_name: mikefranusich/spiral-package-fftx

max_issues_repo_name: mikefranusich/spiral-package-fftx

max_forks_repo_name: mikefranusich/spiral-package-fftx


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

_B_dct1 := n -> let(j:=Ind(n), VStack(
    RowVec(diagDirsum(FList(TReal, [2,1]), fConst(TReal, n-2, 0.0))),
    HStack(ColVec(fConst(TReal, n-3, 0.0)),
           RowTensor(n-3, 1, Mat([[1,1]])),
           ColVec(fConst(TReal, n-3, 0.0))),
    RowVec(diagDirsum(fConst(TReal, n-2, 0.0), FList(TReal, [1, 2]))),
    RowVec(Lambda(j, cond(neq(imod(j, 2),0), -1, 1)))
));

RulesFor(DCT1, rec(
    #F DCT1_Base2: (base case for size 2)
    #F
    #F   DCT1_2 = F_2        (scaled)
    #F
    DCT1_Base2 := rec(
	info             := "DCT_1 -> F_2",
	forTransposition := false,
	isApplicable     := P -> P[1]=2,
	allChildren      := P -> [[]],
	rule             := (P, C) -> F(2)
    ),

    #F DCT1_Base4: DCT1_4 = (M tensor F_2) ^ P
    #F
    DCT1_Base4 := rec(
	info             := "DCT_1 -> (M tensor F_2)^P",
	forTransposition := false,
	isApplicable     := P -> P[1]=4,
	allChildren      := P -> [[]],
	rule := (P, C) ->
	    Tensor(Mat([[1, 1], [1, -1/2]]), F(2)) ^ Perm((2,4), 4)
    ),

    #F DCT1_DCT1and3: 1984
    #F
    #F   DCT1_(n+1) = perm * (DCT1_(n/2+1) dirsum DCT3_n/2 ^ perm) *
    #F             (1 tensor DFT_2)^perm
    #F
    #F   Wang:
    #F     Fast Algorithms for the Discrete W Transform and the
    #F     Discrete Fourier Transform.
    #F     IEEE Trans. on ASSP, 1984, pp. 803--814
    #F
    DCT1_DCT1and3 := rec (
	info             := "DCT1_n --> DCT1_n/2, DCT3_n/2",
	isApplicable     := P -> (P[1] - 1) mod 2 = 0,
	allChildren := P ->
	    When(P[1]=3, [[ DCT1(2) ]],
		      [[ DCT1((P[1]+1)/2), DCT3((P[1]-1)/2) ]]),

	rule := (P, C) ->
	    When(P[1]=3,
		perm2(3) *
		DirectSum(C[1], I(1)) *
		DirectSum(I(1), F(2)) ^ perm4(3),

		perm2(P[1]) *
		DirectSum(C[1], C[2] ^ J((P[1]-1)/2)) *
		DirectSum(I(1), Tensor(I((P[1]-1)/2), F(2))) ^ perm4(P[1]))
    ),

    DCT1_toDCT2 := rec(
        isApplicable := P -> P[1] > 5,
        allChildren := P -> [[DCT2(P[1]-1)]],
        rule := (P, C) ->  let(
            n := P[1], j := Ind(n-1), nn := n-1,
            DirectSum(Diag(Lambda(j, 1/(2*cospi(fdiv(j, 2*nn))))) *
                      C[1], I(1)) *
            _B_dct1(n)
        )
    )
));

ext: gi


text:
%Terminals
    DollarSign ::= '$'
    Percent ::= '%'
    _
    a b c d e f g h i j k l m n o p q r s t u v w x y z
    A B C D E F G H I J K L M N O P Q R S T U V W X Y Z
%End

%Headers
    /.
        inline static int tokenKind[128] ={};

        static  bool init()
        {
            tokenKind['$'] = $sym_type::$prefix$DollarSign$suffix$;
            tokenKind['%'] = $sym_type::$prefix$Percent$suffix$;
            tokenKind['_'] = $sym_type::$prefix$_$suffix$;

            tokenKind['a'] = $sym_type::$prefix$a$suffix$;
            tokenKind['b'] = $sym_type::$prefix$b$suffix$;
            tokenKind['c'] = $sym_type::$prefix$c$suffix$;
            tokenKind['d'] = $sym_type::$prefix$d$suffix$;
            tokenKind['e'] = $sym_type::$prefix$e$suffix$;
            tokenKind['f'] = $sym_type::$prefix$f$suffix$;
            tokenKind['g'] = $sym_type::$prefix$g$suffix$;
            tokenKind['h'] = $sym_type::$prefix$h$suffix$;
            tokenKind['i'] = $sym_type::$prefix$i$suffix$;
            tokenKind['j'] = $sym_type::$prefix$j$suffix$;
            tokenKind['k'] = $sym_type::$prefix$k$suffix$;
            tokenKind['l'] = $sym_type::$prefix$l$suffix$;
            tokenKind['m'] = $sym_type::$prefix$m$suffix$;
            tokenKind['n'] = $sym_type::$prefix$n$suffix$;
            tokenKind['o'] = $sym_type::$prefix$o$suffix$;
            tokenKind['p'] = $sym_type::$prefix$p$suffix$;
            tokenKind['q'] = $sym_type::$prefix$q$suffix$;
            tokenKind['r'] = $sym_type::$prefix$r$suffix$;
            tokenKind['s'] = $sym_type::$prefix$s$suffix$;
            tokenKind['t'] = $sym_type::$prefix$t$suffix$;
            tokenKind['u'] = $sym_type::$prefix$u$suffix$;
            tokenKind['v'] = $sym_type::$prefix$v$suffix$;
            tokenKind['w'] = $sym_type::$prefix$w$suffix$;
            tokenKind['x'] = $sym_type::$prefix$x$suffix$;
            tokenKind['y'] = $sym_type::$prefix$y$suffix$;
            tokenKind['z'] = $sym_type::$prefix$z$suffix$;

            tokenKind['A'] = $sym_type::$prefix$A$suffix$;
            tokenKind['B'] = $sym_type::$prefix$B$suffix$;
            tokenKind['C'] = $sym_type::$prefix$C$suffix$;
            tokenKind['D'] = $sym_type::$prefix$D$suffix$;
            tokenKind['E'] = $sym_type::$prefix$E$suffix$;
            tokenKind['F'] = $sym_type::$prefix$F$suffix$;
            tokenKind['G'] = $sym_type::$prefix$G$suffix$;
            tokenKind['H'] = $sym_type::$prefix$H$suffix$;
            tokenKind['I'] = $sym_type::$prefix$I$suffix$;
            tokenKind['J'] = $sym_type::$prefix$J$suffix$;
            tokenKind['K'] = $sym_type::$prefix$K$suffix$;
            tokenKind['L'] = $sym_type::$prefix$L$suffix$;
            tokenKind['M'] = $sym_type::$prefix$M$suffix$;
            tokenKind['N'] = $sym_type::$prefix$N$suffix$;
            tokenKind['O'] = $sym_type::$prefix$O$suffix$;
            tokenKind['P'] = $sym_type::$prefix$P$suffix$;
            tokenKind['Q'] = $sym_type::$prefix$Q$suffix$;
            tokenKind['R'] = $sym_type::$prefix$R$suffix$;
            tokenKind['S'] = $sym_type::$prefix$S$suffix$;
            tokenKind['T'] = $sym_type::$prefix$T$suffix$;
            tokenKind['U'] = $sym_type::$prefix$U$suffix$;
            tokenKind['V'] = $sym_type::$prefix$V$suffix$;
            tokenKind['W'] = $sym_type::$prefix$W$suffix$;
            tokenKind['X'] = $sym_type::$prefix$X$suffix$;
            tokenKind['Y'] = $sym_type::$prefix$Y$suffix$;
            tokenKind['Z'] = $sym_type::$prefix$Z$suffix$;
            return true;
        };
     inline const static bool dd = init();
         int getKind(int c)
        {
            return (((c & 0xFFFFFF80) == 0) /* 0 <= c < 128? */ ? tokenKind[c] : 0);
        }

    ./
%End

max_stars_repo_name: kuafuwang/LPG-VScode

max_issues_repo_name: kuafuwang/LPG-VScode

max_forks_repo_name: kuafuwang/LPG-VScode


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

IsTag := x -> IsRec(x) and IsBound(x.isTag) and x.isTag;

Declare(ANoTag);

####################################################################
#F AGenericTag - base class for tags, uses RewritableObject features
#F    which provides
#F      __call__
#F      print
#F      rSetChild
#F      rChildren
#F      lessThan
#F      equal
#F
#F Constructor saves all parameters into .params.
#F  Please read Doc(RewritableObject).
#F
Class(AGenericTag, RewritableObject, rec(
    isTag := true,
    kind := self >> ObjId(self),

    # YSV: this should be called when we transpose the parent structure
    transpose := self >> self,

    isSticky  := false,
    right     := self >> When(IsBound(self.isRight), self.isRight, false),
    left      := self >> When(IsBound(self.isRight), not self.isRight, false),
    bisided   := self >> not IsBound(self.isRight),

    leftChild  := self >> When(IsBound(self.isLeftChild),  self.isLeftChild,  false),
    rightChild := self >> When(IsBound(self.isRightChild), self.isRightChild, false),

    edge       := self >> When(IsBound(self.isEdge),  self.isEdge,  false),

    leftEdge   := self >> self.edge() and self.leftChild(),
    rightEdge  := self >> self.edge() and self.rightChild(),

    # distCompose(<chlist>, <opts>)
    #   Distributes tag over Compose(<chlist>).
    #   Returns list of tags to assign to each composition member.
    #   ANoTag is thrown out later during normalization.
    distCompose := (self, chlist, opts) >>
        Cond( self.left(),
                  [self] :: Replicate(Length(chlist)-1, ANoTag),
              self.right(),
                  Replicate(Length(chlist)-1, ANoTag) :: [self],
              # else bisided
                  Replicate(Length(chlist), self) ),

));

Class(ASingletonTag, AGenericTag, rec(
    params := [],
    from_rChildren := (self, rch) >> self,
    kind := self >> self
));

#F ANoRecurse() -- tag that prevents further recursion.
#F
#F Currently only used in autolib.
#F
Class(ANoRecurse, AGenericTag);

#F ANoTag - denotes absence of tags
#F
Class(ANoTag, ASingletonTag, rec());

#F ATopLevel - used to denote that special top level processing should be used,
#F             generally param[1] denotes at which 'level' the tag should be
#F             dropped.
#F
#F             Marek uses it in his WHT/bit-perms expansions. This tag will
#F             be moved elsewhere soon.
Class(ATopLevel, AGenericTag, rec());

####################################################################
#F Mixin class for taggable objects
#F
#F This object/class will add tagging functionality to whatever object inherits it.
#F Like all mixin style classes, it is not meant to be the sole parent class.
#F
Class(TaggedObjectMixin, rec(
    # this is the actual list of tags. This array should not be accessed directy, rather
    # the methods should be used to parse it.
    tags := [],

    # returns the first tag or ANoTag
    firstTag := self >> When(Length(self.tags) >= 1, self.tags[1], ANoTag),

    # returns true if first tag is has 'tag' object id
    firstTagIs := (self, tag) >> self.firstTag().kind()=tag,

    # returns the tags
    getTags := self >> self.tags,

    # returns true if the object has tags
    hasTags := self >> self.tags <> [],

    # returns a copy of the object with 'tags' added on to any existing tags.
    withTags := (self, tags) >> Checked(IsList(tags), self.setTags(self.tags :: tags)),

    # returns a copy of the object with tags given by 'tags'
    setTags := (self, tags) >> Checked(IsList(tags), CopyFields(self, rec(tags := tags))),

    # return a copy of the whole object without the first tag
    withoutFirstTag := self >> self.setTags(Drop(self.tags, 1)),

    # return a copy of the whole object without the tags given by object id 'tag'
    withoutTag := (self, tag) >> Checked(IsTag(tag), self.setTags(Filtered(self.tags, e -> e.kind() <> tag))),

    # returns true when the object id of at least one of the tags is 'tag'
    hasTag := (self, tag) >> Checked(IsTag(tag), ForAny(self.tags, e -> e.kind() = tag)),

    hasAnyTag := (self, tags) >> Checked(IsList(tags), ForAny(self.tags, e -> e.kind() in tags)),

    # returns true when tag number 'n' has the object id given by 't'
    isTag := (self, n, t) >> When(Length(self.tags) >= n, self.tags[n].kind() = t, false),

    numTags := (self) >> Length(self.tags),

    dropTags := (self) >> CopyFields(self, rec(tags := [])),

    #
    #F getTag
    #
    # get the tag by name or number, function has 3 forms:
    #   getTag(2)           returns 2nd tag
    #   getTag(ASomeTag)    if there is only 1 ASomeTag, returns it, if >1, returns array
    #   getTag(ASomeTag, 2) returns 2nd ASomeTag
    #
    # in case of error, getTag returns false
    #
    getTag := meth(arg)
        local s, t, n;

        Constraint(Length(arg) = 2 or Length(arg) = 3);

        s := arg[1];

        # arg[2] is a number.
        if IsInt(arg[2]) then
            return When(Length(s.tags) >= arg[2],
                s.tags[arg[2]],
                false);

        # arg[2] must be a tag id.
        else
            t := Filtered(arg[1].tags, e -> e.kind() = arg[2]);

            # if tag and number are specified
            if Length(arg) = 3 then
                return When(Length(t) >= arg[3],
                    t[arg[3]],
                    false
                );

            else
                return When(t = [],
                    false,
                    When(Length(t) = 1,
                        t[1],
                        t
                    )
                );
            fi;
        fi;
    end,

    getAnyTag := (self, tags) >> First(self.tags, x->x.kind() in tags)

));

####################################################################
#F Tagged non-terminal. This should be used for all new taggable
#F nonterminals.
#F
#F This supercedes the old system which puts tags into .params field
#F
#F Tags on nonterminals contain information that is beyond the mathematical
#F details of the object. Often tags contain information about hardware
#F which is used to direct the breakdown.
#F
Class(TaggedNonTerminal, TaggedObjectMixin, NonTerminal, rec(
    _short_print := true,
    #--------- Transformation rules support --------------------------------
    from_rChildren := (self, rch) >> let(
        len := Length(rch),
        transposed := rch[len-1],
        tags := rch[len],
        t := ApplyFunc(ObjId(self), rch{[1..len-2]}),
        tt := When(transposed, t.transpose(), t),
        attrTakeA(tt.withTags(tags), self)
    ),

    rChildren := self >>
        Concatenation(self.params, [self.transposed, self.tags]),

    rSetChild := meth(self, n, newChild)
        local l;
        l := Length(self.params);
        if n <= l then
            self.params[n] := newChild;
        elif n = l+1 then
            self.transposed := newChild;
        elif n = l+2 then
            self.tags := newChild;
        else Error("<n> must be in [1..", l+2, "]");
        fi;
        # self.canonizeParams(); ??
        self.dimensions := self.dims();
    end,

    # this print method OVERRIDES the print method in NonTerminal.
    # here, we print the TAG INFO.
    print := (self, i, is) >> Print(
	Inherited(i, is),
        When(IsBound(self.tags) and IsList(self.tags) and self.tags<>[],
             Print(".withTags(", self.tags, ")"))),

    takeTA := (self, o) >> self.setTags(o.getTags()).takeAobj(o),
));


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

#F qST( n, i ) - S Gate non-terminal
#F Definition: (2^n x 2^n)-matrix  that applies an n-point S Transform to the target qubits
#F Note:       qST(1) denotes the matrix[[1, 0], [0, e^(-ipi/2)]], or [[1,0], [0, -i]]
#F             qZT(_) is symmetric.
#F
#F qST(n , i) -> an n-qubit S transform, or Sdag transform if i = -1
Class(qST, TaggedNonTerminal, rec(
  abbrevs   := [ (n, i) -> Checked(IsPosInt(n), [n, When(i <= -1, -1, 1)]) ],
  dims      := self >> let(size := 2^self.params[1], [size, size]),
  terminate := self >> When(self.params[2] = -1, Tensor(Replicate(self.params[1], qSdag())), Tensor(Replicate(self.params[1], qS()))),
  isReal    := self >> false,
  rChildren := self >> self.params,
  from_rChildren := (self, rch) >> self.__bases__[1](rch[1], rch[2]),
  groups    := self >> [Replicate(self.params[1], 1)],
  recursive_def := (self, arch) >> self.__bases__[1](self.params[1], self.params[2]),
  SmallRandom := () -> Random([2..5]),
  LargeRandom := () -> Random([6..15]),
  normalizedArithCost := self >> Error("ArithCost not implemented"),
  TType := T_Complex(64),
));

NewRulesFor(qST, rec(

    # qST_BinSplit rule
    # qST_BinSplit qST_(k) -> (qST_(k1) tensor I) (I tensor qST_(k2))
    # k1 + k2 = k
    qST_BinSplit := rec (
        forTransposition := false,
        minSize          := 2,
        applicable       := (self, nt) >> nt.params[1] > 1,
        children         := nt -> List( [1..nt.params[1] - 1],  i -> [ Tensor(qST(i, nt.params[2]), qST(nt.params[1]-i, nt.params[2])).withTags(nt.getTags()) ]  ),
        apply            := (nt, c, cnt) -> c[1],
        switch := true,
    ),

    #F qST_Base: qST(1, _) = qS() or qS().inverse() SPL object
    #F Directly represent as an implementable gate
    qST_Base := rec(
        info             := "qST_(1, _) -> qS()",
        forTransposition := false,
        applicable       := (self, nt) >> nt.params[1]=1,
        apply            := (nt, c, cnt) -> When(nt.params[2] = -1, qSdag(), qS()),
    )

));

max_stars_repo_name: spiral-software/spiral-package-quantum

max_issues_repo_name: spiral-software/spiral-package-quantum

max_forks_repo_name: spiral-software/spiral-package-quantum


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

Declare(simtBlockIdxX, simtBlockIdxY, simtBlockIdxZ);
Declare(simtThreadIdxX, simtThreadIdxY, simtThreadIdxZ);
Declare(simt_syncgrid, simt_syncblock, simt_synccluster);

Declare(ASIMTListDim);

# The basic assumption is that threads are organized in a 6-dimensional structure (a 3D grid of 3D blocks).

# A dimension tag is a list of (thread) ranges over a one or more Grid/Block dimensions.
# e.g., ASIMTGridDimX(32) refers to 32 blocks (indexed from 0 to 31) over the X grid dimension.

# A range can also be a subset of the whole set of allocated elements:
# e.g., ASIMTGridDimX(32, [4,8]) refers to 5 blocks (indexed from 4 to 8) over the X grid dimension.

# More than one ranges can also be expressed over the same dimension:
# e.g., ASIMTGridDimX(32, [0,3], [4,8])
# As well as over different ones:
# e.g., ASIMTListDim( ASIMTGridDimX(32, [4,8]), ASIMTGridDimY(96) )

# SIMT Grid/BLock Dimension tags
Class(ASIMTDim, ASingletonTag, rec(
    isSIMTTag := true,
    # When lowering to icode a dimension tag instantiates a dimension index
    _idx      := Ignore,
    _build_gen_idx := (self, Idx, plist) >> ApplyFunc(Idx, plist),
    simt_idx    := (arg) >> let(self := arg[1], idx := Copy(self._idx), Cond(idx=Ignore, idx, Length(arg)>1, ApplyFunc(idx.set_rng, self.simt_range(arg[2])), ApplyFunc(idx.set_rng, self.simt_range(1)) ) ),
    need_cond := (self, i) >> let(L := Length(self.params), i < L and L > 1 and (self.params[i+1] <> [0,self.params[1]-1]) ),
    simt_cond   := (self, i) >> let( rng := self.simt_range(i), When(
                            self.need_cond(i), Cond(
                                rng[1] = 0, leq(self._idx, V(Last(rng))),
                                Last(rng) = self.params[1]-1, geq(self._idx, V(rng[1])),
                                rng[1] = Last(rng), eq(self._idx, V(rng[1])),
                                logic_and(geq(self._idx, V(rng[1])), leq(self._idx, V(Last(rng))))
                                ),
                            V(true)
                            )),
    update_simt_range := meth(self, i, new_rng)
        local pre, post, dim;

        dim := self;
        if ObjId(new_rng) = ObjId(dim) then
            if Length(dim.params) < i then
            dim.params := dim.params :: List([Length(dim.params)..i], p -> dim.simt_range(p));
            fi;
            pre  := Take(dim.params, i);
            post := When(Length(dim.params)>i+1, Drop(dim.params, i+1), []);
            dim.params := pre::new_rng.simt_all_ranges()::post;
        elif Length(dim.params) <= 2 and i = 1 then
            dim := new_rng;
        else
            if Length(dim.params) < i then
            dim.params := dim.params :: List([Length(dim.params)..i], p -> dim.simt_range(p));
            fi;
            pre  := [ApplyFunc( ObjId(dim), Take(dim.params, i) ) ];
            post := When(Length(dim.params)>i+1, [ApplyFunc( ObjId(dim), Drop(dim.params, i+1) )], []);
            dim := ApplyFunc(ASIMTListDim, pre::[new_rng]::post);
        fi;
        return dim;
    end,
    simt_range  := (self, i) >> When(Length(self.params) > i, self.params[i+1], [0, self.params[1]-1]),
    simt_all_ranges := (self) >> When(Length(self.params) > 1, Drop(self.params,1), [ [0, self.params[1]-1] ]),
    simt_sync   := (self) >> skip(),
    num_ranges := (self) >> Maximum(1, Length(self.params)-1)
    ));

Class(ASIMTGridDim,  ASIMTDim, rec(
    simt_sync := (self) >> simt_syncgrid()
    ));

Class(ASIMTGridDimX,  ASIMTGridDim, rec(
    updateParams := meth(self)
        self._idx := self._build_gen_idx(simtBlockIdxX, self.params{[1]});
        return self;
    end,
    ));
Class(ASIMTGridDimY,  ASIMTGridDim, rec(
    updateParams := meth(self)
        self._idx := self._build_gen_idx(simtBlockIdxY, self.params{[1]});
        return self;
    end,
    ));
Class(ASIMTGridDimZ,  ASIMTGridDim, rec(
    updateParams := meth(self)
        self._idx := self._build_gen_idx(simtBlockIdxZ, self.params{[1]});
        return self;
    end,
    ));

Class(ASIMTBlockDim, ASIMTDim, rec(
    simt_sync := (self) >> simt_syncblock()
    ));

Class(ASIMTBlockDimX,  ASIMTBlockDim, rec(
    simt_sync := (self) >> simt_synccluster(),
    updateParams := meth(self)
        self._idx := self._build_gen_idx(simtThreadIdxX, self.params{[1]});
        return self;
    end,
    ));
Class(ASIMTBlockDimY,  ASIMTBlockDim, rec(
    updateParams := meth(self)
        self._idx := self._build_gen_idx(simtThreadIdxY, self.params{[1]});
        return self;
    end,
    ));
Class(ASIMTBlockDimZ,  ASIMTBlockDim, rec(
    updateParams := meth(self)
        self._idx := self._build_gen_idx(simtThreadIdxZ, self.params{[1]});
        return self;
    end,
    ));

Class(ASIMTListDim, ASIMTDim, rec(

    _point_to_rng := meth(self, i)
        local p, curr, curr_nr;

        curr := self.params[1];
        curr_nr := curr.num_ranges();
        p := 2;
        while i > curr_nr and p <= Length(self.params) do
            i := i - curr_nr;
            curr := self.params[p];
            curr_nr := curr.num_ranges();
            p := p + 1;
        od;
        #[<inner-dim containing requested rng>, <pos of rng within inner-dim>, <pos of inner-dim in params>]
        return [curr, i, p-1];
    end,
    simt_idx    := (arg) >> let(self := arg[1], nump := Length(self.params),
                                When(Length(arg)>1, let( d_i := self._point_to_rng(arg[2]), d_i[1].simt_idx(d_i[2]) ),
                                                    List(self.params, p -> p.simt_idx())
                                ) ),
    need_cond := (self, i) >> let( d_i := self._point_to_rng(i), d_i[1].need_cond(d_i[2]) ),
    simt_cond   := (self, i) >> let( d_i := self._point_to_rng(i), d_i[1].simt_cond(d_i[2]) ),
    update_simt_range := meth(self, i, new_rng)
        local d_i;

        d_i := self._point_to_rng(i);
        self.params[d_i[3]] := d_i[1].update_simt_range(d_i[2], new_rng);

        return self;
    end,
    num_ranges := (self) >> Sum(List(self.params, i -> i.num_ranges())),
    simt_range  := (self, i) >> let( d_i := self._point_to_rng(i), d_i[1].simt_range(d_i[2]) ),
    simt_all_ranges := (self) >> let( numi := self.num_ranges(), List([1..numi], i -> self.simt_range(i)) ),
    simt_sync   := (self) >> Maximum( List(self.params, p -> p.simt_sync()) ),
    simt_dim    := (self, i) >> self.params[i],
    simt_dim_at := (self, i) >> let( d_i := self._point_to_rng(i), d_i[1])
    ));

Class(ASIMTLoopDim,  ASIMTDim, rec(
    simt_range  := (self, i) >> [],
    simt_all_ranges := (self) >> [],
    num_ranges := (self) >> 0,
    update_simt_range := meth(self, i, new_rng)
        local dim;

        dim := self;
        if ObjId(new_rng) <> ObjId(dim) then
            dim := When(i=1, new_rng, ApplyFunc(ASIMTListDim, Replicate(i-1, dim)::[new_rng]) );
        fi;
        return dim;
    end,
    ));

Class(ASIMTFlag, ASIMTDim, rec(
    updateParams := meth(self)
        self._simt_dim := self.params[1];
    end,
    simt_dim    := (self) >> self._simt_dim,
    simt_idx    := (arg) >> let(self := arg[1], ApplyFunc(self._simt_dim.simt_idx, arg{[2..Length(arg)]}) ),
    need_cond := (self, i) >> self._simt_dim.need_cond(i),
    simt_cond   := (self, i) >> self._simt_dim.simt_cond(i),
    update_simt_range := meth(self, i, new_rng)
        self._simt_dim := self._simt_dim.update_simt_range(i, new_rng);
        return self;
    end,
    simt_range  := (self, i) >> self._simt_dim.simt_range(i),
    simt_all_ranges := (self) >> self._simt_dim.simt_all_ranges(),
    simt_sync   := (self) >> self._simt_dim.simt_sync(),
    num_ranges := (self) >> self._simt_dim.num_ranges()
    ));

Class(ASIMTKernelFlag, ASIMTFlag);

Class(ASIMTFissionFlag, ASIMTFlag);

Class(ASIMTTileFlag, ASIMTFlag, rec(
    updateParams := meth(self)
        Inherited();
        self._tsize  := self.params[2];
    end,

    tsize := (self) >> self._tsize
    ));

IsASIMTDim := (tag) -> ObjId(tag) in ASIMTDim._all_descendants();

##### Wrapping ######

Class(PlaceHolderMat, BaseMat, rec(

    abbrevs := [ (id, ttype, dims) -> [ rec(
                            id := id,
                            TType := Copy(ttype),
                            dimensions := Copy(dims) ) ] ],

    new := (self, spl) >> let( is_spl := IsSPL(spl), SPL(WithBases(self, rec(
                            id := When(is_spl, var._id(""), spl.id),
                            TType := When(is_spl, Copy(spl.TType), spl.TType),
                            dimensions := When(is_spl, Copy(spl.dimensions), spl.dimensions) )))
                        ),

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

    rChildren := (self) >> [ self.id, self.TType, self.dimensions ],
    rSetChild := rSetChildFields("id", "TType", "dimensions"),

    rng := meth(self) local d;
        d := [self.dimensions[1]];
        if IsBound(self.a.t_out) then
            return List(TransposedMat([self.a.t_out, d]), e -> TArray(e[1], e[2]));
        else
            return List(d, e -> TArray(self.TType, e));
        fi;
    end,

    dmn := meth(self) local d;
        d := [self.dimensions[2]];
        if IsBound(self.a.t_in) then
            return List(TransposedMat([self.a.t_in, d]), e -> TArray(e[1], e[2]));
        else
            return List(d, e -> TArray(self.TType, e));
        fi;
    end,

    ));

wrap_rule_apply := function(rules, wrapper)
    local rule;
    rules := Flat([rules]);
    for rule in rules do
        if IsNewRule(rule) then
            rule._apply := rule.apply;
            rule.apply := wrapper(rule);
        else
            rule._rule := rule.rule;
            rule.rule  := wrapper(rule);
        fi;
    od;

end;

unwrap_rule_apply := function(rules)
    local rule;
    rules := Flat([rules]);
    for rule in rules do
        if IsNewRule(rule) then
            rule.apply := rule._apply;
        else rule.rule := rule._rule;
        fi;
    od;

end;

_original_apply := (R, nt, children, nonterms) -> Checked(IsRule(R), IsSPL(nt),
    When(IsNewRule(R),
     R._apply(nt, children, nonterms),
     When(NumGenArgs(R._rule)=2,
          R._rule(nt.params, children),
          R._rule(nt.params, children, nonterms)))
);

dummy_wrap := (R) -> ( (nt, c, cnt) -> let(PrintLine("Using wrapped rule."), _apply(R, nt, c, cnt)) );

_strl := (n, c) -> StringJoin("", Replicate(n, c));

TaggableOps := [Tensor, SUM, ISum, DirectSum, IDirSum, VStack, IterVStack, HStack1, IterHStack1];

wrap_apply := (R) -> function(arg)
                        local nt, c, cnt, simt_exp, simt_dims, phs, ss;

                        nt := arg[1]; c := arg[2]; cnt := arg[3];

                        # Use placeholders for children to prevent Collect from going beyond rule.apply's boundaries
                        phs := List(c, spl -> PlaceHolderMat(spl));
                        simt_exp := _original_apply(R, nt, phs, cnt);  # This dispatches to right call: _rule (for Rules) or _apply (for NewRules)
                        ss := Collect(simt_exp, @(1, TaggableOps) );

                        simt_dims := let( cd := nt.getA("simt_dim", [ASIMTLoopDim()]),
                                        Cond(cd=Ignore, [],
                                             IsList(cd) and Length(cd) = 1, Replicate(Length(ss), cd[1]),
                                             cd)
                                    );
                        if Length(simt_dims) <> Length(ss) then
                            Error("Number of tags (", Length(simt_dims), ") <> number of taggable operators (", Length(ss), ")");
                        fi;

                        # PrintLine("@@@@@@@@@@@@@ Beg ", R, "@@@@@@@@@@@@@@@");
                        # PrintLine("NT      >> ", nt);
                        # PrintLine("simt_dims >> ", simt_dims);
                        # PrintLine("simt_expr >>", simt_exp);
                        # PrintLine("phs     >>", phs);
                        # PrintLine("ch      >>", c);
                        # PrintLine("@@@@@@@@@@@@@ End ", R, "@@@@@@@@@@@@@@@");

                        for i in [1..Length(ss)] do
                            simt_exp := SubstTopDownNR(simt_exp, @(1, ObjId(ss[i])).cond(e->e=ss[i]),
                                        e -> ApplyFunc(AllClasses.("SIMT"::ObjId(ss[i]).name), [simt_dims[i], ss[i]]));
                        od;

                        for i in [1..Length(c)] do
                            simt_exp := SubstTopDownNR(simt_exp, @(1).cond(e->e=phs[i]), e -> c[i] );
                        od;

                        return simt_exp;
                    end;

# Node of a TagTree
tag_rec := (tag, tag_list) -> rec(simt_dim := tag, ch_tags := tag_list);

# Merge RuleTree and TagTree
tmp_tag_rt := function(rt, dims, opts)
    if dims = Ignore then
        rt.node := rt.node.setA(["simt_dim", [ASIMTLoopDim()] ]);
        DoForAll([1..Length(rt.children)], i -> tmp_tag_rt(rt.children[i], dims, opts) );
    elif IsASIMTDim(dims) then
        rt.node := rt.node.setA(["simt_dim", [dims] ]);
        DoForAll([1..Length(rt.children)], i -> tmp_tag_rt(rt.children[i], Ignore, opts) );
    elif IsList(dims) then
        rt.node := rt.node.setA(["simt_dim", [ASIMTLoopDim()] ]);
        DoForAll([1..Length(rt.children)], i -> tmp_tag_rt(rt.children[i], dims[i], opts) );
    else
        rt.node := rt.node.setA(["simt_dim", dims.simt_dim]);
        DoForAll([1..Length(rt.children)], i -> tmp_tag_rt(rt.children[i], dims.ch_tags[i], opts) );
    fi;
end;

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:
for i in [10, 9 .. 0] do
    Print(i, "\n");
od;

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

Class(NEONUnparser, CMacroUnparserProg, rec(

    T_Real := (self, t, vars, i, is) >> Cond(
	t.params[1] = 32, Print("float32_t ", self.infix(vars, ", ", i+is)),
	Inherited(t, vars, i, is)),

    # composite instructions, normally they should go away before unparser
    vrev_neon  := (self, o, i, is) >> self.printf("vrev64q_f32($1)", o.args),
    vrev_half  := (self, o, i, is) >> self.printf("vrev64_f32($1)", o.args),

#    vunpacklo2_neon    := (self, o, i, is) >> self.printf("vzipq_f32(vtrnq_f32($1, $2).val[0], vtrnq_f32($1, $2).val[1]).val[0]", [o.args[1], o.args[2]]),
#    vunpackhi2_neon    := (self, o, i, is) >> self.printf("vzipq_f32(vtrnq_f32($1, $2).val[0], vtrnq_f32($1, $2).val[1]).val[1]", [o.args[1], o.args[2]]),

    vunpacklo_half    := (self, o, i, is) >> self.printf("vtrn_f32($1,$2).val[0]", o.args),
    vunpackhi_half    := (self, o, i, is) >> self.printf("vtrn_f32($1,$2).val[1]", o.args),

    vunpacklolo2_neon := (self, o, i, is) >> self.printf(
	"vcombine_f32(vget_low_f32($1),vget_low_f32($2))", o.args),
    vunpacklohi2_neon := (self, o, i, is) >> self.printf(
	"vcombine_f32(vget_low_f32($1),vget_high_f32($2))", o.args),
    vunpackhilo2_neon := (self, o, i, is) >> self.printf(
	"vcombine_f32(vget_high_f32($1),vget_low_f32($2))", o.args),
    vunpackhihi2_neon := (self, o, i, is) >> self.printf(
	"vcombine_f32(vget_high_f32($1),vget_high_f32($2))", o.args),

    # native instructions
    vuzpq_32f  := (self, o, i, is) >> self.printf("vuzpq_f32($1, $2)", o.args),
    vzipq_32f  := (self, o, i, is) >> self.printf("vzipq_f32($1, $2)", o.args),
    vtrnq_32f  := (self, o, i, is) >> self.printf("vtrnq_f32($1, $2)", o.args),
    vtrnq_half := (self, o, i, is) >> self.printf("vtrn_f32($1, $2)", o.args),

    vextract_neon_4x32f  := (self, o, i, is) >> self.printf(
	"$1.val[$2]", [o.args[1], o.args[2].p[1]]),
    vextract_half        := (self, o, i, is) >> self.printf(
	"$1.val[$2]", [o.args[1], o.args[2].p[1]]),

    vext_half  := (self, o, i, is) >> self.prefix("vext_f32", o.args),
    vext_neon  := (self, o, i, is) >> self.prefix("vextq_f32", o.args),

    fdiv := (self, o, i, is) >> self.printf("(((float)$1) / $2)", o.args),
    div  := (self, o, i, is) >> self.printf("($1 / $2)", o.args),
    idiv := (self, o, i, is) >> self.printf("($1 / $2)", o.args),
    imod := (self, o, i, is) >> self.printf("($1 % $2)", o.args),

    Value := (self, o, i, is) >> let(zero := "0", Cond(
	    (o.t = TVect(T_Real(32), 4) or o.t = TVect(TReal, 4))  and Length(Set(o.v)) = 1,
	        Print("vdupq_n_f32(", self(o.v[1], i, is), ")"),
	    (o.t = TVect(T_Real(32), 2) or o.t = TVect(TReal, 2))  and Length(Set(o.v)) = 1,
	        Print("vdup_n_f32(", self(o.v[1], i, is), ")"),
	    o.t = T_Real(32) or o.t = TReal,
	        When(o.v < 0, Print("(", o.v, ")"), Print(o.v)),
	    IsIntT(o.t) or IsUIntT(o.t),
	        When(o.v < 0, Print("(", o.v, ")"), Print(o.v)),
	    o.t = TVect(T_Real(32), 4) or o.t = TVect(TReal, 4),
	        Print("(float32x4_t){",o.v[1].v,", ",o.v[2].v,", ",o.v[3].v,", ",o.v[4].v,"}"),
	    o.t = TVect(T_Real(32), 2) or o.t = TVect(TReal, 2),
	        Print("(float32x2_t){",o.v[1].v,", ",o.v[2].v,"}"),
	    o.t = TString,
	        Print(o.v),
	    Inherited(o, i, is))),

    nth := (self, o, i, is) >> self.printf("$1[$2]", [o.loc, o.idx]),

    TReal := (self, t, vars, i, is) >> Print("float ", self.infix(vars, ", ")),

    TVect := (self, t, vars, i, is) >> Print(
        Cond( t = NEON.t,           "float32x4_t ",
              t = NEON_HALF.t,      "float32x2_t ",
              t = TVect(TReal, 4), "float32x4_t ",
              t = TVect(T_Real(32), 2), "float32x2_t ",
              t = TVect(TReal, 2), "float32x2_t ",
              t = TVect(NEON.t, 2), "float32x4x2_t ",
              t = TVect(NEON_HALF.t, 2), "float32x2x2_t ",
              Error("unexpected vector type")),
        self.infix(vars, ", ")),

    vcombine_neon   := (self, o, i, is) >> self.prefix("vcombine_f32", o.args),
    vload_half_neon := (self, o, i, is) >> self.printf("vld1_f32((float32_t*)$1)", o.args),

    vload1_neon  := (self, o, i, is) >> self.prefix("vld1q_lane_f32",
	[o.args[1], o.args[2], o.args[3].p-1]),
    vload1_half  := (self, o, i, is) >> self.prefix("vld1_lane_f32",
	[o.args[1], o.args[2], o.args[3].p-1]),
    vstore1_neon  := (self, o, i, is) >> Print(Blanks(i),
	self.prefix("vst1q_lane_f32", [o.args[1], o.args[2], o.args[3]-1]), ";\n"),
    vstore1_half  := (self, o, i, is) >> Print(Blanks(i),
	self.prefix("vst1_lane_f32", [o.args[1], o.args[2], o.args[3]-1]), ";\n"),

    vstore2lo_neon  := (self, o, i, is) >> Print(Blanks(i), self.printf(
	    "vst1_f32($1,vget_low_f32($2));\n", o.args)),
    vstore2hi_neon  := (self, o, i, is) >> Print(Blanks(i), self.printf(
	    "vst1_f32($1,vget_high_f32($2));\n", o.args)),

    data := (self, o, i, is) >> Cond(
	o.var.t=TVect(T_Real(32), 4) or o.var.t=TVect(TReal, 4),
	    Print(Blanks(i), self.printf(
              "const float32x4_t $1 = {(float32_t)$2,(float32_t)$3,(float32_t)$4,(float32_t)$5}",
	      [o.var] :: List(o.value.v, x -> x.v)), ";\n", self(o.cmd, i, is)),

	o.var.t=TVect(T_Real(32), 2) or o.var.t=TVect(TReal, 2),
	    Print(Blanks(i), self.printf(
	      "const float32x2_t $1 = {(float32_t)$2, (float32_t)$3}",
	      [o.var] :: List(o.value.v, x->x.v)), ";\n", self(o.cmd, i, is)),

	o.var.t = T_Real(32) or o.var.t = TReal,
	    Print(Blanks(i), self.printf(
	      "float32_t $1 = (float32_t)$2", [o.var, o.value.v]), ";\n", self(o.cmd, i, is)),

        Inherited(o, i, is)),

    vpack  := (self, o, i, is) >> Cond(
	o.t=TVect(T_Real(32), 4) or o.t=TVect(TReal, 4),
	self.printf("(float32x4_t){(float32_t)$1, (float32_t)$2, (float32_t)$3, (float32_t)$4}", o.args),
	o.t=TVect(T_Real(32), 2) or o.t=TVect(TReal, 2),
	self.printf("(float32x2_t){(float32_t)$1, (float32_t)$2}", o.args),

	Inherited(o, i, is)),

    vdup := (self, o, i, is) >> Cond(
	o.t=TVect(T_Real(32), 4) or o.t=TVect(TReal, 4),
	    self.printf("vdupq_n_f32($1)", o.args),
	o.t=TVect(T_Real(32), 2) or o.t=TVect(TReal, 2),
	    self.printf("vdup_n_f32($1)", o.args),
	Inherited(o, i, is)),

    vdup_lane_neon := (self, o, i, is) >> Cond(
	o.args[2] = 1 or o.args[2] = 2,
	    self.printf("vdupq_lane_f32(vget_low_f32($1), $2)", [o.args[1], o.args[2]-1]),
	o.args[2] = 3 or o.args[2] = 4,
	    self.printf("vdupq_lane_f32(vget_high_f32($1), $2)", [o.args[1], o.args[2]-3])
    ),

    vdup_lane_half := (self, o, i, is) >> self.prefix("vdup_lane_f32",
	[o.args[1], o.args[2]-1]),

    vswapcx_neon := (self, o, i, is) >> self.printf("vrev64q_f32($1)", o.args),

    vswapcx_half := (self, o, i, is) >> self.printf("vrev64_f32($1)", o.args),

    mul := (self, o, i, is) >> Cond(ObjId(o.t) = TVect and o.t.size=4,
	CondPat( o,
	    [mul, @TVect, @TVect],
	      self.printf("vmulq_f32($1, $2)", o.args),
	    [mul, @TReal, @TVect],
	      self.printf("vmulq_f32(vdupq_n_f32($1), $2)", o.args),
	    [mul, @TReal, @TReal],
	      self.printf("($1 * $2)", o.args),
	    Inherited(o, i, is)),

	CondPat( o,
	    [mul, @TVect, @TVect],
	      self.printf("vmul_f32($1, $2)", o.args),
	    [mul, @TReal, @TVect],
	      self.printf("vmul_f32(vdup_n_f32($1), $2)", o.args),
	    [mul, @TReal, @TReal],
	      self.printf("($1 * $2)", o.args),
	    Inherited(o, i, is))),

    mul  := (self, o, i, is) >> CondPat( o,
	[mul, @(1).cond(x -> x=TVect(T_Real(32), 4) or x=TVect(TReal,4)),
	      @(2).cond(x -> x=TVect(T_Real(32), 4) or x=TVect(TReal,4))],
	    self.printf("vmulq_f32($1, $2)", o.args),
	[mul, @(1).cond(x -> x=TVect(T_Real(32), 2) or x=TVect(TReal,2)),
	      @(2).cond(x -> x=TVect(T_Real(32), 2) or x=TVect(TReal, 2))],
	    self.printf("vmul_f32($1, $2)", o.args),
	[mul, @TReal, @(1).cond(x -> x=TVect(T_Real(32), 4) or x=TVect(TReal,4))],
	    self.printf("vmulq_f32(vdupq_n_f32($1), $2)", o.args),
	[mul, @TReal, @(1).cond(x -> x=TVect(T_Real(32), 2) or x=TVect(TReal,2))],
	    self.printf("vmul_f32(vdup_n_f32($1), $2)", o.args),
	[mul, @TReal, @TReal],
	    self.printf("($1 * $2)", o.args),

        Inherited(o.ooo, i, is)),

    add  := (self, o, i, is) >> Cond(
	not IsVecT(o.t),
	  self.pinfix(o.args, " + "),
	Length(o.args)=2 and o.t.size=4,
	  self.printf("vaddq_f32($1, $2)", o.args),
	Length(o.args)=2 and o.t.size=2, self.printf("vadd_f32($1, $2)", o.args),
	let(rem := ApplyFunc(add, Drop(o.args, 1)),
	    self.printf(Cond(o.t.size=4,"vaddq_f32($1,$2)","vadd_f32($1,$2)"), [o.args[1], rem]))),

    sub  := (self, o, i, is) >> Cond(
	not IsVecT(o.t),                 self.pinfix(o.args, " - "),
	Length(o.args)=2 and o.t.size=4, self.printf("vsubq_f32($1, $2)", o.args),
	Length(o.args)=2 and o.t.size=2, self.printf("vsub_f32($1, $2)", o.args),
	let(rem := ApplyFunc(add, Drop(o.args, 1)),
	    self.printf(Cond(o.t.size=4,"vsubq_f32($1,$2)","vsub_f32($1,$2)"),
		[o.args[1], rem]))),

		neg  := (self, o, i, is) >> Cond(
		    IsVecT(o.t),
		        self.printf(Cond(o.t.size=4,"vnegq_f32($1)","vneg_f32($1)"), o.args),
		    o.t=TReal,
		        self.printf("(-$1)", o.args),
		    Inherited(o, i, is)),

		fma  := (self, o, i, is) >> CondPat(o,
		    [fma, @TVect, @TVect, @TVect],
		       self.printf(Cond(o.t.size=4,"vmlaq_f32($1, $2, $3)","vmla_f32($1, $2, $3)"), o.args),
		    [fma, @TVect, @TReal, @TVect],
		       self.printf(Cond(o.t.size=4,"vmlaq_f32($1, vdupq_n_f32($2), $3)","vmla_f32($1, vdup_n_f32($2), $3)"), o.args),
		    [fma, @TReal, @TReal, @TReal],
		       self.printf("($1 + ($2 * $3))", o.args),
		    Inherited(o, i, is)),

		fms  := (self, o, i, is) >> CondPat(o,
		    [fms, @TVect, @TVect, @TVect],
		       self.printf(Cond(o.t.size=4,"vmlsq_f32($1, $2, $3)","vmls_f32($1, $2, $3)"), o.args),
		    [fms, @TVect, @TReal, @TVect],
		       self.printf(Cond(o.t.size=4,"vmlsq_f32($1, vdupq_n_f32($2), $3)","vmls_f32($1, vdup_n_f32($2), $3)"), o.args),
		    [fms, @TReal, @TReal, @TReal],
		       self.printf("($1 - ($2 * $3)", o.args),
		    Inherited(o, i, is)),

		nfma  := (self, o, i, is) >> CondPat(o,
		    [nfma, @TVect, @TVect, @TVect],
		       self.printf(Cond(o.t.size=4,"vmlaq_f32(vnegq_f32($1), $2, $3)","vmla_f32(vneg_f32($1), $2, $3)"), o.args),
		    [nfma, @TVect, @TReal, @TVect],
		       self.printf(Cond(o.t.size=4,"vmlaq_f32(vnegq_f32($1), vdupq_n_f32($2), $3)","vmla_f32(vneg_f32($1), vdup_n_f32($2), $3)"), o.args),
		    [nfma, @TReal, @TReal, @TReal],
		       self.printf("(($2 * $3) - $1)", o.args),
		    Inherited(o, i, is)),

		_const := (self,o) >> Cond(
		    o.t = T_Real(32) or o.t = TReal,
		        ["float32_t", [o.v]],
		    IsVecT(o.t), Cond(
			o.t.t = T_Real(32) and o.t.size=4,
			    ["float32x4_t", List(o.v, x->x.v)],
			o.t.t = T_Real(32) and o.t.size=2,
			    ["float32x2_t", List(o.v, x->x.v)],
			Inherited(o)
		    ),
#				o.t = TString,                     ["", [o.v]],
		    Inherited(o)
		)
));


text:
Read("SchreierSims.gi");

#######################################################################
# Description: The deterministic verification, parallelized w.r.t.    #
#              base points, described as method 1 in section 5.1.     #
# Input: - BSGS data structure returned by RandomSchreierSims         #
# Output: - boolean, true if BSGS is correct, else false              #
#######################################################################

Method1 := function(BSGS)
    local l, tasks, results;

    Print("Start verification.\n");
    Print("Apply parallel computing method w.r.t. base points.\n");

    tasks := [1..Length(BSGS[1])];  results := [1..Length(BSGS[1])];

    for l in [1..Length(BSGS[1])] do
        tasks[l] := RunTask(SimpleVerify, BSGS, l);
    od;
    for l in [1..Length(BSGS[1])] do
        results[l] := TaskResult(tasks[l]);
    od;

    Print("\nVerification done.");

    if false in results then
        Print("\nBSGS is wrong!\n");
    else
        Print("\nCorrectness of BSGS verified!\n");
    fi;
    Print("*********************************************************\n");
end;

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

#######################################################################
# Description: The deterministic verification, parallelized w.r.t.    #
#              orbit points, described as method 2 in section 5.1.    #
# Input: - BSGS data structure returned by RandomSchreierSims         #
# Output: - boolean, true if BSGS is correct, else false              #
#######################################################################

Method2 := function(BSGS)
    local l, tasks, results;

    Print("Start verification.\n");
    Print("Apply parallel computing method w.r.t. orbit points.\n");

    tasks := [1..Length(BSGS[1])];  results := [1..Length(BSGS[1])];

    for l in Reversed([1..Length(BSGS[1])]) do
        results[l] := ParallelVerify(BSGS, l);
        if false in results then
            break;
        fi;
    od;

    Print("\nVerification done.");

    if false in results then
        Print("\nBSGS is wrong!\n");
    else
        Print("\nCorrectness of BSGS verified!\n");
    fi;
    Print("*********************************************************\n");
end;

max_stars_repo_name: AlexanderKlemps/parallel-schreier-sims

max_issues_repo_name: AlexanderKlemps/parallel-schreier-sims

max_forks_repo_name: AlexanderKlemps/parallel-schreier-sims


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

# Convention: N - range
#             n - domain

## H(<N>,<n>,<base>,<stride>)
##    Stride index mapping: i -> base + stride * i
#
Class(H, FuncClass, rec(
    def := (N,n,base,stride) -> rec(),
    domain := self >> self.params[2],
    range := self >> self.params[1],
    lambda := self >> let(
	base := self.params[3], str := self.params[4], i := Ind(self.params[2]),
	Lambda(i, base + str * i))
));

# HZ(<N>,<n>,<base>,<stride>)
##    Mod-stride index mapping: i -> (base + stride * i) mod N
##
Class(HZ, FuncClass, rec(
    isCyclic := true,

    def := (N,n,base,str) -> rec(),
    domain := self >> self.params[2],
    range := self >> self.params[1],

    lambda := self >> let(
	n := self.params[2], N := self.params[1], base := self.params[3], str := self.params[4],
	i := Ind(n),
	ii := When(str * n = N, no_mod(i), i), # condition to omit the mod is str*n = N
	Lambda(i, imod(base + str * ii, N)))
));

#F RM(<N>,<n>,<phi>,<g>)
##    Exponential index mapping: i -> (phi * g^i) mod N
##    <g> - generating element for integers mod <N>.
##    <phi> - "phase"
##
Class(RM, FuncClass, rec(
    isCyclic := true,
    abbrevs := [ N    -> [N, N-1, 1, PrimitiveRootMod(N)],
	        (N,n) -> [N, n,   1, PrimitiveRootMod(N)] ],

    domain := self >> self.params[2],
    range := self >> self.params[1],
    def := (N,n,phi,g) -> Checked(IsPosIntSym(N), IsPosIntSym(n), n <= N, rec()),

    lambda := self >> let(
	phi := self.params[3], g := self.params[4], i := Ind(self.params[2]),
	Lambda(i, powmod(phi, g, i, self.params[1]))),
));

#F RR(<N>,<phi>,<g>)
##    Rader permutation: i -> {    0,     if i = 0 }
##                            { RM(i-1),  if i > 0 }
##    <phi> and <g> are parameters to RM
##
Class(RR, PermClass, rec(
    isCyclic := true,

    abbrevs := [ N -> [N,1,PrimitiveRootMod(N)] ],

    def := (N,phi,g) -> Checked(IsPosInt(N), IsPrime(N), rec()),

    domain := self >> self.params[1],
    range := self >> self.params[1],

    # alternatively cond(i, imod(phi * pow(g, i-1), N), 0) we do g^-1 for easier poweropt
    lambda  := self >> let(
	N := self.params[1], phi := self.params[2], g := self.params[3], i:= Ind(N),
	When(not self.transposed,
	    Lambda(i, cond(i, powmod(phi*(g^-1 mod N) mod N, g, i, N), 0)),
	    Lambda(i, cond(i, 1+ilogmod(cond(i,i,1)/phi mod N, g, N), 0)))
    )
));


text:
#############################################################################
##
#W  id384e.gap             GAP library of id's             Hans Ulrich Besche
##

ID_GROUP_TREE.next[384].next[137].next[1].next[16]:=
rec(
  fp:= [ 1476, 2636, 6600, 12238, 14104, 14474, 14980, 18057, 19809, 28110,
31099, 32061, 33987, 36179, 44069, 44294, 45675, 46424, 48838, 57531, 64443,
66649, 68297, 69424, 70887, 75849, 78400, 78938, 80096, 81806, 82147, 83246,
85867, 87468, 89093, 94091, 95053, 99213 ],
  next:= [ rec(
    desc:= [ 108003 ],
    fp:= [ 10, 416, 614, 812 ],
    next:= [ 8918, rec(
      desc:= [ 113011 ],
      fp:= [ 22, 814 ],
      next:= [ 8927, 15548 ] ), rec(
      desc:= [ 302009 ],
      fp:= [ 742, 1081 ],
      next:= [ 8915, rec(
        desc:= [ 206002 ],
        fp:= [ 4, 14 ],
        next:= [ 8897, 8941 ] ) ] ), rec(
      desc:= [ 302008 ],
      fp:= [ 64681, 74821 ],
      next:= [ rec(
        desc:= [ 108003 ],
        fp:= [ 8, 414 ],
        next:= [ 8916, 8913 ] ), 15539 ] ) ] ), rec(
    desc:= [ 109003 ],
    fp:= [ 6, 214, 412 ],
    next:= [ rec(
      desc:= [ 302008, 302012, 110003, 210002 ],
      fp:= [ 2, 12 ],
      next:= [ 9565, 9568 ] ), 9574, rec(
      desc:= [ 302014, 211002 ],
      fp:= [ 4, 212 ],
      next:= [ 9564, 9569 ] ) ] ), rec(
    desc:= [ 109003 ],
    fp:= [ 414, 612 ],
    next:= [ 15444, rec(
      desc:= [ 302006, 302009, 211002 ],
      fp:= [ 214, 412 ],
      next:= [ 7365, 7364 ] ) ] ), rec(
    desc:= [ 107003 ],
    fp:= [ 8, 18, 216, 414, 612 ],
    next:= [ rec(
      desc:= [ 302005, 208002 ],
      fp:= [ 8, 18 ],
      next:= [ 8974, 8997 ] ), 8956, 8958, rec(
      desc:= [ 112011 ],
      fp:= [ 12, 418 ],
      next:= [ rec(
        desc:= [ 302007 ],
        fp:= [ 742, 1081 ],
        next:= [ 8990, rec(
          desc:= [ 114007 ],
          fp:= [ 4, 14 ],
          next:= [ 8981, 15565 ] ) ] ), 8977 ] ), 8989 ] ), rec(
    desc:= [ 108003 ],
    fp:= [ 8, 18, 216, 414, 612 ],
    next:= [ rec(
      desc:= [ 302007, 109003, 209003, 303010, 213003 ],
      fp:= [ 8, 18 ],
      next:= [ 14950, 14949 ] ), 14925, 14936, rec(
      desc:= [ 302008 ],
      fp:= [ 742, 1081 ],
      next:= [ 14961, rec(
        desc:= [ 107003 ],
        fp:= [ 4, 14 ],
        next:= [ 14963, 14929 ] ) ] ), rec(
      desc:= [ 107003 ],
      fp:= [ 216, 612 ],
      next:= [ 14937, 14930 ] ) ] ), rec(
    desc:= [ 106003 ],
    fp:= [ 4, 14 ],
    next:= [ 7444, 15465 ] ), rec(
    desc:= [ 105003 ],
    fp:= [ 16, 412 ],
    next:= [ 15568, 8999 ] ), rec(
    desc:= [ 302008 ],
    fp:= [ 44881, 54682, 55021 ],
    next:= [ rec(
      desc:= [ 208002 ],
      fp:= [ 4, 311 ],
      next:= [ 9680, 9677 ] ), rec(
      desc:= [ 303008 ],
      fp:= [ 2, 82, 181, 261 ],
      next:= [ 9676, 9665, 9679, 9678 ] ), rec(
      desc:= [ 109003 ],
      fp:= [ 4, 14 ],
      next:= [ 9675, 9664 ] ) ] ), rec(
    desc:= [ 107003 ],
    fp:= [ 616, 814, 1012 ],
    next:= [ rec(
      desc:= [ 302007 ],
      fp:= [ 643, 74581 ],
      next:= [ 9622, 9618 ] ), rec(
      desc:= [ 302007 ],
      fp:= [ 544, 64681 ],
      next:= [ 9649, rec(
        desc:= [ 302007 ],
        fp:= [ 4, 183, 282, 381, 461 ],
        next:= [ 9659, 9648, 9654, 9646, 9653 ] ) ] ), 9652 ] ), rec(
    desc:= [ 107003 ],
    fp:= [ 4, 14 ],
    next:= [ rec(
      desc:= [ 302006, 113007, 303009, 214003 ],
      fp:= [ 214, 412 ],
      next:= [ 9668, 9667 ] ), 9660 ] ), rec(
    desc:= [ 114007 ],
    fp:= [ 16, 412 ],
    next:= [ 15466, rec(
      desc:= [ 302006, 302008, 216002 ],
      fp:= [ 14, 212 ],
      next:= [ 7447, 7446 ] ) ] ), rec(
    desc:= [ 302008, 110003, 210003 ],
    fp:= [ 4, 212 ],
    next:= [ 9554, 9559 ] ), rec(
    desc:= [ 110003 ],
    fp:= [ 10, 614, 812 ],
    next:= [ 7326, rec(
      desc:= [ 116011 ],
      fp:= [ 8, 414 ],
      next:= [ rec(
        desc:= [ 302011 ],
        fp:= [ 742, 1081 ],
        next:= [ 15439, 15427 ] ), rec(
        desc:= [ 302010 ],
        fp:= [ 643, 74581 ],
        next:= [ rec(
          desc:= [ 208002 ],
          fp:= [ 214, 412 ],
          next:= [ rec(
            desc:= [ 210002 ],
            fp:= [ 4, 14 ],
            next:= [ 7353, 7338 ] ), rec(
            desc:= [ 108003 ],
            fp:= [ 14, 212 ],
            next:= [ 7305, 7324 ] ) ] ), 7310 ] ) ] ), rec(
      desc:= [ 302010 ],
      fp:= [ 544, 64681 ],
      next:= [ rec(
        desc:= [ 302007 ],
        fp:= [ 3, 281 ],
        next:= [ rec(
          desc:= [ 108003, 302010, 304008, 210004, 40307010 ],
          fp:= [ 1, 641 ],
          next:= [ 7349, 7350 ] ), 7325 ] ), 7322 ] ) ] ), rec(
    desc:= [ 108003 ],
    fp:= [ 8, 18, 414, 612 ],
    next:= [ rec(
      desc:= [ 106003 ],
      fp:= [ 4, 14 ],
      next:= [ 8951, 15554 ] ), 15556, rec(
      desc:= [ 105003 ],
      fp:= [ 4, 14 ],
      next:= [ 15553, 8949 ] ), 8952 ] ), rec(
    desc:= [ 302006, 109003, 209002 ],
    fp:= [ 218, 614 ],
    next:= [ 9469, 9468 ] ), rec(
    desc:= [ 112007 ],
    fp:= [ 4, 14 ],
    next:= [ rec(
      desc:= [ 105003, 302005, 210002 ],
      fp:= [ 14, 212 ],
      next:= [ 9002, 9001 ] ), 15570 ] ), rec(
    desc:= [ 107003 ],
    fp:= [ 12, 418, 616, 814, 1012 ],
    next:= [ 9645, 9580, rec(
      desc:= [ 302006 ],
      fp:= [ 384, 48361 ],
      next:= [ 9592, 9586 ] ), 9619, rec(
      desc:= [ 302007 ],
      fp:= [ 54781, 64582, 64921 ],
      next:= [ 9628, 9630, 9629 ] ) ] ), rec(
    desc:= [ 109003 ],
    fp:= [ 8, 414, 612 ],
    next:= [ rec(
      desc:= [ 302006 ],
      fp:= [ 384, 48361 ],
      next:= [ 9517, rec(
        desc:= [ 211002 ],
        fp:= [ 216, 414 ],
        next:= [ 9536, 9513 ] ) ] ), rec(
      desc:= [ 302006 ],
      fp:= [ 384, 48361 ],
      next:= [ rec(
        desc:= [ 208002 ],
        fp:= [ 8, 414 ],
        next:= [ 9508, 9526 ] ), rec(
        desc:= [ 302012 ],
        fp:= [ 4, 442, 541, 781 ],
        next:= [ 9486, 9540, rec(
          desc:= [ 210002 ],
          fp:= [ 8, 414 ],
          next:= [ 9503, 9516 ] ), 9530 ] ) ] ), rec(
      desc:= [ 302006 ],
      fp:= [ 384, 48361 ],
      next:= [ rec(
        desc:= [ 210002 ],
        fp:= [ 214, 412 ],
        next:= [ 9515, 9514 ] ), rec(
        desc:= [ 209002 ],
        fp:= [ 18, 414 ],
        next:= [ 9518, 9537 ] ) ] ) ] ), rec(
    desc:= [ 108003 ],
    fp:= [ 8, 216, 414, 612 ],
    next:= [ rec(
      desc:= [ 302009 ],
      fp:= [ 544, 64681 ],
      next:= [ 7409, rec(
        desc:= [ 302008 ],
        fp:= [ 461, 18361 ],
        next:= [ rec(
          desc:= [ 302011 ],
          fp:= [ 342, 681 ],
          next:= [ 7408, 7438 ] ), 7407 ] ) ] ), 15446, rec(
      desc:= [ 114011 ],
      fp:= [ 8, 18, 414 ],
      next:= [ 15463, 7388, rec(
        desc:= [ 302008 ],
        fp:= [ 742, 1081 ],
        next:= [ rec(
          desc:= [ 207002 ],
          fp:= [ 18, 414 ],
          next:= [ 7387, 7402 ] ), 7400 ] ) ] ), rec(
      desc:= [ 302006 ],
      fp:= [ 48361, 58341 ],
      next:= [ rec(
        desc:= [ 302012 ],
        fp:= [ 343, 442, 781, 44581 ],
        next:= [ rec(
          desc:= [ 212002 ],
          fp:= [ 2, 12 ],
          next:= [ 7420, 7439 ] ), rec(
          desc:= [ 211002 ],
          fp:= [ 214, 412 ],
          next:= [ 7426, 7428 ] ), 7435, 7432 ] ), rec(
        desc:= [ 212002 ],
        fp:= [ 14, 212 ],
        next:= [ 15452, 15459 ] ) ] ) ] ), rec(
    desc:= [ 302010 ],
    fp:= [ 445, 54781, 64582 ],
    next:= [ rec(
      desc:= [ 215002 ],
      fp:= [ 12, 418, 814 ],
      next:= [ rec(
        desc:= [ 115011 ],
        fp:= [ 12, 814 ],
        next:= [ 15542, 8945 ] ), rec(
        desc:= [ 115007 ],
        fp:= [ 4, 14 ],
        next:= [ 8920, 8921 ] ), rec(
        desc:= [ 115011 ],
        fp:= [ 8, 414 ],
        next:= [ 15544, rec(
          desc:= [ 108003, 112003, 303008, 303010, 214003 ],
          fp:= [ 214, 412 ],
          next:= [ 8936, 8934 ] ) ] ) ] ), rec(
      desc:= [ 209003 ],
      fp:= [ 214, 412 ],
      next:= [ 8943, 8919 ] ), rec(
      desc:= [ 209004 ],
      fp:= [ 6, 214, 412 ],
      next:= [ 8944, 8902, 8908 ] ) ] ), rec(
    desc:= [ 302006, 108003, 208003 ],
    fp:= [ 4, 212 ],
    next:= [ 14967, 14968 ] ), rec(
    desc:= [ 108003 ],
    fp:= [ 8, 216, 414, 612 ],
    next:= [ rec(
      desc:= [ 302008 ],
      fp:= [ 64681, 74821 ],
      next:= [ 8988, rec(
        desc:= [ 208003 ],
        fp:= [ 8, 612 ],
        next:= [ 8992, 8983 ] ) ] ), 15558, rec(
      desc:= [ 302010 ],
      fp:= [ 544, 64681 ],
      next:= [ rec(
        desc:= [ 207002 ],
        fp:= [ 8, 414 ],
        next:= [ 8967, 8973 ] ), 8959 ] ), rec(
      desc:= [ 302010 ],
      fp:= [ 544, 64681 ],
      next:= [ rec(
        desc:= [ 303014 ],
        fp:= [ 263, 362, 461 ],
        next:= [ 8980, rec(
          desc:= [ 211003 ],
          fp:= [ 2, 12 ],
          next:= [ 15561, 15564 ] ), 8991 ] ), rec(
        desc:= [ 208002 ],
        fp:= [ 216, 414 ],
        next:= [ 8987, 8979 ] ) ] ) ] ), rec(
    desc:= [ 108003 ],
    fp:= [ 6, 214, 412 ],
    next:= [ 9566, 9553, 9567 ] ), 9535, rec(
    desc:= [ 108003 ],
    fp:= [ 818, 1016, 1214, 1412 ],
    next:= [ 15547, rec(
      desc:= [ 302008 ],
      fp:= [ 64582, 64921 ],
      next:= [ 8914, 8896 ] ), rec(
      desc:= [ 105003 ],
      fp:= [ 8, 414 ],
      next:= [ 8942, 15552 ] ), 8917 ] ), rec(
    desc:= [ 105003 ],
    fp:= [ 6, 214, 412 ],
    next:= [ 9534, rec(
      desc:= [ 107003 ],
      fp:= [ 8, 414 ],
      next:= [ rec(
        desc:= [ 302007 ],
        fp:= [ 742, 1081 ],
        next:= [ 9498, rec(
          desc:= [ 302014, 112003, 212003 ],
          fp:= [ 6, 412 ],
          next:= [ 9495, 9490 ] ) ] ), rec(
        desc:= [ 302010 ],
        fp:= [ 544, 74821 ],
        next:= [ 9473, 9519 ] ) ] ), rec(
      desc:= [ 302009 ],
      fp:= [ 544, 64681 ],
      next:= [ rec(
        desc:= [ 206002 ],
        fp:= [ 4, 14, 212 ],
        next:= [ 9500, 9533, 9507 ] ), rec(
        desc:= [ 302015 ],
        fp:= [ 263, 362, 621, 16621, 36421 ],
        next:= [ 9523, rec(
          desc:= [ 110003 ],
          fp:= [ 8, 414 ],
          next:= [ 9499, 9509 ] ), 9544, 9532, rec(
          desc:= [ 113003 ],
          fp:= [ 8, 414 ],
          next:= [ 9531, rec(
            desc:= [ 305011, 214005 ],
            fp:= [ 4, 212 ],
            next:= [ 9494, 9489 ] ) ] ) ] ) ] ) ] ), rec(
    desc:= [ 109003 ],
    fp:= [ 10, 218, 416, 614, 812 ],
    next:= [ 7270, 7285, rec(
      desc:= [ 302009 ],
      fp:= [ 742, 1081 ],
      next:= [ rec(
        desc:= [ 207002 ],
        fp:= [ 16, 412 ],
        next:= [ 7288, rec(
          desc:= [ 107003 ],
          fp:= [ 4, 14 ],
          next:= [ 7346, 7265 ] ) ] ), 7286 ] ), rec(
      desc:= [ 302009 ],
      fp:= [ 643, 74581 ],
      next:= [ rec(
        desc:= [ 108003 ],
        fp:= [ 4, 14 ],
        next:= [ 7339, 7269 ] ), 7287 ] ), rec(
      desc:= [ 302006 ],
      fp:= [ 841, 58261 ],
      next:= [ rec(
        desc:= [ 108003 ],
        fp:= [ 4, 14 ],
        next:= [ 15420, 15432 ] ), rec(
        desc:= [ 108003 ],
        fp:= [ 16, 412 ],
        next:= [ 7345, 7314 ] ) ] ) ] ), rec(
    desc:= [ 109003 ],
    fp:= [ 10, 218, 416, 614, 812 ],
    next:= [ 8895, 8899, rec(
      desc:= [ 302010 ],
      fp:= [ 643, 74581 ],
      next:= [ rec(
        desc:= [ 107003 ],
        fp:= [ 4, 14 ],
        next:= [ 8939, 8925 ] ), rec(
        desc:= [ 108003 ],
        fp:= [ 4, 14 ],
        next:= [ 8900, 8898 ] ) ] ), rec(
      desc:= [ 302009 ],
      fp:= [ 1321, 74581 ],
      next:= [ 8940, rec(
        desc:= [ 109003 ],
        fp:= [ 4, 14 ],
        next:= [ 8901, 8922 ] ) ] ), rec(
      desc:= [ 114011 ],
      fp:= [ 22, 814 ],
      next:= [ 8938, rec(
        desc:= [ 107003 ],
        fp:= [ 4, 14 ],
        next:= [ 15538, rec(
          desc:= [ 116007 ],
          fp:= [ 4, 14 ],
          next:= [ 8937, 15549 ] ) ] ) ] ) ] ), rec(
    desc:= [ 302006, 108003, 208003 ],
    fp:= [ 4, 212 ],
    next:= [ 14924, 14922 ] ), rec(
    desc:= [ 302009 ],
    fp:= [ 445, 54781, 64582 ],
    next:= [ rec(
      desc:= [ 110003 ],
      fp:= [ 10, 614, 812 ],
      next:= [ 7343, rec(
        desc:= [ 116011 ],
        fp:= [ 8, 414 ],
        next:= [ rec(
          desc:= [ 206002 ],
          fp:= [ 4, 14 ],
          next:= [ 15430, 15436 ] ), rec(
          desc:= [ 206002 ],
          fp:= [ 4, 14 ],
          next:= [ 7308, 7330 ] ) ] ), rec(
        desc:= [ 209002 ],
        fp:= [ 2, 111 ],
        next:= [ 7341, 7340 ] ) ] ), rec(
      desc:= [ 111003, 211003 ],
      fp:= [ 216, 414 ],
      next:= [ 7342, 7344 ] ), rec(
      desc:= [ 302007 ],
      fp:= [ 4, 381 ],
      next:= [ 7355, rec(
        desc:= [ 112003, 212004 ],
        fp:= [ 6, 412 ],
        next:= [ 7299, 7309 ] ) ] ) ] ), rec(
    desc:= [ 107003 ],
    fp:= [ 8, 18, 216, 414, 612 ],
    next:= [ rec(
      desc:= [ 105003 ],
      fp:= [ 18, 414 ],
      next:= [ rec(
        desc:= [ 302006, 210002 ],
        fp:= [ 4, 14 ],
        next:= [ 14960, 14958 ] ), rec(
        desc:= [ 302007 ],
        fp:= [ 483, 58261 ],
        next:= [ 14948, 14959 ] ) ] ), 14955, 14938, rec(
      desc:= [ 105003 ],
      fp:= [ 18, 414 ],
      next:= [ 14952, rec(
        desc:= [ 302007 ],
        fp:= [ 483, 58261 ],
        next:= [ 14927, rec(
          desc:= [ 302007 ],
          fp:= [ 5, 481 ],
          next:= [ rec(
            desc:= [ 210002 ],
            fp:= [ 4, 14 ],
            next:= [ 14964, 14962 ] ), 14951 ] ) ] ) ] ), rec(
      desc:= [ 302006 ],
      fp:= [ 483, 58261 ],
      next:= [ rec(
        desc:= [ 208002 ],
        fp:= [ 214, 412 ],
        next:= [ 14942, 14945 ] ), rec(
        desc:= [ 211002 ],
        fp:= [ 14, 113, 212 ],
        next:= [ 14935, 14947, 14939 ] ) ] ) ] ), rec(
    desc:= [ 106003 ],
    fp:= [ 8, 414, 612 ],
    next:= [ 8998, rec(
      desc:= [ 111011 ],
      fp:= [ 12, 418 ],
      next:= [ 15566, 8975 ] ), 8994 ] ), rec(
    desc:= [ 108003 ],
    fp:= [ 8, 18, 414 ],
    next:= [ rec(
      desc:= [ 106003 ],
      fp:= [ 4, 14 ],
      next:= [ 7362, 15442 ] ), 15440, 7358 ] ), rec(
    desc:= [ 109003 ],
    fp:= [ 10, 416, 614, 812 ],
    next:= [ 7319, rec(
      desc:= [ 115011 ],
      fp:= [ 8, 18 ],
      next:= [ 15418, 7267 ] ), rec(
      desc:= [ 302009 ],
      fp:= [ 643, 74581 ],
      next:= [ rec(
        desc:= [ 302007 ],
        fp:= [ 282, 461 ],
        next:= [ 7312, 7317 ] ), 7316 ] ), rec(
      desc:= [ 302009 ],
      fp:= [ 544, 64681 ],
      next:= [ rec(
        desc:= [ 107003 ],
        fp:= [ 8, 414 ],
        next:= [ 7321, 15433 ] ), 7313 ] ) ] ), rec(
    desc:= [ 107003 ],
    fp:= [ 8, 18, 216, 414, 612 ],
    next:= [ rec(
      desc:= [ 302007 ],
      fp:= [ 64681, 74821 ],
      next:= [ rec(
        desc:= [ 207002 ],
        fp:= [ 16, 214 ],
        next:= [ rec(
          desc:= [ 107003 ],
          fp:= [ 16, 412 ],
          next:= [ 7414, 7410 ] ), rec(
          desc:= [ 107003 ],
          fp:= [ 2, 12 ],
          next:= [ 7412, 7401 ] ) ] ), 7419 ] ), 7371, rec(
      desc:= [ 302008 ],
      fp:= [ 643, 1321, 74581 ],
      next:= [ 15447, 7379, 7406 ] ), rec(
      desc:= [ 113011 ],
      fp:= [ 8, 18, 414 ],
      next:= [ rec(
        desc:= [ 302005 ],
        fp:= [ 841, 58261 ],
        next:= [ rec(
          desc:= [ 208002 ],
          fp:= [ 4, 14 ],
          next:= [ 15464, 15462 ] ), rec(
          desc:= [ 209002 ],
          fp:= [ 4, 14, 212 ],
          next:= [ 7442, 7415, 7434 ] ) ] ), 7378, rec(
        desc:= [ 302007 ],
        fp:= [ 742, 1081 ],
        next:= [ rec(
          desc:= [ 206002 ],
          fp:= [ 16, 412 ],
          next:= [ 7373, 7403 ] ), 7382 ] ) ] ), rec(
      desc:= [ 302009 ],
      fp:= [ 544, 64681, 74821 ],
      next:= [ rec(
        desc:= [ 302008 ],
        fp:= [ 342, 441, 681 ],
        next:= [ 7413, rec(
          desc:= [ 208002 ],
          fp:= [ 2, 111 ],
          next:= [ 15453, 15450 ] ), 7436 ] ), rec(
        desc:= [ 207002 ],
        fp:= [ 6, 16 ],
        next:= [ 7411, 7418 ] ), rec(
        desc:= [ 207002 ],
        fp:= [ 16, 412 ],
        next:= [ 7437, 7433 ] ) ] ) ] ), 7323, rec(
    desc:= [ 108003 ],
    fp:= [ 12, 418, 616, 814, 1012 ],
    next:= [ 9638, 9578, rec(
      desc:= [ 302008 ],
      fp:= [ 643, 1321, 74581 ],
      next:= [ rec(
        desc:= [ 302010 ],
        fp:= [ 442, 541, 781, 1021 ],
        next:= [ 9591, rec(
          desc:= [ 215004 ],
          fp:= [ 218, 614, 812 ],
          next:= [ 9603, 9614, 9595 ] ), 9599, 9598 ] ), 9616, 9615 ] ), rec(
      desc:= [ 302008 ],
      fp:= [ 544, 64681 ],
      next:= [ rec(
        desc:= [ 208002 ],
        fp:= [ 3, 112, 211 ],
        next:= [ 9626, 9644, 9642 ] ), rec(
        desc:= [ 108003 ],
        fp:= [ 4, 14 ],
        next:= [ 9627, 9617 ] ) ] ), rec(
      desc:= [ 302008 ],
      fp:= [ 445, 54781, 64582, 64921 ],
      next:= [ 9624, rec(
        desc:= [ 302008 ],
        fp:= [ 3, 182, 281, 361 ],
        next:= [ 9641, 9625, rec(
          desc:= [ 210002 ],
          fp:= [ 3, 112 ],
          next:= [ 9647, 9643 ] ), 9639 ] ), 9631, 9650 ] ) ] ), rec(
    desc:= [ 109003 ],
    fp:= [ 818, 1016, 1214, 1412 ],
    next:= [ 15419, rec(
      desc:= [ 302009 ],
      fp:= [ 445, 54781 ],
      next:= [ 7311, 7315 ] ), 15434, 7320 ] ) ] );

max_stars_repo_name: hulpke/smallgrp

max_issues_repo_name: hulpke/smallgrp

max_forks_repo_name: gap-packages/smallgrp


text:
# As is, will not work if val is a String
Assign := function(var, val)
	Read(InputTextString(Concatenation(var, " := ", String(val), ";")));
end;


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

Class(RulesId, RuleSet);

RewriteRules(RulesId, rec(

 IdISumRight := ARule( Compose,
       [ @(1, Id), @(2, ISum, canReorder) ],
  e -> [ ISum(@2.val.var, @2.val.domain, @1.val * @2.val.child(1)).attrs(@(2).val) ]),

 IdISumLeft := ARule( Compose,
       [ @(1, ISum, canReorder), @(2, Id) ],
  e -> [ ISum(@1.val.var, @1.val.domain, @1.val.child(1) * @2.val).attrs(@(1).val) ]),

 IdIComposeRight := ARule( Compose,
       [ @(1, Id), @(2, ICompose, canReorder) ],
  e -> [ ICompose(@2.val.var, @2.val.domain, @1.val * @2.val.child(1)).attrs(@(2).val) ]),

 IdIComposeLeft := ARule( Compose,
       [ @(1, ICompose, canReorder), @(2, Id) ],
  e -> [ ICompose(@1.val.var, @1.val.domain, @1.val.child(1) * @2.val).attrs(@(1).val) ]),

 # Gath * Id
 CommuteGathId := ARule( Compose,
       [ @(1, Gath), @(2, Id) ], # o 1-> 2->
  e -> [ Id(fCompose(@2.val.element, @1.val.func)).attrs(@(2).val), @1.val ]),

 # Prm * Id
 CommutePrmId := ARule( Compose,
       [ @(1, Prm), @(2, Id) ], # o 1-> 2->
  e -> [ Id(fCompose(@2.val.element, @1.val.func)).attrs(@(2).val), @1.val ]),

 # Id * Scat
 CommuteIdScat := ARule( Compose,
       [ @(1, Id), @(2, Scat) ], # <-1 <-2 o
  e -> [ @2.val, Id(fCompose(@1.val.element, @2.val.func)).attrs(@(1).val) ]),

  #(A X B) o (C X D) -> (A o C) X (B o D)
  CommuteIdTensorfTensor := Rule([fCompose, [idTensor, @(1), @(2)],
      [@(5, fTensor).cond(e->Length(e.children()) = 2), @(3).cond(e->e.range() = @1.val.domain()), @(4).cond(e->e.range() = @2.val.domain())]],
      e -> idTensor(fCompose(@(1).val, @(3).val), fCompose(@(2).val, @(4).val))
  ),

 IdNoDiagPullinRight := ARule( Compose,
       [ @(1, Id), @(2, NoDiagPullinRight) ],
  e -> [ NoDiagPullinRight(@(1).val * @2.val.child(1)) ]),

 IdNoDiagPullinLeft := ARule( Compose,
       [ @(1, NoDiagPullinLeft), @(2, Id) ],
  e -> [ NoDiagPullinLeft(@1.val.child(1) * @(2).val) ]),

#(A X B) o ((C X D) X E) -> (A o (C X D)) X (B o E)
  CommuteIdTensorfTensor2 := Rule([fCompose, [idTensor, @(1), @(2)],
      [@(3, fTensor).cond(e->
           Length(e.children()) = 3 and
       fTensor(e.child(1), e.child(2)).range() = @1.val.domain() and
       e.child(3).range() = @2.val.domain()
       ),
       @(4), @(5), @(6)
      ]], e->

      idTensor(fCompose(@(1).val, fTensor(@(3).val.child(1), @(3).val.child(2))),
               fCompose(@(2).val, @(3).val.child(3)))
      ),

  # i o (j)n = (j)1->N, where i is an idId
  IdComposeConst := Rule([fCompose, @(1, idId), @(2, fBase)],
                    e -> idConst(1, @2.val.params[2])
  ),

  idTensorId1 := ARule(idTensor, [@(1), [@(2,[idId]), 1]], e -> [@(1).val]),
  idTensorId2 := ARule(idTensor, [[@(2,[idId,J]), 1], @(1)], e -> [@(1).val]),

  # This is a special-puprose modification of the IdComposeConst rule, above.
  # It's possible that this is not the best way to do this.  I'd be willing
  # to change it later (Peter).
  # i o ((j)n X (k)m ) = (j)1->N X (k)1->N
#  IdComposeTensConst := Rule([fCompose,
#     @(1, idId),
#     @(2, fTensor).cond(e->
#         Length(e.children()) = 2  and
#         ObjId(e.child(1)) = fBase and
#         ObjId(e.child(2)) = fBase)],
#      e -> @(2).val
#  ),

  # id(a x b) = id(c) where c is derived from merging a and b and
  # computing corresponding domains and ranges.
  FoldConstants := ARule(idTensor, [@(1, idConst), @(2, idConst)],
      e -> [idConst(@1.val.domain()*@2.val.domain(), @1.val.params[2]*@2.val.params[2])]
  )

));


text:
Class(ACubeRot_ZXY, ASingletonTag);
Class(ACubeRot_YZX, ASingletonTag);
Class(ACubeRot_XYZ, ASingletonTag);
Class(ACubeRot_ZYX, ASingletonTag);

Class(TTensorII, Tagged_tSPL_Container, rec(
    abbrevs :=  [ (nt, s, l, r) -> Checked(
        IsSPL(nt),
	[nt, s, l, r])
    ],

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

    SPLtSPL := (self, nt, P) >> Error("not implemented"),

    terminate := self >> Error("not implemented"),

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

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

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

    doNotMeasure := true,

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

max_stars_repo_name: spiral-software/spiral-package-mpi

max_issues_repo_name: spiral-software/spiral-package-mpi

max_forks_repo_name: spiral-software/spiral-package-mpi


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

# Option Records
# ==============

#F SPL Options Records
#F -------------------

#F An spl options record is used to collect options passed to the
#F spl compiler by functions that communicate with the spl compiler,
#F such as measure, verify, and search functions.
#F
#F There is a system wide options record, SpiralDefailts (set in config.g),
#F which sets the default values.
#F
#F An spl options record is a record, which contains a subset of the
#F following fields. Adding a new fields requires to change the
#F functions PrintSpecSPLOptionsRecord, CheckSPLOptionsRecord, and
#F MakeSPLCallSPLOptionsRecord.
#F
#F rec(
#F   customDataType  = "int_cplx" | "int_fpN" | "Ipp16sc" (XScale) | "Ipp32sc" (XScale)
#F                                     where N - number of fractional bits, i.e. int_fp8
#F   zeroBits        = zero | <positive int>
#F   dataType        = "real" | "complex",
#F   precision       = "single" | "double" | "extended"
#F   subName         = <string>
#F   schedule        = <integer>
#F   globalUnrolling = <positive int> | "none" | "full",
#F   language        = "fortran" | "c"    # see config.g for languages
#F   compiler        = not to be used
#F   compflags       = <flags for compiler as string>
#F   splflags        = <flags for spl compiler as string>
#F )
#F
#F Note: switching language automatically switches compilers
#F and flags.
#F
#F spl options records should be used as follows:
#F - create your desired spl options record R
#F - merge with defaults, R1 := MergeSPLOptionsRecord(R)
#F - create spl prog for external operations with ProgSPL(SPL, R)
#F

#F CheckSPLOptionsRecord ( <spl-options-record> )
#F   checks whether <spl-options-record> is a valid spl options record
#F   with valid spl options set. If a field name or a field value
#F   is invalid, then an error is signaled, otherwise true is
#F   returned.
#F
CheckSPLOptionsRecord := function ( R )
  local r;
  if not IsRec(R) then Error("<R> must be an spl options record"); fi;

  for r in RecFields(R) do
    Cond(r = "dataType",
        Constraint(R.dataType in ["no default", "real", "complex"]),
         r = "customDataType",
        Constraint(IsString(R.customDataType)),
         r = "customReal",
        Constraint(IsString(R.customReal)),
         r = "customComplex",
        Constraint(IsString(R.customComplex)),
     r = "zeroBits",
        Constraint(IsInt(R.zeroBits) and R.zeroBits >= 0),
         r = "precision",
            Constraint(R.precision in  ["single", "double", "extended"]),
         r = "subName",
        Constraint(IsString(R.subName)),
         r = "file",
        Constraint(IsString(R.file)),
         r = "schedule",
        Constraint(IsInt(R.schedule) and R.schedule > 0),
         r = "globalUnrolling",
        Constraint(
        (IsInt(R.globalUnrolling) and R.globalUnrolling > 0)
        or R.globalUnrolling in ["none", "full"]),
         r = "compiler",
        Constraint(IsString(R.compiler)),
         r = "compflags",
        Constraint(IsString(R.compflags)),
         r = "dmpcompflags",
        Constraint(IsString(R.dmpcompflags)),
       r = "faultTolerant",
        Constraint(IsBool(R.faultTolerant)),
     r = "cgen",
        Constraint(IsFunc(R.cgen)),
     r = "x" or IsSystemRecField(r),
        Ignore(),
     0); # do nothing if field is unrecognized
  od;

  return true;
end;

#F MergeSPLOptionsRecord ( <spl-options-record> )
#F   returns the option record obtained by starting with SpiralDefaults
#F   (config.g, set at installation) and merging or overwriting with the
#F   spl options given by <spl-options-record>
#F
MergeSPLOptionsRecord := R -> Checked(CheckSPLOptionsRecord(R),
    CopyFields(SpiralDefaults, R)
);