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

# ==========================================================================
# SPLSymbol(<name>, <params>)  - obsolete function for creating spl symbols
#
# Note: in this revision all symbol names can be used as functions.
#       SPLSymbol("F", 2) is now F(2)
#
# ==========================================================================
SPLSymbol := function(s, params)
    local sym, perm;
    if not IsList(params) then params := [params]; fi;
    sym  := IsBound(spl.(s));
    if not sym then
	Error("Symbol '", s, "' does not exist. See Dir(spl).");
    else
	return ApplyFunc(spl.(s), params);
    fi;
end;

#F IsSPLSym(<x>) - returns true for parametrized symbols (F, I, ...)
#F
IsSPLSym := x -> IsRec(x) and IsBound(x._sym) and x._sym;

# ==========================================================================
# Sym - base class for SPL symbols
# ==========================================================================
Class(Sym, BaseMat, rec(
    _sym := true,
    transposed := false,
    type := "Sym",

    new := (self, s, obj) >>
        SPL(Inherit(self, rec( symbol := s,
                   obj := obj,
                   dimensions := obj.dims()))),

   # -------- Transformation rules support ---------------------------------
    rChildren := self >> When(IsList(self.params), self.params, [self.params]),
    rSetChild := meth(self, n, newChild)
        if IsList(self.params) then
            self.params[n] := newChild;
        else
            if n<>1 then Error("<n> must be 1"); fi;
            self.params := newChild;
        fi;
        # self.canonizeParams(); ??
        self.dimensions := self.dims();
    end,

    from_rChildren := (self, rch) >> let(
	res := ApplyFunc(ObjId(self), rch).appendAobj(self),
	When(self.transposed, res.transpose(), res)),
    # ----------------------------------------------------------------------

    checkParams := meth(self, params)
        local nargs, nump;
	nargs := NumArgs(self.def);
	nump := When(IsList(params), Length(params), 1);
        if nargs <> -1 and nargs <> nump then
	    Error("Symbol '", self.name, "' needs ", NumArgs(self.def), " parameters: ",
		ParamsFunc(self.def), "\n");
	fi;
    end,

    canonizeParams := meth(self, params)
        local A, nump;
	nump := Length(params);
        if IsBound(self.abbrevs) then
	    for A in self.abbrevs do
                if NumArgs(A) = -1 or NumArgs(A) = nump then
		    return ApplyFunc(A, params);
		fi;
	    od;
	    return params;
	else
	    return params;
	fi;
    end,

    fromDef := meth(arg)
        local result, self, params, h, lkup;
	self := arg[1];
	params := arg{[2..Length(arg)]};
	params := self.canonizeParams(params);
        self.checkParams(params);
	params := When(IsList(params), params, [params]);

	h := self.hash;
	if h<>false then
	    lkup := h.objLookup(self, params);
	    if lkup[1] <> false then return lkup[1]; fi;
	fi;

	result := SPL(WithBases(self, rec(params := params, transposed := false)));
	result.obj := ApplyFunc(result.def, params);
	result.dimensions := Dimensions(result.obj);

	if h<>false then return h.objAdd(result, lkup[2]);
	else return result;
	fi;
    end,

    __call__ := ~.fromDef,

    getObj := self >> When(self.transposed, self.obj.transpose(), self.obj),

    dims          := self >> self.getObj().dims(),
    isPermutation := self >> self.obj.isPermutation(),
    isReal        := self >> self.obj.isReal(),
    toAMat        := self >> self.getObj().toAMat(),
    free          := self >> Union(List(self.params, FreeVars)),

    #-----------------------------------------------------------------------
    transpose := self >>
        Inherit(self, rec(transposed := not self.transposed,
                  dimensions := Reversed(self.dimensions))),
    #-----------------------------------------------------------------------
    equals := meth(self, other)
        return ObjId(other) = ObjId(self)
           and other.params = self.params
           and self.transposed = other.transposed
	   and self.a = other.a;
    end,
    #-----------------------------------------------------------------------
    print := (self, i, is) >> let(
        params := When(IsList(self.params), self.params, [self.params]),
        Print(self.__name__, "(", PrintCS(params), ")",
              When(self.transposed, ".transpose()", ""),
              self.printA())),
    #-----------------------------------------------------------------------
    arithmeticCost := (self, costMul, costAddMul) >>
        self.getObj().arithmeticCost(costMul, costAddMul)
));

#F SymFunc(<name>, <func>)
#F      Convert a function (which returns SPL) to SPL symbol.
#F      Resulting symbol is automatically added to spiral.spl.symbols namespace.
SymFunc := function(name, func)
    local sym;
    sym := rec(notSupported := true, def := func);
    sym.name := name;
    sym.operations := ClassOps;
    sym.__bases__ := [ Sym ];
    CantCopy(sym);
    return sym;
end;


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

# ================================================================
# Diagonals
# ================================================================
Class(RulesDiag, RuleSet);

LeftPull  := [Diag, DiagCpxSplit, RCDiag];
RightPull := [Diag, DiagCpxSplit, RCDiag];
PullCont  := [SUM, BB, Buf, Inplace, Grp];
PullContL  := [SUM, BB, Buf, Inplace, Grp, NoDiagPullinRight];
PullContR  := [SUM, BB, Buf, Inplace, Grp, NoDiagPullinLeft];

RewriteRules(RulesDiag, rec(
 DiagPullInRight := ARule( Compose,
       [ @(1, LeftPull), @(2, PullContL) ],
  e -> [ ApplyFunc(ObjId(@2.val), List(@2.val.children(), c -> @1.val * c)) ]),

 DiagPullInLeft  := ARule( Compose,
       [ @(1, PullContR), @(2, RightPull) ],
  e -> [ ApplyFunc(ObjId(@1.val), List(@1.val.children(), c -> c * @2.val)) ]),

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

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

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

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

 DiagDataRight := ARule( Compose,
       [ @(1, LeftPull), @(2, Data, canReorder) ],
  e -> [ Data(@2.val.var, @2.val.value, @1.val * @2.val.child(1)).attrs(@(2).val) ]),

 DiagDataLeft := ARule( Compose,
       [ @(1, Data, canReorder), @(2, RightPull) ],
  e -> [ Data(@1.val.var, @1.val.value, @1.val.child(1) * @2.val).attrs(@(1).val) ]),

 DiagRecursStepRight := ARule( Compose,
       [ @(1, LeftPull), @(2, RecursStep) ],
  e -> [ RecursStep(@2.val.yofs, @2.val.xofs, @1.val * @2.val.child(1)) ]),

 DiagRecursStepLeft := ARule( Compose,
       [ @(1, RecursStep), @(2, RightPull) ],
  e -> [ RecursStep(@1.val.yofs, @1.val.xofs, @1.val.child(1) * @2.val) ])
));

Class(RulesDiagStandalone, RuleSet);
RewriteRules(RulesDiagStandalone, rec(
 # Gath * Diag
 CommuteGathDiag := ARule( Compose,
       [ @(1, Gath), @(2, Diag) ], # o 1-> 2->
  e -> [ Diag(fCompose(@2.val.element, @1.val.func)).attrs(@(2).val), @1.val ]),

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

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

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

 # Gath * RCDiag
 CommuteGathRCDiag := ARule( Compose,
       [ [@(1, Gath), [@(0,fTensor), ..., [fId,@(2).cond(IsEvenInt)]]],
      @(4, RCDiag) ],
  e -> [ RCDiag(fCompose(@(4).val.element, @(0).val), @(4).val.post),
         @(1).val ]),

 # RCDiag * Scat
 CommuteRCDiagScat := ARule( Compose,
       [ @(4, RCDiag),
     [@(1, Scat), [@(0,fTensor), ..., [fId,@(2).cond(IsEvenInt)]]] ],
  e -> [ @(1).val,
         RCDiag(fCompose(@(4).val.element, @(0).val), @(4).val.post) ]),

 # Gath * IDirSum
 CommuteGathIDirSum := ARule( Compose,
       [ [@(1, Gath), [@(3,fTensor), ..., [fId,@(2)]]],
     [@(0, [PushL, PushLR]), @(4, IDirSum, e -> (@(2).val mod Rows(e.child(1))) = 0)] ],
  e -> let(ch := DropLast(@(3).val.children(), 1),
           ratio := @(2).val / Rows(@(4).val.child(1)),
           f := When(ratio=1, fTensor(ch), fTensor(Concatenation(ch, [fId(ratio)]))),
           j := @(4).val.var,
           jj := Ind(f.domain()),

           [ ObjId(@(0).val)(IDirSum(jj,      #SubstVars(@(4).val.child(1), rec((j.id) := f.at(jj)))),
                   Data(j, f.at(jj), SubstTopDown(@(4).val.child(1), @(5,fBase,e->e.params[2]=j),
                                       e -> fCompose(f, fBase(jj)))))),
             @(1).val ])
 ),

 # IDirSum * Scat
 CommuteIDirSumScat := ARule( Compose,
       [ [@(0,[PushR, PushLR]), @(4, IDirSum)],
         [@(1, Scat), [@(3,fTensor), ..., [fId,@(2).cond(e -> (e mod Cols(@(4).val.child(1))) = 0)]]]],

  e -> let(ch := DropLast(@(3).val.children(), 1),
           ratio := @(2).val / Cols(@(4).val.child(1)),
           f := When(ratio=1, fTensor(ch), fTensor(Concatenation(ch, [fId(ratio)]))),
           j := @(4).val.var,
           jj := Ind(f.domain()),

           [ @(1).val,
             ObjId(@(0).val)(IDirSum(jj,      #SubstVars(@(4).val.child(1), rec((j.id) := f.at(jj)))),
                   Data(j, f.at(jj), SubstTopDown(@(4).val.child(1), @(5,fBase,e->e.params[2]=j),
                                       e -> fCompose(f, fBase(jj)))))) ])
  ),

));

#--------------------------------------------------
#   new stuff for II

Class(RulesDiagII, RuleSet);
RewriteRules(RulesDiagII, rec(
    fCompose_diagXX := Rule([fCompose, @(1, [diagAdd, diagMul]), @(2)],
        e -> ApplyFunc(ObjId(@(1).val), List(@(1).val.children(), c->(fCompose(c, @(2).val ))))),

    simp1 := Rule([fCompose, @(1, II, e->e.params[2]=0), [fTensor, @(2, fId), @(3, fBase, e->@(1).val.params[3] = e.range() )] ],
        e -> II(@(2).val.domain(), 0, 1)),

    fBase_cond := Rule([fCompose, @(1, II), @(2, fBase)],
        e -> let(ii := @(1).val.params,
            diagCond(ivElOf(ivII(ii[1], ii[2], ii[3]), @(2).val.params[2]), fConst(1,1), fConst(1,0)))),

    tensor_cond := ARule(diagTensor, [@(1), @(2, diagCond)],
        e -> [diagCond(@(2).val.child(1), diagTensor(@(1).val, @(2).val.child(2)), diagTensor(@(1).val, @(2).val.child(3)))]),

    ivElOf_fold_val := Rule([ivElOf, @(1, ivII), @(2, Value)],
        e-> When(@2.val.v >= @1.val.params[2] and @2.val.v < @1.val.params[3], V(1), V(0))),

    diagCond_fold := Rule([diagCond, @(1, Value), @(2), @(3)],
        e-> When(@1.val.v = 1, @(2).val, @(3).val)),

    diagAdd_combine1 := ARule(diagAdd, [@(1, II), @(2, II, e->@(1).val.params[3]=e.params[2])],
        e -> [II(@(1).val.params[1], @(1).val.params[2], @(2).val.params[3])]),

    diagAdd_combine2 := ARule(diagAdd, [@(1, II), @(2, II, e->@(1).val.params[2]=e.params[3])],
        e -> [II(@(1).val.params[1], @(2).val.params[2], @(1).val.params[3])]),

    II_to_fConst := Rule(@(1, II, e->e.params[1]=e.params[3] and e.params[2]=0), e->fConst(@(1).val.params[1], 1)),

    diagMul0r := ARule(diagMul, [@(1), @(2, fConst, e->e.params[2] = 0)], e -> [@(2).val]),
    diagMul0l := ARule(diagMul, [@(1, fConst, e->e.params[2] = 0), @(2)], e -> [@(1).val]),

    diagAdd0r := ARule(diagAdd, [@(1), @(2, fConst, e->e.params[2] = 0)], e -> [@(1).val]),
    diagAdd0l := ARule(diagAdd, [@(1, fConst, e->e.params[2] = 0), @(2)], e -> [@(2).val]),

    Diag_fConst_to_I := Rule([@(1, Diag), @(2, fConst, e->e.params[2] = 1)],
        e->Prm(fId(@(2).val.params[1]))),

    adjust_ivII := Rule ([ivElOf, @(1, ivII), [@(2, add), @(3, var), @(4, Value)]],
        e -> let(iv := @(1).val.params, shift := @(4).val.v,
            ivElOf(ivII(iv[1], Maximum(0, iv[2]-shift), Maximum(0, iv[3]-shift)), @(3).val))),

    ivElOf_fold_var1 := Rule([ivElOf, @(1, ivII, e->e.params[2]=0),
        @(2, var, e->IsBound(e.range) and e.range <= @(1).val.params[3])],
        e-> V(1)),

    ivElOf_fold_var2 := Rule([ivElOf, @(1, ivII, e->e.params[3]=0), @(2, var)],
        e-> V(0)),
));

