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

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

Class(HockneyMLC_SIMTRules, RuleSet);

# Set of rules for applying Tile and Fission SIMT flags
RewriteRules(HockneyMLC_SIMTRules, rec(
	#Enable synchronization when different thread clusters are mapped to innermost sums
	# SUM(A*[S,IS], B*[S,IS]) => SUM(A, B)*SUM([S,IS],[S,IS])
    collect := Rule([SIMTSUM, [@(1, Compose),
                                @(2).cond(e -> not ObjId(e) in [Scat]), @(3).cond(s -> IsBound(s.simt_dim))
                            ],
                            [@(6, Compose),
                                @(7).cond(e -> not ObjId(e) in [Scat]), @(8).cond(s -> IsBound(s.simt_dim))
                            ]
                    ], e -> let(dsum1 := @(3).val.dimensions[1], dsum2 := @(8).val.dimensions[1],
                                Compose(    SIMTSUM(Copy(e.simt_dim),   @(2).val*Gath(H(dsum1+dsum2, dsum1, 0, 1)),
                                                                    @(7).val*Gath(H(dsum1+dsum2, dsum2, dsum1, 1)) ),
                                            SIMTSUM(Copy(e.simt_dim),   Scat(H(dsum1+dsum2, dsum1, 0, 1))*@(3).val,
                                                                    Scat(H(dsum1+dsum2, dsum2, dsum1, 1))*@(8).val )
                                    )
                            ),
                        "CollectCompose"
    ),

    tile    := Rule(@(1, SIMTISum, e->ObjId(e.simt_dim)=ASIMTTileFlag ),
                    e -> let(d1 := e.domain/e.simt_dim.tsize(),
                             v1 := var.fresh(e.var.id{[1]}, e.var.t, d1),
                             d2 := e.simt_dim.tsize(),
                             v2 := var.fresh(e.var.id{[1]}, e.var.t, d2),
                             newb := SubstVars(e.children()[1], rec( (e.var.id) := v1*d2 + v2 ) ),
                             simt_dlst := e.simt_dim.simt_dim(),
                             SIMTISum(Copy(simt_dlst.simt_dim(1)), v1, d1,
                                SIMTISum(Copy(simt_dlst.simt_dim(2)), v2, d2, newb )
                                )
                        ),
                    "tile_SIMTISum"
    ),

    fission_SIMTISum := Rule([@(1, SIMTISum, e->ObjId(e.simt_dim)=ASIMTFissionFlag ),
                            @(2, Compose, e-> not ForAny(e.children(), cch -> ObjId(cch) in [Gath, Scat]) )
                        ], function(e)
                            local inner_comp_list, ch, cch, cchlen, pos, simt_dlst, b, newb, v, sl, gl;

                            ch := @(2).val;
                            inner_comp_list := [];
                            cchlen :=  Length(ch.children());
                            for cch in ch.children() do
                                b := Copy( cch );
                                v := e.var.clone();
                                newb := SubstVars(b, rec( (e.var.id) := v ));
                                sl := newb.dimensions[1];
                                gl := newb.dimensions[2];
                                pos := Length(inner_comp_list);
                                simt_dlst := e.simt_dim.simt_dim();
                                Add(inner_comp_list,
                                    Cond(
                                        pos = 0,        SIMTISum(Copy(simt_dlst.simt_dim(pos+1)), v, e.domain, newb*Gath(H(e.domain*gl, gl, gl*v, 1)) ),
                                        pos = cchlen-1, SIMTISum(Copy(simt_dlst.simt_dim(pos+1)), v, e.domain, Scat(H(e.domain*sl, sl, sl*v, 1))*newb ),
                                        SIMTISum(Copy(simt_dlst.simt_dim(pos+1)), v, e.domain, Scat(H(e.domain*sl, sl, sl*v, 1))*newb*Gath(H(e.domain*gl, gl, gl*v, 1)) )
                                    )
                                );
                            od;
                            e := Compose( inner_comp_list );
                            return e;
                        end,
                        "fission_SIMTISum"
    ),

    # If outer loop is mapped to cluster (warp) dimension and it contains other block-parallel-loops then loop must be fissioned
    # to enable synchronization.
    xy_fission := Rule([@(1, SIMTISum, e->ObjId(e.simt_dim)=ASIMTBlockDimX ),
                            @(2, Compose, e->   not ForAny(e.children(), cch -> ObjId(cch) in [Gath, Scat]) and
                                                ForAny(List(Collect(e, @(3,[SIMTISum, SIMTSUM]) ), s -> s.simt_dim.simt_sync()),
                                                        s -> s > @(1).val.simt_dim.simt_sync())
                            )
                        ], function(e)
                            local   inner_comp_list, sync_list,
                                    ch, cch, cchlen, b, newb, v, sl, gl;

                            inner_comp_list := Collect(e.children(),
                                                    @@(4, Compose,
                                                        (expr, cx) -> IsBound(cx.SIMTSUM) and cx.SIMTSUM <> []
                                                                    or
                                                                    IsBound(cx.SIMTISum) and cx.SIMTISum <> []
                                                                    or
                                                                    IsBound(cx.Compose) and cx.Compose <> []
                                                    )
                                                );
                            sync_list := List(Collect(inner_comp_list, @(5,[SIMTISum, SIMTSUM]) ), s -> s.simt_dim.simt_sync());
                            if ForAny(sync_list, s -> s > e.simt_dim.simt_sync()) then
                                Error("Cannot properly syncronize within following structure: ", e);
                            fi;

                            ch := @(2).val;
                            inner_comp_list := [];
                            cchlen :=  Length(ch.children());
                            for cch in ch.children() do
                                b := Copy( cch );
                                v := e.var.clone();
                                newb := SubstVars(b, rec( (e.var.id) := v ));
                                sl := newb.dimensions[1];
                                gl := newb.dimensions[2];
                                Add(inner_comp_list,
                                    Cond(
                                        Length(inner_comp_list) = 0,        SIMTISum(Copy(e.simt_dim), v, e.domain, newb*Gath(H(e.domain*gl, gl, gl*v, 1)) ),
                                        Length(inner_comp_list) = cchlen-1, SIMTISum(Copy(e.simt_dim), v, e.domain, Scat(H(e.domain*sl, sl, sl*v, 1))*newb ),
                                        SIMTISum(Copy(e.simt_dim), v, e.domain, Scat(H(e.domain*sl, sl, sl*v, 1))*newb*Gath(H(e.domain*gl, gl, gl*v, 1)) )
                                    )
                                );
                            od;
                            e := Compose( inner_comp_list );
                            return e;
                        end,
                        "xy_fission_SIMTISum"
    ),
));

HockneyMLC_fftxGen := function(t, opts)
    local sb, cb, rtb, symbl, tt, tb,
    nd3, nfreq, n2, ns3, nfreq, ns2,
    tag_prdft_pd_loop1, rtp_prdft_pd_loop1, tag_dft_pd_loop1, rtp_dft_pd_loop1, tag_prdft_pd_loop2, rtp_prdft_pd_loop2, tag_dft_pd_loop2, rtp_dft_pd_loop2,
    tag_prdft_rec1, rtp_prdft_rec1, tag_prdft_rec2, rtp_prdft_rec2, tag_prdft_rec_z, rtp_prdft_rec_z, tag_dft_rec_y, rtp_dft_rec_y, BX, tag_idft_rec_x, rtp_idft_rec_x,
    tag_dft_rec_x, rtp_dft_rec_x, tag_conv_x, rtp_conv_x, tag_pridft_rec_z, rtp_pridft_rec_z, tag_idft_rec_y, rtp_idft_rec_y, TS_stg1, TS_stg2, TS_stg543, trstg5to1;

    # promote t
    tt := opts.preProcess(t);

    # parse parameters
    symbl := tt.params[2].params[1];
    tb := tt.params[1];
    nd3 := Length(tt.params[1].params[4][1]);
    nfreq := tt.params[1].params[1][1]/2+1;
    n2 := tt.params[1].params[1][1];
    ns3 := Length(tt.params[1].params[6][1]);
    ns2 := Length(tt.params[1].params[6][1]);
    opts.symbol := [symbl];

	if IsBound(tt.params[2].fname) then
	   opts.cudasubName := tt.params[2].fname;
	fi;

    # code gen script
    rtb := RuleTreeMid(tb, opts);

    ## Pruned Real DFT vvvvv
    # rec1: PRDFT_PD_loop presents 4 parallelizable objects
    tag_prdft_pd_loop1 := [ASIMTBlockDimY(32,[0,0],[0,31]), ASIMTLoopDim(), ASIMTBlockDimY(32, [1,6]), ASIMTLoopDim()];
    # rec1: PRDFT_PD_loop produces 0 subrules
    rtp_prdft_pd_loop1 := tag_rec(tag_prdft_pd_loop1, Ignore);

    # rec1: DFT_PD_loop presents 3 parallelizable objects
    tag_dft_pd_loop1 := [ASIMTBlockDimY(32,[7,7],[0,31]), ASIMTBlockDimY(32, [8,14]), ASIMTLoopDim()];
    # rec1: DFT_PD_loop produces 0 subrules
    rtp_dft_pd_loop1 := tag_rec(tag_dft_pd_loop1, Ignore);

    # rec2: PRDFT_PD_loop presents 4 parallelizable objects
    tag_prdft_pd_loop2 := [ASIMTBlockDimY(32,[15,15],[0,31]), ASIMTLoopDim(), ASIMTBlockDimY(32, [16,21]), ASIMTLoopDim()];
    # rec2: PRDFT_PD_loop produces 0 subrules
    rtp_prdft_pd_loop2 := tag_rec(tag_prdft_pd_loop2, Ignore);

    # rec2: DFT_PD_loop presents 3 parallelizable objects
    tag_dft_pd_loop2 := [ASIMTBlockDimY(32,[22,22],[0,31]), ASIMTBlockDimY(32, [23,29]), ASIMTLoopDim()];
    # rec2: DFT_PD_loop produces 0 subrules
    rtp_dft_pd_loop2 := tag_rec(tag_dft_pd_loop2, Ignore);

    tag_prdft_rec1 := [ASIMTBlockDimY(32), ASIMTLoopDim(), ASIMTBlockDimY(32), ASIMTBlockDimY(32, [0,3]), ASIMTBlockDimY(32, [4,12])];
    rtp_prdft_rec1 := tag_rec(tag_prdft_rec1, [rtp_prdft_pd_loop1, rtp_dft_pd_loop1] :: Replicate(2, Ignore));

    tag_prdft_rec2 := [ASIMTBlockDimY(32), ASIMTLoopDim(), ASIMTBlockDimY(32), ASIMTBlockDimY(32, [13,15]), ASIMTBlockDimY(32, [16,25])];
    rtp_prdft_rec2 := tag_rec(tag_prdft_rec2, [rtp_prdft_pd_loop2, rtp_dft_pd_loop2] :: Replicate(2, Ignore));

    #Top prdft rec rule presents 3 parallelizable objects
    tag_prdft_rec_z := [ASIMTBlockDimX(16,[0,0],[0,15]), ASIMTBlockDimX(16), ASIMTBlockDimY(32)];
    #Top prdft rec rule produces 4 subrules
    rtp_prdft_rec_z := tag_rec(tag_prdft_rec_z, [Ignore, Ignore, rtp_prdft_rec1, rtp_prdft_rec2]);

    ## Pruned Real DFT ^^^^^

    ## Pruned DFT vvvvv

    #Top dft rec rule presents 4 parallelizable objects
    tag_dft_rec_y := [ASIMTBlockDimY(13, [0,9]), ASIMTBlockDimY(13), ASIMTBlockDimY(13, [0,6]), ASIMTBlockDimY(13, [7,12])];
    #Top dft rec rule produces 3 subrules
    rtp_dft_rec_y := tag_rec(tag_dft_rec_y, Replicate(3, Ignore));

    ## Pruned DFT ^^^^^

    ## Pruned Conv vvvvv
    BX := 16;
    #Top idft rec rule presents 4 parallelizable objects
    tag_idft_rec_x := [ASIMTBlockDimX(BX), ASIMTBlockDimX(BX, [0,7]), ASIMTBlockDimX(BX, [8,12]), ASIMTBlockDimX(BX, [0,9])];
    #Top idft rec rule produces 3 subrules
    rtp_idft_rec_x := tag_rec(tag_idft_rec_x, Replicate(3, Ignore));

    #Top dft rec rule presents 4 parallelizable objects
    tag_dft_rec_x := [ASIMTBlockDimX(BX, [0,9]), ASIMTBlockDimX(BX), ASIMTBlockDimX(BX, [0,6]), ASIMTBlockDimX(BX, [7,12])];
    #Top dft rec rule produces 3 subrules
    rtp_dft_rec_x := tag_rec(tag_dft_rec_x, Replicate(3, Ignore));

    #Top IOPrunedConv rule presents no parallelizable objects
    tag_conv_x := Ignore;
    #Top IOPrunedConv rule produces 2 subrules
    rtp_conv_x := tag_rec(tag_conv_x, [rtp_idft_rec_x, rtp_dft_rec_x]);

    ## Pruned IDFT YZ vvvvv

    BX := 32;
    #Z-dim pridft rec rule presents 3 parallelizable objects
    tag_pridft_rec_z := [ASIMTBlockDimX(BX, [0,1]), ASIMTBlockDimX(BX, [0,0], [0,31]), ASIMTBlockDimX(BX)];
    #Z-dim pridft rec rule produces 3 subrules
    rtp_pridft_rec_z := tag_rec(tag_pridft_rec_z, Replicate(3, Ignore));

    #Y-dim idft rec rule presents 4 parallelizable objects
    tag_idft_rec_y := [ASIMTBlockDimX(BX), ASIMTBlockDimX(BX, [0,7]), ASIMTBlockDimX(BX, [8,12]), ASIMTBlockDimX(BX, [0,9])];
    #Y-dim idft rec rule produces 3 subrules
    rtp_idft_rec_y := tag_rec(tag_idft_rec_y, Replicate(3, Ignore));

    ## Pruned IDFT YZ ^^^^^

    TS_stg1 := 11;
    TS_stg2 := 11;
    TS_stg543 := 4;

    # Define TagTree
    trstg5to1 := tag_rec(
                        # Taggable operators (Tensors, Sums, etc.) as they appear in Collect( IOPrunedMDRConv.apply(...), @([Tag. ops]) )
                        [ ASIMTLoopDim(), ASIMTKernelFlag( ASIMTGridDimX(nd3) ), ASIMTBlockDimY(6), ASIMTBlockDimY(6), # Stages 7-6: 4 tag operators
                            ASIMTLoopDim(), ASIMTLoopDim(), ASIMTTileFlag(
                                                            ASIMTListDim(
                                                                ASIMTKernelFlag(ASIMTGridDimX(nfreq*n2/TS_stg543)),
                                                                ASIMTBlockDimY(TS_stg543)
                                                            ), TS_stg543), ASIMTLoopDim(), ASIMTLoopDim(), # Stages 5-4-3: 5 tag operators
                            ASIMTKernelFlag(ASIMTGridDimX(ns3)), ASIMTTileFlag(
                                                                ASIMTListDim(
                                                                    ASIMTGridDimY(nfreq/TS_stg2),
                                                                    ASIMTBlockDimX(11) # << 11 OK because perfectly nests Inner loop mapped to Y (so block sync is used at compose boundaries)
                                                                ), TS_stg2), # Stage 2: 2 tag operators
                            ASIMTKernelFlag( ASIMTGridDimX(ns2) ),
                                    ASIMTTileFlag(
                                            ASIMTListDim(
                                                ASIMTGridDimY(ns3/TS_stg1),
                                                ASIMTFissionFlag( ASIMTListDim(ASIMTBlockDimY(32), ASIMTBlockDimX(16)) )
                                            ), TS_stg1), ASIMTLoopDim() ], # Stage 1: 3 tag operators
                        # Children as they appear in IOPrunedMDRConv.children(...)
                        [ rtp_prdft_rec_z, rtp_dft_rec_y, rtp_conv_x, rtp_idft_rec_y, rtp_pridft_rec_z, Ignore ]);

    tmp_tag_rt(rtb, trstg5to1, opts); # RuleTree + TagTree = TaggedRuleTree

    sb := SumsRuleTree(rtb, opts);
    cb := CodeSums(sb, opts);
    cb.transform := t;

    return cb;
