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

RewriteRules(RulesRC, rec(
    RC_ISumAcc := Rule([RC, @(1, ISumAcc)], e -> ISumAcc(@(1).val.var, @(1).val.domain, RC(@(1).val.child(1)))),
    RC_COND := Rule([RC, @(1, COND)], e -> ApplyFunc(COND, [ @(1).val.cond ]::List(@(1).val._children, i->RC(i)))),
));

RewriteRules(RulesFuncSimp, rec(
    H_fTensor_fBase_fId := ARule(fCompose, [@(1, H, e->e.domain()+1 = e.range() and e.params[4] = 1),
        [@(0, fTensor), @(2, fBase), @(3, fId)]],
        e-> [fAdd(@(1).val.range(), @(3).val.domain(), @(1).val.params[3] + @(0).val.at(0))]
    ),
    fTensor_fBase_fId_H := Rule(@(1, fTensor, e->(
            ForAll(e.children(), c->ObjId(c) in [fBase, fId, H] and When(ObjId(c) = H, c.domain() = 1, true))
            and ForAny(e.children(), c->ObjId(c) = H))),
        e->fTensor(List(@(1).val.children(), c->When(ObjId(c) = H, fBase(c.params[1], c.params[3]), c)))
    )
));

RewriteRules(RulesSums, rec(
    Scat_ISumAcc := ARule(Compose,
       [ @(1, Scat), @(2, ISumAcc) ],
        e -> [ CopyFields(@(2).val,
        rec(_children :=  List(@(2).val._children, c -> ScatAcc(@(1).val.func) * c), dimensions := [Rows(@(1).val), Cols(@(2).val)] )) ]),

    PropagateBB := Rule([@(1,BB), @(2,Compose, e->ForAny(e.children(), i->ObjId(i)=COND))],
        e -> Compose(List(@(2).val.children(), i->Cond(ObjId(i) = COND, ApplyFunc(COND, [i.cond]::List(i.children(), j->BB(j))), BB(i))))),

    MergeBB := ARule(Compose, [@(1, BB), @(2, BB)], e-> [ BB(@(1).val.child(1) * @(2).val.child(1)) ]),

    COND_true := Rule(@(1, COND, e->e.cond = V(true)), e -> @(1).val.child(1)),

    eq_false := Rule([@(1, eq), [add, @(3, var, e->e.t=TInt), @(2, Value, e->e.v >0) ], @(0, Value, e->e.v=0)], e-> V(false)),

    eq_false2 := Rule([@(1, eq), @(2, add, e->ForAny(e.args, j->IsValue(j) and j.v > 0)), @(0, Value, e->e.v=0)], e-> V(false)),

    COND_false := Rule(@(1, COND, e->e.cond = V(false)), e -> @(1).val.child(2)),

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

    COND_Diag :=  ARule( Compose, [ @(1, COND), @(2, [Diag, RCDiag]) ],
        e -> [ CopyFields(@(1).val, rec(_children := List(@(1).val._children, c -> c * @(2).val),
                dimensions := [Rows(@(1).val), Cols(@(2).val)] )) ])
));

RewriteRules(RulesStrengthReduce, rec(
    mul_cond := Rule([@(1, mul), [cond, @(2), @(3, Value), @(4, Value)], @(5)],
        e->cond(@(2).val, @(5).val*@(3).val, @(5).val*@(4).val)),

    eq_add := Rule([eq, [add, @(1,Value), @(2)], @(3,Value)], e->eq(@(2).val, @(3).val-@(1).val)),

    eq_false := Rule([@(1, eq), [add, @(3, var, e->e.t=TInt), @(2, Value, e->e.v >0) ], @(0, Value, e->e.v=0)], e-> V(false)),

    eq_false2 := Rule([@(1, eq), @(2, add, e->ForAny(e.args, j->IsValue(j) and j.v > 0)), @(0, Value, e->e.v=0)], e-> V(false)),

    fix_cond := Rule([eq, [mul, @(0,Value, e->e.t=TInt and e.v>0), @(1,var,IsLoopIndex)], @(2, Value, e->e.t=TInt and e.v<0)], e->V(false)),
));

fixISumAccPDLoop := function(s, opts)
    local ss, vr, scat, dim, krnl, lvar, it0, itn, sct, isum, its, sx, scand, scats;

    while Length(Collect(s, ISumAcc)) >0  do
        scand := Collect(s, ISumAcc);
        ss := scand[1];
        vr := ss.var;
        scats := Collect(ss, ScatAcc);
        if Length(scats) > 1 then Error("There should be only one ScatAcc in an ISumAcc..."); fi;
        scat := scats[1];
        dim := Cols(scat);

        krnl := ss.child(1);
        lvar := Ind(vr.range-1);

        it0 := RulesSums(SubstTopDown(SubstVars(Copy(krnl), rec((vr.id) := V(0))), ScatAcc, e-> Scat(fId(dim))));
        itn := RulesSums(RulesSums(SubstTopDown(SubstVars(Copy(krnl), rec((vr.id) := lvar+1)), ScatAcc, e-> ScatAcc(fAdd(dim, dim, 0)))));
        sct := Scat(Collect(krnl, ScatAcc)[1].func);
        isum := ISum(lvar, itn);
        isum.doNotMarkBB := true;
        its := BB(sct) * SUM(it0, isum);

        s := SubstTopDown(s, ss, e->its);
        s := MergedRuleSet(RulesSums, RulesDiag, RulesDiagStandalone, RulesStrengthReduce, RulesRC)(s);
    od;

    return s;
end;

fixScatter := function(s, opts)
    local tags, scts, sct, f, fl, vars, fll, rng, it, itspace, vals, vby2, vby2d, vby2d2, accf2, flnew, sctnew, srec;

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

    # very brittle and danger of infinite loop
    scts := Collect(s, @(1, Scat, e->Cols(e) = 4 and Length(Collect(e, BH)) > 0 and ForAny(e.free(), i->i.t <> TInt)));
    while Length(scts) > 0  do
        sct := scts[1];

        f := sct.func;
        fl := f.lambda();
        vars := Filtered(fl.free(), i->i.t = TInt);
        fll := fl.tolist();
        rng := List(vars, i->i.range);
        itspace := Cartesian(List(rng, i->[1..i]));

        vals := [];
        for it in itspace do
            srec := rec();
            for i in [1..Length(vars)] do
               srec := CopyFields(srec, rec((vars[i].id):=V(it[i]-1)));
            od;
            vals := vals :: List(Copy(fll), e->EvalScalar(SubstVars(e, srec)));
        od;

        vby2 := List([1..Length(vals)/2], i->vals[2*i-1]/2);
        vby2d := fTensor(FData(List(vby2, i->V(i))), fId(2));

        vby2d2 := FData(List(vby2, i->V(i)));
        accf2 := fTensor(List(vars, fBase) :: [fId(fl.vars[1].range/2)]);
        flnew := fTensor(fCompose(vby2d2, accf2).lambda().setRange(Rows(sct)/2), fId(2));

        sctnew := Scat(flnew);
        s := SubstTopDownNR(s, sct, e->sctnew);
        # very brittle and danger of infinite loop
        scts := Collect(s, @(1, Scat, e->Cols(e) = 4 and Length(Collect(e, BH)) > 0 and ForAny(e.free(), i->i.t <> TInt)));
    od;
    s := CopyFields(tags, s);

    return s;
end;

fixGather := function(s, opts)
    local tags, gaths, gath, f, fl, vars, fll, rng, it, itspace, vals, vby2, vby2d, vby2d2, accf2, flnew, gathnew, srec;

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

    # very brittle and danger of infinite loop
    gaths := Collect(s, @(1, Gath, e->Rows(e) = 4 and Length(Collect(e, BH)) > 0 and ForAny(e.free(), i->i.t <> TInt)));
    while Length(gaths) > 0  do
        gath := gaths[1];

        f := gath.func;
        fl := f.lambda();
        vars := Filtered(fl.free(), i->i.t = TInt);
        fll := fl.tolist();
        rng := List(vars, i->i.range);
        itspace := Cartesian(List(rng, i->[1..i]));

        vals := [];
        for it in itspace do
            srec := rec();
            for i in [1..Length(vars)] do
               srec := CopyFields(srec, rec((vars[i].id):=V(it[i]-1)));
            od;
            vals := vals :: List(Copy(fll), e->EvalScalar(SubstVars(e, srec)));
        od;

        vby2 := List([1..Length(vals)/2], i->vals[2*i-1]/2);
        vby2d := fTensor(FData(List(vby2, i->V(i))), fId(2));

        vby2d2 := FData(List(vby2, i->V(i)));
        accf2 := fTensor(List(vars, fBase) :: [fId(fl.vars[1].range/2)]);
        flnew := fTensor(fCompose(vby2d2, accf2).lambda().setRange(Cols(gath)/2), fId(2));

        gathnew := Gath(flnew);
        s := SubstTopDownNR(s, gath, e->gathnew);
        # very brittle and danger of infinite loop
        gaths := Collect(s, @(1, Gath, e->Rows(e) = 4 and Length(Collect(e, BH)) > 0 and ForAny(e.free(), i->i.t <> TInt)));
    od;
    s := CopyFields(tags, s);
    return s;
end;

fixInitFunc := function(c, opts)
    local cis, ci, ci2, fc, freei, offending, myd, tags;

    tags := rec();
    if IsBound(c.dimensions) then tags.dimensions  := c.dimensions; fi;
    if IsBound(c.ruletree) then tags.ruletree  := c.ruletree; fi;

    cis := Collect(c, @(1, func, e->e.id="init"));

    for ci in cis do
        fc := c.free();
        freei := ci.free();
        offending := Filtered(fc, i->i in freei);
        for myd in offending do
           ci2 := func(ci.ret, ci.id, ci.params,
               data(myd, myd.value, ci.cmd));
           c := SubstTopDown(c, @(1, func, e->e = ci), e->ci2);
        od;
    od;

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

fixUpSigmaSPL := function(s, opts)
    s := fixISumAccPDLoop(s, opts);
    s := fixScatter(s, opts);
    s := fixGather(s, opts);

    s := SubstTopDown(s, [@(1,RCDiag, e->e.dims()=[2,2]),
           @(2, FDataOfs, e->IsBound(e.var.value) and e.var.value = V(Flat(List([1..e.var.value.t.size / 2], i->[V(1.0), V(0.0)])))), @(3, I)], e->I(2));

    s := SubstTopDown(s, [@(1,RCDiag, e->e.dims()=[4,4]), @(2, FDataOfs,
            e->IsBound(e.var.value) and ForAll([1..e.var.value.t.size/4], i->e.var.value.v[4*i-3] = V(1.0)) and ForAll([1..e.var.value.t.size/4],
            i->e.var.value.v[4*i- 2] = V(0.0))), @(3, I)],
        e->let(vals := V(Flat(List([1..@(2).val.var.value.t.size/4], i->[@(2).val.var.value.v[4*i-1], @(2).val.var.value.v[4*i]]))),
            dt := Dat(vals.t).setValue(vals),
            fdo := FDataOfs(dt, 2, V(2) * @(2).val.ofs.free()[1]),
            rcdnew  := RCDiag(fdo, I(2)),
            SumsSPL(DirectSum(I(2), rcdnew), opts))
    );
    return s;
end;

fixReplicatedData := function(c, opts)
    local datas, d, replicas, cc, r, tags;

    tags := rec();
    if IsBound(c.dimensions) then tags.dimensions  := c.dimensions; fi;
    if IsBound(c.ruletree) then tags.ruletree  := c.ruletree; fi;

    c  := SubstBottomUp(c, @(1, data, e->IsArray(e.value.t) and ObjId(e.cmd)=data and IsArray(e.cmd.value.t) and e.value.t.size < e.cmd.value.t.size),
        e->data(@(1).val.cmd.var, @(1).val.cmd.value, data(@(1).val.var, @(1).val.value, @(1).val.cmd.cmd)));

    datas := Collect(c, data);
    for d in datas do
        replicas := Collect(d.cmd, @(1, data,
            e->e.value.t.size <= d.value.t.size and e.value.v = d.value.v{[1..e.value.t.size]}));
        for r in replicas do
            cc := r.cmd;
            cc := SubstVars(cc, rec((r.var.id) := d.var));
            c := SubstTopDown(c, @(1, data, e->e.var = r.var), e->cc);
        od;
    od;
    c := CopyFields(tags, c);
    return c;