RulesDiagAll := MergedRuleSet(RulesDiag, RulesDiagStandalone);

ext: gi


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

ImportAll(paradigms.smp);
ImportAll(paradigms.vector);

Class(IAGlobals, rec(
    getOpts := meth(arg)
        local opts, self, smpopts, isa, sseopts, optrec, tid;

        #   set up defaults
        optrec := rec(
            cpu := LocalConfig.cpuinfo,
            useSIMD := true,
            useSMP := true,
            useNewSettingsForICC := false,
            dataType := T_Real(64),
            globalUnrolling := 128,
            useArea := false,
            processIntTables := false,
            use64bit := LocalConfig.osinfo.is64bit()
        );

        smpopts := rec(
                numproc := LocalConfig.cpuinfo.cores,
                api := "OpenMP"
        );

        sseopts := rec(
            svct:=true,
            splitL:=false,
            oddSizes := false,
            stdTTensor := true,
            tsplPFA := false
        );

        self := arg[1];
        arg := Flat(Drop(arg, 1));

        if IsList(arg[1]) then arg := arg[1]; fi;
        if Length(arg) >= 1 then optrec := CopyFields(optrec, arg[1]); fi;
        if Length(arg) >= 2 then smpopts := CopyFields(smpopts, arg[2]); fi;
        if Length(arg) >= 3 then sseopts := CopyFields(sseopts, arg[3]); fi;

        if optrec.useSIMD then
            #   handle SSE
            #   let SSE work the unrolling magic
            sseopts.globalUnrolling := optrec.globalUnrolling;
            sseopts.useArea := optrec.useArea;
            sseopts.processIntTables := optrec.processIntTables;
            if IsBound(optrec.mode) then sseopts.mode := optrec.mode; fi;

            isa := optrec.cpu.getSimdIsa(optrec.dataType);
            opts := CopyFields(SIMDGlobals.getOpts(isa, sseopts), rec(
                IAconf := rec(
                    optrec := optrec,
                    smpopts := smpopts,
                    sseopts := sseopts
                ),

                unparser := When(smpopts.api = "OpenMP",
                    When(isa in [AVX_8x32f, AVX_4x64f],
                        When(IsBound(smpopts.OmpMode) and smpopts.OmpMode = "for",
                            spiral.libgen.OpenMP_AVXUnparser_ParFor,
                            spiral.libgen.OpenMP_AVXUnparser),
                        When(IsBound(smpopts.OmpMode) and smpopts.OmpMode = "for",
                            spiral.libgen.OpenMP_SSEUnparser_ParFor,
                            spiral.libgen.OpenMP_SSEUnparser)),
                    spiral.libgen.SMP_SSEUnparser),
                codegen := spiral.libgen.VecRecCodegen
            ));

            if optrec.useNewSettingsForICC then
                if IsBound(opts.language) then Unbind(opts.language); fi;
                if LocalConfig.osinfo.isLinux() then
                    if optrec.use64bit then
                        opts.profile := default_profiles.linux_x64_icc;
                    else
                        opts.profile := default_profiles.linux_x86_icc;
                    fi;
                elif LocalConfig.osinfo.isDarwin() then
                    Error("No Darwin SMP profiles defined yet");
                elif LocalConfig.osinfo.isWindows() then
                    if optrec.use64bit then
                        opts.profile := default_profiles.win_x64_icc;
                    else
                        opts.profile :=default_profiles.win_x86_icc;
                    fi;
                else
                    Error("Unknow OS");
                fi;
            fi;

        else
            if optrec.useSMP then Error("scalar SMP not implemented"); fi;

            opts := Copy(SpiralDefaults);

            opts.globalUnrolling := optrec.globalUnrolling;
            opts.useArea := optrec.useArea;
            opts.processIntTables := optrec.processIntTables;

            if optrec.dataType = T_Real(32) then
                opts := InitDataType(opts, "f32re");
            elif optrec.dataType = T_Real(64) then
                opts := InitDataType(opts, "f64re");
            else
                Error("unknown data type");
            fi;

            if optrec.useNewSettingsForICC then
                if IsBound(opts.language) then Unbind(opts.language); fi;
                if LocalConfig.osinfo.isLinux() then
                    if optrec.use64bit then
                        opts.profile := default_profiles.linux_x64_icc;
                    else
                        opts.profile := default_profiles.linux_x86_icc;
                    fi;
                elif LocalConfig.osinfo.isDarwin() then
                    Error("No Darwin SMP profiles defined yet");
                elif LocalConfig.osinfo.isWindows() then
                    if optrec.use64bit then
                        opts.profile := default_profiles.win_x64_icc;
                    else
                        opts.profile :=default_profiles.win_x86_icc;
                    fi;
                else
                    Error("Unknow OS");
                fi;
            fi;

            if optrec.useArea then
                opts.globalUnrolling := 2 * opts.globalUnrolling; # * Log2Int(opts.globalUnrolling);
                opts.markBlock := MarkBlocksAreaSums;
            fi;

            return opts;
        fi;

        # handle the SMP pthreads/OpenMP case
        if optrec.useSMP then
            opts.breakdownRules.GT := [ GT_Base, GT_NthLoop,
                CopyFields(GT_Par, rec(parEntireLoop := false, splitLoop := true)), GT_Par_odd,
                GT_Vec_AxI, GT_Vec_IxA, GT_Vec_IxA_Push, GT_Vec_IxA_L, GT_Vec_L_IxA
            ];
            opts.breakdownRules.TTensorI := Concat([
            CopyFields(TTensorI_toGT, rec(
                applicable := (self, t) >> t.hasTags() and ObjId(t.getTags()[1])=AParSMP ))],
                opts.breakdownRules.TTensorI);

            opts.breakdownRules.TTensorInd :=
                Concat([dsA_base_smp, dsA_smp, L_dsA_L_base_smp, L_dsA_L_smp],
                    opts.breakdownRules.TTensorInd);

            tid := When(smpopts.api = "OpenMP", threadId(), var("tid", TInt));
            opts.tags := Concat([ AParSMP(smpopts.numproc, tid)  ], opts.tags);

            if optrec.useNewSettingsForICC then
                if IsBound(opts.language) then Unbind(opts.language); fi;
                if LocalConfig.osinfo.isLinux() then
                    if optrec.use64bit then
                        opts.profile := When(smpopts.api = "OpenMP", default_profiles.linux_x64_icc_openmp, default_profiles.linux_x64_threads);
                    else
                        opts.profile := When(smpopts.api = "OpenMP", default_profiles.linux_x86_icc_openmp, default_profiles.linux_x86_threads);
                    fi;
                elif LocalConfig.osinfo.isDarwin() then
                    Error("No Darwin SMP profiles defined yet");
                elif LocalConfig.osinfo.isWindows() then
                    if optrec.use64bit then
                        opts.profile := When(smpopts.api = "OpenMP", default_profiles.win_x64_icc_openmp, default_profiles.win_x64_icc_threads);
                    else
                        opts.profile := When(smpopts.api = "OpenMP", default_profiles.win_x86_icc_openmp, default_profiles.win_x86_icc_threads);
                    fi;
                else
                    Error("Unknow OS");
                fi;
            else
                opts.language := When(smpopts.api = "OpenMP", optrec.cpu.OpenMP_lang, optrec.cpu.smp_lang);
            fi;

            if smpopts.api = "threads" then opts.subParams := [var("num_threads", TInt), var("tid", TInt)]; fi;
            opts.smp := smpopts;
        fi;

        Add(opts.includes, "<include/mm_malloc.h>");
        if not IsBound(opts.globalUnrolling) then opts.globalUnrolling := optrec.globalUnrolling; fi;
        opts.operations := rec(Print := (s) -> Print("<IA options>"));
        return opts;
    end
));


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

#-----------------------------------------------------------------------------
# Cell SMP Paradigm (preferred)
#-----------------------------------------------------------------------------

#F SymSPL = Symmetric SPL. In this particular case, we know about the symmetry
#F of this SymSPL, so we can remove it.
Class(KillSymSPL, RuleSet);
RewriteRules(KillSymSPL, rec(
     kill_symspl := Rule([@(1,SymSPL), ISum], e -> @(1).val.child(1)),
));