end;

Class(HockneyMlcCUDADeviceOpts, TitanVDefaults, rec(
    operations := rec(Print := s -> Print("<Hockney MLC CUDA Device options record>")),
    tagIt := (self, t) >> t.withTags(self.tags)
));

hockneyMlcCUDADeviceOpts := function(arg)
    local opts;

    opts := Copy(HockneyMlcCUDADeviceOpts);
    opts.includes := [];

    opts.breakdownRules.PRDFT := List([ PRDFT1_Base2, CopyFields(PRDFT_PD, rec(maxSize := 7)), PRDFT_PD_loop], _noT);
    opts.breakdownRules.IPRDFT := List([ IPRDFT1_Base1, IPRDFT1_Base2, IPRDFT_PD_loop, CopyFields(IPRDFT_PD, rec(maxSize := 7))], _noT);
    opts.breakdownRules.IPRDFT2 := List([ IPRDFT2_Base1, IPRDFT2_Base2, IPRDFT2_CT], _noT);
    opts.breakdownRules.PRDFT3 := [ ];
    opts.breakdownRules.URDFT := List([ URDFT1_Base1, URDFT1_Base2, URDFT1_Base4, URDFT1_CT ], _noT);
    opts.breakdownRules.DFT := List([ DFT_Base, CopyFields(DFT_PD, rec(maxSize := 7)),  DFT_PD_loop], _noT);
    opts.breakdownRules.PrunedPRDFT := List([ CopyFields(PrunedPRDFT_base, rec(maxSize := 5)), PrunedPRDFT_CT_rec_block ], _noT);
    opts.breakdownRules.PrunedIPRDFT := List([ CopyFields(PrunedIPRDFT_base, rec(maxSize := 5)), PrunedIPRDFT_CT_rec_block ], _noT);
    opts.breakdownRules.PrunedDFT := List([ PrunedDFT_base, PrunedDFT_CT_rec_block ], _noT);
    opts.breakdownRules.IOPrunedMDRConv := List([ IOPrunedMDRConv_3D_2trip_zyx_freqdata ], _noT);
    opts.breakdownRules.MDRConv := [ MDRConv_3D_2trip_zyx_freqdata ]; # [MDRConv_Base]

    wrap_rule_apply(opts.breakdownRules.PRDFT, wrap_apply);
    wrap_rule_apply(opts.breakdownRules.IPRDFT, wrap_apply);
    wrap_rule_apply(opts.breakdownRules.IPRDFT2, wrap_apply);
    wrap_rule_apply(opts.breakdownRules.URDFT, wrap_apply);
    wrap_rule_apply(opts.breakdownRules.DFT, wrap_apply);
    wrap_rule_apply(opts.breakdownRules.PrunedPRDFT, wrap_apply);
    wrap_rule_apply(opts.breakdownRules.PrunedIPRDFT, wrap_apply);
    wrap_rule_apply(opts.breakdownRules.PrunedDFT, wrap_apply);
    wrap_rule_apply(opts.breakdownRules.PrunedIDFT, wrap_apply);
    wrap_rule_apply(opts.breakdownRules.TTensorI{[2..3]}, wrap_apply);
    wrap_rule_apply(opts.breakdownRules.IOPrunedConv, wrap_apply);
    wrap_rule_apply(opts.breakdownRules.IOPrunedMDRConv, wrap_apply);

    opts.codegen._Formula := Copy(opts.codegen.Formula);
    opts.codegen.Formula := meth ( self, o, y, x, opts )
        local prog, prog1, s;
        s := Copy(o);
        s := fixUpSigmaSPL(s, opts);

        opts.("postCompileStrategy") := function(icode, opts)
            local icode1;
            icode1 := skip();
            # This reduces the amount of constant matrices but might prevent to store
            # some of them in constant mem. E.g., two identical matrices but only one has broadcast-type access.
            # while icode1 <> icode do
            #     icode1 := icode;
            #     icode := fixReplicatedData(icode, opts);
            # od;
            icode := scalarizeAccumulator(icode, opts);
            return icode;
        end;

        prog := opts.codegen._Formula(s, y, x, opts);

        return prog;
    end;

    opts.formulaStrategies.postProcess := opts.formulaStrategies.postProcess
        :: opts.formulaStrategies.sigmaSpl :: [MergedRuleSet(HockneyMLC_SIMTRules, RulesSums)];

    opts.fftxGen := (self, t) >> HockneyMLC_fftxGen(t, self);
    opts.tags := [];

    opts.preProcess := (self, t) >> let(t1 := RulesFFTXPromoteNT(Copy(t)), RulesFFTXPromoteNT_Cleanup(t1));
    opts.prettyPrint := meth(self, c)
        local name;
        name := "fftx_generated";
        if (IsBound(c.transform) and IsBound(c.transform.params[2].fname)) then
            name := c.transform.params[2].fname;
        fi;
        PrintCode(name, c, self);
    end;

    return opts;
end;

confHockneyMlcCUDADevice := rec(
    defaultName := "confHockneyMlcCUDADevice",
    defaultOpts := (arg) >> rec(useHockneyMlcCUDADevice := true),
    confHandler := hockneyMlcCUDADeviceOpts
);

fftx.FFTXGlobals.registerConf(confHockneyMlcCUDADevice);

max_stars_repo_name: franzfranchetti/spiral-package-fftx

max_issues_repo_name: franzfranchetti/spiral-package-fftx

max_forks_repo_name: franzfranchetti/spiral-package-fftx


text:
#
# JupyterZMQ: UUID function
#

#! @Description
#!   Create a new zero UUID
InstallGlobalFunction(NewUUID,
function()
    local uuid;
    uuid := rec( bits := BlistList([1..128], []));
    return Objectify(UUIDType, uuid);
end);

#! @Description
#!   Generate a random UUID according to RFC4122
InstallGlobalFunction(RandomUUID,
function()
    local puuid, uuid;

    puuid := HexStringInt(Random(0,2^128-1));
    if Length(puuid) < 32 then
        puuid := Concatenation(RepeatedString("0", 32 - Length(puuid)), puuid);
    fi;
    uuid := BlistStringDecode(puuid);
    # Set version to 4
    uuid{[49..52]} := [false, true, false, false];
    # Set variant to RFC4122
    uuid{[70..72]} := [true, false, false];

    return Objectify( UUIDType, rec( bits := uuid ) );
end);

InstallGlobalFunction(StringUUID,
function(uuid)
    local hex;
    hex := LowercaseString(HexStringBlist(uuid!.bits));
    return JoinStringsWithSeparator([
                hex{[1..8]}, hex{[9..12]}, hex{[13..16]},
                hex{[17..20]}, hex{[20..32]}],
                "-");
end);

InstallGlobalFunction(HexStringUUID,
uuid -> HexStringBlist(uuid!.bits));

InstallMethod(ViewString, "for UUID", [IsUUID and IsUUIDBlistRep],
uuid -> Concatenation("<uuid ", String(uuid), ">"));

InstallMethod(String, "for UUID", [IsUUID and IsUUIDBlistRep],
    StringUUID);

ext: gi

max_stars_repo_name: mtorpey/uuid

max_issues_repo_name: mtorpey/uuid

max_forks_repo_name: mtorpey/uuid


text:
InstallMethod( JupyterRender, "default fallback"
             , [ IsObject ],
function(obj)
    local str;
    str := ViewString(obj);
    RemoveCharacters(str, "\<\>\n");
    return Objectify( JupyterRenderableType
                    , rec( data := rec( text\/plain := str )
                         , metadata := rec( text\/plain := "") ) );
end);

InstallMethod( JupyterRender, "default fallback"
               , [ IsJupyterRenderableRep ],
               IdFunc);

InstallMethod( JupyterRenderableData, "for a JupyterRenderable"
               , [  IsJupyterRenderableRep ]
               , x -> x!.data );

InstallMethod( JupyterRenderableMetadata, "for a JupyterRenderable"
               , [  IsJupyterRenderableRep ]
               , x -> x!.metadata );

InstallMethod( ViewString, "for a JupyterRenderable"
               , [  IsJupyterRenderableRep ]
               , x -> "<jupyter renderable>" );

InstallMethod( JupyterRenderable, "for a record and a record"
               , [ IsObject, IsObject ],
function(data, metadata)
    return Objectify( JupyterRenderableType
                    , rec( data := data, metadata := metadata ) );
end);

max_stars_repo_name: isuruf/GapJupyterKernel

max_issues_repo_name: isuruf/GapJupyterKernel

max_forks_repo_name: isuruf/GapJupyterKernel


text:
InstallMethod(ParabolicSystem,
              "Roots of standard parabolic",
              [IsRootSystem,IsList],
              function(R,I)
              local P,graph,order,i_levels;
              P:=Objectify(NewType(NewFamily("ParabolicSystemFamily"),
                                   IsAttributeStoringRep and
                                   IsParabolicSystem),
                           rec());
              SetUnderlyingRootSystem(P,R);
              SetSupportRoots(P,I);
              SetRadicalRoots(P,Filtered(positiveRoots(R),i->InRadical(P,i)));
              SetiRadicalRoots(P,List(RadicalRoots(P),i->Position(positiveRoots(R),i)));
              SetHeight(P,Maximum(List([1..Length(positiveRoots(R))],i->RootHeight(P,i))));
              SetLevels(P,List([1..Height(P)],i->Filtered(positiveRoots(R),j->RootHeight(P,j)=i)));
              i_levels:=List([1..Height(P)],i->Filtered([1..Length(positiveRoots(R))],j->RootHeight(P,positiveRoots(R)[j])=i));
              SetiLevels(P,i_levels);
              i_levels:=Concatenation(i_levels);
              Set0Level(P,Filtered([1..Length(positiveRoots(R))],i->not i in i_levels));
              SetShapes(P,Set(List(positiveRoots(R),i->Shape(P,i))));
              SetLeviModules(P,List(Levels(P),i->Set(List(i,j->LeviModule(P,j)))));
              SetiLeviModules(P,List(iLevels(P),i->Set(List(i,j->iLeviModule(P,j)))));
              SetName(P,Concatenation("<parabolic of ",
                                      SemiSimpleType(UnderlyingLieAlgebra(R)),
                                      " root system with support in ",
                                      String(I),">"));

              order:=function(P)
                  local levels,result,lung,i;
                  levels:=Concatenation(iLevels(P));
                  lung:=Length(levels);
                  result:=List([1..lung],i->0);
                  for i in [1..lung] do
                      result[levels[i]]:=i;
                  od;
                  return result;
              end;

              SetOrdering(P,order(P));

#              SetRootGraph(P,graph(P));

              return P;
end);

InstallMethod(Shape,
              "Shape of root",
              [IsParabolicSystem,IsList],
              function(P,root_coeffs)
              local I;
              if not root_coeffs in positiveRoots(UnderlyingRootSystem(P)) then Error("Vector not in root system. "); fi;
              I:=SupportRoots(P);
              return root_coeffs{Filtered([1..Length(root_coeffs)],i->not i in I)};
end);

InstallMethod(Shape,
              "Shape of root",
              [IsParabolicSystem,IsPosInt],
              function(P,root_index)
              if not root_index in [1..Length(positiveRoots(UnderlyingRootSystem(P)))] then Error("Not a root index. "); fi;
              return Shape(P,positiveRoots(UnderlyingRootSystem(P))[root_index]);
end);

InstallMethod(LeviModule,
              "Levi module (vectors) for root vector",
              [IsParabolicSystem,IsList],
              function(P,root_coeffs)
              local I,s;
              if not root_coeffs in positiveRoots(UnderlyingRootSystem(P)) then Error("Vector not in root system. "); fi;
              I:=SupportRoots(P);
              s:=Shape(P,root_coeffs);
              return Filtered(positiveRoots(UnderlyingRootSystem(P)),i->Shape(P,i)=s);
end);

InstallMethod(LeviModule,
              "Levi module (vectors) for root index",
              [IsParabolicSystem,IsPosInt],
              function(P,root_index)
              if not root_index in [1..Length(positiveRoots(UnderlyingRootSystem(P)))] then Error("Not a root index. "); fi;
              return LeviModule(P,positiveRoots(UnderlyingRootSystem(P)[root_index]));
end);

InstallMethod(iLeviModule,
              "Levi module (vectors) for root vector",
              [IsParabolicSystem,IsList],
              function(P,root_coeffs)
              local I,s;
              if not root_coeffs in positiveRoots(UnderlyingRootSystem(P)) then Error("Vector not in root system. "); fi;
              I:=SupportRoots(P);
              s:=Shape(P,root_coeffs);
              return Filtered([1..Length(positiveRoots(UnderlyingRootSystem(P)))],i->Shape(P,i)=s);
end);

InstallMethod(iLeviModule,
              "Levi module (indices) for root vector",
              [IsParabolicSystem,IsPosInt],
              function(P,root_index)
              if not root_index in [1..Length(positiveRoots(UnderlyingRootSystem(P)))] then Error("Not a root index. "); fi;
              return iLeviModule(P,positiveRoots(UnderlyingRootSystem(P))[root_index]);
end);

InstallMethod(LeviModules,
              "Levi modules at a given height",
              [IsParabolicSystem,IsPosInt],
              function(P,i)
              if i > Height(P) then return []; fi;
              return Filtered(LeviModules(P),j->RootHeight(P,j[1])=i);
end);

InstallMethod(iLeviModules,
              "Levi modules at a given height",
              [IsParabolicSystem,IsPosInt],
              function(P,i)
              if i > Height(P) then return []; fi;
              return Filtered(iLeviModules(P),j->RootHeight(P,j[1])=i);
end);

InstallMethod(RootHeight,
              "the level of a root",
              [IsParabolicSystem,IsList],
              function(P,root_coeffs)
              return Sum(Shape(P,root_coeffs));
end);

InstallMethod(RootHeight,
              "the level of a root",
              [IsParabolicSystem,IsPosInt],
              function(P,root_index)
              return Sum(Shape(P,root_index));
end);