end;

scalarizeAccumulator := function(c, opts)
    local scand, svr, svars, tags;

    tags := rec();
    if IsBound(c.dimensions) then tags.dimensions  := c.dimensions; fi;
    if IsBound(c.ruletree) then tags.ruletree  := c.ruletree; fi;

    scand := Filtered(Set(Flat(List(Collect(c, decl), i->i.vars))), e->IsArray(e.t) and ForAll(Collect(c, @(1, nth, f-> f.loc = e)), g->IsValue(g.idx)));
    for svr in scand do
        svars := List([1..svr.t.size], e->var.fresh_t("q", svr.t.t));
        c := SubstTopDown(c, @(1, decl, e->svr in e.vars), e-> decl(svars :: Filtered(@(1).val.vars, i-> i <>svr), @(1).val.cmd));
        c := SubstTopDown(c, @(1, nth, e->e.loc = svr), e->svars[e.idx.v+1]);
    od;
    c := CopyFields(tags, c);
    return c;
end;

max_stars_repo_name: broderickpt/spiral-package-fftx

max_issues_repo_name: broderickpt/spiral-package-fftx

max_forks_repo_name: broderickpt/spiral-package-fftx


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

#############################################################################
##
#M  FrancyMessage( <messageType>, <string>, <string> )  . .  create a Label
##
InstallMethod(FrancyMessage,
  "message type, a title, a value",
  true,
  [IsFrancyMessageType,
   IsString,
   IsString],
  0,
function(messageType, title, value)
local id;
id := GenerateID();
  return Objectify(FrancyMessageObjectType, rec(
    id    := id,
    type  := messageType!.value,
    title := title,
    text  := value
  ));
end);

InstallOtherMethod(FrancyMessage,
  "a title, a value",
  true,
  [IsString,
   IsString],
  0,
function(title, value)
  return FrancyMessage(FrancyMessageType.DEFAULT, title, value);
end);

InstallOtherMethod(FrancyMessage,
  "message type, a value",
  true,
  [IsFrancyMessageType,
   IsString],
  0,
function(messageType, value)
  return FrancyMessage(messageType, "", value);
end);

InstallOtherMethod(FrancyMessage,
  "a value",
  true,
  [IsString],
  0,
function(value)
  return FrancyMessage(FrancyMessageType.DEFAULT, "", value);
end);

ext: gi

max_stars_repo_name: LaGuer/francy

max_issues_repo_name: LaGuer/francy

max_forks_repo_name: LaGuer/francy


text:
# The idiomatic way to compute with polynomials

x := Indeterminate(Rationals, "x");

# This is a value in a polynomial ring, not a function
p := 6*x^3 - 4*x^2 + 7*x - 19;

Value(p, 3);
# 128

u := CoefficientsOfUnivariatePolynomial(p);
# [ -19, 7, -4, 6 ]

# One may also create the polynomial from coefficients
q := UnivariatePolynomial(Rationals, [-19, 7, -4, 6], x);
# 6*x^3-4*x^2+7*x-19

p = q;
# true

# Now a Horner implementation
Horner := function(coef, x)
	local v, c;
	v := 0;
	for c in Reversed(coef) do
		v := x*v + c;
	od;
	return v;
end;

Horner(u, 3);
# 128

ext: gap

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


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

Class(PkRDFT12_Base, PRDFT_Base, rec(
    dims := self >> [self.params[1], self.params[1]],
    terminate := self >> let(N := self.params[1], k := self.params[2],
	rr := When(self.transposed, Cols(self), Rows(self)),
	mat := Cond(IsEvenInt(N),
            Mat(Concatenation(
                    [List([0..N-1], c-> When(IsEvenInt(c), 1,1))],
                    [List([0..N-1], c-> When(IsEvenInt(c), 1,-1))],
                    List([2..rr-1], r -> When(r mod 2 = 0,
                        List([0..N-1], c -> self.projRe(self.omega(N,k,Int(r/2),c))),
                        List([0..N-1], c -> self.projIm(self.omega(N,k,Int(r/2),c))))))),
            Mat(Concatenation(
                    [List([0..N-1], c -> 1)],
                    List([2..rr], r -> When(r mod 2 = 0,
                        List([0..N-1], c -> self.projRe(self.omega(N,k,Int(r/2),c))),
                        List([0..N-1], c -> self.projIm(self.omega(N,k,Int(r/2),c)))))))),
        When(self.transposed, mat.transpose(), mat)),
));

Class(PkRDFT1, PkRDFT12_Base, rec(
    abbrevs := [ (n) -> Checked(IsPosInt(n),  [n, 1]),
                 (n,k) -> Checked(IsPosIntSym(n), IsIntSym(k),
                                  IsSymbolic(n) or IsSymbolic(k) or Gcd(n,k)=1, [n, k mod n]) ],
    omega := (N,k,r,c) -> E(N)^(k*r*c)
));
Class(PkRDFT2, PkRDFT12_Base, rec(
    abbrevs := [ (n) -> Checked(IsPosIntSym(n),  [n, 1]),
                 (n,k) -> Checked(IsPosIntSym(n), IsIntSym(k),
                                  IsSymbolic(n) or IsSymbolic(k) or Gcd(n,k)=1, [n, k mod (2*n)]) ],
    omega := (N,k,r,c) -> E(2*N)^(k*r*(2*c+1))
));

Class(PkDHT1, PkRDFT1, rec(
    terminate := self >> let(n := self.params[1], k := self.params[2], i := When(IsEvenInt(n), 2, 1),
        DirectSum(I(i), Tensor(I(Int((n-1)/2)), F(2))) *
        PkRDFT1(n, k).terminate()
    )
));

Class(PkDHT2, PkRDFT2, rec(
    terminate := self >> let(n := self.params[1], k := self.params[2], i := When(IsEvenInt(n), 2, 1),
        DirectSum(I(i), Tensor(I(Int((n-1)/2)), Diag(1,-1)*F(2))) *
        PkRDFT2(n, k).terminate()
    )
));

Class(URealDFT_Base, NonTerminal, rec(
    abbrevs := [ (n) -> Checked(IsPosInt(n),  [n, 1]),
                 (n,k) -> Checked(IsPosInt(n), IsInt(k), Gcd(n,k)=1, [n, k mod n]) ],
    terminate := self >> let(N := self.params[1], k := self.params[2],
	rr := When(self.transposed, Cols(self), Rows(self)),
	Mat(Concatenation(
	    [List([0..N-1], c-> When(IsEvenInt(c), 1,0))],
	    [List([0..N-1], c-> When(IsOddInt(c), 1,0))],
	    List([2..rr-1], r -> When(r mod 2 = 0,
		 List([0..N-1], c -> self.projRe(self.omega(N,k,Int(r/2),c))),
		 List([0..N-1], c -> self.projIm(self.omega(N,k,Int(r/2),c)))))))),
    isReal := True,
    SmallRandom := () -> Random([2..16]),
    LargeRandom := () -> 2 ^ Random([6..15])
));

Class(URDFT_Base, URealDFT_Base, rec(
    projRe := Re,
    projIm := Im,
    hashAs := self >> self
));

Class(URDFT, URDFT_Base, rec(
    dims := self >> let(n:=self.params[1], [ 2*(Int((n+1)/2)), n ]),
    omega := (N,k,r,c) -> E(N)^(k*r*c),
));

URDFT1:=URDFT;

# Regularized Rules
# R1 -> (R1, C1) (R1')
# R2 -> (R2, C2) (R1')
PRF12_CTReg_Children := (N,k,PRFt,DFTt,PRF1prime, maxRadix) -> Map2(
    Filtered(DivisorPairs(N),d->d[1] <= maxRadix/2 and d[2]<>2),
    (m,n) -> When(IsEvenInt(n),
	[ PRFt(2*m,k), DFTt(m,k), PRF1prime(n,k) ],
	[ PRFt(m,k),   DFTt(m,k), PRF1prime(n,k) ] )
);

PRF12_CTReg_Rule := (N,k,C,Conj,Tw) -> let(mm:=Cols(C[1]), m:=When(IsEvenInt(N), mm/2, mm),
    n:=Cols(C[3]), Nf:=Int(N/2), Nc:=Int((N+1)/2),
    nf:=Int(n/2), nc:=Int((n+1)/2), mf:=Int(m/2), mc:=Int((m+1)/2), j:=Ind(nc-1),

    SUM(
	BB(RC(Scat(H(Nc,m,0,n/2))) * C[1] * L(mm,2) * RC(Gath(H(nc*m, m, 0, nc)))),

	When(nc=1, [],
	ISum(j, BB(
	     RC(Scat(BH(Nc,N,m,j+1,n))) *
	     Conj * RC(C[2]) * Tw(j) *
	     RC(Gath(H(nc*m, m, j+1, nc))))))
    ) *
    Tensor(I(m), C[3]) * L(N,m)
);

PRF12_CTReg2_Rule := (N,k,C,Conj,Tw) -> let(mm:=Cols(C[1]), m:=When(IsEvenInt(N), mm/2, mm),
    n:=Cols(C[3]), Nf:=Int(N/2), Nc:=Int((N+1)/2),
    nf:=Int(n/2), nc:=Int((n+1)/2), mf:=Int(m/2), mc:=Int((m+1)/2), j:=Ind(nc-1),

    SUM(
	RC(Scat(H(Nc,m,0,n/2))) * C[1] * RC(Gath(H(nc*m, m, 0, 1))),

	When(nc=1, [],
	ISum(j,
	     RC(Scat(BH(Nc,N,m,j+1,n))) *
	     Conj * RC(C[2]) * Tw(j) * L(2*m, m) *
	     RC(Gath(H(nc*m, m, m*(j+1), 1)))))
    ) *
    Tensor(C[3], I(m))
);

PRF12_CTReg3_Rule := (N,k,C,Conj,Tw) -> let(mm:=Cols(C[1]), m:=When(IsEvenInt(N), mm/2, mm),
    n:=Cols(C[3]), Nf:=Int(N/2), Nc:=Int((N+1)/2),
    nf:=Int(n/2), nc:=Int((n+1)/2), mf:=Int(m/2), mc:=Int((m+1)/2), j:=Ind(nc-1),

    SUM(
	RC(Scat(H(Nc,m,0,n/2))) * C[1] * L(2*m, 2) * Gath(H(N, 2*m, 0, nc)),

	When(nc=1, [],
	ISum(j,
	     RC(Scat(BH(Nc,N,m,j+1,n))) *
	     Conj * RC(C[2]) * Tw(j) *
	     Gath(H(N, 2*m, j+1, nc))))
    ) *
    Tensor(I(m), L(n,2)*C[3]) * L(N,m)
);

PRF12_CTReg4_Rule := (N,k,C,Conj,Tw) -> let(mm:=Cols(C[1]), m:=When(IsEvenInt(N), mm/2, mm),
    n:=Cols(C[3]), Nf:=Int(N/2), Nc:=Int((N+1)/2),
    nf:=Int(n/2), nc:=Int((n+1)/2), mf:=Int(m/2), mc:=Int((m+1)/2), j:=Ind(nc-1),

    RC(Tensor(I(m/2), DirectSum(I(n/2+1), J(n/2-1)))) *
    RC(L(N/2, m)) *
    DirectSum(
        C[1] * L(2*m, 2),
	IterDirectSum(j, RC(M(m,m/2)) * Conj * RC(C[2]) * Tw(j))
    ) *
    L(N, nc) *
    Tensor(I(m), L(n,2)*C[3]) * L(N,m)
);