NewRulesFor(GT, rec(
#F For cases where p doesn't divide n nicely. (eg: Parallelize (I6 x A) across 4 processors
    GT_Cell_nonmultiple := rec(
        maxSize       := false,
        minSize       := false,

        requiredFirstTag := ParCell,

        applicable := (self, t) >> let(
                rank   := Length(t.params[4]),
                procs  := t.firstTag().params[1],
                rank = 1
                  and let(its := t.params[4][1],
                    PatternMatch(t, [GT, @(1), @(2,XChain), @(3,XChain), ...], empty_cx())
                    # Only do IxA for now
                    and t.params[2] = GTPar and t.params[3] = GTPar
                    and IsPosInt(idiv(its, procs).v)
                    and IsPosInt(imod(its, procs).v)
                  )
        ),
        children := (self, t) >> let(
            spl     := t.params[1],
            g       := t.params[2],
            s       := t.params[3],
            its     := t.params[4][1],
            procs   := t.firstTag().params[1],
            pksize  := t.firstTag().params[2],
            remits  := imod(its, procs).v,
            parits  := its - remits,
            origtags:= t.getTags(),
            partag  := t.firstTag(),
            tags    := Drop(t.getTags(), 1),
            When(remits > 1,
                [ [ GT(spl, g, s, [parits]).withTags(origtags), GT(spl, g, s, [remits]).setTags(Concatenation(ParCell(remits, pksize), tags)) ] ],
                [ [ GT(spl, g, s, [parits]).withTags(origtags), spl.withTags(tags) ] ]
            )
        ),

        apply := (self, t, C, Nonterms) >> let(
                 DirectSum(C[1], C[2])
        ) #apply
    ), #GT_Cell_nonmultiple

    GT_Cell := rec(
        maxSize       := false,
        minSize       := false,

        requiredFirstTag := ParCell,

        applicable := (self, t) >> let(
                rank   := Length(t.params[4]),
                procs  := t.firstTag().params[1],
                pkSize := t.firstTag().params[2],
                rank = 1
                  and let(its := t.params[4][1],
                    PatternMatch(t, [GT, @(1), @(2,XChain), @(3,XChain), ...], empty_cx())
                    # The (IxA) construct doesn't need the pkSize constraint that
                    # the (AxI) construct (and (IxA)L, L(IxA)) does:
                    and When(t.params[2] = GTPar and t.params[3] = GTPar,
                                IsPosInt(its/(procs)),
                                IsPosInt(its/(procs*pkSize))
                        )
                  )
        ),

        children := (self, t) >> let(
            spl    := t.params[1],
            g      := t.params[2],
            s      := t.params[3],
            its    := t.params[4][1],
            procs  := t.firstTag().params[1],
            remits := its/procs,
            tags   := Drop(t.getTags(), 1), # Eat up tag at topmost level. NOTE: Nested parallelism not considered for now
            When(remits > 1,
            [ [ GT(spl, g, s, [its / procs]).withTags(tags), ] ],
            [ [ spl.withTags(tags) ] ])
        ),

        apply := (self, t, C, Nonterms) >> let(
            spl     := t.params[1],
            g       := t.params[2],
            s       := t.params[3],
            its     := t.params[4][1],
            gg      := When(g=GTVec, XChain([0,1,2]), XChain([1,2,0])),
            ss      := When(s=GTVec, XChain([0,1,2]), XChain([1,2,0])),
            procs   := t.firstTag().params[1],
            pkSize  := t.firstTag().params[2],
            i       := var("spuid", TInt, procs),
            #z := Error("Breakpoint\n"),

            # NOTE: There is a degree of freedom in allocating packet sizes for
            # ScatDist and GathDist. Here's how we handle things:

            # When we do on-chip exchanges, we "combine" ScatSends with
            # GathDists (and vice versa). So we want the packet size of the
            # ScatSends and GathDists (and vice versa) to match across multiple
            # factors of the transform (like the DFT). So we issue a ScatDist
            # here with a matching packet size.

            # However, when we do off-chip stuff, we compose ScatDists with
            # ScatMems instead. Here, we want the highest packet size possible
            # for our ScatDists and GathDists.

            # We need to find a universal solution that works neatly for
            # everything. Perhaps issue the largest available pkSize for
            # ScatDist here, and then step it down if needed in the rewrite
            # rules (since we know it's a ScatDist and that this is possible).

            cscatter:= When(s = GTPar,
                        #ScatDist(ss.part(1, i, Rows(spl), [procs, its/(procs*pkSize)]).range(), pkSize, procs, i),
                        ScatDist(procs, its*Rows(spl)/procs, procs, i),
                        ScatSend(ss.part(1, i, Rows(spl), [procs, its/(procs*pkSize)]),         pkSize, procs, i)),
            cgather := When(g = GTPar,
                        #GathDist(gg.part(1, i, Cols(spl), [procs, its/(procs*pkSize)]).range(), pkSize, procs, i),
                        GathDist(procs, its*Rows(spl)/procs, procs, i),
                        GathRecv(gg.part(1, i, Cols(spl), [procs, its/(procs*pkSize)]),         pkSize, procs, i)),

            DistSum(procs, i, procs, cscatter * C[1] * cgather)
        ) #apply
    ), #GT_Cell

    GT_Cell_auto := rec(
        minPkSize := 4, # Hardcoded for real single precision (4 real elements == 16 bytes)
        requiredFirstTag := ParCell_auto,

        applicable := (self, t) >> let(
                rank   := Length(t.params[4]),
                procs  := t.firstTag().params[1],
                rank = 1
                  and let(its := t.params[4][1],
                    PatternMatch(t, [GT, @(1), @(2,XChain), @(3,XChain), ...], empty_cx())
                    # The (IxA) construct doesn't need the pkSize constraint that
                    # the (AxI) construct (and (IxA)L, L(IxA)) does:
                    and When(t.params[2] = GTPar and t.params[3] = GTPar,
                                IsPosInt(its/(procs)),
                                IsPosInt(its/(procs * self.minPkSize)) #Need a packet size of at least 4
                        )
                  )
        ),

        children := (self, t) >> let(
            spl    := t.params[1],
            g      := t.params[2],
            s      := t.params[3],
            n      := t.params[4][1],
            p      := t.firstTag().params[1],
            r      := n/p,
            tags   := Drop(t.getTags(), 1), # Eat up tag at topmost level. NOTE: Nested parallelism not considered for now
            When(r > 1,
            [ [ GT(spl, g, s, [n/p]).withTags(tags) ] ],
            [ [ spl.withTags(tags) ] ])
        ),

        apply := (self, t, C, Nonterms) >> let(
#           z := Error("Breakpoint\n"),
            spl     := t.params[1],
            g       := t.params[2],
            s       := t.params[3],
            n       := t.params[4][1],
            gg      := When(g.params[1]=[0,1], XChain([0,1,2]), XChain([1,2,0])),
            ss      := When(s.params[1]=[0,1], XChain([0,1,2]), XChain([1,2,0])),
            p       := t.firstTag().params[1],
            pkSize  := n/p,       # Automatically use up the remaining stuff for pkSize
                       # NOTE: pkSize is WRONG! Will be different for Par and Vec!
            i       := var("spuid", TInt, p),
            cscatter:= When(s.params[1]=[1,0],
                        ScatDist(ss.part(1, i, Rows(spl), [p, n/(p*pkSize)]).range(), pkSize, p, i),
                        ScatSend(ss.part(1, i, Rows(spl), [p, n/(p*pkSize)]),         pkSize, p, i)),
            cgather := When(g.params[1]=[1,0],
                        GathDist(gg.part(1, i, Cols(spl), [p, n/(p*pkSize)]).range(), pkSize, p, i),
                        GathRecv(gg.part(1, i, Cols(spl), [p, n/(p*pkSize)]),         pkSize, p, i)),

            DistSum(p, i, p, cscatter * C[1] * cgather)
        ) #apply
    ), #GT_Cell
));

#-----------------------------------------------------------------------------
# Cell DMP Paradigm
#-----------------------------------------------------------------------------
NewRulesFor(GT, rec(
    # Converts (In x A) to A DistSum (parallel loop)
    GT_CellDMP_base_old := rec(
        maxSize       := false,
        minSize       := false,

        requiredFirstTag := ParCellDMP_old,

        applicable := (self, t) >> let(
                rank   := Length(t.params[4]),
                procs  := t.firstTag().params[1],
                rank = 1
                  and let(its := t.params[4][1],
                    PatternMatch(t, [GT, @(1), @(2,XChain), @(3,XChain), ...], empty_cx())
                    and IsPosInt(its/(procs)))
                  and (t.params[2]=GTPar and t.params[3]=GTPar)

        ),

        children := (self, t) >> let(
            spl    := t.params[1],
            g      := t.params[2],
            s      := t.params[3],
            its    := t.params[4][1],
            procs  := t.firstTag().params[1],
            tags  := Drop(t.getTags(), 1),  # Children don't need tags since we parallelize only at the topmost level
            [ [ GT(spl, g, s, [its / procs]).withTags(tags), InfoNt(procs) ] ]
        ),

        apply := (self, t, C, Nonterms) >> let(
            spl     := t.params[1],
            g       := t.params[2],
            s       := t.params[3],
            its     := t.params[4][1],
            gg      := When(g.params[1]=[0,1], XChain([0,1,2]), XChain([1,2,0])),
            ss      := When(s.params[1]=[0,1], XChain([0,1,2]), XChain([1,2,0])),
            procs   := t.firstTag().params[1],
            i       := var("spuid", TInt, procs),
            cscatter:= ScatDist(ss.part(1, i, Rows(spl), [procs, its/(procs)]).range(), 1, procs, i),
            cgather := GathDist(gg.part(1, i, Cols(spl), [procs, its/(procs)]).range(), 1, procs, i),

            DistSum(procs, i, procs, cscatter * C[1] * cgather)

        ) #apply
    ), #GT_CellDMP_Base

    # Converts all other TTensor constructs (IxA)L, L(IxA), (AxI) into L's and (IxA)s
    # The L's are factorized to PTensor * Comm_Cell * PTensor
    GT_CellDMP_gen_old := rec(
        maxSize    := false,
        minSize    := false,
        forTransposition := false,

        requiredFirstTag := ParCellDMP_old,

        applicable := (self, t) >> let(
                rank   := Length(t.params[4]),
                procs  := t.firstTag().params[1],
                rank = 1
                  and let(its := t.params[4][1],
                    PatternMatch(t, [GT, @(1), @(2,XChain), @(3,XChain), ...], empty_cx())
                    and IsPosInt(its/(procs))
                    )
                  and not(t.params[2]=GTPar and t.params[3]=GTPar)
             ),

        children := (self, t) >> let(
            spl    := t.params[1],
            g      := t.params[2],
            s      := t.params[3],
            its    := t.params[4][1],
            procs  := t.firstTag().params[1],
            tags   := t.getTags(),
            m      := Rows(spl),
            n      := its,

            Cond(g = GTVec and s = GTPar, # L(IxA) case
                  [ [
                      TCompose([
                         GT(spl, GTPar, GTPar, t.params[4]),
                         TL(m*n, n)
                      ]).withTags(tags),
                  ] ],

                  g = GTPar and s = GTVec, # (IxA)L case
                  [ [
                      TCompose([
                         TL(m*n, m),
                         GT(spl, GTPar, GTPar, t.params[4]),
                      ]).withTags(tags),
                  ] ],

                  g = GTVec and s = GTVec, # AxI case
                  [ [
                      TCompose([
                         TL(m*n, m),
                         GT(spl, GTPar, GTPar, t.params[4]),
                         TL(m*n, n)
                      ]).withTags(tags),
                  ] ]
            )
        ),

        apply := (self, t, C, Nonterms) >> C[1]
    ) #GT_CellDMP_gen

));

NewRulesFor(TL, rec(
   TL_CellDMP_old := rec(
       maxSize       := false,
       minSize       := false,
       forTransposition := false,

       requiredFirstTag := ParCellDMP_old,

       applicable := (self, t) >> let(
           mn := t.params[1],
           m  := t.params[2],
           p  := t.firstTag().params[1],
           n  := mn/m,
           p2 := p*p,
           IsPosInt(m/p) and IsPosInt(mn/p2) and IsPosInt(n) and (n>p)
       ),

       # Note: TLs with a ParCellDMP tag should never have either of the Ix set (Should always be I1 x L x I1)
       children := (self, t) >> let(
           mn := t.params[1],
           m  := t.params[2],
           p  := t.firstTag().params[1],
           n  := mn/m,
           p2 := p*p,

           # Eat up the parallel (ParCellDMP) tag.
           tags  := Drop(t.getTags(), 1),

           [ [ TL(mn/p, m/p).withTags(tags), TL(n, p, 1, m/p).withTags(tags), InfoNt(m,n,p) ] ]
       ),

       apply := (self, t, C, Nonterms) >> let(
           m   := Nonterms[3].params[1],
           n   := Nonterms[3].params[2],
           p   := Nonterms[3].params[3],
           p2  := p*p,
           mn  := m*n,

           PTensor(C[1], p) * Comm_Cell(p, mn/p2) * PTensor(C[2], p)
       )
   ) #TL_CellDMP
));

setStandAlone := function(tags)
   local i, retval;
   retval := Copy(tags);

   for i in retval do
      i.isEdge := true;
   od;
   return(retval);
end;