InstallMethod(Level,
              "roots at level",
              [IsParabolicSystem,IsPosInt],
              function(P,i)
              if i > Height(P) then return []; fi;
              return Levels(P)[i];
end);

InstallMethod(InRadical,
              "Check for root in radical",
              [IsParabolicSystem,IsList],
              function(P,root_coeffs)
              if not root_coeffs in positiveRoots(UnderlyingRootSystem(P)) then Error("Vector not in root system. "); fi;
              return Length(Filtered([1..Length(root_coeffs)],i->not i in SupportRoots(P) and root_coeffs[i]<>0))>0;
end);


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

ImportAll(spiral);

SkewRDFT := BSkewPRDFT;
BRDFT := PolyBDFT;

alphai := (n,i,a) -> Cond(i mod 2 = 0, (a + Int(i/2))/n, (1-a + Int(i/2))/n);

time_r := a -> [[1, -CosPi(2*a)/SinPi(2*a)],
                [0, 1/SinPi(2*a)]];

r_time := a -> [[1, CosPi(2*a)],
                [0, SinPi(2*a)]];

B1_time_r := n -> Cond(IsOddInt(n),
    DirectSum(I(1), List([1..(n-1)/2], i-> Mat(time_r(i/n)))),
    DirectSum(I(1), List([1..n/2-1], i-> Mat(time_r(i/n))), I(1)));

RDFT1 := SymFunc("RDFT1", n -> Cond(
	n=2,F(2),
	IsEvenInt(n), VStack(Gath(H(n+2,1,0,1)), Gath(H(n+2,n-2,2,1)), Gath(H(n+2,1,n,1))) * PRDFT1(n),
	IsOddInt(n), VStack(Gath(H(n+1,1,0,1)), Gath(H(n+1,n-1,2,1))) * PRDFT1(n)));

RDFT3 := SymFunc("RDFT3", n -> Cond(
	n=2,I(2),
	IsEvenInt(n), PRDFT3(n),
	IsOddInt(n), VStack(Gath(H(n+1,n-1,0,1)), Gath(H(n+1,1,n-1,1))) * PRDFT3(n)));

BRDFT1 := SymFunc("BRDFT1", n -> B1_time_r(n) * RDFT1(n));

rdft1 := (k,m) ->
   Perm(PermFunc(rP(2*k*m, 2*k), 2*k*m), 2*k*m).transpose() *
   DirectSum(RDFT1(m), List([1..k-1], i->Mat(MatrDFT(2*m,i/(2*k)))), RDFT3(m)) *
   Tensor(RDFT1(2*k), I(m));

rdft1odd := (k,m) ->
   Perm(PermFunc(rPhat((2*k-1)*m, 2*k-1), (2*k-1)*m), (2*k-1)*m).transpose() *
   DirectSum(RDFT1(m), List([1..k-1], i->Mat(MatrDFT(2*m,i/(2*k-1))))) *
   Tensor(RDFT1(2*k-1), I(m));

test_rdft1 := (k,m) -> PrintLine(inf_norm(MatSPL(rdft1(k,m)) - MatSPL(RDFT1(2*k*m))));
test_rdft1odd := (k,m) -> PrintLine(inf_norm(MatSPL(rdft1odd(k,m)) - MatSPL(RDFT1((2*k-1)*m))));
do_test_rdft1 := function()
   test_rdft1(4,4);
   test_rdft1(4,8);
   test_rdft1(6,8);
   test_rdft1(6,2);
   test_rdft1(6,5);
   test_rdft1(11,5);
end;
do_test_rdft1odd := function()
   test_rdft1odd(4,4);
   test_rdft1odd(4,8);
   test_rdft1odd(6,8);
   test_rdft1odd(6,2);
   test_rdft1odd(6,5);
   test_rdft1odd(5,5);
   test_rdft1odd(11,5);
end;

rdft3 := (k,m) ->
   RC(K(k*m, m)) *
   DirectSum(List([0..k-1], i->Mat(MatrDFT(2*m,(i+1/2)/(2*k))))) * #alphai(k,i,1/4))))) *
   Tensor(RDFT3(2*k),I(m));

rdft3odd := (k,m) ->
   Perm(PermFunc(rQhat((2*k-1)*m, 2*k-1), (2*k-1)*m), (2*k-1)*m).transpose() *
   DirectSum(List([0..k-2], i->Mat(MatrDFT(2*m,(i+1/2)/(2*k-1)))), RDFT3(m)) * #alphai(k,i,1/4))))) *
   Tensor(RDFT3(2*k-1),I(m));

test_rdft3 := (k,m) -> PrintLine(inf_norm(MatSPL(rdft3(k,m)) - MatSPL(RDFT3(2*k*m))));
test_rdft3odd := (k,m) -> PrintLine(inf_norm(MatSPL(rdft3odd(k,m)) - MatSPL(RDFT3((2*k-1)*m))));
do_test_rdft3 := function()
   test_rdft3(4,4);
   test_rdft3(4,8);
   test_rdft3(6,8);
   test_rdft3(6,2);
   test_rdft3(6,5);
   test_rdft3(11,5);
end;
do_test_rdft3odd := function()
   test_rdft3odd(4,4);
   test_rdft3odd(4,8);
   test_rdft3odd(6,8);
   test_rdft3odd(6,2);
   test_rdft3odd(6,5);
   test_rdft3odd(5,5);
   test_rdft3odd(11,5);
end;

rrdft := (k,m,a) ->
   RC(K(k*m, m)) *
   DirectSum(List([0..k-1], i->Mat(MatrDFT(2*m,(i+a)/(2*k))))) *
   Tensor(Mat(MatrDFT(2*k,a)),I(m));

brdft := (k,m,a) ->
   RC(K(k*m, m)) *
   DirectSum(List([0..k-1], i->BRDFT(2*m,alphai(k,i,a)))) *
   Tensor(BRDFT(2*k, a), I(m));

bruun_rdft3 := (k,m) ->
   RC(K(k*m, m)) *
   DirectSum(List([0..k-1], i->SkewRDFT(2*m,alphai(k,i,a)))) *
   Tensor(BRDFT(2*k,1/4),I(m));

bruun_rdft1 := (k,m) ->
   DirectSum(RDFT1(m), List([1..k-1], i->SkewRDFT(2*m,i/(2*k))), RDFT3(m)) *
   Tensor(BRDFT1(2*k), I(m));

rcred := mat -> List(mat{[0..Length(mat)/2-1]*2+1}, row -> row{[0..Length(row)/2-1]*2+1});


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

Class(FixCodeAVX, RuleSet);
RewriteRules(FixCodeAVX, rec(
    fix_noneExp := Rule( noneExp, e -> e.t.zero()),

    vpermf128_8x32f_to_vextract := Rule( [assign, [deref, @(1)], [vpermf128_8x32f, @(2), @(3), @(4)]],
        e -> let( p := @(4).val.p,
                  a := [[@(2).val, [0]], [@(2).val, [1]], [@(3).val, [0]], [@(3).val, [1]]],
                  dst := tcast(TPtr(TVect(T_Real(32), 4)), @(1).val),
                  chain(
                      assign( deref(dst  ), ApplyFunc(vextract_4l_8x32f, a[p[1]])),
                      assign( deref(dst+1), ApplyFunc(vextract_4l_8x32f, a[p[2]]))
                  )
        )),

    flatten_chains := ARule( chain, [@(1, chain)], e -> [@(1).val.cmds] ),

    addsub_4x64f_to_mul := Rule( [addsub_4x64f, _0, @(1)],
        e -> mul(e.t.value([-1,1,-1,1]), @(1).val)),
    addsub_8x32f_to_mul := Rule( [addsub_8x32f, _0, @(1)],
        e -> mul(e.t.value([-1,1,-1,1,-1,1,-1,1]), @(1).val)),

    # replacing 32 byte load by 16 byte load which is faster for out of L1 problems.
    deref_to_doubleload := Rule(
        @@(1, deref, (x, cx) -> _avxT(x.t, cx.opts) and cx.isInside(assign) and Last(cx.assign).loc<>x and
                                cx.isInside(rdepth_marker) and Last(cx.rdepth_marker).depth=1),
        (e, cx) -> Cond(
            e.t = TVect(T_Real(64), 4) or e.t = TVect(TReal, 4),
                let( p := tcast(TPtr(TVect(T_Real(64), 2)), e.loc),
                     vinsert_2l_4x64f( vinsert_2l_4x64f(e.t.zero(), deref(p), [0]), deref(p + 1), [1])),
            e.t = TVect(T_Real(32), 8) or e.t = TVect(TReal, 8),
                let( p := tcast(TPtr(TVect(T_Real(32), 4)), e.loc),
                     vinsert_4l_8x32f( vinsert_4l_8x32f(e.t.zero(), deref(p), [0]), deref(p + 1), [1])),
            Error("unexpected type"))
    ),

    vloadu_8x32f_to_doubleload := Rule( @@(1, vloadu_8x32f, (x, cx) -> cx.isInside(rdepth_marker) and Last(cx.rdepth_marker).depth=1),
        (e, cx) -> vinsert_4l_8x32f( vinsert_4l_8x32f(e.t.zero(), vloadu_4x32f(e.args[1]), [0]), vloadu_4x32f(e.args[1] + 4), [1])
    ),

    vloadu_4x64f_to_doubleload := Rule( @@(1, vloadu_4x64f, (x, cx) -> cx.isInside(rdepth_marker) and Last(cx.rdepth_marker).depth=1),
        (e, cx) -> vinsert_2l_4x64f( vinsert_2l_4x64f(e.t.zero(), vloadu_2x64f(e.args[1]), [0]), vloadu_2x64f(e.args[1] + 2), [1])
    ),

));