PRF12_CTReg5_Rule := (N,k,C,Conj,Tw) -> let(mm:=Cols(C[1]), m:=When(IsEvenInt(N), mm/2, mm),
    n:=Cols(C[3]), Nf:=Int(N/2), Nc:=Int((N+1)/2),
    nf:=Int(n/2), nc:=Int((n+1)/2), mf:=Int(m/2), mc:=Int((m+1)/2), j:=Ind(nc-1),

    RC(Tensor(I(m/2), DirectSum(I(n/2+1), J(n/2-1)))) *
    RC(L(N/2, m)) *
    DirectSum(
        C[1],
	IterDirectSum(j, RC(M(m,m/2)) * Conj * RC(C[2]) * Tw(j) * L(2*m, m))
    ) *
    Tensor(C[3], I(m))
);

unrc := m -> let(r := Dimensions(m)[1], c := Dimensions(m)[2],
    MatSPL(Gath(H(r,r/2,0,2))) * m * MatSPL(Scat(H(c,c/2,0,2))));

RulesFor(URDFT1, rec(
    URDFT1_Base1 := BaseRule(URDFT1, [1, @]),
    URDFT1_Base2 := BaseRule(URDFT1, [2, @]),
    URDFT1_Base4 := BaseRule(URDFT1, [4, @]),

    URDFT1_CT := rec(
        maxRadix := 32,
	isApplicable := P -> not IsPrime(P[1]) and not (IsEvenInt(P[1]) and IsPrime(P[1]/2)),
	allChildren  := (self, P) >> PRF12_CTReg_Children(P[1], P[2], URDFT1, DFT1, URDFT1, self.maxRadix),
	rule := (P,C) -> let(N:=P[1], k:=P[2], m:=When(IsEvenInt(N), Cols(C[1])/2, Cols(C[1])),
	    PRF12_CTReg_Rule(N, k, C, Diag(BHD(m,1,-1)),
                j -> RC(DirectSum(I(1), Diag(fPrecompute(fCompose(Twid(N,m,k,0,0,j+1), fAdd(m, m-1, 1))))))))),

    # slightly different variant of same rule,
    URDFT1_CT2 := rec(
	isApplicable := P -> not IsPrime(P[1]) and not (IsEvenInt(P[1]) and IsPrime(P[1]/2)),
	allChildren  := P -> PRF12_CTReg_Children(P[1], P[2], URDFT1, DFT1, URDFT1),
	rule := (P,C) -> let(N:=P[1], k:=P[2], m:=When(IsEvenInt(N), Cols(C[1])/2, Cols(C[1])),
	    PRF12_CTReg2_Rule(N, k, C, Diag(BHD(m,1,-1)), j->RC(Diag(fPrecompute(Twid(N,m,k,0,0,j+1))))))),

    URDFT1_CT3 := rec(
	isApplicable := P -> not IsPrime(P[1]) and not (IsEvenInt(P[1]) and IsPrime(P[1]/2)),
	allChildren  := P -> PRF12_CTReg_Children(P[1], P[2], URDFT1, DFT1, URDFT1),
	rule := (P,C) -> let(N:=P[1], k:=P[2], m:=When(IsEvenInt(N), Cols(C[1])/2, Cols(C[1])),
	    PRF12_CTReg3_Rule(N, k, C, Diag(BHD(m,1,-1)), j->RC(Diag(fPrecompute(Twid(N,m,k,0,0,j+1))))))),

    URDFT1_CT4 := rec(
	isApplicable := P -> not IsPrime(P[1]) and not (IsEvenInt(P[1]) and IsPrime(P[1]/2)),
	allChildren  := P -> PRF12_CTReg_Children(P[1], P[2], URDFT1, DFT1, URDFT1),
	rule := (P,C) -> let(N:=P[1], k:=P[2], m:=When(IsEvenInt(N), Cols(C[1])/2, Cols(C[1])),
	    PRF12_CTReg4_Rule(N, k, C, Diag(BHD(m,1,-1)), j->RC(Diag(fPrecompute(Twid(N,m,k,0,0,j+1))))))),

    URDFT1_CT5 := rec(
	isApplicable := P -> not IsPrime(P[1]) and not (IsEvenInt(P[1]) and IsPrime(P[1]/2)),
	allChildren  := P -> PRF12_CTReg_Children(P[1], P[2], URDFT1, DFT1, URDFT1),
	rule := (P,C) -> let(N:=P[1], k:=P[2], m:=When(IsEvenInt(N), Cols(C[1])/2, Cols(C[1])),
	    PRF12_CTReg5_Rule(N, k, C, Diag(BHD(m,1,-1)), j->RC(Diag(fPrecompute(Twid(N,m,k,0,0,j+1))))))),

));


text:
%Headers
/.
        //
        // These functions are only needed in order to use a lexer as
        // a special purpose "diff" program. To do so, one can invoke
        // the main program below directly from the command line.
        //
        private static int LINES = 0,
                           TOKENS = 1;
        private static int differ_mode = LINES; // default
        private static String extension = "";
        private static int changeCount = 0,
                           insertCount = 0,
                           deleteCount = 0,
                           moveCount = 0;

        private static void compareFiles(String old_file, String new_file)
        {
            try
            {
                $action_type old_lexer, new_lexer;
                if (old_file.equals(""))
                {
                    char [] input_chars = new char[0];
                    old_lexer = new $action_type(input_chars, "null_file");
                }
                else old_lexer = new $action_type(old_file);

                if (new_file.equals(""))
                {
                    char [] input_chars = new char[0];
                    new_lexer = new $action_type(input_chars, "null_file");
                }
                else new_lexer = new $action_type(new_file);

                PrsStream old_stream = new PrsStream(old_lexer.getILexStream());
                old_lexer.lexer(old_stream);

                PrsStream new_stream = new PrsStream(new_lexer.getILexStream());
                new_lexer.lexer(new_stream);

                Differ diff = (differ_mode == TOKENS ? (Differ) new DifferTokens(old_stream, new_stream)
                                                     : (Differ) new DifferLines(old_stream, new_stream));
                diff.compare();

                if (diff.getChangeCount() > 0)
                {
                    diff.outputChanges();

                    changeCount += diff.getChangeCount();
                    insertCount += (diff.getInsertCount() + diff.getReplaceInsertCount());
                    deleteCount += (diff.getDeleteCount() + diff.getReplaceDeleteCount());
                    moveCount += diff.getMoveCount();
                }
            }
            catch (Exception e)
            {
                System.err.println(e.getMessage());
                e.printStackTrace();
            }
        }

        private static void compareDirectories(java.io.File old_dir, java.io.File new_dir)
        {
            try
            {
                java.io.File old_file[] = old_dir.listFiles(),
                             new_file[] = new_dir.listFiles();
                java.util.HashMap old_map = new java.util.HashMap();
                for (int i = 0; i < old_file.length; i++)
                {
                    String name = old_file[i].getName();
                    if (old_file[i].isDirectory() || name.endsWith(extension))
                        old_map.put(name, old_file[i]);
                }

                for (int i = 0; i < new_file.length; i++)
                {
                    java.io.File file = (java.io.File) old_map.get(new_file[i].getName());
                    if (file != null)
                    {
                        old_map.remove(new_file[i].getName());

                        if (file.isDirectory() && new_file[i].isDirectory())
                             compareDirectories(file, new_file[i]);
                        else compareFiles(file.getPath(), new_file[i].getPath());
                    }
                    else if (new_file[i].isDirectory() ||
                             new_file[i].getName().endsWith(extension))
                    {
                        String s = new_file[i].getName() +
                                   " found in directory " +
                                   new_dir.getPath() +
                                   " does not exist in directory " +
                                   old_dir.getPath();
                        System.err.println("*Warning: " + s);

                        if (! new_file[i].isDirectory())
                            compareFiles("", new_file[i].getPath());
                    }
                }

                for (java.util.Iterator i = old_map.entrySet().iterator(); i.hasNext(); )
                {
                    java.util.Map.Entry e = (java.util.Map.Entry) i.next();
                    java.io.File file = (java.io.File) e.getValue();

                    String s = file.getName() +
                               " not found in directory " +
                               new_dir.getPath();
                    System.err.println("*Warning: " + s);

                    if (! file.isDirectory())
                        compareFiles(file.getPath(), "");
                }
            }
            catch (Exception e)
            {
                System.err.println(e.getMessage());
                e.printStackTrace();
            }
        }

        public static void main(String[] args)
        {
            String new_file = null,
                   old_file = null;
            boolean help = false;

            int i;
            for (i = 0; i < args.length; i++)
            {
                if (args[i].charAt(0) == '-')
                {
                    if (args[i].equals("-ext"))
                         extension = (i + 1 < args.length ? args[++i] : "");
                    else if (args[i].equals("-h"))
                         help = true;
                    else if (args[i].equals("-l"))
                         differ_mode = LINES;
                    else if (args[i].equals("-t"))
                         differ_mode = TOKENS;
                }
                else break;
            }
            if (i < args.length)
            {
                new_file = args[i++];
                old_file = new_file; // assume only one file is specified
            }
            if (i < args.length)
                old_file = args[i++];
            for (; i < args.length; i++)
                System.err.println("Invalid argument: " + args[i]);

            if (help || (new_file == null &&  old_file == null))
            {
                System.out.println();
                System.out.println("Usage: diff [OPTION]... file1 [file2]");
                System.out.println("Compute stats for file1 or compare file1 to file2 statement by statement.");
                System.out.println();
                System.out.println("  -ext s -- if file1 and file2 are directories, compare only files that end\n" +
                                   "            with the extension (suffix) s.");
                System.out.println("  -h     -- print this help message");
                System.out.println("  -l     -- compare line by line instead of statement by statement");
                System.out.println("  -t     -- compare token by token instead of statement by statement");
            }
            else if (old_file.equals(new_file))
            {
                java.io.File old_dir = new java.io.File(old_file);
                // if (old_dir.isDirectory())
                //     computeStats(old_dir);
                // else computeStats(old_file);

                System.out.println("*** No difference ***");
                // System.out.println("    Number of files: " + fileCount);
                // System.out.println("    Number of lines: " + lineCount);
                // System.out.println("    Number of types (classes/interfaces): " + (classCount + interfaceCount) + " (" + classCount + "/" + interfaceCount + ")");
                // System.out.println("    Number of statements: " + statementCount);
                // System.out.println("    Number of braces (left/right): (" + leftBraceCount + "/" + rightBraceCount + ")");
            }
            else
            {
                java.io.File old_dir = new java.io.File(old_file),
                     new_dir = new java.io.File(new_file);
                if (old_dir.isDirectory() && new_dir.isDirectory())
                     compareDirectories(old_dir, new_dir);
                else compareFiles(old_file, new_file);

                if (changeCount == 0)
                    System.out.println("***** No difference *****");
                else
                {
                    System.out.println("***** " +
                                       changeCount +
                                       " different " +
                                       (changeCount == 1 ? "section" : "sections") + " *****");
                    System.out.println("    " + moveCount    + " statements moved");
                    System.out.println("    " + insertCount  + " statements inserted");
                    System.out.println("    " + deleteCount  + " statements deleted");
                }
            }

            return;
        }
./
%End

max_stars_repo_name: kuafuwang/LPG-VScode

max_issues_repo_name: kuafuwang/LPG-VScode

max_forks_repo_name: kuafuwang/LPG-VScode


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

Class(SIMDchain, chain);

update_opts_params:= function(opts, new_param)
   opts.params := opts.params::[new_param];
end;

update_opts_includes:= function(opts, new_param)
   opts.includes := opts.includes::[new_param];
end;