#-----------------------------------------------------------------------------#-----------------------------------------------------------------------------#----------------------------------------------------------------------------- #-----------------------------------------------------------------------------
#-----------------------------------------------------------------------------
# Cell DMP Paradigm - new
#-----------------------------------------------------------------------------
NewRulesFor(GT, rec(
    GT_StickyL := rec(
        requiredFirstTag := StickyL,
        applicable := (self, t) >> true,
        children   := (self, t) >> [[ t.withoutFirstTag() ]],
        apply := (self, t, C, Nonterms) >> C[1]
    ),

    # Converts (In x A) to A DistSum (parallel loop)
    GT_CellDMP_base := rec(
        maxSize       := false,
        minSize       := false,

        requiredFirstTag := ParCellDMP,

        applicable := (self, t) >> let(
                rank   := Length(t.params[4]),
                procs  := t.firstTag().params[1],
                rank = 1
                  and let(its := t.params[4][1],
                    PatternMatch(t, [GT, @(1), @(2,XChain), @(3,XChain), ...], empty_cx())
                    and IsPosInt(its/(procs)))
                  and (t.params[2]=GTPar and t.params[3]=GTPar)

        ),

        children := (self, t) >> let(
            spl    := t.params[1],
            g      := t.params[2],
            s      := t.params[3],
            its    := t.params[4][1],
            procs  := t.firstTag().params[1],
            tags  := Drop(t.getTags(), 1),  # No nested parallelism for now
            [ [ GT(spl, g, s, [its / procs]).withTags(tags), InfoNt(procs) ] ]
        ),

        apply := (self, t, C, Nonterms) >> let(
            spl     := t.params[1],
            g       := t.params[2],
            s       := t.params[3],
            its     := t.params[4][1],
            gg      := When(g.params[1]=[0,1], XChain([0,1,2]), XChain([1,2,0])),
            ss      := When(s.params[1]=[0,1], XChain([0,1,2]), XChain([1,2,0])),
            procs   := t.firstTag().params[1],
            i       := var("spuid", TInt, procs),
            cscatter:= ScatDist(ss.part(1, i, Rows(spl), [procs, its/(procs)]).range(), 1, procs, i),
            cgather := GathDist(gg.part(1, i, Cols(spl), [procs, its/(procs)]).range(), 1, procs, i),

            DistSum(procs, i, procs, cscatter * C[1] * cgather)

        ) #apply
    ), #GT_CellDMP_Base

    # Converts all other TTensor constructs (IxA)L, L(IxA), (AxI) into L's and (IxA)s
    # The L's are factorized to PTensor * Comm_Cell * PTensor
    GT_CellDMP_gen := rec(
        maxSize    := false,
        minSize    := false,
        forTransposition := false,

        requiredFirstTag := ParCellDMP,

        # NOTE: Ensure applicability for vectorized distributed algo
        applicable := (self, t) >> let(
                rank   := Length(t.params[4]),
                procs  := t.firstTag().params[1],
                v      := When(Length(t.firstTag().params)>1, t.firstTag().params[2], 1),
                n      := t.params[4][1],
                rank = 1
                and PatternMatch(t, [GT, @(1), @(2,XChain), @(3,XChain), ...], empty_cx())
                and IsPosInt((n/v) / procs)
                and not(t.params[2]=GTPar and t.params[3]=GTPar)
                #and t.params[2]=GTVec and t.params[3]=GTVec #NOTE: remove this later
             ),

        children := (self, t) >> let(
            spl    := t.params[1],
            g      := t.params[2],
            s      := t.params[3],
            n      := t.params[4][1],
            procs  := t.firstTag().params[1],
            v      := When(Length(t.firstTag().params)>1, t.firstTag().params[2], 1),
            tags   := t.getTags(),

            ltags  := When(tags[1].leftChild(),  setStandAlone(tags), tags),
            rtags  := When(tags[1].rightChild(), setStandAlone(tags), tags),

           # Error("BP"),

            m      := Rows(spl),
            #splnew := When(v=1, spl, GT(spl.withTags([StickyL(v)]), g, s, [v])),
            splnew := When(v=1, spl, GT(spl, g, s, [v])),

            Cond(s = GTPar and g = GTVec, # (IxA)L case
                  [ [
                         GT(splnew, GTPar, GTPar, [n/v]).withTags(tags),
                         TL(m*n/v, n/v, 1, v).withTags(rtags)
                  ] ],

                  # TEMP: commenting this out. We shouldn't see this case
                  s = GTVec and g = GTPar, # L(IxA) case
                  Error("WTF?"),
                  #[ [
                  #       TL(m*n/v, m, 1, v).withTags(ltags),
                  #       GT(splnew, GTPar, GTPar, [n/v]).withTags(tags)
                  #] ],

                  s = GTVec and g = GTVec, # AxI case
                  [ [
                         TL(m*n/v, m, 1, v).withTags(ltags),
                         GT(splnew, GTPar, GTPar, [n/v]).withTags(tags),
                         TL(m*n/v, n/v, 1, v).withTags(rtags)
                  ] ]
            )
        ),

        apply := (self, t, C, Nonterms) >> When(Length(C)=2, C[1]*C[2], C[1]*C[2]*C[3])
    ) #GT_CellDMP_gen

));

NewRulesFor(TL, rec(
   TL_CellDMP := rec(
       maxSize       := false,
       minSize       := false,
       forTransposition := false,

       requiredFirstTag := ParCellDMP,

       applicable := (self, t) >> let(
           mn := t.params[1],
           m  := t.params[2],
           p  := t.firstTag().params[1],
           n  := mn/m,
           p2 := p*p,
           IsPosInt(m/p) and IsPosInt(mn/p2) and IsPosInt(n) and (n>p)
       ),

       children := (self, t) >> let(
           mnByv := t.params[1],
           m     := t.params[2],
           v     := t.params[4],
           p     := t.firstTag().params[1],
           nByv  := mnByv/m,
           p2    := p*p,

           # Eat up the parallel (ParCellDMP) tag.
           tags  := Drop(t.getTags(), 1),
           #When(IsBound(tags[1].isEdgeTL), Error("Edge"), Print("")),

           [ [ TL(mnByv/p, m/p, 1, v).withTags(tags), TL(nByv, p, 1, (m*v)/p).withTags(tags), InfoNt(m,nByv,p,v) ] ]
       ),

       apply := (self, t, C, Nonterms) >> let(
           m      := Nonterms[3].params[1],
           nByv   := Nonterms[3].params[2],
           p      := Nonterms[3].params[3],
           v      := Nonterms[3].params[4],
           p2     := p*p,
           tags   := t.getTags(),

           # Hacking RulesTerm inside here was done by ff.
           #fl     := PTensor(paradigms.vector.rewrite.RulesTerm(C[1]), p),
           #fr     := PTensor(paradigms.vector.rewrite.RulesTerm(C[2]), p),

           #fl     := When(tags[1].leftEdge(),   LeftEdge(PTensor(C[1], p)), PTensor(C[1], p)),
           #fr     := When(tags[1].rightEdge(), RightEdge(PTensor(C[2], p)), PTensor(C[2], p)),

           # Works, but need to terminate blockvperm
           #Cond(tags[1].leftEdge(),  (LeftEdge(fl) * Comm_Cell(p, (m*nByv*v)/p2) * fr),
           #     tags[1].rightEdge(), (    fl  * Comm_Cell(p, (m*nByv*v)/p2) * RightEdge(fr)),
           #                          (    fl  * Comm_Cell(p, (m*nByv*v)/p2) *      fr ))

           frules := MergedRuleSet(RulesSums, RulesFuncSimp, KillSymSPL, RulesVec, RulesTerm, RulesPropagate),

           csums := frules(C[2].sums()),
           sg := ScatGath(fTensor(fId(csums.dimensions[1]/csums.v), fId(csums.v)), fTensor(csums.func, fId(csums.v)) ),
           sgnew := frules(sg.toloop(csums.v)),

           csums1 := frules(C[1].sums()),
           sg1 := ScatGath(fTensor(fId(csums1.dimensions[1]/csums1.v), fId(csums1.v)), fTensor(csums1.func, fId(csums1.v)) ),
           sgnew1 := frules(sg1.toloop(csums1.v)),

           fl     := When(tags[1].leftEdge(),
                     #PTensor(csums1, p),
                     PTensor(sgnew1, p),
                     PTensor(C[1], p)),

           fr     := When(tags[1].rightEdge(),
                     #PTensor(csums, p),
                     PTensor(sgnew, p),
                     PTensor(C[2], p)),

           #Error("BP0"),

           fl * Comm_Cell(p, (m*nByv*v)/p2) * fr

       )
   ) #TL_CellDMP
));


text:
################################################################################
##
#W get_duality.gi              DwG Package
##
#W Paul Bruillard, Cesar Galindo, Siu-Hung Ng, Julia Plavnik, Eric Rowell,
#W Zhenghan Wang
##
## Installation file for get_duality functions of the DwG Package
##
#Y Copyright (C) 2016, Battelle Memorial Institute
##
################################################################################

################################################################################
##
#F get_duality(<list>). . . . a helper function for computing duality
##
InstallGlobalFunction(get_duality, function(N)
  local dual_pairs, unpaired, x, x_dual,r;
  r:=Size(N);
  unpaired:=[1..r];
  while Size(unpaired) > 0 do
    x:=Remove(unpaired,1);
    Display(x);
    if Size(unpaired) > 0 then
      for x_dual in unpaired do
        if N[x][x_dual][1] = 1 then
          break;
        fi;
      od;
    else
      x_dual := x;
    fi;
    if N[x][x_dual][1] <> 1 then
      Error("corrupt fusion matrices.");
    fi;
    Add(dual_pairs,[x,x_dual]);
  od;
  return dual_pairs;
end);

#E get_duality.gi . . . . . . . . . . . . . . . . . . . . . . . ends here

max_stars_repo_name: pnnl/GroupTheoretical

max_issues_repo_name: pnnl/GroupTheoretical

max_forks_repo_name: pnnl/GroupTheoretical


text:
#
# JupyterMsg: Jupyter kernel using ZeroMQ
#
# Implementations
#
# TODO: Check signature
InstallGlobalFunction( JupyterMsgDecode,
function(kernel, raw)
    local result, bindIfBound, sl, ids, tmp;

    bindIfBound := function(name, pos)
        if IsBound(raw[pos]) then
            result.(name) := JsonStringToGap(raw[pos]);
        fi;
    end;

    result := rec();

    ids := [];
    sl := 1;
    while raw[sl] <> "<IDS|MSG>" do
        Add(ids, raw[sl]);
        sl := sl + 1;
    od;
    result.hmac := raw[sl + 1];
    result.remainder := raw;

    bindIfBound("header", sl + 2);
    bindIfBound("parent_header", sl + 3);
    bindIfBound("metadata", sl + 4);
    bindIfBound("content", sl + 5);

    tmp := CRYPTING_SHA256_HMAC( kernel!.SessionKey
                               , Concatenation( raw[sl + 2]
                                              , raw[sl + 3]
                                              , raw[sl + 4]
                                              , raw[sl + 5] ) );
    tmp := List(tmp, CRYPTING_HexStringIntPad8);
    tmp := LowercaseString(Concatenation(tmp));
    result.hmac_verify := tmp;

    if result.hmac <> result.hmac_verify then
        PrintTo( "*errout*", "HMAC verification for message failed: "
                 , result.hmac, " <> ", result.hmac_verify, "\n" );
    fi;

    return result;
end);

InstallGlobalFunction( JupyterMsgEncode,
function(kernel, msg)
    local raw, k, bindIfBound, tmp;

    bindIfBound := function(pos, name)
        if IsBound(msg.(name)) then
            raw[pos] := GapToJsonString(msg.(name));
        fi;
    end;

    raw := [];
    # TODO: What is the correct behaviour here?
    if IsBound(msg.uuid) then
        raw[1] := msg.uuid;
    else
        raw[1] := "";
    fi;
    raw[2] := "<IDS|MSG>";
    raw[3] := msg.hmac;
    bindIfBound(4, "header");
    bindIfBound(5, "parent_header");
    bindIfBound(6, "metadata");
    bindIfBound(7, "content");

    # TODO: Ugly
    if Length(raw) > 3 then
        tmp := CRYPTING_SHA256_HMAC(msg.key,
                                    Concatenation( raw[4]
                                                 , raw[5]
                                                 , raw[6]
                                                 , raw[7]));
        tmp := List(tmp, CRYPTING_HexStringIntPad8);
        tmp := LowercaseString(Concatenation(tmp));
        raw[3] := LowercaseString(tmp);
    fi;
    return raw;
end);

InstallGlobalFunction(JupyterMsgRecv,
function(kernel, sock)
    local raw;
    raw := ZmqReceiveList(sock);
    if IsBound(kernel!.ProtocolLog) then
        AppendTo(kernel!.ProtocolLog, raw);
        AppendTo(kernel!.ProtocolLog, "\n");
    fi;
    return JupyterMsgDecode(kernel, raw);
end);

InstallGlobalFunction(JupyterMsgSend,
function(kernel, sock, msg)
    local msg2;
    if IsBound(kernel!.ProtocolLog) then
        msg2 := JupyterMsgEncode(kernel, msg);
        AppendTo(kernel!.ProtocolLog, msg2);
        AppendTo(kernel!.ProtocolLog, "\n");
    fi;
    ZmqSend(sock, JupyterMsgEncode(kernel, msg));
end);

# Create a message template with the necessasry fields filled
InstallGlobalFunction(JupyterMsg,
function(kernel, msg_type, parent_header, content, metadata)
    return rec( uuid := kernel!.ZmqIdentity
              , sep := "<IDS|MSG>"                 # This could be in JupyterEncode
              , hmac := ""
              , header := rec( username := kernel!.Username
                             , session := kernel!.SessionID
                             , msg_type := msg_type
                             , version := kernel!.ProtocolVersion
                             , date := ISO8601Stamp()
                             , msg_id := HexStringUUID(RandomUUID())
                             )
              , parent_header := parent_header
              , metadata := metadata
              , content := content
              , key := kernel!.SessionKey
                             # This shouldn't be here as all
                             # messaging functions
                             # should just be running in kernel context
              );
end);

max_stars_repo_name: ZachNewbery/JupyterKernel

max_issues_repo_name: ZachNewbery/JupyterKernel

max_forks_repo_name: ZachNewbery/JupyterKernel


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

same_params := (pat1, targ, pat2) -> pat1.target(targ).cond(e->e.params = pat2.val.params);
one_by_one := x -> x.range() = 1 and x.domain() = 1;

# does not work with .cond (!!!) NOTE
Enum(100, x->@(x), @z, @n, @m, @k, @l, @j, @g, @phi, @r, @s, @ss, @b, @bb);
Enum(200, x->@(x), @M, @N, @R, @S, @W, @F, @G, @D);