RewriteRules(RulesStrengthReduce, rec(
    # a dupe
    #tcast_fold := Rule([tcast, @(1), [tcast, @(2), @(3)]], e->tcast(@(1).val, @(3).val)),

    avx_add_addsub_vzero := Rule([add, @(1), [@(2, [addsub_4x64f, addsub_8x32f]), _0, @(3)]],
        e -> ObjId(@(2).val)(@(1).val, @(3).val)),

    avx_addsub_vzero := Rule([@(1, [addsub_4x64f, addsub_8x32f]), @(2), _0], e -> @(2).val),

    # is it valid?
    avx_addsub_noneExp := Rule([@(1, [addsub_4x64f, addsub_8x32f]), @(2, noneExp), @(3)],
	e -> ObjId(@(1).val)(@(2).val.t.zero(), @(3).val)),

    avx_addsub_mul_mul_values := Rule([@(0, [add, sub, addsub_4x64f, addsub_8x32f]), [mul, @(1, Value), @(2)], [mul, @@(3, Value, (x, cx)->x=@(1).val and _avxT(@(0).val.t, cx.opts)), @(4)]],
        e -> mul(@(1).val, ObjId(@(0).val)(@(2).val, @(4).val))),

    # replace by mul const?
    avx_addsub_zero_neg := Rule( [@(1, [addsub_4x64f, addsub_8x32f]), _0, [neg, @(2)]],
        e -> neg(ObjId(@(1).val)(e.t.zero(), @(2).val))),

    avx_vunpack_4x64f_neg := Rule([@(1, [vunpacklo_4x64f, vunpackhi_4x64f]), [neg, @(2)], [neg, @(3)]],
       e -> neg(ObjId(e)(@(2).val, @(3).val))),

    #NOTE: this is real hack, unpacking from vdup(TVect(TReal,2),2) is simply a vdup(TReal, 4), how to do this right?
    # need the same for AVX_8x32f
    avx_vunpacklo_4x64f_vdup_2 := Rule([vunpacklo_4x64f, [@(1, vdup), nth, _2],
           [@(2, vdup, x -> x=@(1).val), @(3, nth, x -> x.loc _is var and not IsBound(x.loc.value)), _2]],
       e -> let(n := @(3).val, vdup(deref(nth(n.loc, 0).toPtr(e.t.t) + 2*n.idx), 4))),
    avx_vunpackhi_4x64f_vdup_2 := Rule([vunpackhi_4x64f, [@(1, vdup), nth, _2],
           [@(2, vdup, x -> x=@(1).val), @(3, nth, x -> x.loc _is var and not IsBound(x.loc.value)), _2]],
       e -> let(n := @(3).val, vdup(deref(nth(n.loc, 0).toPtr(e.t.t) + 2*n.idx+1), 4))),

    avx_vperm_neg := Rule([@(0, [vperm_4x64f, vperm_8x32f]), [neg, @(1)], @(2)],
        e -> neg(ObjId(e)(@(1).val, @(2).val))),

    avx_vshuffle_neg := Rule([@(0, [vshuffle_4x64f, vshuffle_8x32f]), [neg, @(1)], [neg, @(2)], @(3)],
       e -> neg(ObjId(@(0).val)(@(1).val, @(2).val, @(3).val))),

    avx_addsub_neg_neg := Rule([@(0, [addsub_4x64f, addsub_8x32f]), [neg, @(1)], [neg, @(2)]],
        e -> neg(ObjId(e)(@(1).val, @(2).val))),

    avx_fold_vshuffle := Rule(
	[@(1, [vshuffle_4x64f, vshuffle_8x32f]), @(2, Value, e->IsVecT(e.t)),
	                                         @(3, Value, e->IsVecT(e.t)), @(4)],
        e -> Value(@(2).val.t, @(1).val.semantic(@(2).val.v, @(3).val.v, @(4).val.p))),

    # op(v1*a, v2*b) where v1*v2=0 turns into op(v1,v2)*vblend(a,b)
    avx_addsub_to_blend := Rule([@(1, [add, sub, addsub_4x64f, addsub_8x32f]), [mul, @(2, Value), @(3)], [mul, @@(4, Value, (x, cx) -> _avxT(x.t, cx.opts) and x*@(2).val=x.t.zero()), @(5)]],
        e -> let( c1 := @(2).val,
                  c2 := @@(4).val,
                  c  := ObjId(e)(c1, c2),
                  b  := Cond( e.t.size = 4, vblend_4x64f, vblend_8x32f ),
                  c*b(@(3).val, @(5).val, List(c1.v, x -> Cond(x=0, 2, 1))))),

    avx_addsub_ushuffles := Rule([@(0, [add, sub]),
            [@(1, [vshuffle_4x64f, vshuffle_8x32f]), @(2), @(3).cond(x -> x=@(2).val), @(4)],
            [@(5, [vshuffle_4x64f, vshuffle_8x32f]), @(6), @(7).cond(x -> x=@(6).val), @(8).cond(x -> x = @(4).val and ObjId(@(5).val)=ObjId(@(1).val))]],
        e -> ObjId(@(1).val)(ObjId(e)(@(2).val, @(6).val), ObjId(e)(@(2).val, @(6).val), @(4).val)),
    avx_addsub_vperm := Rule([@(0, [add, sub]),
            [vperm_4x64f, @(1), @(2)], [vperm_4x64f, @(3), @(4).cond(x -> x = @(2).val)]],
        e -> vperm_4x64f(ObjId(e)(@(1).val, @(3).val), @(2).val)),

    avx_fmaddsub_4x64f := Rule(
	[@@(1, addsub_4x64f, (e,cx)->IsBound(cx.opts.vector) and cx.opts.vector.SIMD = "AVX3"),
	    @(2, mul, e->Length(e.args)=2), @(3)],
        e-> fmaddsub_4x64f(@(2).val.args[1], @(2).val.args[2], @(3).val)),

    avx_fmaddsub_8x32f := Rule(
	[@@(1, addsub_8x32f, (e,cx)->IsBound(cx.opts.vector) and cx.opts.vector.SIMD = "AVX3"),
	    @(2, mul, e->Length(e.args)=2), @(3)],
        e-> fmaddsub_8x32f(@(2).val.args[1], @(2).val.args[2], @(3).val)),

    avx_add_mul_to_addsub := Rule([add, [mul, @(1, Value), @(2)],
                                        [mul, @@(3, Value, (x, cx) -> _avxT(x.t, cx.opts) and
                                                ForAll([1..x.t.size], i -> x.v[i]*(1 - 2*(i mod 2))=@(1).val.v[i])), @(4)]],
        e -> mul(@(1).val, Cond(e.t.size=4, addsub_4x64f, addsub_8x32f)(@(2).val, @(4).val))),
    avx_sub_mul_to_addsub := Rule([sub, [mul, @(1, Value), @(2)],
                                        [mul, @@(3, Value, (x, cx) -> _avxT(x.t, cx.opts) and
                                                ForAll([1..x.t.size], i -> x.v[i]*(-1 + 2*(i mod 2))=@(1).val.v[i])), @(4)]],
        e -> mul(@(1).val, Cond(e.t.size=4, addsub_4x64f, addsub_8x32f)(@(2).val, @(4).val))),

    avx_drop_insert_mulmask := Rule( [mul, @(1, Value), [@(2, vinsert_2l_4x64f), _0, @,
            @(3).cond( x -> (x.p[1]=0 and @(1).val=@(2).val.t.value([1, 1, 0, 0])) or
                            (x.p[1]=1 and @(1).val=@(2).val.t.value([0, 0, 1, 1])))]],
        e -> @(2).val),

    avx_drop_store_mulmask_a := Rule( [@(1, vextract_2l_4x64f), [@(2, mul), @(3, Value), @(4)],
            @(5).cond( x -> (x.p[1]=0 and @(3).val=@(2).val.t.value([1, 1, 0, 0])) or
                            (x.p[1]=1 and @(3).val=@(2).val.t.value([0, 0, 1, 1])))],
        e -> vextract_2l_4x64f(@(4).val, @(5).val)),

    avx_drop_store_mulmask_b := Rule( [@(1, vextract_2l_4x64f), [@(2, [add, sub]), @(3), [mul, @(4, Value), @(5)]],
            @(6).cond( x -> (x.p[1]=0 and @(4).val=@(2).val.t.value([1, 1, 0, 0])) or
                            (x.p[1]=1 and @(4).val=@(2).val.t.value([0, 0, 1, 1])))],
        e -> vextract_2l_4x64f(ObjId(@(2).val)(@(3).val, @(5).val), @(6).val)),

    avx_drop_mulmask32f_a := Rule( [mul, @(1, Value), [@(2, vpermf128_8x32f), @(3), _0,
            @(4).cond( x -> x.p[2]>=3 and @(1).val=@(2).val.t.value([1,1,1,1,0,0,0,0]) )]],
        e -> @(2).val),

    avx_drop_mulmask32f_b := Rule( [vextract_4l_8x32f, [@(1, [add, sub]), @(2), [@(3, vperm_8x32f), [mul, @(4, Value), @(5)], @(6)]],
            @(7).cond( x -> (x.p[1]=0 and @(4).val = @(1).val.t.value([1,1,1,1,0,0,0,0])) or
                            (x.p[1]=1 and @(4).val = @(1).val.t.value([0,0,0,0,1,1,1,1])))],
        e -> vextract_4l_8x32f(ObjId(@(1).val)(@(2).val, vperm_8x32f(@(5).val, @(6).val)), @(7).val)),

    avx_drop_mulmask32f_c := Rule( [mul, @(1, Value), [@(2, vinsert_4l_8x32f), _0, [vload_2l_4x32f, _0, @(3)],
            @(4).cond( x -> x.p[1]=0 and @(1).val = @(2).val.t.value([1,1,0,0,0,0,0,0]) )]],
        e -> @(2).val ),

    avx_term_cxtr_8x32f := Rule( vcxtr_8x32f,
        e -> vblend_8x32f(e.args[1], vperm_8x32f(vupermf128_8x32f(e.args[1], [ 2, 1 ]), [3,4,1,2]), [1,1,2,2,2,2,1,1])),

    # complex vectorization loop rules
    # 4-way AVX
    AVX_4x64f_load_TComplex := Rule([@@(1, vpack, (e,cx)-> IsBound(cx.opts) and IsBound(cx.opts.vector) and IsBound(cx.opts.vector.isa) and
        cx.opts.vector.isa = AVX_4x64f),
        [@(2,re), nth], @(3, im, e->e.args[1]=@(2).val.args[1]), @(4,re, e->e.args[1]=@(2).val.args[1]), @(5,im, e->e.args[1]=@(2).val.args[1])],

        e->let(dptr := @(2).val.args[1].toPtr(TVect(T_Real(64), 2)), src := vinsert_2l_4x64f(vzero_4x64f(),deref(dptr),[0]),
            vpermf128_4x64f(src, src, [1,1]))),

    AVX_4x64f_fold_shuffle := Rule([@(1, vshuffle_4x64f), [@(2, vpermf128_4x64f), @(3, vinsert_2l_4x64f), @(6), @(7)],
        @(4).cond(e->e=@(2).val), @(5, vparam, e->Length(Set(e.p))=1)],
        e->vdup(deref(@(3).val.args[2].toPtr(T_Real(64)) + (@(5).val.p[1]-1)), 4)),

    AVX_4x64f_load_TVect := Rule([@@(1, vpack, (e,cx)-> IsBound(cx.opts) and IsBound(cx.opts.vector) and IsBound(cx.opts.vector.isa) and
        cx.opts.vector.isa = AVX_4x64f),
        @(2,nth,e->IsBound(e.loc.init)), @(3, nth, e->e.loc=@(2).val.loc), @(4,nth, e->e.loc=@(2).val.loc), @(5,nth, e->e.loc=@(2).val.loc)],
        e->deref(@(2).val.toPtr(TVect(T_Real(64), 4)))),

#    # 8-way AVX

    AVX_8x32f_broadcast_re := Rule([@(0, vshuffle_8x32f), [@@(1, vpack, (e,cx)-> IsBound(cx.opts) and IsBound(cx.opts.vector) and IsBound(cx.opts.vector.isa) and cx.opts.vector.isa = AVX_8x32f),
                [@(2,re), nth], @(3, im, e->e.args[1]=@(2).val.args[1]),
                @(4,re, e->e.args[1]=@(2).val.args[1]), @(5,im, e->e.args[1]=@(2).val.args[1]),
                @(6,re, e->e.args[1]=@(2).val.args[1]), @(7,im, e->e.args[1]=@(2).val.args[1]),
                @(8,re, e->e.args[1]=@(2).val.args[1]), @(9,im, e->e.args[1]=@(2).val.args[1])],
            @(10).cond(e->e=@@(1).val), @(11,vparam, e->e.p=[1,1,3,3])],
        e->vdup(@(2).val, 8)),

    AVX_8x32f_broadcast_im := Rule([@(0, vshuffle_8x32f), [@@(1, vpack, (e,cx)-> IsBound(cx.opts) and IsBound(cx.opts.vector) and IsBound(cx.opts.vector.isa) and cx.opts.vector.isa = AVX_8x32f),
                [@(2,re), nth], @(3, im, e->e.args[1]=@(2).val.args[1]),
                @(4,re, e->e.args[1]=@(2).val.args[1]), @(5,im, e->e.args[1]=@(2).val.args[1]),
                @(6,re, e->e.args[1]=@(2).val.args[1]), @(7,im, e->e.args[1]=@(2).val.args[1]),
                @(8,re, e->e.args[1]=@(2).val.args[1]), @(9,im, e->e.args[1]=@(2).val.args[1])],
            @(10).cond(e->e=@@(1).val), @(11,vparam, e->e.p=[2,2,4,4])],
        e->vdup(@(3).val, 8)),

    AVX_8x32f_load_TVect := Rule([@@(1, vpack, (e,cx)-> IsBound(cx.opts) and IsBound(cx.opts.vector) and IsBound(cx.opts.vector.isa) and
        cx.opts.vector.isa = AVX_8x32f),
        @(2,nth,e->IsBound(e.loc.init)), @(3, nth, e->e.loc=@(2).val.loc), @(4,nth, e->e.loc=@(2).val.loc), @(5,nth, e->e.loc=@(2).val.loc),
        @(6,nth,e->IsBound(e.loc.init)), @(7, nth, e->e.loc=@(2).val.loc), @(8,nth, e->e.loc=@(2).val.loc), @(9,nth, e->e.loc=@(2).val.loc)],
        e->deref(@(2).val.toPtr(TVect(T_Real(64), 8))))
));

AVXFixProblems := function(c, opts)
    local zerovar, emptyvar, cc;

#    emptyvar := var.fresh_t("ZERO", TVect(TReal, opts.vector.isa.v));
#    zerovar := var.fresh_t("ZERO", TVect(TReal, opts.vector.isa.v));
#
#    cc := SubstTopDown(c, [@(1, neg), @(2)], e->sub(zerovar, @(2).val));
#    cc := SubstTopDown(c, @(1, var, e->IsBound(e.value) and IsList(e.value.v) and ForAll(e.value.v, IsZero)), e->zerovar);
#
#    c := decl([zerovar, emptyvar], chain(assign(zerovar, bin_xor(emptyvar, emptyvar)), cc));
#
#    if Collect(c, vpack) <> [] then Error("Caught in FixProblems"); fi;
    c := UntilDone(c, MergedRuleSet(RulesStrengthReduce, FixCodeAVX), opts);

    return c;
end;

AVX_4x64f.fixProblems := AVXFixProblems;
AVX_8x32f.fixProblems := AVXFixProblems;


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

#
# AVX conversion bridges
#

ISA_Bridge.add(Class(CVT_SSE_4x32f_AVX_8x32f, ISA_Bridge_I, rec(
    isa_from    := AVX_8x32f,
    isa_to      := SSE_4x32f(T_Real(32)),
    props       := [],
    code := (self, y, x, opts) >> let( t := var.fresh_t("c", self.isa_from.t),
        decl( [t], chain(
            assign(t, self._x(x,0)),
            assign(self._y(y,0), vextract_4l_8x32f(t, [0])),
            assign(self._y(y,1), vextract_4l_8x32f(t, [1]))
        )))
)));

ISA_Bridge.add(Class(CVT_AVX_8x32f_AVX_4x64f, ISA_Bridge_I, rec(
    isa_from    := AVX_4x64f,
    isa_to      := AVX_8x32f,
    props       := [],
    code := (self, y, x, opts) >>
        assign(self._y(y,0), vpermf128_8x32f(vcvt_8x32f_4x64f(self._x(x,0)), vcvt_8x32f_4x64f(self._x(x,1)), [1,3]))
)));

ISA_Bridge.add(Class(CVT_AVX_4x64f_SSE_4x32f, ISA_Bridge_I, rec(
    isa_from    := SSE_4x32f(T_Real(32)),
    isa_to      := AVX_4x64f,
    props       := [],
    code := (self, y, x, opts) >>
        assign(self._y(y,0), vcvt_4x64f_4x32f(self._x(x,0)))
)));

ISA_Bridge.add(Class(CVT_AVX_8x32f_SSE_4x32i, ISA_Bridge_I, rec(
    isa_from    := SSE_4x32f(T_Int(32)),
    isa_to      := AVX_8x32f,
    props       := [],
    code := (self, y, x, opts) >>
        assign(self._y(y,0), vcvt_8x32f_8x32i(vinsert_4l_8x32f(vinsert_4l_8x32f(self.isa_to.t.zero(), self._x(x,0), [0]), self._x(x,1), [1])))
)));

ISA_Bridge.add(Class(CVT_SSE_4x32i_AVX_8x32f_round, ISA_Bridge_I, rec(
    isa_from    := AVX_8x32f,
    isa_to      := SSE_4x32f(T_Int(32)),
    props       := ["round"],
    code := (self, y, x, opts) >> let( t := var.fresh_t("c", TVect(T_Int(32), 8)),
        decl([t], chain(
            assign(t, vcvt_8x32i_8x32f(self._x(x,0))),
            assign(self._y(y,0), vextract_4l_8x32f(t, [0])),
            assign(self._y(y,1), vextract_4l_8x32f(t, [1]))
        )))
)));

ISA_Bridge.add(Class(CVT_SSE_4x32i_AVX_8x32f_trunc, ISA_Bridge_I, rec(
    isa_from    := AVX_8x32f,
    isa_to      := SSE_4x32f(T_Int(32)),
    props       := ["trunc"],
    code := (self, y, x, opts) >> let( t := var.fresh_t("c", TVect(T_Int(32), 8)),
        decl([t], chain(
            assign(t, vcvtt_8x32i_8x32f(self._x(x,0))),
            assign(self._y(y,0), vextract_4l_8x32f(t, [0])),
            assign(self._y(y,1), vextract_4l_8x32f(t, [1]))
        )))
)));

ISA_Bridge.add(Class(CVT_AVX_4x64f_SSE_4x32i, ISA_Bridge_I, rec(
    isa_from    := SSE_4x32f(T_Int(32)),
    isa_to      := AVX_4x64f,
    props       := [],
    code := (self, y, x, opts) >>
        assign(self._y(y,0), vcvt_4x64f_4x32i(self._x(x,0)))
)));

ISA_Bridge.add(Class(CVT_SSE_4x32i_AVX_4x64f_round, ISA_Bridge_I, rec(
    isa_from    := AVX_4x64f,
    isa_to      := SSE_4x32f(T_Int(32)),
    props       := ["round"],
    code := (self, y, x, opts) >>
        assign(self._y(y,0), vcvt_4x32i_4x64f(self._x(x,0))),
)));

ISA_Bridge.add(Class(CVT_SSE_4x32i_AVX_4x64f_trunc, ISA_Bridge_I, rec(
    isa_from    := AVX_4x64f,
    isa_to      := SSE_4x32f(T_Int(32)),
    props       := ["trunc"],
    code := (self, y, x, opts) >>
        assign(self._y(y,0), vcvtt_4x32i_4x64f(self._x(x,0))),
)));

ISA_Bridge.add(Class(CVT_AVX_8x32f_8x32f_clip32i, ISA_Bridge_I, rec(
    isa_from    := AVX_8x32f,
    isa_to      := AVX_8x32f,
    props       := ["saturation"],
    range       := (self) >> RangeT(self.clip.min, self.clip.max, T_Real(32).range().eps),
    clip        := rec( min := T_Int(32).range().min, max := T_Int(32).range().max ),
    code := (self, y, x, opts) >> assign( self._y(y,0), min(max(self._x(x,0), self.isa_from.t.value(self.clip.min)), self.isa_from.t.value(self.clip.max)) ),
)));