Class(HCOLFuncGen, rec(
    generate := (self, ss, opts) >>  let(
          err_chk := When(opts.errorCheck,
            let( x:= var.fresh_t("x", T_Real(32)), j := Ind(1),
                 ch := List(Filtered(opts.params, v->v.t in [T_Real(32), T_Real(64)]),
                                        vv-> PointWise(1, Lambda([x, j], cond(isinf(vv), 1, cond(isnan(vv), 2, 0))))),
                 n := Length(ch),
                 t :=
                 OLCompose(
                    VStack(
                        Reduction(n, (a, b)->add(a, b), V(0), False),
                        I(n)
                    ),
                    VStack(    ch )
                 ),
                 tt := Rewrite(SumsSPL(t, opts), [RulesSumsHCOL, RulesTerminateReductionHCOL, RulesSumsHCOL], opts),
                 rtn_var := var("err", TPtr(TBool)), rtn := var("err_rtn", TBool),
                 update_opts_params(opts, rtn_var),
                 update_opts_includes(opts, "<math.h>"),
                 decl([rtn],
                    chain(
                      SubstBottomUp(opts.codegen.OLCompose(tt, rtn_var, self.getX(opts), opts),
                                    [nth, @(4, var, e->e.id = rtn_var.id), @(5, Value, e->e.v=0)], e->rtn),
                      creturnCond(neq(rtn, 0), -2)
                 ))),
            chain([])),

          code := opts.codegen.(ss.name)(ss, self.getY(opts), self.getX(opts), opts),

          finalCode := When(opts.useCReturn,
                let( w:= var.fresh_t("w", self.getRetType(opts)),
                     decl([w],
                        chain(
                            When(opts.errorCheck, err_chk, chain([])),
                            SIMDchain(
                                SubstBottomUp(code,
                                    [nth, @(2, var, e->e.id = self.getY(opts).id),
                                          @(3, Value, e->e.v=0)], e->w),
                               creturn(w))))),
             SIMDchain(err_chk, code)),

        func(self.getRetType(opts), "transform", self.getParams(finalCode, opts), finalCode)) ,

        getParams := (self, fc, opts) >> let(

            y := self.getY(opts), x := self.getX(opts),
            usedX := Length(Collect(fc, @(1, var, e->e.id=x.id)))>0,
            Cond(
                opts.useCReturn and not usedX,       opts.params,
                not opts.useCReturn and not usedX,   [y]::opts.params,
                opts.useCReturn and usedX,           [x]::opts.params,
                not opts.useCReturn and usedX,       [y, x]::opts.params
        )
    ),
    getRetType := (self, opts) >> When(opts.useCReturn, opts.YType.t, TVoid),
    getX :=  (self, opts) >> CopyFields(X, rec(t:=opts.XType)),
    getY :=  (self, opts) >> CopyFields(Y, rec(t:=opts.YType))
));

max_stars_repo_name: spiral-software/spiral-package-hcol

max_issues_repo_name: spiral-software/spiral-package-hcol

max_forks_repo_name: spiral-software/spiral-package-hcol


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

RandInput := () -> Complex(
    RandomList([1..2^27])/(1.0*2^27),
    RandomList([1..2^27])/(1.0*2^27));

RandValExp := function(exp)
    local e, rargs;
    if IsBound(exp.randVal) then
	return exp.randVal;

    elif IsValue(exp) then
	exp.randVal := exp.v;
	return exp.randVal;

    elif IsLoc(exp) then
	exp.randVal := RandInput();
	#Print("RandInput : ", exp, " = ", exp.randVal, "\n");
	return exp.randVal;
    else
	rargs := List(exp.args, RandValExp);
	e := ApplyFunc(ObjId(exp), rargs);
	e := e.ev();
	Constraint(not IsExp(e) and not IsLoc(e));
	exp.randVal := e;
	return exp.randVal;
    fi;
end;

RandValAssign := function(cmd)
    local r;
    if IsBound(cmd.loc.randVal) then return cmd.loc.randVal;
    elif IsBound(cmd.exp.randVal) then
	cmd.loc.randVal := cmd.exp.randVal;
	return cmd.loc.randVal;
    else
	r := RandValExp(cmd.exp);
	cmd.exp.randVal := r;
	cmd.loc.randVal := r;
	return r;
    fi;
end;

Class(RCSE, rec(
    cseLookup := (self, randVal) >> HashLookup(self.csetab, V(randVal)),

    cseAdd := meth(self, cmd)
        local op, ent, v;
	if not IsBound(cmd.exp.isExpComposite) then return; fi;
	HashAdd(self.csetab, V(cmd.loc.randVal), cmd.loc);
	HashAdd(self.csetab, V(-cmd.loc.randVal), neg(cmd.loc));
    end,

    cseCmd := meth(self, cmd)
        local lkup;
	RandValAssign(cmd);
	lkup := self.cseLookup(cmd.loc.randVal);
	if lkup <> false then
	    cmd.exp := lkup;
	else
	    self.cseAdd(cmd);
	fi;
    end,

    remove_node := (self, exp) >> Chain(
	HashDelete(self.csetab, V(exp.randVal)),
	HashDelete(self.csetab, V(-exp.randVal))),

    node := meth(self, t, exp)
        local lkup, v;
	if IsValue(exp) or IsLoc(exp) then
	    RandValExp(exp);
	    return exp;
	fi;

	exp := ApplyFunc(ObjId(exp), List(exp.args, a->self.node(t, a)));
	RandValExp(exp);

	lkup := self.cseLookup(exp.randVal);
	if lkup <> false then
	    return lkup;
	else
	    v := var.fresh_t("t", t);
	    v.randVal := exp.randVal;
	    v.def := assign(v, exp);
	    self.cseAdd(v.def);
	    return v;
	fi;
    end,

    flush := meth(self)
        self.csetab := HashTable( (val,size) -> val.t.hash(val.v, size) );
    end,

    __call__ := meth(self, code)
        self.csetab := HashTable( (val,size) -> val.t.hash(val.v, size) );
	DoForAll(Collect(code, assign), cmd -> self.cseCmd(cmd));
	return code;
    end
));


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

# ==========================================================================
# Mon(<perm-func>, <value-diag>, <structural-diag>) = perm * diag
# ==========================================================================
Class(Mon, BaseMat, SumsBase, rec(
    new := (self, perm, vdiag, sdiag) >> Checked(
       domain(perm) = domain(vdiag) and domain(perm) = domain(sdiag),
       SPL(WithBases(self,
		rec( perm := perm,
		     vdiag := vdiag,
		     sdiag := sdiag ))).setDims()),

    #-----------------------------------------------------------------------
    dims := self >> [self.perm.domain(), self.perm.range()],
    setDims := meth(self) self.dimensions := self.dims(); return self; end,
    #-----------------------------------------------------------------------
    sums := self >> self,
    rChildren := self >> [self.perm, self.vdiag, self.sdiag],
    rSetChild := rSetChildFields("perm", "vdiag", "sdiag"),

    children := ~.rChildren,
    setChild := ~.rSetChild,
    child := (self, n) >> Cond(n=1, self.perm, n=2, self.vdiag, n=3, self.sdiag,
	Error("<n> must be in [1..3]")),
    #-----------------------------------------------------------------------
    isPermutation := False,
    #-----------------------------------------------------------------------
    isReal := self >> self.vdiag.isReal(),
    #-----------------------------------------------------------------------
    toAMat := self >> AMatSPL( Gath(self.perm) * Diag(self.vdiag) * Diag(self.sdiag) ),
    #-----------------------------------------------------------------------
    transpose := self >>   # we use inherit to copy all fields of self
        Inherit(self, rec(
		perm := self.perm.transpose(),
		vdiag := fCompose(self.vdiag, self.perm),
		sdiag := fCompose(self.sdiag, self.perm),
		dimensions := Reversed(self.dimensions))),
    #-----------------------------------------------------------------------
    print := (self, i, is) >>
       Print(self.name, "(", self.perm, ", ", self.vdiag, ", ", self.sdiag, ")"),
    #-----------------------------------------------------------------------
    arithmeticCost := (self, costMul, costAddMul) >> Sum(self.vdiag.tolist(), costMul)
));

Declare(SMon, GMon, SGMon);

# ==========================================================================
# GMon(<gath-func>, <value-diag>, <structural-diag>) = diag * gath
# ==========================================================================
Class(GMon, Mon, rec(
    #-----------------------------------------------------------------------
    dims := self >> [domain(self.perm), range(self.perm)],
    #-----------------------------------------------------------------------
    toAMat := self >> AMatSPL( self.vdiag * self.sdiag * Gath(self.perm) ),
    #-----------------------------------------------------------------------
    transpose := self >> SMon(self.perm, self.vdiag, self.sdiag)
    #-----------------------------------------------------------------------
));

# ==========================================================================
# SMon(<scat-func>, <value-diag>, <structural-diag>) = scat * diag
# ==========================================================================
Class(SMon, Mon, rec(
    #-----------------------------------------------------------------------
    dims := self >> [range(self.perm), domain(self.perm)],
    #-----------------------------------------------------------------------
    toAMat := self >> AMatSPL( Scat(self.perm) * self.vdiag * self.sdiag ),
    #-----------------------------------------------------------------------
    transpose := self >> GMon(self.perm, self.vdiag, self.sdiag)
    #-----------------------------------------------------------------------
));

# ==========================================================================
# SGMon(<scat-func>, <gath-func>, <value-diag>, <structural-diag>) = scat * diag * gath
# ==========================================================================
Class(SGMon, Mon, rec(
    new := (self, scat, gath, vdiag, sdiag) >> Checked(
       domain(scat) = domain(vdiag) and domain(vdiag) = domain(sdiag) and
       domain(sdiag) = domain(gath),
       SPL(WithBases(self,
		rec( scat := scat,
		     gath := gath,
		     vdiag := vdiag,
		     sdiag := sdiag ))).setDims()),

    rChildren := self >> [self.scat, self.gath, self.vdiag, self.sdiag],
    rSetChild := rSetChildFields("scat", "gath", "vdiag", "sdiag"),
    child := (self, n) >> Cond(n=1, self.scat,  n=2, self.gath,
	                       n=3, self.vdiag, n=4, self.sdiag,
			       Error("<n> must be in [1..3]")),
    #-----------------------------------------------------------------------
    dims := self >> [range(self.scat), range(self.gath)],
    #-----------------------------------------------------------------------
    toAMat := self >> AMatSPL(Scat(self.scat) * self.vdiag * self.sdiag *
	                      Gath(self.gath)),
    #-----------------------------------------------------------------------
    transpose := self >> SGMon(self.gath, self.scat, self.vdiag, self.sdiag),
    #-----------------------------------------------------------------------
    print := (self, i, is) >>
       Print(self.name, "(", self.scat, ", ", self.gath, ", ", self.vdiag, ", ", self.sdiag, ")")
 ));


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