IsIndexMapping := x -> IsBound(x._perm) and x._perm;

# L^mn_n
ltensor_flip := function(match, mn, n)
    local m, i, pos, dom, ldom, lran, ll, rr, res, ldom;
    dom := Product(match, x->x.domain());
    ldom := 1; lran := 1; pos := 1;
    while lran<>n do
        lran:=lran*range(match[pos]);
        ldom:=ldom*domain(match[pos]);
	pos:=pos+1;
    od;
    ll := match{[1..pos-1]};
    rr := match{[pos..Length(match)]};
    res := fTensor(Concatenation(rr, ll));

    if ldom=1 or ldom=dom then
    return res;
    else
    return fCompose(res, L(dom, ldom));
    fi;
end;

# for given diagDirsum looking for child which relative domain start <= offset and relative domain end >= (offset + domain)
# returns [child index, offset relative to child]
_findDirsumChild := function( dirsum, offset, domain)
    local i, ch, d;
    ch := dirsum.children();
    d := 0;
    for i in [1..Length(ch)] do
        if offset >= d then
            d := d + ch[i].domain();
            if offset + domain <= d then
                return [i, offset - d + ch[i].domain()];
            fi;
        else
            return [0, 0];
        fi;
    od;
    return [0, 0];
end;

# *******************************************************************
Class(RulesFuncSimp, RuleSet);
RewriteRules(RulesFuncSimp, rec(
 # ===================================================================
 # Operator flattening and identity
 # ===================================================================
 FlattenTensor  := ARule(fTensor,  [@(1, fTensor) ], e -> @(1).val.children()),
 FlattenCompose := ARule(fCompose, [@(1, fCompose)], e -> @(1).val.children()),

 ComposeId1 := ARule(fCompose, [@(1), fId ], e -> [@(1).val]),
 ComposeId2 := ARule(fCompose, [fId,  @(1)], e -> [@(1).val]),

 fBase_1x1 := Rule(@(1, fBase, one_by_one), e->fId(1)),

 fId_drop_Value := Rule([@(1, fId), @(2, Value)], e->fId(@(2).val.v)),

 Const_fbase := ARule(fCompose, [@(1), @(2, fBase, e -> IsInt(e.params[2]))],
     e -> [ When(IsInt(@(1).val.range()),
                 fBase(@(1).val.range(), @(1).val.at(@(2).val.params[2])),
                 fConst(@(1).val.range(), 1, @(1).val.at(@(2).val.params[2]))) ]),

 Compose1x1_1 := ARule(fCompose, [@(1).cond(one_by_one), @(2)],  e -> [@(2).val]),
 Compose1x1_2 := ARule(fCompose, [@(2), @(1).cond(one_by_one) ], e -> [@(2).val]),

 TensorId1 := ARule(fTensor, [@(1), [@(2,[fId,J]), 1]], e -> [@(1).val]),
 TensorId2 := ARule(fTensor, [[@(2,[fId,J]), 1], @(1)], e -> [@(1).val]),

 TensorId1V := ARule(fTensor, [@(1), [@(2,[fId,J]), @(3, Value, e->e.v=1)]], e -> [@(1).val]),
 TensorId2V := ARule(fTensor, [[@(2,[fId,J]), @(3, Value, e->e.v=1)], @(1)], e -> [@(1).val]),

 diagTensorId1  := ARule(diagTensor, [@(1), [fConst, @(2), 1, 1]], e -> [@(1).val]),
 diagTensorIdV1 := ARule(diagTensor, [@(1), [fConst, @(2), 1, _1]], e -> [@(1).val]),
 diagTensorId2  := ARule(diagTensor, [[fConst, @(2), 1, 1], @(1)], e -> [@(1).val]),
 diagTensorIdV2 := ARule(diagTensor, [[fConst, @(2), 1, _1], @(1)], e -> [@(1).val]),

 diagTensorZero1  := ARule(diagTensor, [@(1), [fConst, @(2), 1, 0]], e -> [fConst(@(1).val.domain(), 0)]),
 diagTensorZeroV1 := ARule(diagTensor, [@(1), [fConst, @(2), 1, _0]], e -> [fConst(@(1).val.domain(), 0)]),
 diagTensorZero2  := ARule(diagTensor, [[fConst, @(2), 1, 0], @(1)], e -> [fConst(@(1).val.domain(), 0)]),
 diagTensorZeroV2 := ARule(diagTensor, [[fConst, @(2), 1, _0], @(1)], e -> [fConst(@(1).val.domain(), 0)]),

 DropL := Rule(@(1, [L,OS], e -> e.params[2] = 1 or e.params[2] = e.params[1]), e -> fId(@1.val.params[1])),
 DropTr := Rule(@(1, Tr, e -> e.params[1] = 1 or e.params[2] = 1), e -> fId(@1.val.domain())),

 IP_toI  := Rule([IP, @(2), fId], e -> fId(@(2).val)),
 J_1   := Rule([J, 1], e->fId(1)),
 OS_2  := Rule([OS, 2, @], e->fId(2)),

 OS_fAdd := ARule(fCompose, [ [OS, @(1), @(2).cond(e -> e in [@(1).val-1, -1])],
                              [fAdd, @(3), @(4).cond(e -> e = @(3).val - 1), 1] ],
     e -> [ fAdd(@(3).val, @(4).val, 1), J(@(4).val) ]),

 Lfold := ARule(fCompose, [@(1,L), @(2,L),e->let(
             prod := @(1).val.params[2] * @(2).val.params[2],
             a:=@(1).val.params[1],
             When(prod < a, prod in FactorsInt(a), (prod / a) in FactorsInt(a)))],
     e -> let(
         prod := @(1).val.params[2] * @(2).val.params[2], a:=@(1).val.params[1],
         [ L(a, When(prod < a, prod, prod / a)) ])),

 Hfold :=  ARule(fCompose, [@(1,H), @(2,H)],
     e-> [H(@(1).val.params[1],
             @(2).val.params[2],
             @(1).val.params[3]+@(2).val.params[3]*@(1).val.params[4],
             @(1).val.params[4]*@(2).val.params[4])]),

 # ===================================================================
 # Diagonal Functions
 # ===================================================================
 ComposePrecompute := Rule([fCompose, ..., [fPrecompute, @(1)], ...],
     e -> fPrecompute(ApplyFunc(fCompose, List(e.children(),
         x->When(ObjId(x)=fPrecompute, x.child(1), x))))),

 PrecomputePrecompute := Rule([fPrecompute, [fPrecompute, @(1)]], e->fPrecompute(@(1).val)),

 # drop diagDirsum when H refers to only one its child
 drop_diagDirsum_H := ARule( fCompose, [ @(1, diagDirsum), @(2, H, x -> x.params[4] = 1 and _findDirsumChild(@(1).val, x.params[3], x.params[2])[1]>0)],
     e -> let( f := _findDirsumChild(@(1).val, @(2).val.params[3], @(2).val.params[2]),
               obj := @(1).val.children()[f[1]],
               [ obj, H(obj.domain(), @(2).val.params[2], f[2], 1) ] )),

 RCData_fCompose_TReal := Rule([RCData, [@(1,fCompose), @(2).cond(e->e.range = TReal), ...]],
     e -> fCompose(RCData(@(2).val), diagTensor(fCompose(Drop(@(1).val.children(),1)), fId(2)))),

 RCData_fCompose_TInt := Rule([RCData, [@(1,fCompose), @(2).cond(e->e.range = TInt), ...]],
     e -> fCompose(RCData(@(2).val), fTensor(fCompose(Drop(@(1).val.children(),1)), fId(2)))),

# RCData_fCompose := Rule([RCData, [@(1,fCompose), @(2), ...]],
#     e -> fCompose(RCData(@(2).val), fTensor(fCompose(Drop(@(1).val.children(),1)), fId(2)))),

 RCData_Precompute := Rule([RCData, [fPrecompute, @(1)]], e->fPrecompute(RCData(@(1).val))),
 RCData_CRData     := Rule([RCData, [CRData, @(1)]], e -> @(1).val),
 CRData_RCData     := Rule([CRData, [RCData, @(1)]], e -> @(1).val),

 FConj_Precompute := Rule([@(0,[FConj,FRConj]), [fPrecompute, @(1)]], e->fPrecompute(ObjId(@(0).val)(@(1).val))),

 ComposeConstX := ARule(fCompose, [@(1,fConst), @(2)],
     e -> [fConst(@(1).val.params[1], @(2).val.domain(), @(1).val.params[3])]),

 fAdd0 := Rule([fAdd, @(1), @(2).cond(e->e=@(1).val), 0], e->fId(@(1).val)),

 fAddPair := ARule(fCompose, [@(1, fAdd), @(2, fAdd)],
     e -> [ fAdd(@(1).val.range(), @(2).val.domain(),
             @(1).val.params[3] + @(2).val.params[3]) ]),

 fTensor_fAdd_H := ARule(fTensor, [ [@(1, fAdd),@, _1, @], @(2, H)],
     e -> [ H(@(1).val.params[1] * @(2).val.params[1],
              @(2).val.params[2],
              @(2).val.params[1] * @(1).val.params[3] + @(2).val.params[3],
              @(2).val.params[4]) ]),

 # ===================================================================
 # fTensor
 # ===================================================================
 # fId(m) (X) fId(n) -> fId(m*n)
  TensorIdId := ARule(fTensor, [@(1,fId), @(2,fId)],
    e -> [ fId(@(1).val.params[1] * @(2).val.params[1]) ]),

  TensorId1 := ARule(fTensor, [@(1, fId, x->x.params[1]=1)],
    e -> []),

 # L(mn,m) o (fBase(m,j) (X) f) -> f (X) fBase(j,m)
 LTensorFlip := ARule(fCompose,
      [ @(1,L), [ @(3,fTensor), @(2).cond(e->range(e) = @(1).val.params[2] and domain(e)=1), ...] ],
 	e -> [ fTensor(Copy(Drop(@(3).val.children(), 1)), Copy(@(2).val)) ] ),

 # L(mn,n) o f (X) fBase(j,m) -> (fBase(m,j) (X) f)
 LTensorFlip1 := ARule(fCompose,
      [ @(1,L), [ @(3,fTensor), ...,
             @(2, fBase, e->range(e) = @(1).val.params[1]/@(1).val.params[2]) ]],
    e -> [fTensor(Copy(Last(@(3).val.children())), Copy(DropLast(@(3).val.children(), 1)))] ),

 L_H := ARule(fCompose, [ @(1, L), [ @(2, fTensor), @(3, fBase), @(4, fId, e->IsInt(@(1).val.params[1]/(@(1).val.params[2] * e.params[1])))] ],
    e-> [ H(@(1).val.params[1], @(4).val.params[1],
            @(4).val.params[1] * @(1).val.params[2] * imod(@(3).val.params[2], @(1).val.params[1] / (@(4).val.params[1] * @(1).val.params[2])) +
            idiv(@(3).val.params[2], @(1).val.params[1] / (@(4).val.params[1] * @(1).val.params[2])),
            @(1).val.params[2]) ]),

  Refl0_u_H0 := ARule(fCompose, [@(1, Refl0_u), @(2, H, e -> e.params[3]=0 and e.params[4]=1)], e -> [ let(
     k := @(1).val.params[1],
     H(@(1).val.range(), @(2).val.domain(), 0, k)) ]),

  Refl0_u_Hrest := ARule(fCompose, [
          @(1, Refl0_u),
          @(2, H, e -> e.params[3]=@(1).val.params[2] and e.params[4]=1),
          [fTensor, @(3,fBase), @(4, fId, e->e.params[1]=@(1).val.params[2])]],
      e -> [ let(k := @(1).val.params[1],
              BH(@(1).val.range(), 2*@(1).val.range(), @(4).val.domain(), 1+@(3).val.params[2], 2*k)) ]),

  Refl0_u_Hrest_1it := ARule(fCompose, [
          @(1, Refl0_u),
          @(2, H, e -> e.params[3]=@(1).val.params[2] and e.params[4]=1 and e.params[3]=e.params[2]) ],
      e -> [ let(k := @(1).val.params[1],
              BH(@(1).val.range(), 2*@(1).val.range(), @(2).val.domain(), 1, 2*k)) ]),

  Refl1_H := ARule(fCompose, [@(1, Refl1), [fTensor, @(2,fBase), @(3, fId, e->e.params[1]=@(1).val.params[2])]],
      e -> let(
          k := @(1).val.params[1],
          [ BH(@(1).val.range(), 2*@(1).val.range() - 1, @(3).val.domain(), @(2).val.params[2], 2*k) ])),

  fTensor_X_fId_H := ARule(fCompose, [[@(1, fTensor), ..., [fId, @(2)]], # NOTE: assumes d|base
        [@(3, H), @(4).cond(e->_divides(@(2).val,e)), @(5).cond(e->_divides(@(2).val,e)),
                  @(6).cond(e->_divides(@(2).val,e)), _1]], e ->
       let(p := @(3).val.params, d := @(2).val, base := p[3],
           [ @(1).val, fTensor(H(p[1]/d, p[2]/d, base/d, 1), fId(d)) ])),

 TrTensor := ARule(fCompose, [@(1, Tr), [ fTensor, @(2), @(3).cond(e->e.range()=@(1).val.params[1]) ]],
     e -> [ fCompose(fTensor(@(3).val, @(2).val), Tr(@(3).val.domain(), @(2).val.domain())) ]
 ),

 LTensorGeneral := ARule(
fCompose,
[@(1,[L,Tr]),
 @(2,fTensor,
 e -> let(
    p := @(1).val.params,
    a := When(ObjId(@(1).val)=Tr, p[1], p[1]/p[2]),
    b := p[2],
    ForAll(e.children(), x->AnySyms(x.domain(), x.range()) or x.domain()>1 and x.range()>1) and # workaround old fId(1) bug
    full_merge_tensor_chains(@(3), [b, a], e.children(), (x,y)->y,
 Product, fTensor,
                             fId, x->x, x->x, range) <> false))],
     e -> [ ltensor_flip(@(3).val, @(1).val.params[1], @(1).val.params[2]) ]
 ),

 # fTensor o fTensor -> fTensor( ..o.., ..o.. , ...)
 TensorMerge := ARule(fCompose,
      [ @(1,fTensor), @(2,fTensor,e->full_merge_tensor_chains(
          @(3), @(1).val.children(), e.children(),
                fCompose, fTensor, fTensor, x->x, x->x, domain, range) <> false) ],
 e -> [ fTensor(@(3).val) ]
 ),

 TensorMerge_fIdInt := ARule(fCompose,
      [ @(1,fTensor,e->let(l1:=Last(e.children()), ObjId(l1)=fId and IsInt(l1.domain()))),
        @(2,fTensor, e->let(l1:=Last(@(1).val.children()), l2 := Last(e.children()), IsInt(l2.domain()) and ObjId(l2)=fId and Gcd(l1.domain(), l2.domain())>1))],
        e -> let(l1 := Last(@(1).val.children()).domain(), l2 := Last(@(2).val.children()).domain(), gcd := Gcd(l1, l2),
            [ fTensor(fCompose(fTensor(DropLast(@(1).val.children(), 1), fId(l1/gcd)), fTensor(DropLast(@(2).val.children(), 1), fId(l2/gcd))), fId(gcd)) ])
    ),

 # diagTensor o fTensor -> diagTensor( ..o.., ..o.. , ...)
 diagTensorMerge := ARule(fCompose,
      [ @(1,diagTensor), @(2,fTensor,e->full_merge_tensor_chains(
          @(3), @(1).val.children(), e.children(),
                fCompose, diagTensor, fTensor, x->x, x->fConst(range(x), 1), domain, range) <> false) ],
 e -> [ diagTensor(@(3).val) ]
 ),

 # fTensor(..., Y, ...) o X -> fTensor(..., Y o X, ...)
 TensorComposeMergeRight := ARule(fCompose,
      [ @(1,fTensor), @(2).cond(e -> compat_tensor_chains(
                                @(1).val.children(), [e], domain, range)) ],
 e -> [ fTensor(merge_tensor_chains(
         @(1).val.children(), [@(2).val], fCompose, x->x, x->x, domain, range)) ]
 ),

 # X o fTensor(..., Y, ...) -> fTensor(..., X o Y, ...)
 TensorComposeMergeLeft := ARule(fCompose,
      [ @(1), @(2,fTensor,e-> IsInt(@(1).val.range()) and compat_tensor_chains(  # NOTE: BAD HACK!!! Cannot do that with diag functions...
                             [@(1).val], e.children(), domain, range)) ],
 e -> [ fTensor(merge_tensor_chains(
         [@(1).val], @(2).val.children(), fCompose, x->x, x->x, domain, range)) ]
 ),

 # [fBase,fId] X (fB o fAny) -> fB2 o (fI X fAny)
 TensorB := ARule(fTensor, [
        @(1, [fBase, fId], e -> Is2Power(e.range())),
        @(2, fCompose, e -> ObjId(e.child(1)) = fB)
    ],
    e -> [fCompose(
        let(i := @(1).val,
            b := @(2).val.child(1),
            fB(i.range() * b.range(), List(b.params[3], e -> Log2Int(i.range()) + e))
        ),
        fTensor(
            @(1).val,
            @(2).val.child(2)
        )
    )]
 ),

 # ===================================================================
 # Other
 # ===================================================================

 # primitive constant folding for functions
 Compose_FList_const := ARule(fCompose, [@(1, FList), @(2).cond(x->x.free()=[])],
     e -> let(
	 lst := @(1).val,  idx := List(@(2).val.tolist(), EvalScalar),
	 [ FList(lst.t, lst.list{1+idx}) ])),

 # fCond o f  (note: fCond o fCond should only be handled by Compose_fCond_lft)
 Compose_fCond_rt := ARule(fCompose, [[fCond, @(1), @(2), @(3)], @(4).cond(e->ObjId(e)<>fCond)],
     e -> [ fCond(fCompose(@(1).val, @(4).val), fCompose(@(2).val, @(4).val), fCompose(@(3).val, @(4).val)) ]),

 Compose_fCond_lft := ARule(fCompose, [@(4), [fCond, @(1), @(2), @(3)]],
     e -> [ fCond(@(1).val, fCompose(@(4).val, @(2).val), fCompose(@(4).val, @(3).val)) ]),

 fCond_fold0  := Rule([fCond, [fConst, @, @,  0], @(2), @(3)], e -> @(3).val),
 fCond_foldV0 := Rule([fCond, [fConst, @, @, _0], @(2), @(3)], e -> @(3).val),

 fCond_fold1  := Rule([fCond, [fConst, @, @,  1], @(2), @(3)], e -> @(2).val),
 fCond_foldV1 := Rule([fCond, [fConst, @, @, _1], @(2), @(3)], e -> @(2).val),

 # B2 o fTensor

    fCompose_B2_fTensor := ARule(fCompose, [@(1,B2), @(2,fTensor)],
        e -> [fB2(@(1).val.getES(), @(2).val.children())]
    ),

 # L(abc,b) o L(abc,c) -> L(abc,bc)
 fCompose_Tr := ARule( fCompose, [ [@(1, [Tr, L]), @, @(2)], [@(3, [Tr, L]), @, @(4).cond( x -> let( abc := @(1).val.range(), @(3).val.domain() = abc and (abc mod (@(2).val*x) = 0)))]],
     e -> let( abc := @(1).val.range(), bc := @(2).val*@(4).val,
         When( ObjId(@(1).val) = L, [L(abc, bc)], [Tr(abc/bc, bc)] ))),

 fDirsum_fIds := ARule( fDirsum, [@(1, fId), @(2, fId)], e -> [fId(@(1).val.range()+@(2).val.range())]),

));