ISA_Bridge.add(Class(CVT_AVX_8x32f_8x32f_clip32ui, CVT_AVX_8x32f_8x32f_clip32i, rec(
    clip        := rec( min := T_UInt(32).range().min, max := T_UInt(32).range().max ),
)));

ISA_Bridge.add(Class(CVT_AVX_8x32f_8x32f_clip16i, CVT_AVX_8x32f_8x32f_clip32i, rec(
    clip        := rec( min := T_Int(16).range().min,  max := T_Int(16).range().max ),
)));

ISA_Bridge.add(Class(CVT_AVX_8x32f_8x32f_clip16ui, CVT_AVX_8x32f_8x32f_clip32i, rec(
    clip        := rec( min := T_UInt(16).range().min, max := T_UInt(16).range().max ),
)));

ISA_Bridge.add(Class(CVT_AVX_8x32f_8x32f_clip8i, CVT_AVX_8x32f_8x32f_clip32i, rec(
    clip        := rec( min := T_Int(8).range().min,  max := T_Int(8).range().max ),
)));

ISA_Bridge.add(Class(CVT_AVX_8x32f_8x32f_clip8ui, CVT_AVX_8x32f_8x32f_clip32i, rec(
    clip        := rec( min := T_UInt(8).range().min, max := T_UInt(8).range().max ),
)));

ISA_Bridge.add(Class(CVT_AVX_4x64f_4x64f_clip32i, ISA_Bridge_I, rec(
    isa_from    := AVX_4x64f,
    isa_to      := AVX_4x64f,
    props       := ["saturation"],
    range       := (self) >> RangeT(self.clip.min, self.clip.max, T_Real(32).range().eps),
    clip        := rec( min := T_Int(32).range().min, max := T_Int(32).range().max ),
    code := (self, y, x, opts) >> assign( self._y(y,0), min(max(self._x(x,0), self.isa_from.t.value(self.clip.min)), self.isa_from.t.value(self.clip.max)) ),
)));

ISA_Bridge.add(Class(CVT_AVX_4x64f_4x64f_clip32ui, CVT_AVX_4x64f_4x64f_clip32i, rec(
    clip        := rec( min := T_UInt(32).range().min, max := T_UInt(32).range().max ),
)));


text:
TanPlus := function(a, b)
    return (a + b) / (1 - a * b);
end;

TanTimes := function(n, a)
    local x;
    x := 0;
    while n > 0 do
        if IsOddInt(n) then
            x := TanPlus(x, a);
        fi;
        a := TanPlus(a, a);
        n := QuoInt(n, 2);
    od;
    return x;
end;

Check := function(a)
    local x, p;
    x := 0;
    for p in a do
        x := TanPlus(x, SignInt(p[1]) * TanTimes(AbsInt(p[1]), p[2]));
    od;
    return x = 1;
end;

ForAll([
    [[1, 1/2], [1, 1/3]],
    [[2, 1/3], [1, 1/7]],
    [[4, 1/5], [-1, 1/239]],
    [[5, 1/7], [2, 3/79]],
    [[5, 29/278], [7, 3/79]],
    [[1, 1/2], [1, 1/5], [1, 1/8]],
    [[5, 1/7], [4, 1/53], [2, 1/4443]],
    [[6, 1/8], [2, 1/57], [1, 1/239]],
    [[8, 1/10], [-1, 1/239], [-4, 1/515]],
    [[12, 1/18], [8, 1/57], [-5, 1/239]],
    [[16, 1/21], [3, 1/239], [4, 3/1042]],
    [[22, 1/28], [2, 1/443], [-5, 1/1393], [-10, 1/11018]],
    [[22, 1/38], [17, 7/601], [10, 7/8149]],
    [[44, 1/57], [7, 1/239], [-12, 1/682], [24, 1/12943]],
    [[88, 1/172], [51, 1/239], [32, 1/682], [44, 1/5357], [68, 1/12943]]], Check);

Check([[88, 1/172], [51, 1/239], [32, 1/682], [44, 1/5357], [68, 1/12944]]);

ext: gap

max_stars_repo_name: mullikine/RosettaCodeData

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(Rules_fAdd, RuleSet);
RewriteRules(Rules_fAdd, rec(
    fTensor := Rule([@(1, fTensor), @(2, fAdd, e->ForAll(e.params, i->IsInt(i) or IsValue(i))), @(3, fId, e-> IsInt(e.params[1]) or IsValue(e.params[1]))],
        e-> let(fadd := @(2).val, fid := @(3).val.domain(), fAdd(fadd.params[1]*fid, fadd.params[2]*fid, fadd.params[3]*fid))),

    fTensorVScatVGath :=ARule(fTensor, [ @@(2, fAdd, (e,cx)->(IsBound(cx.VScat) and Length(cx.VScat)>=1) or (IsBound(cx.VGath) and Length(cx.VGath)>=1)), @(3, fId) ],
        e-> let(fadd := @@(2).val, fid := @(3).val.domain(), [ fAdd(fadd.params[1]*fid, fadd.params[2]*fid, fadd.params[3]*fid) ] )),

    fAddFBase := ARule(fCompose, [@@(2, fAdd, (e,cx)->(IsBound(cx.VScat) and Length(cx.VScat)>=1) or (IsBound(cx.VGath) and Length(cx.VGath)>=1)), @(3, fBase)],
        e-> let(fadd := @(2).val, fbase := @(3).val, [ fAdd(fadd.params[1], 1, fadd.params[3]+fbase.params[2]) ] ))

));


text:
#############################################################################
##
#W  vspc.gi                     GAP library                     Thomas Breuer
##
##
#Y  Copyright (C)  1997,  Lehrstuhl D für Mathematik,  RWTH Aachen,  Germany
#Y  (C) 1998 School Math and Comp. Sci., University of St Andrews, Scotland
#Y  Copyright (C) 2002 The GAP Group
##
##  This file contains generic methods for vector spaces.
##

#############################################################################
##
#M  SetLeftActingDomain( <extL>, <D> )
##
##  check whether the left acting domain <D> of the external left set <extL>
##  knows that it is a division ring.
##  This is used, e.g.,  to tell a free module over a division ring
##  that it is a vector space.
##
InstallOtherMethod( SetLeftActingDomain,
    "method to set also 'IsLeftActedOnByDivisionRing'",
    [ IsAttributeStoringRep and IsLeftActedOnByRing, IsObject ],0,
    function( extL, D )
    if HasIsDivisionRing( D ) and IsDivisionRing( D ) then
      SetIsLeftActedOnByDivisionRing( extL, true );
    fi;
    TryNextMethod();
    end );

#############################################################################
##
#M  IsLeftActedOnByDivisionRing( <M> )
##
InstallMethod( IsLeftActedOnByDivisionRing,
    "method for external left set that is left acted on by a ring",
    [ IsExtLSet and IsLeftActedOnByRing ],
    function( M )
    if IsIdenticalObj( M, LeftActingDomain( M ) ) then
      TryNextMethod();
    else
      return IsDivisionRing( LeftActingDomain( M ) );
    fi;
    end );

#############################################################################
##
#F  VectorSpace( <F>, <gens>[, <zero>][, "basis"] )
##
##  The only difference between `VectorSpace' and `FreeLeftModule' shall be
##  that the left acting domain of a vector space must be a division ring.
##
InstallGlobalFunction( VectorSpace, function( arg )
    if Length( arg ) = 0 or not IsDivisionRing( arg[1] ) then
      Error( "usage: VectorSpace( <F>, <gens>[, <zero>][, \"basis\"] )" );
    fi;
    return CallFuncList( FreeLeftModule, arg );
    end );

#############################################################################
##
#M  AsSubspace( <V>, <C> )  . . . . . . . for a vector space and a collection
##
InstallMethod( AsSubspace,
    "for a vector space and a collection",
    [ IsVectorSpace, IsCollection ],
    function( V, C )
    local newC;

    if not IsSubset( V, C ) then
      return fail;
    fi;
    newC:= AsVectorSpace( LeftActingDomain( V ), C );
    if newC = fail then
      return fail;
    fi;
    SetParent( newC, V );
    UseIsomorphismRelation( C, newC );
    UseSubsetRelation( C, newC );

    return newC;
    end );

#############################################################################
##
#M  AsLeftModule( <F>, <V> )  . . . . . .  for division ring and vector space
##
##  View the vector space <V> as a vector space over the division ring <F>.
##
InstallMethod( AsLeftModule,
    "method for a division ring and a vector space",
    [ IsDivisionRing, IsVectorSpace ],
    function( F, V )

    local W,        # the space, result
          base,     # basis vectors of field extension
          gen,      # loop over generators of 'V'
          b,        # loop over 'base'
          gens,     # generators of 'V'
          newgens;  # extended list of generators

    if Characteristic( F ) <> Characteristic( LeftActingDomain( V ) ) then

      # This is impossible.
      return fail;

    elif F = LeftActingDomain( V ) then

      # No change of the left acting domain is necessary.
      return V;

    elif IsSubset( F, LeftActingDomain( V ) ) then

      # Check whether 'V' is really a space over the bigger field,
      # that is, whether the set of elements does not change.
      base:= BasisVectors( Basis( AsField( LeftActingDomain( V ), F ) ) );
      for gen in GeneratorsOfLeftModule( V ) do
        for b in base do
          if not b * gen in V then

            # The field extension would change the set of elements.
            return fail;

          fi;
        od;
      od;

      # Construct the space.
      W:= LeftModuleByGenerators( F, GeneratorsOfLeftModule(V), Zero(V) );

    elif IsSubset( LeftActingDomain( V ), F ) then

      # View 'V' as a space over a smaller field.
      # For that, the list of generators must be extended.
      gens:= GeneratorsOfLeftModule( V );
      if IsEmpty( gens ) then
        W:= LeftModuleByGenerators( F, [], Zero( V ) );
      else

        base:= BasisVectors( Basis( AsField( F, LeftActingDomain( V ) ) ) );
        newgens:= [];
        for b in base do
          for gen in gens do
            Add( newgens, b * gen );
          od;
        od;
        W:= LeftModuleByGenerators( F, newgens );

      fi;

    else

      # View 'V' first as space over the intersection of fields,
      # and then over the desired field.
      return AsLeftModule( F,
                 AsLeftModule( Intersection( F,
                     LeftActingDomain( V ) ), V ) );

    fi;

    UseIsomorphismRelation( V, W );
    UseSubsetRelation( V, W );
    return W;
    end );

#############################################################################
##
#M  ViewObj( <V> )  . . . . . . . . . . . . . . . . . . . view a vector space
##
##  print left acting domain, if known also dimension or no. of generators
##
InstallMethod( ViewObj,
    "for vector space with known generators",
    [ IsVectorSpace and HasGeneratorsOfLeftModule ],
    function( V )
    Print( "<vector space over ", LeftActingDomain( V ), ", with ",
           Length( GeneratorsOfLeftModule( V ) ), " generators>" );
    end );

InstallMethod( ViewObj,
    "for vector space with known dimension",
    [ IsVectorSpace and HasDimension ],
    1, # override method for known generators
    function( V )
    Print( "<vector space of dimension ", Dimension( V ),
           " over ", LeftActingDomain( V ), ">" );
    end );

InstallMethod( ViewObj,
    "for vector space",
    [ IsVectorSpace ],
    function( V )
    Print( "<vector space over ", LeftActingDomain( V ), ">" );
    end );

#############################################################################
##
#M  PrintObj( <V> ) . . . . . . . . . . . . . . . . . . .  for a vector space
##
InstallMethod( PrintObj,
    "method for vector space with left module generators",
    [ IsVectorSpace and HasGeneratorsOfLeftModule ],
    function( V )
    Print( "VectorSpace( ", LeftActingDomain( V ), ", ",
           GeneratorsOfLeftModule( V ) );
    if IsEmpty( GeneratorsOfLeftModule( V ) ) and HasZero( V ) then
      Print( ", ", Zero( V ), " )" );
    else
      Print( " )" );
    fi;
    end );

InstallMethod( PrintObj,
    "method for vector space",
    [ IsVectorSpace ],
    function( V )
    Print( "VectorSpace( ", LeftActingDomain( V ), ", ... )" );
    end );

#############################################################################
##
#M  \/( <V>, <W> )  . . . . . . . . .  factor of a vector space by a subspace
#M  \/( <V>, <vectors> )  . . . . . .  factor of a vector space by a subspace
##
InstallOtherMethod( \/,
    "method for vector space and collection",
    IsIdenticalObj,
    [ IsVectorSpace, IsCollection ],
    function( V, vectors )
    if IsVectorSpace( vectors ) then
      TryNextMethod();
    else
      return V / Subspace( V, vectors );
    fi;
    end );

InstallOtherMethod( \/,
    "generic method for two vector spaces",
    IsIdenticalObj,
    [ IsVectorSpace, IsVectorSpace ],
    function( V, W )
    return ImagesSource( NaturalHomomorphismBySubspace( V, W ) );
    end );

#############################################################################
##
#M  Intersection2Spaces( <AsStruct>, <Substruct>, <Struct> )
##
InstallGlobalFunction( Intersection2Spaces,
    function( AsStructure, Substructure, Structure )
    return function( V, W )
    local inters,  # intersection, result
          F,       # coefficients field
          gensV,   # list of generators of 'V'
          gensW,   # list of generators of 'W'
          VW,      # sum of 'V' and 'W'
          B;       # basis of 'VW'

    if LeftActingDomain( V ) <> LeftActingDomain( W ) then

      # Compute the intersection as vector space over the intersection
      # of the coefficients fields.
      # (Note that the characteristic is the same.)
      F:= Intersection2( LeftActingDomain( V ), LeftActingDomain( W ) );
      return Intersection2( AsStructure( F, V ), AsStructure( F, W ) );

    elif IsFiniteDimensional( V ) and IsFiniteDimensional( W ) then

      # Compute the intersection of two spaces over the same field.
      gensV:= GeneratorsOfLeftModule( V );
      gensW:= GeneratorsOfLeftModule( W );
      if   IsEmpty( gensV ) then
        if Zero( V ) in W then
          inters:= V;
        else
          inters:= [];
        fi;
      elif IsEmpty( gensW ) then
        if Zero( V ) in W then
          inters:= W;
        else
          inters:= [];
        fi;
      else
        # Compute a common coefficient space.
        VW:= LeftModuleByGenerators( LeftActingDomain( V ),
                                     Concatenation( gensV, gensW ) );
        B:= Basis( VW );

        # Construct the coefficient subspaces corresponding to 'V' and 'W'.
        gensV:= List( gensV, x -> Coefficients( B, x ) );
        gensW:= List( gensW, x -> Coefficients( B, x ) );

        # Construct the intersection of row spaces, and carry back to VW.
        inters:= List( SumIntersectionMat( gensV, gensW )[2],
                       x -> LinearCombination( B, x ) );

        # Construct the intersection space, if possible with a parent.
        if     HasParent( V ) and HasParent( W )
           and IsIdenticalObj( Parent( V ), Parent( W ) ) then
          inters:= Substructure( Parent( V ), inters, "basis" );
        elif IsEmpty( inters ) then
          inters:= Substructure( V, inters, "basis" );
          SetIsTrivial( inters, true );
        else
          inters:= Structure( LeftActingDomain( V ), inters, "basis" );
        fi;

        # Run implications by the subset relation.
        UseSubsetRelation( V, inters );
        UseSubsetRelation( W, inters );
      fi;

      # Return the result.
      return inters;

    else
      TryNextMethod();
    fi;
    end;
end );