NewRulesFor(GT, rec(

# Why must we do parallelization and streaming in one big rule? Because they
# must both be done together. If done separately, we'll get two sets of outer
# perms (LxI)s, and we don't want that. In addition, the parallel breakdown
# rule has to tell a special mbuf rule how many iterations to multibuffer for,
# and a) this is messy, and b) this is trivial to handle in one single rule,
# and therefore best done within a single rule.

#F GT_ParStreamCore
#F Params: (p, u). p=# of procs, u=minimum packet size
    GT_ParStreamCore := rec(
        kmin          := 16, #128  # Minimum kernel size. DMA packet size = 4*k bytes for single precision real
        kmax          := 2048,  # Based on Max packet size for single precision (2048*4*2 = 16 Kbytes) NOTE: HARDCODED!

        requiredFirstTag := ParStream,

        # Variable naming scheme: either (In x Am) or (Am x In)

        applicable := (self, t) >> let(
            rank := Length(t.params[4]),     # Eg.: (Am x In)
            g    := t.params[2],
            s    := t.params[3],

            m    := Rows(t.params[1]),       # (splsize)
            n    := t.params[4][1],          # Available iterations (n)
            p    := t.firstTag().params[1],  # procs
            u    := t.firstTag().params[2],  # Specified minimum packet size

            b    := Maximum((m*n)/(self.kmax*p), 4), # Number of multibuffered iterations
            r    := n/(b*p),                 # Remaining iterations.

            # Conditions:
            rank = 1
            and PatternMatch(t, [GT, @(1), @(2,XChain), @(3,XChain), ...], empty_cx())
            and r >= 1       # If not, we have too few iterations for mbuf+par
            and When(g=GTPar and s=GTPar, r*m >= u,  r >= u)  # Respect specified packet size
        ),

        children := (self, t) >> let(
            g    := t.params[2],
            s    := t.params[3],

            m    := Rows(t.params[1]),
            n    := t.params[4][1],
            p    := t.firstTag().params[1],  # procs

            b    := Maximum((m*n)/(self.kmax*p), 4), # Number of multibuffered iterations
            r    := n/(b*p),                         # Remaining iterations.

            pkSizeG := When(g=GTPar, r*m, r),
            pkSizeS := When(s=GTPar, r*m, r),

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

            When(r > 1,
              [ [ GT(t.params[1], g, s, [r]).withTags(tags),
                InfoNt(b, pkSizeG, pkSizeS) ] ],
              [ [ t.params[1].withTags(tags), InfoNt(b, pkSizeG, pkSizeS) ] ]
            )
        ),

        apply := (self, t, C, Nonterms) >> let(
            dft := Nonterms[1],
            g   := t.params[2],
            s   := t.params[3],
            n   := t.params[4][1],
            m   := Rows(t.params[1]),
            mc  := Cols(t.params[1]),
            p   := t.firstTag().params[1],
            b   := Nonterms[2].params[1],
            pkSizeG := Nonterms[2].params[2],
            pkSizeS := Nonterms[2].params[3],
            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])),

            i      := var("spuid", TInt, p),
            j      := Ind(b),

            sp     := When(s=GTPar,
                        ScatDist(p, n*m/p, p, i),
                        ScatSend(ss.part(1, i, m, [p, 1]), n/p, p, i)),

            gp     := When(g=GTPar,
                        GathDist(p, n*m/p, p, i),
                        GathRecv(gg.part(1, i, mc, [p, 1]), n/p, p, i)),

            sm     := When(s=GTPar,
                        ScatMem(s.part(1, j, 1, [b]), pkSizeS),
                        ScatMem(s.part(1, j, m,  [b]), pkSizeS)),

            gm     := When(g=GTPar,
                        GathMem(g.part(1, j, 1, [b]), pkSizeG),
                        GathMem(g.part(1, j, mc, [b]), pkSizeG)),

            DistSum(p, i, p,
                sp * MultiBufISum(j, b, sm, C[1], gm) * gp
            )

        ) #apply (let)
    ) #GT_ParStreamCore
));

NewRulesFor(GT, rec(
#F GT_StreamParChip
#F Params: (p, um, uc=u/p). p=# of procs, um=min pksize (mem), uc=min pksize (on-chip)
    GT_StreamParChip := rec(
        kmin          := 16, #128  # Minimum kernel size. DMA packet size = 4*k bytes for single precision real
        kmax          := 2048,  # Based on Max packet size for single precision (2048*4*2 = 16 Kbytes) NOTE: HARDCODED!

        requiredFirstTag := ParStream,

        # Variable naming scheme: either (In x Am) or (Am x In)
        applicable := (self, t) >> let(
            rank := Length(t.params[4]),     # Eg.: (Am x In)
            g    := t.params[2],
            s    := t.params[3],

            m    := Rows(t.params[1]),       # (splsize)
            n    := t.params[4][1],          # Available iterations (n)
            p    := t.firstTag().params[1],  # procs
            u    := t.firstTag().params[2],  # Specified minimum packet size (mem)

            b    := Maximum((m*n)/(self.kmax*p), 4), # Number of multibuffered iterations
            r    := n/b,                     # Remaining iterations.

            # Conditions:
            rank = 1
            and PatternMatch(t, [GT, @(1), @(2,XChain), @(3,XChain), ...], empty_cx())
            and r >= 1       # If not, we have too few iterations for mbuf+par
            and When(g=GTPar and s=GTPar, r*m >= u,  r >= u)  # Respect specified packet size
            #and not GT_ParStreamCore.applicable(t) # For now, prioritize the GT_ParStreamCore rule
        ),

        children := (self, t) >> let(

            g    := t.params[2],
            s    := t.params[3],

            m    := Rows(t.params[1]),
            n    := t.params[4][1],
            p    := t.firstTag().params[1],  # procs
            u    := t.firstTag().params[2],  # Specified minimum packet size (mem)
            uc   := u/p,                     # Minimum packet size (on-chip)

            b    := Maximum((m*n)/(self.kmax*p), 4), # Number of multibuffered iterations
            r    := n/b,                     # Remaining iterations.

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

            # Add a parallel tag to child
            When(r > 1,
               [[ GT(t.params[1], g, s, [r]).setTags(Concatenation([ParCell(p, uc)], tags)) ]],
               [[ t.params[1].setTags(Concatenation([ParCell(p, uc)], tags)) ]]
            )
        ),

        apply := (self, t, C, Nonterms) >> let(
            dft     := Nonterms[1],
            g       := t.params[2],
            s       := t.params[3],

            m    := Rows(t.params[1]),
            n    := t.params[4][1],
            p    := t.firstTag().params[1],  # procs
            u    := t.firstTag().params[2],  # Specified minimum packet size (mem)

            b    := Maximum((m*n)/(self.kmax*p), 4), # Number of multibuffered iterations
            r    := n/b,                             # Remaining iterations.

            pkSizeS := When(s=GTPar, r*m, r),
            pkSizeG := When(g=GTPar, r*m, r),

            i      := Ind(b),
            sm     := When(s=GTPar,
                        ScatMem(s.part(1, i, 1, [b]), pkSizeS),
                        ScatMem(s.part(1, i, Rows(dft)/pkSizeS, [b/pkSizeS]), pkSizeS)),
            gm     := When(g=GTPar,
                        GathMem(g.part(1, i, 1, [b]), pkSizeG),
                        GathMem(g.part(1, i, Cols(dft)/pkSizeG, [b/pkSizeG]), pkSizeG)),

            MultiBufISum(i, b,
                sm,
                C[1],
                gm
            )

        ) #let/apply
    ), #GT rule

));

# NOTE: Works if the specified u is the max u. Else, fails.
# NOTE: Should convert this rule so that u is automatically picked.
NewRulesFor(GT, rec(
#F This is the ParStream version of the MBufCell_spec vector recursion rule
#F Automatically choses k,n, subject to a specified minimum
    GT_DFT_ParStreamCore_VecRecur := rec(
        kmin          := 16, #128  # Minimum kernel size. DMA packet size = 4*k bytes for single precision real
        kmax          := 2048,  # Maximum kernel size on a single SPE

        requiredFirstTag := ParStream,

        # Variable naming scheme: either (In x Am) or (Am x In)

        #applicable := false, # THIS RULE IS INCORRECT!

        applicable := (self, t) >> let(
            rank := Length(t.params[4]),     # Eg.: (Am x In)
            g    := t.params[2],
            s    := t.params[3],

            k    := t.params[4][1],          # Available iterations
            mn   := Rows(t.params[1]),       # (splsize)
            p    := t.firstTag().params[1],  # procs
            u    := t.firstTag().params[2],  # Specified minimum packet size

            mns  := Filtered(DivisorPairs(mn),
                       i->let(m:=i[1], n:=i[2],
                          n >= u                # n is a packet size in second factor
                          and m >= 4*p          # multibuffer and parallelize around m iterations in second factor
                          and m*n >= self.kmin  # First factor kernel size is big enough
                          and m*n <= self.kmax  # First  factor fits on SPEs
                          and k*n <= self.kmax) # Second factor fits on SPEs
                    ),
#           Error("BP--applicable"),
            rank = 1
#           and not GT_ParStreamCore.applicable(t)  # Regular rule is not applicable
            and g=GTVec and s=GTPar             # Only consider (IxD)L
            and PatternMatch(t, [GT, DFT, XChain, XChain, @, @, @], empty_cx()) # Only consider GTs with DFT for a kernel
            and k >= 4*p                  # First loop is multibuffered and parallelized around k
            and k >= u                    # k is also the one of the packet sizes (2nd factor)
            and Length(mns) >= 1          # We can break the DFT based on all our constraints
        ),

        children := (self, t) >> let(
            g    := t.params[2],
            s    := t.params[3],
            k    := t.params[4][1],          # Available iterations (n)
            mn   := Rows(t.params[1]),       # (splsize)
            p    := t.firstTag().params[1],  # procs
            u    := t.firstTag().params[2],  # Specified minimum packet size

            tags := Drop(t.getTags(), 1),
            #oldtags := Drop(t.getTags(), 1),
            #tags := Concatenation([ParStream(p, u)], oldtags),

            mns  := Filtered(DivisorPairs(mn),
                       i->let(m:=i[1], n:=i[2],
                          n >= u                # n is a packet size in second factor
                          and m >= 4*p          # multibuffer around m iterations in second factor
                          and m*n >= self.kmin  # First factor kernel size is big enough
                          and m*n <= self.kmax  # First  factor fits on SPEs
                          and k*n <= self.kmax) # Second factor fits on SPEs
                    ),

            Map2(mns, (m,n)-> [
              GT(DFT(m,1), GTVec, GTVec, [n]).withTags(tags),
              GT(DFT(n,1), GTVec, GTPar, [k]).withTags(tags),
              InfoNt(m, n)
              ])
        ),

        apply := (self, t, C, Nonterms) >> let(
            k    := t.params[4][1],
            m    := Nonterms[3].params[1],
            n    := Nonterms[3].params[2],
            p    := t.firstTag().params[1],  # procs

            j1    := Ind(k/p),
            j2    := Ind(m/p),
            pvar := var("spuid", TInt, p),
            pkSize1 := Rows(C[1]),

            #Error("BP----"),

            # WORKS:
            #MultiBufISum(j1, k,
            #      ScatMem(GTPar.part(1, j1, 1, [k]), pkSize1),
            #      C[1] * Diag(fPrecompute(Tw1(m*n, n, 1))),
            #      GathMem(GTPar.part(1, j1, 1, [k]), pkSize1)
            #)
            #*
            #MultiBufISum(j2, m,
            #     ScatMem(fTensor(fId(k), fBase(j2)), n),
            #     C[2],
            #     GathMem(fTensor(fId(n), fBase(j2)), k)
            #)

            perms := L(k*m, m),
            ixfs  := fTensor(fBase(pvar), fId(m*k/p)),

            permg := L(m*n, m),
            ixfg  := fTensor(fBase(pvar), fId(m*n/p)),

            DistSum(p, pvar, p,
                ScatDist(p, k*n*m/p, p, pvar) *
                MultiBufISum(j1, k/p,
                  ScatMem(GTPar.part(1, j1, 1, [k]), pkSize1),
                  C[1] * Diag(fPrecompute(Tw1(m*n, n, 1))),
                  GathMem(GTPar.part(1, j1, 1, [k]), pkSize1)
                ) * GathDist(p, k*n*m/p, p, pvar)
            ) *
            DistSum(p, pvar, p,
               ScatSend(fCompose(perms, ixfs), n, p, pvar) *

               MultiBufISum(j2, m/p,
                 ScatMem(fTensor(fBase(j2), fId(k)), n),
                  C[2],
                 GathMem(fTensor(fBase(j2), fId(n)), k)
               ) *

               GathRecv(fCompose(permg, ixfg), k, p, pvar)
            )

        ) #let/apply
    )
));