text:
# run this file with
#   gap -b -q < /dev/null group.gap | perl -pe 's/\\\n//s' | indent -kr

Print("/* ----- data generated by group.gap begins ----- */\n\n");
Print("struct group {\n    unsigned long autosize;\n");
Print("    int order, ngens;\n    const char *gens;\n};\n");
Print("struct groups {\n    int ngroups;\n");
Print("    const struct group *groups;\n};\n\n");
Print("static const struct group groupdata[] = {\n");
offsets := [0];
offset := 0;
for n in [2..26] do
  Print("    /* order ", n, " */\n");
  for G in AllSmallGroups(n) do

    # Construct a representation of the group G as a subgroup
    # of a permutation group, and find its generators in that
    # group.

    # GAP has the 'IsomorphismPermGroup' function, but I don't want
    # to use it because it doesn't guarantee that the permutation
    # representation of the group forms a Cayley table. For example,
    # C_4 could be represented as a subgroup of S_4 in many ways,
    # and not all of them work: the group generated by (12) and (34)
    # is clearly isomorphic to C_4 but its four elements do not form
    # a Cayley table. The group generated by (12)(34) and (13)(24)
    # is OK, though.
    #
    # Hence I construct the permutation representation _as_ the
    # Cayley table, and then pick generators of that. This
    # guarantees that when we rebuild the full group by BFS in
    # group.c, we will end up with the right thing.

    ge := Elements(G);
    gi := [];
    for g in ge do
      gr := [];
      for h in ge do
        k := g*h;
	for i in [1..n] do
	  if k = ge[i] then
	    Add(gr, i);
	  fi;
	od;
      od;
      Add(gi, PermList(gr));
    od;

    # GAP has the 'GeneratorsOfGroup' function, but we don't want to
    # use it because it's bad at picking generators - it thinks the
    # generators of C_4 are [ (1,2)(3,4), (1,3,2,4) ] and that those
    # of C_6 are [ (1,2,3)(4,5,6), (1,4)(2,5)(3,6) ] !

    gl := ShallowCopy(Elements(gi));
    Sort(gl, function(v,w) return Order(v) > Order(w); end);

    gens := [];
    for x in gl do
      if gens = [] or not (x in gp) then
        Add(gens, x);
	gp := GroupWithGenerators(gens);
      fi;
    od;

    # Construct the C representation of the group generators.
    s := [];
    for x in gens do
      if Size(s) > 0 then
        Add(s, '"');
        Add(s, ' ');
        Add(s, '"');
      fi;
      sep := "\\0";
      for i in ListPerm(x) do
        chars := "ABCDEFGHIJKLMNOPQRSTUVWXYZ";
        Add(s, chars[i]);
      od;
    od;
    s := JoinStringsWithSeparator(["    {", String(Size(AutomorphismGroup(G))),
                                   "L, ", String(Size(G)),
                                   ", ", String(Size(gens)),
                                   ", \"", s, "\"},\n"],"");
    Print(s);
    offset := offset + 1;
  od;
  Add(offsets, offset);
od;
Print("};\n\nstatic const struct groups groups[] = {\n");
Print("    {0, NULL}, /* trivial case: 0 */\n");
Print("    {0, NULL}, /* trivial case: 1 */\n");
n := 2;
for i in [1..Size(offsets)-1] do
  Print("    {", offsets[i+1] - offsets[i], ", groupdata+",
        offsets[i], "}, /* ", i+1, " */\n");
od;
Print("};\n\n/* ----- data generated by group.gap ends ----- */\n");
quit;

ext: gap


text:
latexPoly:=function(poly,APR,var_name)
    local result,ext_rep,i,j,monom,coeff,
    avars,ext_avars,pos;

    if not IsPolynomialRing(APR) then return String(Int(poly)); fi;

    ext_rep:=ExtRepPolynomialRatFun(poly);
    avars:=IndeterminatesOfPolynomialRing(APR);
    ext_avars:=List(avars,i->ExtRepPolynomialRatFun(i));
    ext_avars:=List(ext_avars,i->i[1][1]);
    result:="";
    for i in 2*[1..Length(ext_rep)/2]-1 do
        monom:="";
        for j in 2*[1..Length(ext_rep[i])/2]-1 do
            pos:=Position(ext_avars,ext_rep[i][j]);
            monom:=Concatenation(monom,String(avars[pos]));
            if ext_rep[i][j+1]<>1 then
                monom:=Concatenation(monom,"^{",String(ext_rep[i][j+1]),"}");
            fi;
        od;
        coeff:=latexPoly(ext_rep[i+1],CoefficientsRing(APR),"a");#String(Int(ext_rep[i+1]));
        if i > 1 then result:=Concatenation(result,"+"); fi;
        if coeff="1" and monom<>"" then
            result:=Concatenation(result,monom);
        else
            result:=Concatenation(result,coeff,monom);
        fi;
    od;
    result:=ReplacedString(result,String(avars[1]){[1]},var_name);
    return result;
end;