#############################################################################
##
#M  Intersection2( <V>, <W> ) . . . . . . . . . . . . . for two vector spaces
##
InstallMethod( Intersection2,
    "method for two vector spaces",
    IsIdenticalObj,
    [ IsVectorSpace, IsVectorSpace ],
    Intersection2Spaces( AsLeftModule, SubspaceNC, VectorSpace ) );

#############################################################################
##
#M  ClosureLeftModule( <V>, <a> ) . . . . . . . . . closure of a vector space
##
InstallMethod( ClosureLeftModule,
    "method for a vector space with basis, and a vector",
    IsCollsElms,
    [ IsVectorSpace and HasBasis, IsVector ],
    function( V, w )
    local   B; # basis of 'V'

    # We can test membership easily.
    B:= Basis( V );
#T why easily?
    if Coefficients( B, w ) = fail then

      # In the case of a vector space, we know a basis of the closure.
      B:= Concatenation( BasisVectors( B ), [ w ] );
      V:= LeftModuleByGenerators( LeftActingDomain( V ), B );
      UseBasis( V, B );

    fi;
    return V;
    end );

#############################################################################
##
##  Methods for collections of subspaces of a vector space
##

#############################################################################
##
#R  IsSubspacesVectorSpaceDefaultRep( <D> )
##
##  is the representation of domains of subspaces of a vector space <V>,
##  with the components 'structure' (with value <V>) and 'dimension'
##  (with value either the dimension of the subspaces in the domain
##  or the string '\"all\"', which means that the domain contains all
##  subspaces of <V>).
##
DeclareRepresentation(
    "IsSubspacesVectorSpaceDefaultRep",
    IsComponentObjectRep,
    [ "dimension", "structure" ] );
#T not IsAttributeStoringRep?

#############################################################################
##
#M  PrintObj( <D> )  . . . . . . . . . . . . . . . . . for a subspaces domain
##
InstallMethod( PrintObj,
    "method for a subspaces domain",
    [ IsSubspacesVectorSpace and IsSubspacesVectorSpaceDefaultRep ],
    function( D )
    if IsInt( D!.dimension ) then
      Print( "Subspaces( ", D!.structure, ", ", D!.dimension, " )" );
    else
      Print( "Subspaces( ", D!.structure, " )" );
    fi;
    end );

#############################################################################
##
#M  Size( <D> ) . . . . . . . . . . . . . . . . . . .  for a subspaces domain
##
##  The number of $k$-dimensional subspaces in a $n$-dimensional space over
##  the field with $q$ elements is
##  $$
##  a(n,k) = \prod_{i=0}^{k-1} \frac{q^n-q^i}{q^k-q^i} =
##           \prod_{i=0}^{k-1} \frac{q^{n-i}-1}{q^{k-i}-1}.
##  $$
##  We have the recursion
##  $$
##  a(n,k+1) = a(n,k) \frac{q^{n-i}-1}{q^{i+1}-1}.
##  $$
##
##  (The number of all subspaces is $\sum_{k=0}^n a(n,k)$.)
##
InstallMethod( Size,
    "method for a subspaces domain",
    [ IsSubspacesVectorSpace and IsSubspacesVectorSpaceDefaultRep ],
    function( D )

    local k,
          n,
          q,
          size,
          qn,
          qd,
          ank,
          i;

    if D!.dimension = "all" then

      # all subspaces of the space
      n:= Dimension( D!.structure );

      q:= Size( LeftActingDomain( D!.structure ) );
      size:= 1;
      qn:= q^n;
      qd:= q;

      # $a(n,0)$
      ank:= 1;

      for k in [ 1 .. Int( (n-1)/2 ) ] do

        # Compute $a(n,k)$.
        ank:= ank * ( qn - 1 ) / ( qd - 1 );
        qn:= qn / q;
        qd:= qd * q;

        size:= size + ank;

      od;

      size:= 2 * size;

      if n mod 2 = 0 then

        # Add the number of spaces of dimension $n/2$.
        size:= size + ank * ( qn - 1 ) / ( qd - 1 );
      fi;

    else

      # number of spaces of dimension 'k' only
      n:= Dimension( D!.structure );
      if   D!.dimension < 0 or
           n < D!.dimension then
        return 0;
      elif n / 2 < D!.dimension then
        k:= n - D!.dimension;
      else
        k:= D!.dimension;
      fi;

      q:= Size( LeftActingDomain( D!.structure ) );
      size:= 1;

      qn:= q^n;
      qd:= q;
      for i in [ 1 .. k ] do
        size:= size * ( qn - 1 ) / ( qd - 1 );
        qn:= qn / q;
        qd:= qd * q;
      od;

    fi;

    # Return the result.
    return size;
    end );

#############################################################################
##
#M  Enumerator( <D> ) . . . . . . . . . . . . . . . .  for a subspaces domain
##
##  Use the iterator to compute the elements list.
#T This is not allowed!
##
InstallMethod( Enumerator,
    "method for a subspaces domain",
    [ IsSubspacesVectorSpace and IsSubspacesVectorSpaceDefaultRep ],
    function( D )
    local iter,    # iterator for 'D'
          elms;    # elements list, result

    iter:= Iterator( D );
    elms:= [];
    while not IsDoneIterator( iter ) do
      Add( elms, NextIterator( iter ) );
    od;
    return elms;
    end );
#T necessary?

#############################################################################
##
#M  Iterator( <D> ) . . . . . . . . . . . . . . . . .  for a subspaces domain
##
##  uses the subspaces iterator for full row spaces and the mechanism of
##  associated row spaces.
##
BindGlobal( "IsDoneIterator_Subspaces",
    iter -> IsDoneIterator( iter!.associatedIterator ) );

BindGlobal( "NextIterator_Subspaces", function( iter )
    local next;
    next:= NextIterator( iter!.associatedIterator );
    next:= List( GeneratorsOfLeftModule( next ),
                 x -> LinearCombination( iter!.basis, x ) );
    return Subspace( iter!.structure, next, "basis" );
    end );

BindGlobal( "ShallowCopy_Subspaces",
    iter -> rec( structure          := iter!.structure,
                 basis              := iter!.basis,
                 associatedIterator := ShallowCopy(
                                           iter!.associatedIterator ) ) );

InstallMethod( Iterator,
    "for a subspaces domain",
    [ IsSubspacesVectorSpace and IsSubspacesVectorSpaceDefaultRep ],
    function( D )
    local V;      # the vector space

    V:= D!.structure;
    return IteratorByFunctions( rec(
               IsDoneIterator     := IsDoneIterator_Subspaces,
               NextIterator       := NextIterator_Subspaces,
               ShallowCopy        := ShallowCopy_Subspaces,
               structure          := V,
               basis              := Basis( V ),
               associatedIterator := Iterator(
                      Subspaces( FullRowSpace( LeftActingDomain( V ),
                                               Dimension( V ) ),
                                 D!.dimension ) ) ) );
    end );

#############################################################################
##
#M  Subspaces( <V>, <dim> )
##
InstallMethod( Subspaces,
    "for a vector space, and an integer",
    [ IsVectorSpace, IsInt ],
    function( V, dim )
    if IsFinite( V ) then
      return Objectify( NewType( CollectionsFamily( FamilyObj( V ) ),
                                     IsSubspacesVectorSpace
                                 and IsSubspacesVectorSpaceDefaultRep ),
                        rec(
                             structure  := V,
                             dimension  := dim
                           )
                      );
    else
      TryNextMethod();
    fi;
    end );

#############################################################################
##
#M  Subspaces( <V> )
##
InstallMethod( Subspaces,
    "for a vector space",
    [ IsVectorSpace ],
    function( V )
    if IsFinite( V ) then
      return Objectify( NewType( CollectionsFamily( FamilyObj( V ) ),
                                     IsSubspacesVectorSpace
                                 and IsSubspacesVectorSpaceDefaultRep ),
                        rec(
                             structure  := V,
                             dimension  := "all"
                           )
                      );
    else
      TryNextMethod();
    fi;
    end );

#############################################################################
##
#F  IsSubspace( <V>, <U> ) . . . . . . . . . . . . . . . . . check <U> <= <V>
##
InstallGlobalFunction( IsSubspace, function( V, U )
    return IsVectorSpace( U ) and IsSubset( V, U );
end );

#############################################################################
##
#M  IsVectorSpaceHomomorphism( <map> )
##
InstallMethod( IsVectorSpaceHomomorphism,
    [ IsGeneralMapping ],
    function( map )
    local S, R, F;
    S:= Source( map );
    if not IsVectorSpace( S ) then
      return false;
    fi;
    R:= Range( map );
    if not IsVectorSpace( R ) then
      return false;
    fi;
    F:= LeftActingDomain( S );
    return ( F = LeftActingDomain( R ) ) and IsLinearMapping( F, map );
    end );

#############################################################################
##
#E

max_stars_repo_name: JavascriptID/sourcerer-app

max_issues_repo_name: JavascriptID/sourcerer-app

max_forks_repo_name: JavascriptID/sourcerer-app


text:
#
# Read("~/Workspace/groupsSB/epi/A2/handle.ext.gi");
#

type:="A";
rank:=2;
nr_pos_roots:=3;

Read("~/Workspace/groupsSB/epi/group.gi");

highest_weight:=[1,1];# corresponding highest weight module is the Lie algebra

Read("~/Workspace/groupsSB/epi/ext.gi");

handleReg:=function()
    local k,ua,ub,uab,u,e,h,h1,h2;
    k:=3; # power of exterior algebra
    ua:=ext_root_group(k,1,1);
    ub:=ext_root_group(k,2,1);
    uab:=ext_root_group(k,3,1/2); #p\neq 2
    u:=ua*ub*uab;
    #
    #uua:=ext_root_group(1,2);
    #uub:=ext_root_group(2,2);
    #uuab:=ext_root_group(3,2);
    #uu:=uua*uub*uuab;
    #
    id_mat:=u^0;
    e:=u-id_mat;

    ese:=Eigenspaces(Rationals,e)[1]; #eigenvalues 0, centralizer
    h1:=ext_e(k,cb0[7]);
    h2:=ext_e(k,cb0[8]);;
    h:=h1*h2;
    esh1:=Eigenspaces(Rationals,h);
    evh1:=Eigenvalues(Rationals,h);

    esh10:=[];
    index_ev0:=Position(evh1,0);
    esh10:=esh1[index_ev0];

    #
    # sum and intersection
    #
    si:=SumIntersectionMat(Basis(ese),Basis(esh10));
    inter:=si[2];
    return inter;
end;

#
# c_2=0
#
#handle2:=function()
#   local k,ua,ub,uab,u,e;
k:=3;
ua:=ext_root_group(k,1,1);
uab:=ext_root_group(k,3,1);
u:=ua*uab;
e:=u-id_mat;

ese:=Eigenspaces(Rationals,e)[1]; #eigenvalues 0, centralizer
h1:=ext_e(k,cb0[7]);;
esh1:=Eigenspaces(Rationals,h1);
evh1:=Eigenvalues(Rationals,h1);

esh10:=[];
index_ev0:=Position(evh1,0);
esh10:=esh1[index_ev0];

#
# sum and intersection
#
si:=SumIntersectionMat(Basis(ese),Basis(esh10));
inter:=si[2];
#    return inter;
#end;

max_stars_repo_name: iuliansimion/groupsSB

max_issues_repo_name: iuliansimion/groupsSB

max_forks_repo_name: iuliansimion/groupsSB


text:
# Built-in

L := List([1 .. 100], n -> Random(1, 10));

MaximumList(L);
# 10


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

# ==========================================================================
# Scale
# ==========================================================================
Class(Scale, BaseContainer, rec(
    new := meth(self, s, spl)
        Constraint(IsSPL(spl));
	if s = 1 then return spl; fi;
        Constraint(IsScalarOrNum(s));
        return SPL(WithBases( self,
	    rec( _children := [spl],
		scalar := s,
		dimensions := spl.dims() )));
    end,

    rChildren := self >> [self.scalar, self._children[1]],
    rSetChild := meth(self, n, what)
        if n=1 then self.scalar := what;
        elif n=2 then self._children[1] := what;
        else Error("<n> must be in [1..2]");
        fi;
    end,
    #from_rChildren  -- inherited
    #-----------------------------------------------------------------------
    isPermutation := False,
    #-----------------------------------------------------------------------
    isReal := self >> IsRealNumber(self.scalar) and IsRealSPL(self._children[1]),
    #-----------------------------------------------------------------------
    toAMat := self >> EvalScalar(self.scalar) * AMatSPL(self._children[1]),
    #-----------------------------------------------------------------------
    arithmeticCost := meth(self, costMul, costAddMul)
        return costMul(self.scalar) * Minimum(self.dimensions)
	       + self._children[1].arithmeticCost(costMul, costAddMul);
    end
));


text:
# Some small pregroups

# The lists of small pregroups were generated by
# Chris Jefferson <caj21@st-andrews.ac.uk> and
# Richard Parker.
ANATPH_small_pregroups[1] :=
   [ [ [ 1 ] ] ];

ANATPH_small_pregroups[2] :=
   [ [ [1, 2],
       [2, 1] ] ];

ANATPH_small_pregroups[3] :=
   [ [ [1, 2, 3],
       [2, 1, 0],
       [3, 0, 1] ] ];
# missing: Cyclic group of order 3

ANATPH_small_pregroups[4] :=
 [ [ [1, 2, 3, 4],
     [2, 1, 0, 0],
     [3, 0, 1, 0],
     [4, 0, 0, 1] ],
   [ [1, 2, 3, 4],
     [2, 1, 4, 3],
     [3, 4, 1, 2],
     [4, 3, 2, 1] ],
   [ [1, 2, 3, 4],
     [2, 1, 0, 0],
     [3, 0, 0, 1],
     [4, 0, 1, 0] ],
   [ [1, 2, 3, 4],
     [2, 1, 0, 0],
     [3, 0, 4, 1],
     [4, 0, 1, 3] ],
   [ [1, 2, 3, 4],
     [2, 1, 4, 3],
     [3, 4, 2, 1],
     [4, 3, 1, 2] ] ];