NewRulesFor(GT, rec(
#F Vector recursion using streaming to entire chip
#F Automatically choses k,n, subject to a specified minimum pksize
    GT_DFT_StreamParChip_VecRecur := rec(
        kmin          := 128,  # Minimum kernel size. DMA packet size = 4*k bytes for single precision real
        kmax          := 2048,  # Maximum kernel size on a single SPE

        requiredFirstTag := ParStream,

        # Variable naming scheme: either (In x Am) or (Am x In)

        applicable := (self, t) >> let(
            rank := Length(t.params[4]),     # Eg.: (Am x In)
            g    := t.params[2],
            s    := t.params[3],

            k    := t.params[4][1],          # Available iterations
            mn   := Rows(t.params[1]),       # (splsize)
            p    := t.firstTag().params[1],  # procs
            u    := t.firstTag().params[2],  # Specified minimum packet size

            mns  := Filtered(DivisorPairs(mn),
                       i->let(m:=i[1], n:=i[2],
                          n >= u                  # n is a packet size in second factor
                          and m >= 4              # multibuffer around m iterations in second factor
                          and m*n >= self.kmin    # First factor kernel size is big enough
                          and m*n <= self.kmax*p  # First  factor fits on SPEs
                          and k*n <= self.kmax*p) # Second factor fits on SPEs
                    ),
            rank = 1
#          and not GT_ParStreamCore.applicable(t)  # None of the 3 regular rules is applicable: GT_StreamParChip and GT_ParStreamCore should also not be applicable
#           and not GT_DFT_ParStreamCore_VecRecur.applicable(t)
           and not GT_StreamParChip.applicable(t)
            and g=GTVec and s=GTPar             # Only consider (IxD)L
            and PatternMatch(t, [GT, DFT, XChain, XChain, @, @, @], empty_cx()) # Only consider GTs with DFT for a kernel
            and k >= 4                    # First loop is multibuffered around k
            and k >= u                    # k is also the one of the packet sizes (2nd factor)
            and Length(mns) >= 1          # We can break the DFT based on all our constraints
        ),

        children := (self, t) >> let(
            g    := t.params[2],
            s    := t.params[3],
            k    := t.params[4][1],          # Available iterations (n)
            mn   := Rows(t.params[1]),       # (splsize)
            p    := t.firstTag().params[1],  # procs
            u    := t.firstTag().params[2],  # Specified minimum packet size (mem)
            uc   := u/p,                     # Minimum packet size (on-chip)

            oldtags := Drop(t.getTags(), 1),
            tags := Concatenation([ParCell(p, uc)], oldtags),

            mns  := Filtered(DivisorPairs(mn),
                       i->let(m:=i[1], n:=i[2],
                          n >= u                  # n is a packet size in second factor
                          and m >= 4              # multibuffer around m iterations in second factor
                          and m*n >= self.kmin    # First factor kernel size is big enough
                          and m*n <= self.kmax*p  # First  factor fits on SPEs
                          and k*n <= self.kmax*p) # Second factor fits on SPEs
                    ),

            Map2(mns, (m,n)-> [
              GT(DFT(m,1), GTVec, GTVec, [n]).withTags(tags),
              GT(DFT(n,1), GTVec, GTPar, [k]).withTags(tags),
              InfoNt(m, n)
              ])
        ),

        apply := (self, t, C, Nonterms) >> let(
            k    := t.params[4][1],
            m    := Nonterms[3].params[1],
            n    := Nonterms[3].params[2],

            i    := Ind(k),
            j    := Ind(m),
            pkSize1 := Rows(C[1]),

            MultiBufISum(i, k,
                ScatMem(GTPar.part(1, i, 1, [k]), pkSize1),
                C[1] * Diag(fPrecompute(Tw1(m*n, n, 1))),
                GathMem(GTPar.part(1, i, 1, [k]), pkSize1)
            )
            *
            MultiBufISum(j, m,
                ScatMem(fTensor(fId(k), fBase(j)), n),
                C[2],
                GathMem(fTensor(fId(n), fBase(j)), k)
            )

        ) #let/apply
    )
));

#NOTES
#-----------------------------------------------------------------------------
# Mon 11 Jan 2010 12:56:58 AM EST

# Outline for cleanup:

# We need to separate and disentangle all these rules: vecrecur should be
# independent of ParStreamCore or StreamParChip. ParStreamCore and
# StreamParChip should simply split the parallelism and let the lower level
# mbuf and par rules handle the rest. SigmaSPL rewrite rules should either
# /naturally/ be able to understand whether something is ParStreamCore or
# StreamParChip and act accordingly, or, if this is not possible naturally,
# this should be marked up somehow at the breakdown-rule level. Ugly ways for
# SigmaSPL rules to 'try' to figure out what's going on is as bad as Clippy.


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

Class(RealDFT_Base, TaggedNonTerminal, rec(
    abbrevs := [ (n) -> Checked(IsPosInt(n),  [n, 1]),
                 (n,k) -> Checked(IsPosInt(n), IsInt(k), Gcd(n,k)=1, [n, k]) ],

    terminate := self >> let(N := self.params[1], k := self.params[2],
	rr := When(self.transposed, Cols(self), Rows(self)),
	mat := Mat(List([0..rr-1], r -> When(r mod 2 = 0,
                    List([0..N-1], c -> self.projRe(self.omega(N,k,Int(r/2),c))),
                    List([0..N-1], c -> self.projIm(self.omega(N,k,Int(r/2),c)))))),
        When(self.transposed, mat.transpose(), mat)),

    toAMat := self >> self.terminate().toAMat(),

    isReal := True,
    SmallRandom := () -> Random([2..16]),
    LargeRandom := () -> 2 ^ Random([6..15]),

    normalizedArithCost := self >> let(n := self.params[1],
       floor(2.5 * n * log(n) / log(2.0))),

    hashAs := self >> let(t:=ObjId(self)(self.params[1], 1).withTags(self.getTags()),
        When(self.transposed, t.transpose(), t))

));

Class(PDHT13_Base, RealDFT_Base, rec(
    projRe := w -> Re(w)+Im(w),
    projIm := w -> Re(w)-Im(w)
));

Class(PDHT24_Base, RealDFT_Base, rec(
    projRe := w -> Re(w)+Im(w),
    projIm := w -> -Re(w)+Im(w)
));

Class(PRDFT_Base, RealDFT_Base, rec(
    projRe := Re,
    projIm := Im
));

Class(IPRDFT_Base, TaggedNonTerminal, rec(
    abbrevs := [ (n) -> Checked(IsPosInt(n),  [n, 1]),
	         (n,k) -> Checked(IsPosInt(n), IsInt(k), Gcd(n,k)=1, [n, k]) ],
    terminate := self >> let(n:=self.params[1], k:=self.params[2],
	self.rdft(n, -k).terminate().transpose() * self.diag(n)),
    isReal := True,
    SmallRandom := () -> Random([2..16]),
    LargeRandom := () -> 2 ^ Random([6..15]),

    diag0 := n -> When(n=2, I(4),
        Diag(diagDirsum(fConst(TReal, 2, 1), fConst(TReal, 2*Int((n-1)/2), 2),
                When(IsEvenInt(n), fConst(TReal, 2, 1), [])))),

    diag1 := n -> Diag(diagDirsum(fConst(TReal, 2*Int(n/2), 2),
                          When(IsOddInt(n), fConst(TReal, 2, 1), []))),

    normalizedArithCost := self >> let(n := self.params[1],
       floor(2.5 * n * log(n) / log(2.0))),

    hashAs := self >> let(t:=ObjId(self)(self.params[1], 1).withTags(self.getTags()),
        When(self.transposed, t.transpose(), t))
));

Declare(IPRDFT, IPRDFT1, IPRDFT2, IPRDFT3, IPRDFT4);

# PRDFT(<N>, <k>) - Packed Real DFT Nonterminal
#   Is a N+2 x N real matrix. N+2 outputs correspond
#   to Floor(N/2)+1 complex outputs of complex DFT
#   the other half of complex DFT outputs are complex
#
Class(PRDFT, PRDFT_Base, rec(
    abbrevs := [ (n) -> Checked(IsPosInt(n),  [n, 1]),
                 (n,k) -> Checked(IsPosInt(n), IsInt(k), Gcd(n,k)=1, [n, k mod n]) ],
    dims := self >> let(n:=self.params[1], d := [2*(Int(n/2)+1), n], When(self.transposed, Reversed(d), d)),
    omega := (N,k,r,c) -> E(N)^(k*r*c),
    inverse := self >> IPRDFT1(self.params[1], -self.params[2]),
    inverseViaDiag := self >> self.transpose() * IPRDFT.diag(self.params[1])
));
PRDFT1 := PRDFT;

Class(PRDFT2, PRDFT_Base, rec(
    abbrevs := [ (n) -> Checked(IsPosInt(n),  [n, 1]),
                 (n,k) -> Checked(IsPosInt(n), IsInt(k), Gcd(n,k)=1, [n, k mod (2*n)]) ],
    dims := self >> let(n:=self.params[1], d := [2*(Int(n/2)+1), n], When(self.transposed, Reversed(d), d)),
    omega := (N,k,r,c) -> E(2*N)^(k*r*(2*c+1)),
    inverse := self >> IPRDFT3(self.params[1], -self.params[2])
));

Class(PRDFT3, PRDFT_Base, rec(
    abbrevs := [ (n) -> Checked(IsPosInt(n),  [n, 1]),
                 (n,k) -> Checked(IsPosInt(n), IsInt(k), Gcd(n,k)=1, [n, k mod (2*n)]) ],
    dims := self >> let(n:=self.params[1], d := [2*(Int((n+1)/2)), n], When(self.transposed, Reversed(d), d)),
    omega := (N,k,r,c) -> E(2*N)^(k*(2*r+1)*c),
    inverse := self >> IPRDFT2(self.params[1], -self.params[2])
));

Class(PRDFT4, PRDFT_Base, rec(
    abbrevs := [ (n) -> Checked(IsPosInt(n),  [n, 1]),
                 (n,k) -> Checked(IsPosInt(n), IsInt(k), Gcd(n,k)=1, [n, k mod (4*n)]) ],
    dims := self >> let(n:=self.params[1], d := [2*(Int((n+1)/2)), n], When(self.transposed, Reversed(d), d)),
    omega := (N,k,r,c) -> E(4*N)^(k*(2*r+1)*(2*c+1)),
    inverse := self >> IPRDFT4(self.params[1], -self.params[2])
));

# Class(SkewRDFT3, PRDFT_Base, rec(
#     abbrevs := [ (n) -> Checked(IsPosInt(n),  [n, 1, 0, 0]),
#       (n,k) -> Checked(IsPosInt(n), IsInt(k), Gcd(n,k)=1, [n, k, 0, 0]),
#       (n,k,lr,rr) -> Checked(IsPosInt(n), IsInt(k), IsRat(lr), IsRat(rr),
# 	  Gcd(n,k)=1, [n, k, lr, rr])
#     ],
#     dims := self >> let(n:=self.params[1], [ 2*Int((n+1)/2), n ]),
#     omega := (self,N,k,r,c) >> let(
# 	lnum:=Numerator(self.params[3]), lden:=Denominator(self.params[3]),
# 	rnum:=Numerator(self.params[4]), rden:=Denominator(self.params[4]),
# 	E(2*N*lden*rden) ^ (k*(lden*(2*r+1)+lnum)*(rden*c+rnum)))
# ));
# Class(SkewRDFT4, SkewRDFT3, rec(
#     dims := self >> let(n:=self.params[1], [ 2*Int((n+1)/2), n ]),
#     omega := (self,N,k,r,c) >> let(
# 	lnum:=Numerator(self.params[3]), lden:=Denominator(self.params[3]),
# 	rnum:=Numerator(self.params[4]), rden:=Denominator(self.params[4]),
# 	E(4*N*lden*rden) ^ (k*(lden*(2*r+1)+lnum)*(rden*(2*c+1)+rnum)))
# ));

#
# Hartley
#
Class(PDHT, PDHT13_Base, rec(
    abbrevs := [ (n) -> Checked(IsPosInt(n),  [n, 1]),
                 (n,k) -> Checked(IsPosInt(n), IsInt(k), Gcd(n,k)=1, [n, k mod (n)]) ],
    dims := self >> let(n:=self.params[1], d := [2*(Int(n/2)+1), n], When(self.transposed, Reversed(d), d)),
    omega := PRDFT.omega
));
PDHT1 := PDHT;

Class(PDHT2, PDHT24_Base, rec(
    abbrevs := [ (n) -> Checked(IsPosInt(n),  [n, 1]),
                 (n,k) -> Checked(IsPosInt(n), IsInt(k), Gcd(n,k)=1, [n, k mod (2*n)]) ],
    dims := self >> let(n:=self.params[1], d := [2*(Int(n/2)+1), n], When(self.transposed, Reversed(d), d)),
    omega := PRDFT2.omega
));