latexUnipotent:=function(u)
    local result,
    coeffs,i,pr,root,poly,APR;

    coeffs:=coefficients(u);
    APR:=ring(chevalleyAdj(u));
    pr:=positiveRoots(chevalleyAdj(u));
    result:="";
    for i in [1..Length(coeffs)] do
        #root:=Concatenation(List(pr[coeffs[i][1]],j->String(j)));
        #poly:=ReplacedString(String(coeffs[i][2]),"*","");
        result:=Concatenation(result,"u_{",String(coeffs[i][1]),"}(",latexPoly(coeffs[i][2],APR,"t"),")");
    od;

    #result:=ReplacedString(result,"a","t");
    return result;
end;

latexUnipotentNice:=function(u)
    local result,
    coeffs,i,pr,root,poly,APR;

    coeffs:=coefficients(u);
    APR:=ring(chevalleyAdj(u));
    pr:=positiveRoots(chevalleyAdj(u));
    result:="";
    for i in [1..Length(coeffs)] do
        if type(chevalleyAdj(u))="F" then root:=Concatenation(List(pr[coeffs[i][1]],j->String(j))){[2,4,3,1]};
        else root:=Concatenation(List(pr[coeffs[i][1]],j->String(j))); fi;
        #poly:=ReplacedString(String(coeffs[i][2]),"*","");
        result:=Concatenation(result,"u_{",root,"}(",latexPoly(coeffs[i][2],APR,"t"),")");
    od;

    #result:=ReplacedString(result,"a","t");
    return result;
end;

latexUnipotentTrimmed:=function(u,step)
    local result,
    coeffs,i,pr,root,poly,APR,index;

    coeffs:=coefficients(u);
    APR:=ring(chevalleyAdj(u));
    pr:=positiveRoots(chevalleyAdj(u));
    result:="";
    index:=1;
    for i in [1..Length(coeffs)] do
        if index > step then
            index:=0;
            result:=Concatenation(result," \\\\ && ");
        fi;
        result:=Concatenation(result,"u_{",String(coeffs[i][1]),"}(",latexPoly(coeffs[i][2],APR,"t"),")\n");
        index:=index+Length(ExtRepPolynomialRatFun(coeffs[i][2]))/2;
    od;

    #result:=ReplacedString(result,"a","t");
    return result;
end;

trimmedString:=function(str,chars,step)
    local result,i;
    result:=[""];
    while Length(str)>0 do
        i:=step;
        while i<Length(str) and not str[i] in chars do
            i:=i+1;
        od;
        if i<Length(str) then
            Add(result,Concatenation(str{[1..i-1]},"\\\\ \n & & "));
            str:=str{[i..Length(str)]};
        else
            Add(result,str);
            str:="";
        fi;
    od;
    return Concatenation(result);
end;

trimmedStringNice:=function(str,chars,step)
    local result,i;
    result:=[""];
    while Length(str)>0 do
        i:=step;
        while i<Length(str) and not str[i] in chars do
            i:=i+1;
        od;
        if i<Length(str) then
            Add(result,Concatenation(str{[1..i-1]},"$\\\\\n&\n$"));
            str:=str{[i..Length(str)]};
        else
            Add(result,str);
            str:="";
        fi;
    od;
    return Concatenation(result);
end;

trimmedString2:=function(str,chars,start,step)
    local result,i,flag,nr_plusuri;
    result:=[""];
    flag:=true;
    nr_plusuri:=0;
    while Length(str)>0 do
        if flag = true then
            i:=step-start;
            flag:=false;
        else i:=step; fi;

        #i:=i+2*Length(Positions(str{[1..Minimum(i,Length(str))]},'+'));

        while i<Length(str) and not str[i] in chars do
            i:=i+1;
            #Print(str[i]);
            #if str[i] in ['+'] then nr_plusuri:=nr_plusuri+1; fi;
        od;
        if i<Length(str) then
            Add(result,Concatenation(str{[1..i-1]},"\\\\ \n & & "));
            str:=str{[i..Length(str)]};
        else
            Add(result,str);
            #Print("\n",nr_plusuri,"\n");
            str:="";
            nr_plusuri:=0;
        fi;
    od;
    return Concatenation(result);
end;

trimmedString2Nice:=function(str,chars,start,step)
    local result,i,flag,nr_plusuri;
    result:=[""];
    flag:=true;
    nr_plusuri:=0;
    while Length(str)>0 do
        if flag = true then
            i:=step-start;
            flag:=false;
        else i:=step; fi;

        #i:=i+2*Length(Positions(str{[1..Minimum(i,Length(str))]},'+'));

        while i<Length(str) and not str[i] in chars do
            i:=i+1;
            #Print(str[i]);
            #if str[i] in ['+'] then nr_plusuri:=nr_plusuri+1; fi;
        od;
        if i<Length(str) then
            Add(result,Concatenation(str{[1..i-1]},"$\\\\\n&\n$"));
            str:=str{[i..Length(str)]};
        else
            Add(result,str);
            #Print("\n",nr_plusuri,"\n");
            str:="";
            nr_plusuri:=0;
        fi;
    od;
    return Concatenation(result);
end;

latexClass:=function(orb,info)
    local fisier,
    u,sys,str,#info,
    step,i;

    #u:=Representative(orb);
    u:=BorelRep(orb);
    sys:=chevalleyAdj(u);
    #info:=handleClass(orb);

    fisier:=Concatenation("tex/data/",type(sys),String(rank(sys)),
                          "/char",String(Characteristic(sys)),
                          "/",Label(orb),".tex");
    fisier:=ReplacedString(fisier,"\\tilde ","s");
    fisier:=Filename(home_dir,fisier);

    step:=150;

    SizeScreen([2000,2000]);
    PrintTo(fisier,"\\textbf{Class $",Label(orb),"$:}");
    AppendTo(fisier,"\\begin{eqnarray*}\n");
    AppendTo(fisier,"&&\\mathbf{u}=");
    AppendTo(fisier,trimmedString(latexUnipotent(u),"u+",step));
#    AppendTo(fisier,"\\end{eqnarray*}\n");
#    AppendTo(fisier,"\\begin{eqnarray*}\n");
    AppendTo(fisier,"\\\\ && \\mathbf{C_U(u)^{\\circ}} =");
    AppendTo(fisier,trimmedString2(latexUnipotent(info[1][1]),"u+",30,step));
#    AppendTo(fisier,"\\end{eqnarray*}\n");
#    AppendTo(fisier,"\\begin{eqnarray*}\n");
    AppendTo(fisier,"\\\\ && \\mathbf{Z(C_U(u)^{\\circ})^{\\circ}} =");
    AppendTo(fisier,trimmedString2(latexUnipotent(info[2]),"u+",40,step));
    AppendTo(fisier,"\\\\ &&=");

    #
    #  short form of decomposition
    #
    str:="";
    for i in info[4][1] do
        if i <> info[4][1][1] then  str:=Concatenation(str," \\mid "); fi;
        str:=Concatenation(str,"W_{",String(Dimension(info[3][i])),"}");
    od;
    AppendTo(fisier,trimmedString(str,"u+",step));
    AppendTo(fisier," \\\|");
    str:="";
    for i in info[4][2] do
        if i <> info[4][2][1] then str:=Concatenation(str," \\mid "); fi;
        str:=Concatenation(str,"W_{",String(Dimension(info[3][i])),"}");
    od;
    AppendTo(fisier,trimmedString(str,"u+",step));
    #
    #  explicite form of decomposition
    #
    AppendTo(fisier,"\\\\ && =");
    str:="";
    for i in info[4][1] do
        if i <> info[4][1][1] then str:=Concatenation(str," \\mid "); fi;
        str:=Concatenation(str,latexUnipotent(info[3][i]));
    od;
    AppendTo(fisier,trimmedString(str,"u+",step));
    AppendTo(fisier," \\\|");
    str:="";
    for i in info[4][2] do
        if i <> info[4][2][1] then str:=Concatenation(str," \\mid "); fi;
        str:=Concatenation(str,latexUnipotent(info[3][i]));
    od;
    AppendTo(fisier,trimmedString(str,"u+",step));

    AppendTo(fisier,"\\end{eqnarray*}\n");
    SizeScreen([156,39]);
end;

latexClasses:=function(orbs,infos)
    local i;
    for i in [1..Length(Classes(orbs))] do
        latexClass(Classes(orbs)[i],infos[i]);
    od;
end;

latexClassShort:=function(orb,info)
    local fisier,
    uu,u,sys,str,#info,
    step,i;

    uu:=Representative(orb);
    u:=BorelRep(orb);
    sys:=chevalleyAdj(u);
    #info:=handleClass(orb);

    fisier:=Concatenation("tex/data/",type(sys),String(rank(sys)),
                          "/char",String(Characteristic(sys)),
                          "/",Label(orb),".tex");
    fisier:=ReplacedString(fisier,"\\tilde ","s");
    fisier:=Filename(home_dir,fisier);

    step:=150;

    SizeScreen([2000,2000]);
    PrintTo(fisier,"\\textbf{Class $",Label(orb),"$:}");
    AppendTo(fisier,"\\begin{eqnarray*}\n");
    AppendTo(fisier,"&&\\mathbf{\\tilde u}=");
    AppendTo(fisier,trimmedString(latexUnipotent(uu),"u+",step));
    AppendTo(fisier,"\\\\ &&\\mathbf{u}=");
    AppendTo(fisier,trimmedString(latexUnipotent(u),"u+",step));
#    AppendTo(fisier,"\\end{eqnarray*}\n");
#    AppendTo(fisier,"\\begin{eqnarray*}\n");
    AppendTo(fisier,"\\\\ && \\mathbf{C_U(u)^{\\circ}} =");
    AppendTo(fisier,trimmedString2(latexUnipotent(info[1]),"u+",30,step));
#    AppendTo(fisier,"\\end{eqnarray*}\n");
#    AppendTo(fisier,"\\begin{eqnarray*}\n");
    AppendTo(fisier,"\\\\ && \\mathbf{Z(C_U(u)^{\\circ})^{\\circ}} =");
    AppendTo(fisier,trimmedString2(latexUnipotent(info[2]),"u+",40,step));
    AppendTo(fisier,"\\end{eqnarray*}\n");
    SizeScreen([156,39]);
end;

latexClassShortNice:=function(orb,info)
    local fisier,
    uu,u,sys,str,#info,
    step,i;

    uu:=Representative(orb);
    u:=BorelRep(orb);
    sys:=chevalleyAdj(u);
    #info:=handleClass(orb);

    fisier:=Concatenation("tex/data/",type(sys),String(rank(sys)),
                          "/char",String(Characteristic(sys)),
                          "/",Label(orb),".tex");
    fisier:=ReplacedString(fisier,"\\tilde ","s");
    fisier:=Filename(home_dir,fisier);

    step:=150;

    SizeScreen([2000,2000]);
#    PrintTo(fisier,"$$\n\\begin{array}{p{3cm} p{20cm}}\n");
#    AppendTo(fisier,"\\textbf{Class $",Label(orb),"$:} &\\\\\n ");
    PrintTo(fisier,"\\textbf{p=",Characteristic(sys),":}");
    AppendTo(fisier,"&$\\mathbf{u}=");
    AppendTo(fisier,trimmedStringNice(latexUnipotentNice(uu),"u+",step));
    AppendTo(fisier,"$\\\\\n&\n$\\mathbf{\\tilde u}=");
    AppendTo(fisier,trimmedStringNice(latexUnipotentNice(u),"u+",step));
    AppendTo(fisier,"$\\\\\n&\n$\\mathbf{C_U(\\tilde u)^{\\circ}} =");
    AppendTo(fisier,trimmedString2Nice(latexUnipotentNice(Unipotent(chevalleyAdj(u),coefficients(info[1]),Ordering(u))),"u+",30,step));
    AppendTo(fisier,"$\\\\\n&\n$\\mathbf{Z(C_G(u))^{\\circ}} =");
    AppendTo(fisier,trimmedString2Nice(latexUnipotentNice(Unipotent(chevalleyAdj(uu),coefficients(FromPositiveBorel(orb,info[2])),Ordering(uu))),"u+",40,step));
    AppendTo(fisier,"$\\\\\n");
#    AppendTo(fisier,"$\n\\end{array}\n$$");
    SizeScreen([156,39]);
end;

latexClassVeryShortNice:=function(orb,info)
    local fisier,
    uu,u,sys,str,#info,
    step,i;

    uu:=Representative(orb);
    u:=BorelRep(orb);
    sys:=chevalleyAdj(u);
    #info:=handleClass(orb);

    fisier:=Concatenation("tex/data/",type(sys),String(rank(sys)),
                          "/char",String(Characteristic(sys)),
                          "/",Label(orb),".tex");
    fisier:=ReplacedString(fisier,"\\tilde ","s");
    fisier:=Filename(home_dir,fisier);

    step:=150;

    SizeScreen([2000,2000]);