ANATPH_small_pregroups[5] :=
   [ [ [1, 2, 3, 4, 5],
       [2, 1, 0, 0, 0],
       [3, 0, 1, 0, 0],
       [4, 0, 0, 1, 0],
       [5, 0, 0, 0, 1] ],
     [ [1, 2, 3, 4, 5],
       [2, 1, 4, 3, 0],
       [3, 4, 1, 2, 0],
       [4, 3, 2, 1, 0],
       [5, 0, 0, 0, 1]],
     [ [1, 2, 3, 4, 5],
       [2, 1, 0, 0, 0],
       [3, 0, 1, 0, 0],
       [4, 0, 0, 0, 1],
       [5, 0, 0, 1, 0]],
     [ [1, 2, 3, 4, 5],
       [2, 1, 0, 0, 0],
       [3, 0, 1, 0, 0],
       [4, 0, 0, 5, 1],
       [5, 0, 0, 1, 4]],
     [ [1, 2, 3, 4, 5],
       [2, 1, 4, 3, 0],
       [3, 4, 2, 1, 0],
       [4, 3, 1, 2, 0],
       [5, 0, 0, 0, 1]],
     [ [1, 2, 3, 4, 5],
       [2, 0, 1, 0, 0],
       [3, 1, 0, 0, 0],
       [4, 0, 0, 0, 1],
       [5, 0, 0, 1, 0]],
     [ [1, 2, 3, 4, 5],
       [2, 3, 1, 0, 0],
       [3, 1, 2, 0, 0],
       [4, 0, 0, 0, 1],
       [5, 0, 0, 1, 0]],
     [ [1, 2, 3, 4, 5],
       [2, 3, 1, 0, 0],
       [3, 1, 2, 0, 0],
       [4, 0, 0, 5, 1],
       [5, 0, 0, 1, 4]],
     [ [1, 2, 3, 4, 5],
       [2, 3, 4, 5, 1],
       [3, 4, 5, 1, 2],
       [4, 5, 1, 2, 3],
       [5, 1, 2, 3, 4]]];

ANATPH_small_pregroups[6] :=
    [ [ [1,2,3,4,5,6],
        [2,1,0,5,4,0],
        [3,6,0,1,0,0],
        [4,0,1,0,0,2],
        [5,0,2,0,0,1],
        [6,3,0,0,1,0] ],
      [ [1,2,3,4,5,6],
        [2,1,5,0,3,0],
        [3,0,0,1,0,2],
        [4,6,1,0,0,0],
        [5,0,0,2,0,1],
        [6,4,0,0,1,0] ],
      [ [1,2,3,4,5,6],
        [2,1,0,6,0,4],
        [3,5,0,1,0,0],
        [4,0,1,0,2,0],
        [5,3,0,0,0,1],
        [6,0,2,0,1,0] ],
      [ [1,2,3,4,5,6],
        [2,1,5,6,3,4],
        [3,5,0,1,0,2],
        [4,6,1,0,2,0],
        [5,3,0,2,0,1],
        [6,4,2,0,1,0] ],
      [ [1,2,3,4,5,6],
        [2,1,6,0,0,3],
        [3,0,0,1,2,0],
        [4,5,1,0,0,0],
        [5,4,0,0,0,1],
        [6,0,0,2,1,0] ],
      [ [1,2,3,4,5,6],
        [2,1,6,5,4,3],
        [3,6,0,1,2,0],
        [4,5,1,0,0,2],
        [5,4,2,0,0,1],
        [6,3,0,2,1,0] ],
      [ [1,2,3,4,5,6],
        [2,1,5,6,3,4],
        [3,6,1,0,0,2],
        [4,5,0,1,2,0],
        [5,4,2,0,0,1],
        [6,3,0,2,1,0] ],
      [ [1,2,3,4,5,6],
        [2,1,0,6,0,4],
        [3,0,1,5,4,0],
        [4,5,6,1,2,3],
        [5,4,0,3,0,1],
        [6,0,4,2,1,0] ],
      [ [1,2,3,4,5,6],
        [2,1,0,5,4,0],
        [3,0,1,6,0,4],
        [4,6,5,1,3,2],
        [5,0,4,2,0,1],
        [6,4,0,3,1,0] ],
      [ [1,2,3,4,5,6],
        [2,1,6,0,0,3],
        [3,5,1,6,2,4],
        [4,0,5,1,3,0],
        [5,3,4,0,0,1],
        [6,0,2,3,1,0] ],
      [ [1,2,3,4,5,6],
        [2,1,4,3,6,5],
        [3,4,1,2,0,0],
        [4,3,2,1,0,0],
        [5,6,0,0,2,1],
        [6,5,0,0,1,2] ],
      [ [1,2,3,4,5,6],
        [2,1,4,3,0,0],
        [3,4,1,2,6,5],
        [4,3,2,1,0,0],
        [5,0,6,0,3,1],
        [6,0,5,0,1,3] ],
      [ [1,2,3,4,5,6],
        [2,1,4,3,0,0],
        [3,4,1,2,0,0],
        [4,3,2,1,6,5],
        [5,0,0,6,4,1],
        [6,0,0,5,1,4] ],
      [ [1,2,3,4,5,6],
        [2,1,6,5,4,3],
        [3,6,1,0,0,2],
        [4,5,0,1,2,0],
        [5,4,0,2,1,0],
        [6,3,2,0,0,1] ],
      [ [1,2,3,4,5,6],
        [2,1,6,0,0,3],
        [3,6,1,5,4,2],
        [4,0,5,1,3,0],
        [5,0,4,3,1,0],
        [6,3,2,0,0,1] ],
      [ [1,2,3,4,5,6],
        [2,1,5,6,3,4],
        [3,5,1,0,2,0],
        [4,6,0,1,0,2],
        [5,3,2,0,1,0],
        [6,4,0,2,0,1] ],
      [ [1,2,3,4,5,6],
        [2,1,0,6,0,4],
        [3,0,1,5,4,0],
        [4,6,5,1,3,2],
        [5,0,4,3,1,0],
        [6,4,0,2,0,1] ],
      [ [1,2,3,4,5,6],
        [2,1,5,0,3,0],
        [3,5,1,6,2,4],
        [4,0,6,1,0,3],
        [5,3,2,0,1,0],
        [6,0,4,3,0,1] ],
      [ [1,2,3,4,5,6],
        [2,1,0,5,4,0],
        [3,0,1,6,0,4],
        [4,5,6,1,2,3],
        [5,4,0,2,1,0],
        [6,0,4,3,0,1] ],
      [ [1,2,3,4,5,6],
        [2,1,4,3,6,5],
        [3,4,1,2,0,0],
        [4,3,2,1,0,0],
        [5,6,0,0,1,2],
        [6,5,0,0,2,1] ],
      [ [1,2,3,4,5,6],
        [2,1,0,0,6,5],
        [3,0,1,5,4,0],
        [4,0,5,1,3,0],
        [5,6,4,3,1,2],
        [6,5,0,0,2,1] ],
      [ [1,2,3,4,5,6],
        [2,1,0,0,6,5],
        [3,0,1,6,0,4],
        [4,0,6,1,0,3],
        [5,6,0,0,1,2],
        [6,5,4,3,2,1] ],
      [ [1,2,3,4,5,6],
        [2,1,4,3,0,0],
        [3,4,1,2,6,5],
        [4,3,2,1,0,0],
        [5,0,6,0,1,3],
        [6,0,5,0,3,1] ],
      [ [1,2,3,4,5,6],
        [2,1,0,5,4,0],
        [3,0,1,0,6,5],
        [4,5,0,1,2,0],
        [5,4,6,2,1,3],
        [6,0,5,0,3,1] ],
      [ [1,2,3,4,5,6],
        [2,1,0,6,0,4],
        [3,0,1,0,6,5],
        [4,6,0,1,0,2],
        [5,0,6,0,1,3],
        [6,4,5,2,3,1] ],
      [ [1,2,3,4,5,6],
        [2,1,4,3,0,0],
        [3,4,1,2,0,0],
        [4,3,2,1,6,5],
        [5,0,0,6,1,4],
        [6,0,0,5,4,1] ],
      [ [1,2,3,4,5,6],
        [2,1,5,0,3,0],
        [3,5,1,0,2,0],
        [4,0,0,1,6,5],
        [5,3,2,6,1,4],
        [6,0,0,5,4,1] ],
      [ [1,2,3,4,5,6],
        [2,1,6,0,0,3],
        [3,6,1,0,0,2],
        [4,0,0,1,6,5],
        [5,0,0,6,1,4],
        [6,3,2,5,4,1] ]
    ];

ANATPH_small_pregroups[7] :=
    [ [ [1, 2, 3, 4, 5, 6, 7],
        [2, 1, 0, 0, 0, 0, 0],
        [3, 0, 1, 0, 0, 0, 0],
        [4, 0, 0, 1, 0, 0, 0],
        [5, 0, 0, 0, 1, 0, 0],
        [6, 0, 0, 0, 0, 1, 0],
        [7, 0, 0, 0, 0, 0, 1] ],
      [ [1, 2, 3, 4, 5, 6, 7],
        [2, 1, 4, 3, 0, 0, 0],
        [3, 4, 1, 2, 0, 0, 0],
        [4, 3, 2, 1, 0, 0, 0],
        [5, 0, 0, 0, 1, 0, 0],
        [6, 0, 0, 0, 0, 1, 0],
        [7, 0, 0, 0, 0, 0, 1] ],
      [ [1, 2, 3, 4, 5, 6, 7],
        [2, 1, 4, 3, 0, 0, 0],
        [3, 4, 1, 2, 0, 0, 0],
        [4, 3, 2, 1, 0, 0, 0],
        [5, 0, 0, 0, 1, 7, 6],
        [6, 0, 0, 0, 7, 1, 5],
        [7, 0, 0, 0, 6, 5, 1] ],
      [ [1, 2, 3, 4, 5, 6, 7],
        [2, 1, 4, 3, 6, 5, 0],
        [3, 4, 1, 2, 0, 0, 0],
        [4, 3, 2, 1, 0, 0, 0],
        [5, 6, 0, 0, 1, 2, 0],
        [6, 5, 0, 0, 2, 1, 0],
        [7, 0, 0, 0, 0, 0, 1] ],
      [ [1, 2, 3, 4, 5, 6, 7],
        [2, 1, 0, 0, 0, 0, 0],
        [3, 0, 1, 0, 0, 0, 0],
        [4, 0, 0, 1, 0, 0, 0],
        [5, 0, 0, 0, 1, 0, 0],
        [6, 0, 0, 0, 0, 0, 1],
        [7, 0, 0, 0, 0, 1, 0] ],
      [ [1, 2, 3, 4, 5, 6, 7],
        [2, 1, 0, 0, 0, 0, 0],
        [3, 0, 1, 0, 0, 0, 0],
        [4, 0, 0, 1, 0, 0, 0],
        [5, 0, 0, 0, 1, 0, 0],
        [6, 0, 0, 0, 0, 7, 1],
        [7, 0, 0, 0, 0, 1, 6] ],
      [ [1, 2, 3, 4, 5, 6, 7],
        [2, 1, 4, 3, 0, 0, 0],
        [3, 4, 1, 2, 0, 0, 0],
        [4, 3, 2, 1, 0, 0, 0],
        [5, 0, 0, 0, 1, 0, 0],
        [6, 0, 0, 0, 0, 0, 1],
        [7, 0, 0, 0, 0, 1, 0] ],
      [ [1, 2, 3, 4, 5, 6, 7],
        [2, 1, 4, 3, 0, 0, 0],
        [3, 4, 1, 2, 0, 0, 0],
        [4, 3, 2, 1, 0, 0, 0],
        [5, 0, 0, 0, 1, 0, 0],
        [6, 0, 0, 0, 0, 7, 1],
        [7, 0, 0, 0, 0, 1, 6] ],
      [ [1, 2, 3, 4, 5, 6, 7],
        [2, 1, 4, 3, 0, 0, 0],
        [3, 4, 1, 2, 0, 0, 0],
        [4, 3, 2, 1, 0, 0, 0],
        [5, 0, 0, 0, 1, 7, 6],
        [6, 0, 0, 0, 7, 5, 1],
        [7, 0, 0, 0, 6, 1, 5] ],
      [ [1, 2, 3, 4, 5, 6, 7],
        [2, 1, 4, 3, 6, 5, 0],
        [3, 4, 1, 2, 0, 0, 0],
        [4, 3, 2, 1, 0, 0, 0],
        [5, 6, 0, 0, 2, 1, 0],
        [6, 5, 0, 0, 1, 2, 0],
        [7, 0, 0, 0, 0, 0, 1] ],
      [ [1, 2, 3, 4, 5, 6, 7],
        [2, 1, 4, 3, 6, 5, 0],
        [3, 5, 1, 0, 2, 0, 0],
        [4, 6, 2, 0, 1, 0, 0],
        [5, 3, 0, 1, 0, 2, 0],
        [6, 4, 0, 2, 0, 1, 0],
        [7, 0, 0, 0, 0, 0, 1] ],
      [ [1, 2, 3, 4, 5, 6, 7],
        [2, 1, 4, 3, 6, 5, 0],
        [3, 5, 1, 6, 2, 4, 0],
        [4, 6, 2, 5, 1, 3, 0],
        [5, 3, 6, 1, 4, 2, 0],
        [6, 4, 5, 2, 3, 1, 0],
        [7, 0, 0, 0, 0, 0, 1] ],
      [ [1, 2, 3, 4, 5, 6, 7],
        [2, 1, 0, 0, 0, 0, 0],
        [3, 0, 1, 0, 0, 0, 0],
        [4, 0, 0, 0, 1, 0, 0],
        [5, 0, 0, 1, 0, 0, 0],
        [6, 0, 0, 0, 0, 0, 1],
        [7, 0, 0, 0, 0, 1, 0] ],
      [ [1, 2, 3, 4, 5, 6, 7],
        [2, 1, 0, 0, 0, 0, 0],
        [3, 0, 1, 0, 0, 0, 0],
        [4, 0, 0, 5, 1, 0, 0],
        [5, 0, 0, 1, 4, 0, 0],
        [6, 0, 0, 0, 0, 0, 1],
        [7, 0, 0, 0, 0, 1, 0] ],
      [ [1, 2, 3, 4, 5, 6, 7],
        [2, 1, 0, 0, 0, 0, 0],
        [3, 0, 1, 0, 0, 0, 0],
        [4, 0, 0, 5, 1, 0, 0],
        [5, 0, 0, 1, 4, 0, 0],
        [6, 0, 0, 0, 0, 7, 1],
        [7, 0, 0, 0, 0, 1, 6] ],
      [ [1, 2, 3, 4, 5, 6, 7],
        [2, 1, 0, 0, 0, 0, 0],
        [3, 0, 1, 0, 0, 0, 0],
        [4, 0, 0, 5, 6, 7, 1],
        [5, 0, 0, 6, 7, 1, 4],
        [6, 0, 0, 7, 1, 4, 5],
        [7, 0, 0, 1, 4, 5, 6] ],
      [ [1, 2, 3, 4, 5, 6, 7],
        [2, 1, 4, 3, 0, 0, 0],
        [3, 4, 2, 1, 0, 0, 0],
        [4, 3, 1, 2, 0, 0, 0],
        [5, 0, 0, 0, 1, 0, 0],
        [6, 0, 0, 0, 0, 1, 0],
        [7, 0, 0, 0, 0, 0, 1] ],
      [ [1, 2, 3, 4, 5, 6, 7],
        [2, 1, 4, 3, 0, 0, 0],
        [3, 4, 2, 1, 0, 0, 0],
        [4, 3, 1, 2, 0, 0, 0],
        [5, 0, 0, 0, 1, 0, 0],
        [6, 0, 0, 0, 0, 0, 1],
        [7, 0, 0, 0, 0, 1, 0] ],
      [ [1, 2, 3, 4, 5, 6, 7],
        [2, 1, 4, 3, 0, 0, 0],
        [3, 4, 2, 1, 0, 0, 0],
        [4, 3, 1, 2, 0, 0, 0],
        [5, 0, 0, 0, 1, 0, 0],
        [6, 0, 0, 0, 0, 7, 1],
        [7, 0, 0, 0, 0, 1, 6] ],
      [ [1, 2, 3, 4, 5, 6, 7],
        [2, 1, 4, 3, 0, 0, 0],
        [3, 4, 2, 1, 0, 0, 0],
        [4, 3, 1, 2, 0, 0, 0],
        [5, 0, 0, 0, 1, 7, 6],
        [6, 0, 0, 0, 7, 5, 1],
        [7, 0, 0, 0, 6, 1, 5] ],
      [ [1, 2, 3, 4, 5, 6, 7],
        [2, 1, 4, 3, 6, 5, 0],
        [3, 4, 2, 1, 0, 0, 0],
        [4, 3, 1, 2, 0, 0, 0],
        [5, 6, 0, 0, 2, 1, 0],
        [6, 5, 0, 0, 1, 2, 0],
        [7, 0, 0, 0, 0, 0, 1] ],
      [ [1, 2, 3, 4, 5, 6, 7],
        [2, 1, 4, 3, 0, 0, 0],
        [3, 0, 0, 0, 1, 2, 0],
        [4, 0, 0, 0, 2, 1, 0],
        [5, 6, 1, 0, 0, 0, 0],
        [6, 5, 0, 1, 0, 0, 0],
        [7, 0, 0, 0, 0, 0, 1] ],
      [ [1, 2, 3, 4, 5, 6, 7],
        [2, 1, 4, 3, 0, 0, 0],
        [3, 0, 0, 0, 1, 2, 4],
        [4, 0, 0, 0, 2, 1, 3],
        [5, 6, 1, 7, 0, 0, 0],
        [6, 5, 7, 1, 0, 0, 0],
        [7, 0, 0, 0, 6, 5, 1] ],
      [ [1, 2, 3, 4, 5, 6, 7],
        [2, 1, 4, 3, 6, 5, 0],
        [3, 4, 0, 0, 1, 2, 0],
        [4, 3, 0, 0, 2, 1, 0],
        [5, 6, 1, 2, 0, 0, 0],
        [6, 5, 2, 1, 0, 0, 0],
        [7, 0, 0, 0, 0, 0, 1] ],
      [ [1, 2, 3, 4, 5, 6, 7],
        [2, 1, 4, 3, 6, 5, 0],
        [3, 4, 5, 6, 1, 2, 0],
        [4, 3, 6, 5, 2, 1, 0],
        [5, 6, 1, 2, 3, 4, 0],
        [6, 5, 2, 1, 4, 3, 0],
        [7, 0, 0, 0, 0, 0, 1] ],
      [ [1, 2, 3, 4, 5, 6, 7],
        [2, 0, 1, 0, 0, 0, 0],
        [3, 1, 0, 0, 0, 0, 0],
        [4, 0, 0, 0, 1, 0, 0],
        [5, 0, 0, 1, 0, 0, 0],
        [6, 0, 0, 0, 0, 0, 1],
        [7, 0, 0, 0, 0, 1, 0] ],
      [ [1, 2, 3, 4, 5, 6, 7],
        [2, 3, 1, 0, 0, 0, 0],
        [3, 1, 2, 0, 0, 0, 0],
        [4, 0, 0, 0, 1, 0, 0],
        [5, 0, 0, 1, 0, 0, 0],
        [6, 0, 0, 0, 0, 0, 1],
        [7, 0, 0, 0, 0, 1, 0] ],
      [ [1, 2, 3, 4, 5, 6, 7],
        [2, 3, 1, 0, 0, 0, 0],
        [3, 1, 2, 0, 0, 0, 0],
        [4, 0, 0, 5, 1, 0, 0],
        [5, 0, 0, 1, 4, 0, 0],
        [6, 0, 0, 0, 0, 0, 1],
        [7, 0, 0, 0, 0, 1, 0] ],
      [ [1, 2, 3, 4, 5, 6, 7],
        [2, 3, 1, 0, 0, 0, 0],
        [3, 1, 2, 0, 0, 0, 0],
        [4, 0, 0, 5, 1, 0, 0],
        [5, 0, 0, 1, 4, 0, 0],
        [6, 0, 0, 0, 0, 7, 1],
        [7, 0, 0, 0, 0, 1, 6] ],
      [ [1, 2, 3, 4, 5, 6, 7],
        [2, 0, 1, 5, 0, 0, 0],
        [3, 1, 0, 0, 4, 0, 0],
        [4, 0, 0, 0, 0, 1, 3],
        [5, 0, 0, 0, 0, 2, 1],
        [6, 0, 7, 1, 0, 0, 0],
        [7, 6, 0, 0, 1, 0, 0] ],
      [ [1, 2, 3, 4, 5, 6, 7],
        [2, 3, 1, 0, 0, 0, 0],
        [3, 1, 2, 0, 0, 0, 0],
        [4, 0, 0, 5, 6, 7, 1],
        [5, 0, 0, 6, 7, 1, 4],
        [6, 0, 0, 7, 1, 4, 5],
        [7, 0, 0, 1, 4, 5, 6] ],
      [ [1, 2, 3, 4, 5, 6, 7],
        [2, 3, 4, 5, 1, 0, 0],
        [3, 4, 5, 1, 2, 0, 0],
        [4, 5, 1, 2, 3, 0, 0],
        [5, 1, 2, 3, 4, 0, 0],
        [6, 0, 0, 0, 0, 0, 1],
        [7, 0, 0, 0, 0, 1, 0] ],
      [ [1, 2, 3, 4, 5, 6, 7],
        [2, 3, 4, 5, 6, 7, 1],
        [3, 4, 5, 6, 7, 1, 2],
        [4, 5, 6, 7, 1, 2, 3],
        [5, 6, 7, 1, 2, 3, 4],
        [6, 7, 1, 2, 3, 4, 5],
        [7, 1, 2, 3, 4, 5, 6] ]
    ];