Class(PDHT3, PDHT13_Base, rec(
    abbrevs := [ (n) -> Checked(IsPosInt(n),  [n, 1]),
                 (n,k) -> Checked(IsPosInt(n), IsInt(k), Gcd(n,k)=1, [n, k mod (2*n)]) ],
    dims := self >> let(n:=self.params[1], d := [2*(Int((n+1)/2)), n], When(self.transposed, Reversed(d), d)),
    omega := PRDFT3.omega
));

Class(PDHT4, PDHT24_Base, rec(
    abbrevs := [ (n) -> Checked(IsPosInt(n),  [n, 1]),
                 (n,k) -> Checked(IsPosInt(n), IsInt(k), Gcd(n,k)=1, [n, k mod (4*n)]) ],
    dims := self >> let(n:=self.params[1], d := [2*(Int((n+1)/2)), n], When(self.transposed, Reversed(d), d)),
    omega := PRDFT4.omega
));

# IPRDFT(<N>, <k>) - Inverse Packed Real DFT Nonterminal
#   Is a N x N+2 real matrix. When applied to N+2 outputs
#   of PRDFT yields the original real input.
Class(IPRDFT, IPRDFT_Base, rec(
    abbrevs := [ (n) -> Checked(IsPosInt(n),  [n, 1]),
                 (n,k) -> Checked(IsPosInt(n), IsInt(k), Gcd(n,k)=1, [n, k mod (n)]) ],
    dims := self >> let(n:=self.params[1], d := [n, 2*(Int(n/2)+1)], When(self.transposed, Reversed(d), d)),
    rdft := PRDFT,
    diag := IPRDFT_Base.diag0
));
IPRDFT1 := IPRDFT;

# IPRDFT2 is inverse if PRDFT3
Class(IPRDFT2, IPRDFT_Base, rec(
    abbrevs := [ (n) -> Checked(IsPosInt(n),  [n, 1]),
                 (n,k) -> Checked(IsPosInt(n), IsInt(k), Gcd(n,k)=1, [n, k mod (2*n)]) ],
    dims := self >> let(n:=self.params[1], d := [n, 2*(Int((n+1)/2))], When(self.transposed, Reversed(d), d)),
    rdft := PRDFT3,
    diag := IPRDFT_Base.diag1
));

# IPRDFT3 is inverse if PRDFT2
Class(IPRDFT3, IPRDFT_Base, rec(
    abbrevs := [ (n) -> Checked(IsPosInt(n),  [n, 1]),
                 (n,k) -> Checked(IsPosInt(n), IsInt(k), Gcd(n,k)=1, [n, k mod (2*n)]) ],
    dims := self >> let(n:=self.params[1], d := [n, 2*(Int(n/2)+1)], When(self.transposed, Reversed(d), d)),
    rdft := PRDFT2,
    diag := IPRDFT_Base.diag0
));

Class(IPRDFT4, IPRDFT_Base, rec(
    abbrevs := [ (n) -> Checked(IsPosInt(n),  [n, 1]),
                 (n,k) -> Checked(IsPosInt(n), IsInt(k), Gcd(n,k)=1, [n, k mod (4*n)]) ],
    dims := self >> let(n:=self.params[1], d := [n, 2*(Int((n+1)/2))], When(self.transposed, Reversed(d), d)),
    rdft := PRDFT4,
    diag := IPRDFT_Base.diag1
));


text:
MEMO_Base64Digits :=
  "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz-_";

InstallGlobalFunction(MEMO_Digits,
function(n, base, args...)
  local minlen, len, digits, str;
  # Adapted from the DigitsNumber function in GAPDoc-1.6.1
  if Length(args) = 0 then
    minlen := 0;
  elif Length(args) = 1 then
    minlen := args[1];
  else
    return fail;
  fi;
  digits := MEMO_Base64Digits;
  str := "";
  while n <> 0 do
    Add(str, digits[(n mod base) + 1]);
    n := QuoInt(n, base);
  od;
  if Length(str) < minlen then
    # Pad with zeroes
    Append(str, ListWithIdenticalEntries(minlen - Length(str), digits[1]));
  fi;
  return Reversed(str);
end);

InstallGlobalFunction(MEMO_Hash,
function(key)
  local str, ints, sum, i;
  IO_ClearPickleCache();  # In case IO_Pickle was recently interrupted
  str := IO_Pickle(key);  # Pickle the key to a string
  ints := SHA256String(str);  # Get the SHA-256 checksum in 32-bit chunks
  sum := 0;  # Bring all 256 bits together into a single integer
  for i in [1..Length(ints)] do
    sum := sum + ints[i] * 2 ^ (32 * (i-1));
  od;
  str := MEMO_Digits(sum, 64, 43);  # Make into a padded base-64 string
  return str;
end);

max_stars_repo_name: gap-packages/Memoisation

max_issues_repo_name: gap-packages/Memoisation

max_forks_repo_name: gap-packages/Memoisation


text: # Comment (till end of line)


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

Class(TInterval, Typ, rec(
     __call__ := (self, size) >>
        WithBases(self, rec(
        size := Checked((IsInt(size) and size >= 0) or IsSymbolic(size), size),
        operations := TypOps)),

    hash := (self, val, size) >> (Sum(val, x->x.t.hash(x.v, size)) mod size)+1,

    equals := (self, o) >> IsBound(o.__bases__) and self.__bases__ = o.__bases__
                           and self.size = o.size,

    rChildren := self >> [self.size],
    rSetChild := rSetChildFields("size"),

    print := self >> When(IsBound(self.size), Print(self.name, "(", self.size, ")"), Print(self.name)),
    isIntervalT := true
));

#F IntervalFunc - base class for interval generating functions
Class(IntervalFunc, Function, Sym, rec(
    isSPL := false,
    range := self >> TInterval(self.size),
    domain := self >> TInterval(self.size),
#    tensor_op := (v,rows,cols,fl,gl) ->
#        Lambda(v, fl.at(idiv(v, cols)) * gl.at(imod(v, cols))),
#
    sums := self >> self,
    rSetChild := meth(self, n, newChild)
        self.params[n] := newChild;
    end,

    print := meth(self, i, is)
        local params;
    params := let(p:=self.params, When(IsList(p), p, [p]));
        Print(self.name, "(",
          DoForAllButLast(params, x -> Print(x, ", ")),
          Last(params), ")");
    end,

    def := arg -> rec(),

    toset := self >> Filtered([0..self.size-1], i ->self.elOf(i).ev() = V(1)),

    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,  Inherit(
        rec(params := params),
        ApplyFunc(self.def,params))));

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

#F ivII - interval function
#F
#F ivII(<N>) - diagonal of N 1's
#F ivII(<N>, <to>) - diagonal with 1's from 0..<to>, and 0's after
#F ivII(<N>, <from>, <to>) - diagonal with 1's in <from>..<to>, 0's elsewhere
#F
Class(ivII, IntervalFunc, rec(
    abbrevs := [ (N)         -> [N,0,N],
             (N,to)      -> [N,0,to],
             (N,from,to) -> [N,from,to] ],
    def := (N, from, to) -> Checked(IsPosInt0(N),IsPosInt0(from),IsPosInt0(to),
                                rec(size:=N)),
    elOf := (self, i) >> let(from:=self.params[2], to:=self.params[3], cond(leq(from, i, to-1), 1, 0))

));

Class(ivElOf, FuncClassOper, rec(
    updateDomain  := self >> TVoid,
    updateRange := self >> TInt,
    eval := self >> self.child(1).elOf(self.child(2))
));

Class(ivUnion, FuncClassOper, rec(
    updateDomain := self >> self.child(1).domain(),
    updateRange  := self >> self.child(1).range(),
    elOf := (self, i) >> When(ForAny(self.children(), j->j.elOf(i).ev() = V(1)), V(1), V(0)),
    toset := self >> Filtered([0..self.child(1).size-1], i ->self.elOf(i).ev() = V(1)),
));

Class(ivIntersect, FuncClassOper, rec(
    updateDomain := self >> self.child(1).domain(),
    updateRange  := self >> self.child(1).range(),
    elOf := (self, i) >> When(ForAll(self.children(), j->j.elOf(i).ev() = V(1)), V(1), V(0)),
    toset := self >> Filtered([0..self.child(1).size-1], i ->self.elOf(i).ev() = V(1)),
));

Class(ivInvert, FuncClassOper, rec(
    updateDomain := self >> self.child(1).domain(),
    updateRange  := self >> self.child(1).range(),
    elOf := (self, i) >> When(self.child(1).elOf(i).ev() = V(1), V(0), V(1)),
    toset := self >> Filtered([0..self.child(1).size-1], i ->self.elOf(i).ev() = V(0)),
));

Class(ivTensor, fTensorBase, rec(
    print := FuncClassOper.print,
    updateDomain := self >> TInterval(Product(List(self.children(), i->i.size))),
    updateRange := self >> TInterval(Product(List(self.children(), i->i.size))),
    split_op := self >> let(
        c := self.children(), v := Ind(self.domain().size),
        c1dom := c[1].domain().size, c2dom := c[2].domain().size,
        Checked(Length(c)=2,
            [ Cond(c1dom=1, 0, c2dom=1, v, idiv(v, c2dom)),
            Cond(c1dom=1, v, c2dom=1, 0, imod(v, c2dom)),
            v ])),
    combine_op := (self, split, F, G) >> F.elOf(split[1]) * G.elOf(split[2]),
    elOf := (self, i) >>
        Cond(self.numChildren() > 2, self.rightBinary().elOf(i),
         let(split := self.split_op(),
            Lambda(split[3],
             self.combine_op(split, self.child(1), self.child(2))).at(i))),
    toset := self >> Filtered([0..self.domain().size-1], i ->self.elOf(i).ev() = V(1)),
));


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