#    PrintTo(fisier,"$$\n\\begin{array}{p{3cm} p{20cm}}\n");
#    AppendTo(fisier,"\\textbf{Class $",Label(orb),"$:} &\\\\\n ");
    PrintTo(fisier,"\\textbf{p=",Characteristic(sys),":}");
    AppendTo(fisier,"&$\\mathbf{u}=");
    AppendTo(fisier,trimmedStringNice(latexUnipotentNice(uu),"u+",step));
    AppendTo(fisier,"$\\\\\n&\n$\\mathbf{\\tilde u}=");
    AppendTo(fisier,trimmedStringNice(latexUnipotentNice(u),"u+",step));
#    AppendTo(fisier,"$\\\\\n&\n$\\mathbf{C_U(\\tilde u)^{\\circ}} =");
#    AppendTo(fisier,trimmedString2Nice(latexUnipotentNice(Unipotent(chevalleyAdj(u),coefficients(info[1]),Ordering(u))),"u+",30,step));
    AppendTo(fisier,"$\\\\\n&\n$\\mathbf{Z(C_G(u))^{\\circ}} =");
    AppendTo(fisier,trimmedString2Nice(latexUnipotentNice(Unipotent(chevalleyAdj(uu),coefficients(FromPositiveBorel(orb,info[2])),Ordering(uu))),"u+",40,step));
    AppendTo(fisier,"$\\\\\n");
#    AppendTo(fisier,"$\n\\end{array}\n$$");
    SizeScreen([156,39]);
end;

latexClassesShort:=function(orbs,infos)
    local i;
    for i in [1..Length(Classes(orbs))] do
        latexClassShort(Classes(orbs)[i],infos[i]);
    od;
end;

latexClassesShortNice:=function(orbs,infos)
    local i;
    for i in [1..Length(Classes(orbs))] do
        latexClassShortNice(Classes(orbs)[i],infos[i]);
    od;
end;

latexClassesVeryShortNice:=function(orbs,infos)
    local i;
    for i in [1..Length(Classes(orbs))] do
        latexClassVeryShortNice(Classes(orbs)[i],infos[i]);
    od;
end;

latexFisiere:=function(orbs)
    local result,
    sys,tip,rang,p,
    fisier,fisi,orb;

    sys:=algebraicU(orbs);
    tip:=type(sys);
    rang:=rank(sys);
    p:=Characteristic(sys);
    fisier:=Concatenation("tex/data/all",tip,String(rang),"char",String(p),".tex");
    fisier:=Filename(home_dir,fisier);

    result:="";
    for orb in Classes(orbs) do
        result:=Concatenation(result,"\\input{data/",tip,String(rang),"/char",String(p),"/",Label(orb),".tex}\n");
        fisi:=Concatenation("tex/data/",tip,String(rang),"/char",String(p),"/",Label(orb),".tex");
        fisi:=ReplacedString(fisi,"\\tilde ","s");
        fisi:=Filename(home_dir,fisi);
        PrintTo(fisi," ");
    od;

    result:=ReplacedString(result,"\\tilde ","s");
    PrintTo(fisier,result);
end;

latexFisiereNice:=function(orbs)
    local result,inainte,dupa,
    sys,tip,rang,p,
    fisier,fisi,orb;

    sys:=algebraicU(orbs);
    tip:=type(sys);
    rang:=rank(sys);
    p:=Characteristic(sys);
    fisier:=Concatenation("tex/data/all",tip,String(rang),"char",String(p),".tex");
    fisier:=Filename(home_dir,fisier);

    dupa:="\\end{array}$";

    result:="";
    for orb in Classes(orbs) do
        inainte:=Concatenation("\\\\\n\\par\\noindent\n$\\begin{array}{p{0.7cm} p{15cm}}\n$\\mathbf{",Label(orb),"}$ &\\\\");
        result:=Concatenation(result,inainte,"\n\\input{data/",tip,String(rang),"/char",String(p),"/",Label(orb),".tex}\n",dupa,"\n");
        fisi:=Concatenation("tex/data/",tip,String(rang),"/char",String(p),"/",Label(orb),".tex");
        fisi:=ReplacedString(fisi,"\\tilde ","s");
        fisi:=Filename(home_dir,fisi);
#        PrintTo(fisi," ");
    od;

    result:=ReplacedString(result,"\\tilde ","s");
    PrintTo(fisier,result);
end;

#
#   SAGE
#

sageAdu:=function(classes)
    local sys,orbs,orbs_len,o,fisier_out,dir,dir_labels,ad,result,i,j,dim;
    sys:=algebraicU(classes);
    orbs:=Classes(classes);
    dir_labels:=Directory(Concatenation("~/workspace/ChevalleyAdj/sage/data/",type(sys),String(rank(sys))));
    dir:=Directory(Concatenation("~/workspace/ChevalleyAdj/sage/data/",
                                 type(sys),String(rank(sys)),"/char",String(Characteristic(sys))));

    SizeScreen([2000,2000]);
    for o in orbs do
        fisier_out:=Concatenation(Label(o),".sage");
        fisier_out:=ReplacedString(fisier_out,"\\tilde ","s");
        fisier_out:=Filename(dir,fisier_out);
        ad:=Descend(Representative(o),chevalleyAdj(chevalleyAdj(Representative(o))));
        ad:=Adu(ad);
        ad:=List(ad,i->List(i,j->Int(j)));
#        result:=List(ad,i->ReplacedString(String(i),", ",","));
#        result:=ReplacedString(result,", ",",");
#        result:=ReplacedString(result,"],","],\n");
#        PrintTo(fisier_out,Concatenation("rep=",result));
        PrintTo(fisier_out,"rep=");
        AppendTo(fisier_out,"[");
        dim:=Length(ad);
        for i in [1..dim] do
            AppendTo(fisier_out,"[");
            for j in [1..dim] do
                AppendTo(fisier_out,String(ad[i][j]));
                if j<>dim then AppendTo(fisier_out,","); fi;
            od;
            AppendTo(fisier_out,"]");
            if i<>dim then AppendTo(fisier_out,",\n"); fi;
        od;
        AppendTo(fisier_out,"]");
    od;

    fisier_out:=Filename(dir_labels,Concatenation("char",String(Characteristic(sys)),"Labels.sage"));
    PrintTo(fisier_out,"labels=[");
    orbs_len:=Length(orbs);
    for i in [1..orbs_len] do
         AppendTo(fisier_out,"\"",ReplacedString(Label(orbs[i]),"\\tilde ","s"),"\"");
         if i<>orbs_len then AppendTo(fisier_out,","); fi;
    od;
    AppendTo(fisier_out,"]");

    SizeScreen([156,39]);
end;

latexJordanBlocks:=function(classes)
    local sys,result,fisier_out,tex_blocks,o;
    sys:=algebraicU(classes);
    fisier_out:=Concatenation("~/workspace/ChevalleyAdj/tex/data/",
                              type(sys),String(rank(sys)),"/char",String(Characteristic(sys)),"blocks.tex");

    tex_blocks:=function(blocks)
        local result;
        result:=Concatenation(List(blocks,b->Concatenation(",",String(b[1]),"^{",String(b[2]),"}")));
        result:=result{[2..Length(result)]};
        return ReplacedString(result,"^{1}","");
    end;

    result:="";
    for o in Classes(classes) do
        result:=Concatenation(result,"$",Label(o),"$ & $",tex_blocks(AduJordanBlocks(Representative(o))),"$\\\\ \n");
    od;

    PrintTo(fisier_out,"\\begin{tabular}{|c|c|}\n \\hline \n");
    AppendTo(fisier_out,"Class & Blocks\\\\ \\hline\n");
    AppendTo(fisier_out,result,"\\hline\n");
    AppendTo(fisier_out,"\\end{tabular}\n");
end;


text:
# Several ways to do it
"Hello world!";

Print("Hello world!\n"); # No EOL appended

Display("Hello world!");

f := OutputTextUser();
WriteLine(f, "Hello world!\n");
CloseStream(f);

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

RulesFor(DCT3, rec(
    #F DCT3_Base2: DCT2_2 = F_2 * diag(1, 1/sqrt(2))
    #F
    DCT3_Base2 := rec (
	info             := "DCT3_2 -> F_2",
	isApplicable     := P -> P[1] = 2,
	rule := (P, C) -> F(2) * Diag(1, Sqrt(1/2))
    ),

    #F DCT3_Base3: DCT3_3 = (I_1 dirsum F_2)^P * (M dirsum c) * Q
    #F well-known, e.g.,
    #F   Pueschel/Moura: Discrete Cosine and Sine Transforms (in preparation)
    DCT3_Base3 := rec (
	info             := "DCT3_3 = (I_1 dirsum F_2)^P * (M dirsum c) * Q",
	isApplicable     := P -> P[1]=3,
	rule := (P, C) ->
            (DirectSum(I(1), F(2)) ^ Perm((1,2), 3)) *
	    DirectSum(Mat([[1, 1/2], [1, -1]]),
		      Diag(CosPi(1/6))) *
	    Perm((2,3), 3)
    ),

    #F DCT3_Base5: (possible base case for size 3)
    #F
    #F   DCT3_5 = P * (I_1 dirsum (SkewDCT3_2's * (pDCT8_2 tensor I_2) * B
    #F
    #F   Pueschel/Moura: Discrete Cosine and Sine Transforms, in preparation
    #F
    DCT3_Base5 := rec (
	info             := "DCT3_5 -> SkewDCT3_2's, pDCT8_2",
	isApplicable     := P -> P[1]=5,
	allChildren := P -> [[ SkewDTT(DCT3(2), 1/5),
		               SkewDTT(DCT3(2), 3/5),
			       PolyDTT(DCT8(2)) ]],
	rule := (P, C) ->
	    Perm((1,2,4,5,3), 5) *
	    DirectSum(I(1),
		      DirectSum(C[1], C[2]) * Tensor(C[3], I(2))) *
	    Mat([ [  1,   0,   -1,   0,    1 ],
		  [  1,   0,  1/2,   0,    0 ],
		  [  0,   1,    0,   0,    0 ],
		  [  0,   0,  1/2,   0,  1/2 ],
		  [  0,   0,    0,   1,    0 ] ])
    ),

    #F DCT3_Orth9: DCT3_9 = hand-derived radix-3 Cooley-Tukey for RDFT
    #F
    DCT3_Orth9 := rec (
	info             := "DCT3_9 -> TRDFT_3's",
	isApplicable     := P -> P[1] = 9,
	allChildren := P -> [ [ RDFT(3).transpose() ] ],
	rule := (P, C) ->
	   Perm( (1,5)(2,7,8,4,9,6), 9 ) *
	   Tensor(C[1],	I(3)) *
	   DirectSum(
	       I(3),
	       (DirectSum(I(2), Rot(2/9), Rot(4/9))) ^ Perm((2,4,5,3), 6)) *
	   Tensor(I(3), C[1]) *
	   DirectSum(
	       I(3),
	       Mat([ [ 1,  0,  0, 0, 0, 0 ],
		     [ 0,  1,  0, 0, 0, 1 ],
		     [ 0,  0,  1, 0, 1, 0 ],
		     [ 0,  0,  0, 1, 0, 0 ],
		     [ 0,  0, -1, 0, 1, 0 ],
		     [ 0, -1,  0, 0, 0, 1 ] ])) *
	   Diag(1, -1, 1, -1, 1, -1, 1, -1, 1) *
	   Perm( (2,7,8)(3,4)(5,9), 9 )
    ),

    #F DCT3_toSkewDCT3:  DCT3_n = DCT3_n(1/2)
    #F
    DCT3_toSkewDCT3 := rec (
	info             := "DCT3_n -> SkewDCT3_n(1/2)",
	forTransposition := false,
	isApplicable     := P -> IsPrimePowerInt(P[1]) and P[1] mod 2 = 0 and P[1] > 2,
	allChildren := P -> [ [ SkewDTT(DCT3(P[1]), 1/2) ] ],
	rule := (P, C) -> C[1]
    )
));


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

myScript := function()
    local k, filename;

    for k in [16, 32, 64, 128] do
        filename := ConcatenationString("fft", String(k));
        HDLSynthesize(streamDFTUnroll(k, 2, 4), 1, 0, 16, 300, 0, filename);
    od;
end;