InstallGlobalFunction("NrSmallPregroups",
function(n)
    if not IsPosInt(n) then
        Error("n has to be a positive integer");
    fi;

    if IsBound(ANATPH_small_pregroups[n]) then
        return Length(ANATPH_small_pregroups[n]);
    fi;
    return fail;
end);

InstallGlobalFunction("SmallPregroup",
function(n,i)
    local names;
    if not IsPosInt(n) then
        Error("n has to be a positive integer");
    fi;

    if IsPosInt(i) and
       IsBound(ANATPH_small_pregroups[n]) and
       IsBound(ANATPH_small_pregroups[n][i]) then
        names := Concatenation( ["1"]
                              , List([2..n], y -> Concatenation("p", String(y-1))));
        return PregroupByTable(names, ANATPH_small_pregroups[n][i]);
    else
        Error("small pregroup of size ", n, " with index ", i, " not available");
    fi;
end);

max_stars_repo_name: RussWoodroofe/walrus

max_issues_repo_name: RussWoodroofe/walrus

max_forks_repo_name: RussWoodroofe/walrus


text:
#
# Read("~/Workspace/Chevalley.gap/init.gi"); Read(Filename(home_dir,"load.gi")); Read(Filename(test_dir,"nilchv.test.gi"));
#
# Read("~/Workspace/Chevalley.gap/init.gi"); Read(Filename(test_dir,"nilchv.test.init.gi")); Read(Filename(test_dir,"nilchv.test.gi"));
#
# Read(Filename(test_dir,"nilchv.test.gi"));
#

sys:=NilpotentChv(ChevalleyAdj("F",4,GF(2)),[[1,1],[2,1]]);
Print("Created NilpotentChv ",sys,"\n");
Print("\t[",IsNilpotentChv(sys)=true,"] IsNilpotentChv(sys)=true\n");


text:
IsBound(a);
# true

Unbind(a);

IsBound(a);
# false


text:
# GAP has an improved floating-point support since version 4.5

Pi := Acos(-1.0);

r := Pi / 5.0;
d := 36;

Deg := x -> x * Pi / 180;

Sin(r);         Asin(last);
Sin(Deg(d));    Asin(last);
Cos(r);         Acos(last);
Cos(Deg(d));    Acos(last);
Tan(r);         Atan(last);
Tan(Deg(d));    Atan(last);

max_stars_repo_name: djgoku/RosettaCodeData

max_issues_repo_name: djgoku/RosettaCodeData

max_forks_repo_name: djgoku/RosettaCodeData


text:
while true do
    a := Random(0, 19);
    Print(a);
    if a = 10 then
        Print("\n");
        break;
    fi;
    a := Random(0, 19);
    Print("\t", a, "\n");
od;

# 11      6
# 5       8
# 1       4
# 5       10
# 1       16
# 10


text:
Median := function(v)
  local n, w;
  w := SortedList(v);
  n := Length(v);
  return (w[QuoInt(n + 1, 2)] + w[QuoInt(n, 2) + 1]) / 2;
end;

a := [41, 56, 72, 17, 93, 44, 32];
b := [41, 72, 17, 93, 44, 32];

Median(a);
# 44
Median(b);
# 85/2


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

# ADram tag, generic class for DRAM data layout optimizations
# params[1] = rb = row-buffer size in words
# params[2] = m = size of local memory in words
Class(ADram, AGenericTag, rec(
	updateParams := meth(self)
		Checked(Length(self.params)=2);
		self.rb := self.params[1];
		self.m := self.params[2];
	end,
));

# ACube tag: Tag meant for cubic memory mapping for DRAM,
# params[1] = k = cube-edge in words (k^3 = row-buffer size)
# params[2] = m = size of local memory in words
Class(ACube, AGenericTag, rec(
	isDMARead := false,
	isDMAWrite := false,
	updateParams := meth(self)
		Checked(Length(self.params)=2);
		self.k := self.params[1];
		self.m := self.params[2];
	end,
));

Class(ACubeRd, AGenericTag, rec(
	isDMARead := true,
	isDMAWrite := false,
	updateParams := meth(self)
		Checked(Length(self.params)=2);
		self.k := self.params[1];
		self.m := self.params[2];
	end,
));

Class(ACubeWr, AGenericTag, rec(
	isDMARead := false,
	isDMAWrite := true,
	updateParams := meth(self)
		Checked(Length(self.params)=2);
		self.k := self.params[1];
		self.m := self.params[2];
	end,
));

# ATile tag: Tag meant for 2D-tiled memory mapping for DRAM,
# params[1] = k = tile-width in words (k^2 = row-buffer size)
# params[2] = m = size of local memory in words
Class(ATile, AGenericTag, rec(
	isDMARead := false,
	isDMAWrite := false,
	updateParams := meth(self)
		Checked(Length(self.params)=2);
		self.k := self.params[1];
		self.m := self.params[2];
	end,
));

Class(ATileRd, ATile, rec(
	isDMARead := true,
	isDMAWrite := false,
	updateParams := meth(self)
		Checked(Length(self.params)=2);
		self.k := self.params[1];
		self.m := self.params[2];
	end,
));

Class(ATileWr, ATile, rec(
	isDMARead := false,
	isDMAWrite := true,
	updateParams := meth(self)
		Checked(Length(self.params)=2);
		self.k := self.params[1];
		self.m := self.params[2];
	end,
));

#  # size = scratchpad segment size
#  # nsgmts = number of segments, defaults to 1.
#  Class(ALStore, AGenericTag, rec(
#      updateParams := meth(self)
#          Checked(Length(self.params)=3);
#          self.size := self.params[1];
#          self.nsgmts := self.params[2];
#          self.linesize := self.params[3];
#      end,
#      isRegCx := false
#  ));
#
#  Class(ALStoreCx, ALStore, rec(
#      updateParams := meth(self)
#          Checked(Length(self.params)=3);
#          self.size := self.params[1]/2;
#          self.nsgmts := self.params[2];
#          self.linesize := self.params[3];
#      end,
#      isRegCx := true
#  ));
#
#  # APad tag.
#  #
#  # APad is a special buffering tag meant for things like scratchpads, hence
#  # the name. It takes four parameters. In order, they are:
#  # b - the block size in number of elements
#  # s - the segment size, in number of blocks
#  # u - the number of segments
#  # n - a string identifier
#  #
#  # The block size is the SMALLEST allowed transfer size when copying data.
#  # The breakdown rules which propagate the APad tag insure that sub-block
#  # numbers of contiguous elements are never moved when the tag is present.
#  #
#  # A segment is a discrete memory separate from other segments, in the case
#  # that the 'u' parameter is >1. In the case of DPA, each segment is
#  # connected to a different compute processor, and we do parallelization
#  # after the APad tag is dropped.
#  #
#  # The string identifier is for the platform writer. You can label things
#  # like "local memory" or "vector register file."
#  #
#  # Also, the software pipelining loop (rather than a standard ISum) is
#  # automatically used with the APad tag.
#
#  Class(APad, AGenericTag, rec(
#      applied := false,
#
#      b := (self) >> self.params[1],
#      s := (self) >> self.params[2],
#      u := (self) >> self.params[3],
#      bs := (self) >> self.params[1] * self.params[2],
#      bsu := (self) >> Product(DropLast(self.params, 1)),
#      n := (self) >> self.params[4],
#
#      apply := self >> CopyFields(self, rec(applied := true)),
#
#      updateParams := meth(self)
#          Checked(Length(self.params)=4);
#      end,
#
#      print := meth(self)
#          Print(self.name, "(", PrintCS([self.b(), self.s(), self.u()]),
#              ", \"", self.n(), "\")");
#      end
#  ));


text:
#
# AoC: Advent of Code solutions in GAP
#
# Implementations
#
AoC.Year2019 := rec();

InstallMethod( Puzzle,
    [IsString, IsString, IsFunction, IsFunction],
    function( Year, Day, PartOne, PartTwo )
        local puzzle;
        puzzle := rec(
            Input := function( )
                return InputTextFile(
                    Filename(
                        DirectoriesPackageLibrary( "AoC", StringFormatted( "input/{1}", Year ) )[1],
                        StringFormatted( "day{1}.txt", Day )
                    )
                );
            end,
            PartOne := PartOne,
            PartTwo := PartTwo
        );
        Objectify( TYPE_PUZZLE, puzzle );
        return puzzle;
    end );
InstallMethod( Input,
  "of a puzzle",
  [ IsPuzzle ],
  puzzle -> puzzle!.Input() );
InstallMethod( PartOne,
  "of a puzzle",
  [ IsPuzzle ],
  puzzle -> puzzle!.PartOne );
InstallMethod( PartTwo,
  "of a puzzle",
  [ IsPuzzle ],
  puzzle -> puzzle!.PartTwo );

max_stars_repo_name: yurrriq/advent-of-code