Class(RollingPointers, rec(
    __call__ := (self, code, opts) >> CopyFields(self, rec(opts := opts)).apply(code),

    _linearPatterns := Set([[], [mul]]),

    _contains    := (v, U) -> CollectNR(U, @.cond(u->u=v))<>[],
    _linear_exp  := (self, e, v) >> Length(Collect(e, @@.cond((x, cx)-> x=v and (List(cx.parents, ObjId) in self._linearPatterns))))=1,

    _basePtr  := (summands) -> let( p := Filtered(summands, e->IsPtrT(e.t) or IsArrayT(e.t)), Checked(Length(p)=1, p[1]) ),

    # _collectPtrs(<l>) -- <l> loop
    _collectPtrs := meth(self, l)
        local free, vptr, pattern, ptradds, ptrs, e, p, d, id, clusters, loopinv, i;

        # Collecting innermost ponter arithmetic additions.
        # Index computation expression expected to be in normalized form.
        pattern := @@(1, [add], (x, cx) -> IsPtrT(x.t) and (Filtered(cx.parents, IsLoop)=[] or Last(Filtered(cx.parents, IsLoop))=l));

        ptradds := Collect(l, pattern);
        ptradds := Filtered(ptradds, x -> Collect(x.args, pattern)=[]);

        ptrs := tab();

        free := Set(l.free() :: [l.var]);
        # extracting base pointer and common subexpression with this pointer.
        # relying on expression print method to get string id.
        for e in ptradds do
            p  := self._basePtr(e.args);
            # we may have buffer defined in the loop referenced by loop variable
            vptr := Collect(p, var);
            if vptr=[] or Intersection(vptr, free)<>[] then
                id := StringPrint(p);
                if IsBound(ptrs.(id)) then
                    ptrs.(id).cmn := Checked(ptrs.(id).ptr = p, Intersection(ptrs.(id).cmn, Set(e.args)));
                else
                    ptrs.(id) := rec( ptr := p, cmn := Set(e.args), offs := [] );
                fi;
            fi;
        od;
        # make sure common expression doesn't have local loop variables
        for id in UserNSFields(ptrs) do
            ptrs.(id).cmn := Filtered(ptrs.(id).cmn, a -> IsSubset(free, Collect(a, var)));
        od;
        # preparing offsets list for each base expression.
        #
        # offset record: rec(
        #    ref  := <full original expression>,
        #    tail := <offset from common base expression, as summands list>
        # )
        for e in ptradds do
            id := StringPrint(self._basePtr(e.args));
            if IsBound(ptrs.(id)) then
                d := Difference(e.args, ptrs.(id).cmn);
                Add(ptrs.(id).offs, rec( ref := e, tail := d ));
            fi;
        od;

        # when loop variable is not in common subexpression 'cmn' but still in 'tail'
        # try to split pointer so that variable goes to 'cmn' in each resulting pointer
        for id in UserNSFields(ptrs) do
            p := ptrs.(id);
            if not self._contains(l.var, p.cmn) then
                [ clusters, loopinv ] := SplitBy(ConcatList(p.offs, o -> o.tail), s -> self._linear_exp(s, l.var));
                if ForAll(loopinv, s -> not self._contains(l.var, s)) then
                    clusters := Set(clusters);
                    # NOTE: hardcoded max number of resulting rolling pointers
                    if Length(clusters)<=4 then
                        for i in [1..Length(clusters)] do
                            ptrs.(id :: "_" :: StringInt(i)) := rec(
                                ptr  := p.ptr,
                                cmn  := Set( p.cmn :: [clusters[i]] ),
                                # each original <offs> can contain only one expression from <clusters>
                                # because previousely summands were grouped by loop variables
                                offs := List(Filtered(p.offs, o -> clusters[i] in o.tail),
                                          o -> rec( ref := o.ref, tail := RemoveList(o.tail, clusters[i]) )),
                            );
                        od;
                        Unbind(ptrs.(id));
                    fi;
                fi;
            fi;
        od;

        return ptrs;
    end,

    _offsetsInit := function(vars, init, coeff)
        vars.linv := vars.linv :: vars.offs;
        init.linv := init.linv :: List( [1..Length(coeff)], i -> assign(vars.offs[i], vars.strides*coeff[i]));
    end,

    # cost is number of rolling pointers + their stride + number of offsets from rolling pointers
    _offsetsCost := (self, k) >> Cond(Length(k.rp_coeff)>4, 100000, self.rpCost*Length(k.rp_coeff) + 1 + self.coeffCost*Length(k.coeff)),

    # default cost of rolling pointer and coefficient
    rpCost    := 3,
    coeffCost := 1,

    # includeValues
    includeValues := false,

    _strideGrp := s -> CondPat( s, @(1, Value),       V(1),
                                   [mul, Value, @],   s.args[2],
                                # else
                                   s ),
    _strideVal := s -> CondPat( s, @(1, Value),       s.v,
                                   [mul, Value, @],   s.args[1].v,
                                # else
                                   1 ),

    _getStrides := meth(self, l, offs)
        local s, strides, offsets, free;

        if self.includeValues then
            s := ConcatList(offs, e -> e.tail);
        else
            s := ConcatList(offs, e -> Filtered(e.tail, x -> not (x _is Value)));
        fi;
        # ignore expressions with loop-local variables
        free := Set(l.free() :: [l.var]);
        s := Filtered(s, a -> IsSubset(free, Collect(a, var)));

        # get strides vector
        strides := Filtered(GroupList(s, self._strideGrp), e -> Length(e[2])>1);
        if strides<>[] then
            strides := TransposedMat(strides)[1];
        fi;

        # get stride coefficients for each offset + leftover expression
        offsets := TransposedMat(List(offs, function(e)
            local c, leftovers, x, p;
            leftovers := [];
            c := Replicate(Length(strides), 0);
            for x in e.tail do
                p := Position(strides, self._strideGrp(x));
                if p<>false then
                    c[p] := self._strideVal(x);
                else
                    Add(leftovers, x);
                fi;
            od;
            return [c, rec(ref := e.ref, leftovers := leftovers)];
        end));

        return [strides, offsets[1], offsets[2]];
    end,

    # reduce <coeff> offsets list to coefficients relative to <k> "evenly" spaced rolling pointers
    _k_func := function(k, coeff)
        local n, rp, offs, i;
        n     := Length(coeff);
        rp    := List([0..k-1], i -> QuoInt(n*i, k)+1) :: [n+1];
        offs := [];
        for i in [1..k] do
            offs := UnionSet(offs, List(coeff{[rp[i]+1..rp[i+1]-1]}, c -> c - coeff[rp[i]]));
        od;
        return rec(
            rp_coeff := List(DropLast(rp, 1), i -> coeff[i]),
            coeff    := offs,
        );
    end,

    # finds best split [<rolling ptrs>, <offsets>] for <coeff> offsets list.
    # returns record:
    #    .rp_coeff - list of coefficients for rolling pointers
    #    .coeff    - reduced <coeff> list so that every original <coeff>
    #                can be expressed as some <rp_coeff> plus some <coeff>
    #                from this list.
    #
    # Brute force, _offsetsCost can be redefined

    _bestSplit := meth(self, coeff)
        local scoeff, min_k, min_cost, i, j, k, o, cost;
        scoeff := Set(coeff);
        min_k := self._k_func(1, scoeff);
        min_cost := self._offsetsCost(min_k);
        for i in [1..Length(scoeff[1])] do
            o := Sort(ShallowCopy(scoeff), (a,b) -> Cond(a[i]=b[i], a<b, a[i]<b[i]));
            for j in [1..Length(Set(List(o, e -> e[i])))] do
                k := self._k_func(j, o);
                cost := self._offsetsCost(k);
                if cost < min_cost then
                    min_k := k;
                    min_cost := cost;
                fi;
            od;
        od;
        return min_k;
    end,

    # find indices in <offsets> of a given <offs> (not including crossed-out in <refs>)
    _offsetsIndices := (offsets, offs, refs) -> Filtered([1..Length(offsets)], i -> offsets[i]=offs and refs[i].ref<>false),

    # add mapping pair to <map> and cross-out index.
    _mapRecord := function(map, rem, summands)
        Add(map, rec(
            ref := rem.ref,
            subst := ApplyFunc(add, summands :: rem.leftovers)
        ));
        rem.ref := false;
    end,

    # split pointer expressions into rolling pointers + offsets
    _splitPtr := meth(self, ptr, loop, vars, init, incr, map)
        local strides, offsets_coeff, offsets_rem, cfg, i, j, k, ld, li, p;

        [strides, offsets_coeff, offsets_rem] := self._getStrides(loop, ptr.offs);

        cfg := self._bestSplit(offsets_coeff);

        # temporary variables lists.
        # Doing CSE is smarter as we may have same stride and offset computation expressions
        # for other pointers. I don't have this situation in DFTs though.
        vars.rp      := List(cfg.rp_coeff, e -> var.fresh_t("rp",   Cond(IsArrayT(ptr.ptr.t), ptr.ptr.t.toPtrType(), ptr.ptr.t)));
        vars.offs    := List(cfg.coeff,    e -> var.fresh_t("offs", TInt));
        vars.strides := List(strides,      e -> var.fresh_t("s",    TInt));
        vars.inc     := var.fresh_t("inc", TInt);

        # vars.offs is not added here as they might be declared in or outside of the loop
        # this is done by initOffsets method later
        vars.linv    := vars.linv :: vars.rp :: vars.strides :: [vars.inc];

        # preparing mapping from original pointer expressions to rolling pointers expressions
        for i in [1..Length(cfg.rp_coeff)] do
            # for each expression which correspond to this rolling pointer
            for k in self._offsetsIndices(offsets_coeff, cfg.rp_coeff[i], offsets_rem) do
                self._mapRecord(map, offsets_rem[k], [vars.rp[i]]);
            od;

            for j in [1..Length(cfg.coeff)] do
                # for each expression which correspond to this offset from current rolling pointer
                for k in self._offsetsIndices(offsets_coeff, cfg.rp_coeff[i]+cfg.coeff[j], offsets_rem) do
                    self._mapRecord(map, offsets_rem[k], [vars.rp[i], vars.offs[j]]);
                od;
            od;
        od;

        Checked(ForAll(offsets_rem, e -> e.ref=false), "paranoid");

        # <ld> loop var dependent summands, <li> - loop invariant summands
        [ld, li] := SplitBy(ptr.cmn, e -> self._contains(loop.var, e));

        # initialization code
        init.strides := List([1..Length(strides)], i -> assign(vars.strides[i], strides[i]));
        init.rp      := List([1..Length(cfg.rp_coeff)], i -> assign(vars.rp[i], SReduce(ApplyFunc(add, li :: SubstVars( Copy(ld), rec( (loop.var.id) := V(0))) :: Cond(vars.strides<>[], [vars.strides*cfg.rp_coeff[i]], [])), self.opts)));
        init.inc     := assign(vars.inc, SReduce(ApplyFunc(add, SubstVars( Copy(ld), rec( (loop.var.id) := V(1)))), self.opts));

        init.linv    := init.linv :: init.strides :: init.rp :: [init.inc];

        for p in vars.rp do
            Add(incr, assign(p, add(p, vars.inc)));
        od;

        self._offsetsInit(vars, init, cfg.coeff);
    end,

    _processPointer := meth(self, ptr, loop, vars, init, incr, map)
        local ld;
        # process pointers with separable base and loop dependent linear offset
        ld := Filtered(ptr.cmn, e -> self._contains(loop.var, e));
        if ld<>[]
           and ForAll(ld, s -> self._linear_exp(s, loop.var))
           and ForAll(ptr.offs, e -> not self._contains(loop.var, e.tail))
        then
            self._splitPtr(ptr, loop, vars, init, incr, map);
        fi;
    end,

    inject := (self, code, vars, prologue, epilogue) >> Cond(
         code _is decl,  decl(code.vars, self.inject(code.cmd, vars, prologue, epilogue)),
         code _is data,  data(code.var, code.value, self.inject(code.cmd, vars, prologue, epilogue)),
         # else
             decl(vars, chain(prologue, code, epilogue))),

    transformLoop := meth(self, l)
        local new_loop, ptrs, field, vars, init, incr, map, idx, oid, cmd, cp, subst;

        # collect references grouped by base pointers and extracted common summands
        ptrs := self._collectPtrs(l);

        # variables updated by _processPointer method
        #   vars.linv - loop invariant (outer) variables list
        #   vars.loop - inner variables list (declared in the loop body)
        #
        # initialization:
        #   init.linv - outer commands
        #   init.loop - inner commands
        #
        # incr - list of commands that increment pointers (loop body footer)
        # map - list of mapping records:
        #   map[i].ref - expression to replace
        #   map[i].subst
        #
        vars := rec( linv := [], loop := [] );
        init := rec( linv := [], loop := [] );
        incr := [];
        map  := [];

        for field in UserNSFields(ptrs) do
            self._processPointer( ptrs.(field), l, vars, init, incr, map );
        od;

        if map<>[] then
            # copy propagation of simple loop invariant expressions
            cp  := tab();
            init.linv := Filtered(init.linv, function(e) local prop;
                    prop := CondPat( e,
                              [assign, var, Value], true,
                              [assign, var, param], true,
                              false);
                    if prop then
                        cp.(e.loc.id) := e.exp;
                    fi;
                    return not prop;
                end);
            vars.linv := Filtered(vars.linv, e -> not IsBound(cp.(e.id)));

            subst := (c) -> SubstVars(c, cp);

            idx := Set(List(map, e -> e.ref));
            oid := Set(List(map, e -> ObjId(e.ref)));
            cmd := SubstTopDownNR_named(l.cmd, @(1, oid, x -> x in idx),
                    e -> subst(map[PositionProperty( map, x -> x.ref=e )].subst), "toRollingPointer");
            new_loop := decl( vars.linv, chain( subst(init.linv) :: [
                                loop( l.var, l.range, self.inject(cmd, subst(vars.loop), subst(init.loop), subst(incr)))
                        ]));
        else
            new_loop := l;
        fi;

        return new_loop;
    end,

    apply := meth(self, code)

        # normalize linear expressions, group summands and pull out loop variables
        code := SubstTopDownNR_named(code, @(1, self.opts.simpIndicesInside), e ->
                   SubstBottomUp(GroupSummandsExp(ESReduce(e, self.opts)), add, a->__groupSummandsVar(a.args)), "rollingPtrNormalize");

        # transform loops, innermost first
        code := SubstBottomUp(code, @(1, [loop, loopn]), e -> self.transformLoop(e));

        return code;
    end,
));