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

Class(PrunedPRDFT, TaggedNonTerminal, rec(
    abbrevs := [
        (n,blk,pat) -> Checked(IsPosIntSym(n), IsPosIntSym(blk), IsList(pat),
            AnySyms(n,blk) or (IsInt(_unwrap(n)/_unwrap(blk)) and ForAll(pat, i->IsInt(i) and i < n/blk)),
            [_unwrap(n), 1, _unwrap(blk), pat]),
        (n,k,blk,pat) -> Checked(IsPosIntSym(n), IsIntSym(k), IsPosIntSym(blk), IsList(pat),
            AnySyms(n,k) or Gcd(_unwrap(n),_unwrap(k))=1,
            AnySyms(n,blk) or (IsInt(_unwrap(n)/_unwrap(blk)) and ForAll(pat, i->IsInt(i) and i < n/blk)),
            [_unwrap(n), When(AnySyms(n,k), k, k mod _unwrap(n)), _unwrap(blk), pat])
        ],

    dims := self >> let(
        size := let(n:=self.params[1], 2*(Int(n/2)+1)),
        d := [size, self.params[3]*Length(self.params[4])],
        When(self.transposed, [d[2], d[1]], d)
    ),

    terminate := self >> let(
        size := self.params[1], blk := self.params[3], pat := self.params[4],
        res := Compose(PRDFT(size, self.params[2]),
               Tensor(Mat(List(pat, i->BasisVec(size/blk, i))).transpose(), I(blk))).terminate(),
        When(self.transposed, res.transpose(), res)
    ),

    isReal    := True,
    normalizedArithCost := self >> let(n := self.params[1], IntDouble(2.5 * n * d_log(n) / d_log(2))),
    TType := TReal
));

_pruned_PRF12_CT_Children := (N,k,PRFt,DFTt,PRFtp,PRF1) -> Map2(DivisorPairs(N),
    (m,n) -> When(IsEvenInt(n),
	[ PRFt(m,k), DFTt(m,k), PRFtp(m,k) ] :: PRF1(m,n,k),
	[ PRFt(m,k), DFTt(m,k) ] :: PRF1(m,n,k))
);

_pruned_PRF12_CT_Stage1 := (N,k,C,Conj,Tw) ->
    let(m:=Cols(C[1]),
        n:=N/m,
        Nf:=Int(N/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(Nf+1,mf+1,0,n))) * C[1] * Gath(H(2*(nf+1)*m, m, 0, 2*(nf+1))),
            When(nc=1, [],
                ISum(j,
                    RC(Scat(BH(Nf+1,N,m,j+1,n))) *
                    Conj * RC(C[2]) * Tw(j) *
                    RC(Gath(H((nf+1)*m, m, j+1, nf+1)))
                )
            ),
            When(IsOddInt(n), [],
                RC(Scat(H(Nf+1,mc,nf,n))) * C[3] * Gath(H(2*(nf+1)*m, m, 2*nf, 2*(nf+1)))
            )
        )
    );

_pruned_IPRF12_CT_Stage2 := (N,k,C,Conj,Tw) ->
    let(m:=Rows(C[1]),
        n:=N/m,
        Nf:=Int(N/2),
        nf:=Int(n/2),
        nc:=Int((n+1)/2),
        mf:=Int(m/2),
        mc:=Int((m+1)/2),
        j:=Ind(nc-1),
        SUM(
            Scat(H(2*(nf+1)*m, m, 0, 2*(nf+1))) * C[1] * RC(Gath(H(Nf+1,mf+1,0,n))),
        	When(nc=1, [],
                ISum(j,
                    RC(Scat(H((nf+1)*m, m, j+1, nf+1))) *
                    Tw(j) * RC(C[2]) * Conj *
    	            RC(Gath(BH(Nf+1,N,m,j+1,n)))
                )
            ),
        	When(IsOddInt(n), [],
    	       Scat(H(2*(nf+1)*m, m, 2*nf, 2*(nf+1))) * C[3] * RC(Gath(H(Nf+1,mc,nf,n))))
        )
    );

NewRulesFor(PrunedPRDFT, rec(
    PrunedPRDFT_base := rec(
       forTransposition := false,
       maxSize := 256,
       applicable := (self, nt) >> not nt.hasTags() and nt.params[1] <= self.maxSize and nt.params[3] = 1,
       children := nt -> [[PRDFT(nt.params[1], nt.params[2])]],
       apply := (nt, C, cnt) -> let(size := nt.params[1], blk := nt.params[3], pat := nt.params[4],
            C[1] * Tensor(Mat(List(pat, i->BasisVec(size/blk, i))).transpose(), I(blk)).terminate())
    ),
    PrunedPRDFT_DFT := rec(
       forTransposition := false,
       applicable := (self, nt) >> nt.params[1] = nt.params[3] and nt.params[4] = [0],
       children := nt -> [[ PRDFT(nt.params[1], nt.params[2]).withTags(nt.getTags()) ]],
       apply := (nt, C, cnt) -> C[1]
    ),
    PrunedPRDFT_CT_rec_block := rec(
       forTransposition := true,
       fPrecompute := f->fPrecompute(f),
       applicable := (self, nt) >> nt.params[1] > 2
            and not nt.hasTags()
            and not IsPrime(nt.params[1])
            and nt.params[3] = 1
            and Last(nt.params[4])+1-nt.params[4][1] = Length(Set(nt.params[4]))
            and ForAny(Map2(DivisorPairs(nt.params[1]), (m,n)->_ctpr_applicable(m, n, nt.params[4])), i->i),
       children := (nt) -> Filtered(_pruned_PRF12_CT_Children(nt.params[1],nt.params[2],PRDFT1, DFT1, PRDFT3,
                                        (m, n, k)->_pruned_children(m, n, nt.params[4], (r, sp) -> PrunedPRDFT(r, k mod n, nt.params[3], sp))),
                              lst -> ForAll(Filtered(lst, i->ObjId(i) = PrunedPRDFT), j->j.params[4] <> [])),
       apply := (self, nt, C, cnt) >> let(N:=nt.params[1], k:=nt.params[2], m:=Cols(C[1]), n := N/m,
            _pruned_PRF12_CT_Stage1(N, k, C, Diag(BHD(m,1,-1)), j->RC(Diag(self.fPrecompute(Twid(N,m,k,0,0,j+1))))) *
            _build_stack(m, n, nt.params[4], let(s1len := When(IsEvenInt(n), 3, 2), C{[s1len+1..Length(C)]}), nt.dims()[2], Gath, IterVStack, VStack, (a,b)->a*b))
    )
));

Class(PrunedIPRDFT, TaggedNonTerminal, rec(
    abbrevs := [
        (n,blk,pat) -> Checked(IsPosIntSym(n), IsPosIntSym(blk), IsList(pat),
            AnySyms(n,blk) or (IsInt(_unwrap(n)/_unwrap(blk)) and ForAll(pat, i->IsInt(i) and i < n/blk)),
            [_unwrap(n), 1, _unwrap(blk), pat]),
        (n,k,blk,pat) -> Checked(IsPosIntSym(n), IsIntSym(k), IsPosIntSym(blk), IsList(pat),
            AnySyms(n,k) or Gcd(_unwrap(n),_unwrap(k))=1,
            AnySyms(n,blk) or (IsInt(_unwrap(n)/_unwrap(blk)) and ForAll(pat, i->IsInt(i) and i < n/blk)),
            [_unwrap(n), When(AnySyms(n,k), k, k mod _unwrap(n)), _unwrap(blk), pat])
        ],

    dims := self >> let(
        size := let(n:=self.params[1], 2*(Int(n/2)+1)),
        d := [self.params[3]*Length(self.params[4]), size],
        When(self.transposed, [d[2], d[1]], d)
    ),

    terminate := self >> let(
        size := self.params[1], blk := self.params[3], pat := self.params[4],
        res := Compose(Tensor(Mat(List(pat, i->BasisVec(size/blk, i))), I(blk)),
                       IPRDFT(size, self.params[2])).terminate(),
        When(self.transposed, res.transpose(), res)
    ),

    isReal    := True,
    normalizedArithCost := self >> let(n := self.params[1], IntDouble(2.5 * n * d_log(n) / d_log(2))),
    TType := TReal
));

NewRulesFor(PrunedIPRDFT, rec(
    PrunedIPRDFT_base := rec(
       forTransposition := false,
       maxSize := 256,
       applicable := (self, nt) >> not nt.hasTags() and nt.params[1] <= self.maxSize and nt.params[3] = 1,
       children := nt -> [[IPRDFT(nt.params[1], nt.params[2])]],
       apply := (nt, C, cnt) -> let(size := nt.params[1], blk := nt.params[3], pat := nt.params[4],
            Tensor(Mat(List(pat, i->BasisVec(size/blk, i))), I(blk)).terminate() * C[1] )
    ),
    PrunedIPRDFT_DFT := rec(
       forTransposition := false,
       applicable := (self, nt) >> nt.params[1] = nt.params[3] and nt.params[4] = [0],
       children := nt -> [[ IPRDFT(nt.params[1], nt.params[2]).withTags(nt.getTags()) ]],
       apply := (nt, C, cnt) -> C[1]
    ),
    PrunedIPRDFT_CT_rec_block := rec(
       forTransposition := true,
       applicable := (self, nt) >> nt.params[1] > 2
            and not nt.hasTags()
            and not IsPrime(nt.params[1])
            and nt.params[3] = 1
            and Last(nt.params[4])+1-nt.params[4][1] = Length(Set(nt.params[4]))
            and ForAny(Map2(DivisorPairs(nt.params[1]), (m,n)->_ctpr_applicable(m, n, nt.params[4])), i->i),
       children := (nt) -> Filtered(_pruned_PRF12_CT_Children(nt.params[1],nt.params[2], IPRDFT1, DFT1, IPRDFT2,
                                        (m, n, k)->_pruned_children(m, n, nt.params[4], (r, sp) -> PrunedIPRDFT(r, k mod n, nt.params[3], sp))),
                              lst -> ForAll(Filtered(lst, i->ObjId(i) = PrunedIPRDFT), j->j.params[4] <> [])),
       apply := (nt, C, cnt) -> let(N:=nt.params[1], k:=nt.params[2], m:=Rows(C[1]), n := N/m,
            _build_stack(m, n, nt.params[4], let(s1len := When(IsEvenInt(n), 3, 2), C{[s1len+1..Length(C)]}), nt.dims()[1], Scat, IterHStack1, HStack1, (a,b)->b*a) *
            _pruned_IPRF12_CT_Stage2(N, k, C, Diag(BHD(m,1,-1)), j->RC(Diag(fPrecompute(Twid(N,m,k,0,0,j+1)))))
       )
    ),
));

Class(IOPrunedConv, TaggedNonTerminal, rec(
    abbrevs := [
        (n,h,oblk,opat,iblk,ipat) -> Checked(IsPosIntSym(n), IsPosIntSym(oblk), IsList(opat), IsPosIntSym(iblk), IsList(ipat),
            AnySyms(n,iblk,oblk) or (IsInt(_unwrap(n)/_unwrap(iblk)) and IsInt(_unwrap(n)/_unwrap(oblk)) and
                ForAll(ipat, i->IsInt(i) and i < n/iblk) and ForAll(opat, i->IsInt(i) and i < n/oblk)),
            [_unwrap(n), h, _unwrap(oblk), opat, _unwrap(iblk), ipat, false]),
        (n,h,oblk,opat,iblk,ipat,isFreqData) -> Checked(IsPosIntSym(n), IsPosIntSym(oblk), IsList(opat), IsPosIntSym(iblk), IsList(ipat),
            AnySyms(n,iblk,oblk) or (IsInt(_unwrap(n)/_unwrap(iblk)) and IsInt(_unwrap(n)/_unwrap(oblk)) and
                ForAll(ipat, i->IsInt(i) and i < n/iblk) and ForAll(opat, i->IsInt(i) and i < n/oblk)),
            [_unwrap(n), h, _unwrap(oblk), opat, _unwrap(iblk), ipat, isFreqData])
        ],
    dims      := self >> [self.params[3]*Length(self.params[4]), self.params[5]*Length(self.params[6])],
    terminate := self >> let(size := self.params[1],
        taps := When(self.params[7], MatSPL(DFT(self.params[1], -1)) *
               Cond(IsList(self.params[2]), self.params[2], IsLambda(self.params[2]), self.params[2].tolist(), IsBound(self.params[2].list), self.params[2].list),
           self.params[2]),
        oblk := self.params[3], opat := self.params[4], iblk := self.params[5], ipat := self.params[6],
        Tensor(Mat(List(opat, i->BasisVec(size/oblk, i))), I(oblk)).terminate() *
        spiral.transforms.filtering.Circulant(size, taps, -size).terminate() *
        Tensor(Mat(List(ipat, i->BasisVec(size/iblk, i))).transpose(), I(iblk)).terminate()),
    isReal    := self >> false,
    normalizedArithCost := self >> let(n := self.params[1], IntDouble(5 * n * d_log(n) / d_log(2))),
    TType := TComplex,
    hashAs := self >> ApplyFunc(ObjId(self), [self.params[1], fUnk(self.params[2].range(), self.params[2].domain())]::Drop(self.params, 2))
));

NewRulesFor(IOPrunedConv, rec(
    IOPrunedConv_PrunedDFT_IPrunedDFT := rec(
       forTransposition := false,
       applicable :=  (self, nt) >> not nt.hasTags(),
       children := nt -> [[ PrunedIDFT(nt.params[1], 1, nt.params[3], nt.params[4]),
                            PrunedDFT(nt.params[1], -1, nt.params[5], nt.params[6]) ]],
       apply := (nt, C, cnt) -> Cond(ObjId(nt.params[2]) = var, # a variable is frequency domain data -> FData
                                    Checked(nt.params[7], C[1] * Diag(FData(nt.params[2])) * C[2]), # if data is provided it needs to be frequency domain
                                    ObjId(nt.params[2]) = Lambda,  # data is inside a larger array
                                    Checked(nt.params[7], C[1] * Diag(nt.params[2]) * C[2]),
                                    nt.params[7],   # if params[7] is true then data is frequency domain data
                                    C[1] * Diag(When(IsList(nt.params[2]), FList(TComplex, nt.params[2]), nt.params[2])) * C[2],
                                    let(cxfftdiag := 1/nt.params[1]*ComplexFFT(List(nt.params[2].tolist(), i->ComplexAny(_unwrap(i)))),
                                        C[1] * Diag(FData(cxfftdiag)) * C[2]))
    ),
));

#===========================================

Class(IOPrunedRConv, TaggedNonTerminal, rec(
    abbrevs := [
        (n,h,oblk,opat,iblk,ipat) -> Checked(IsPosIntSym(n), IsPosIntSym(oblk), IsList(opat), IsPosIntSym(iblk), IsList(ipat),
            AnySyms(n,iblk,oblk) or (IsInt(_unwrap(n)/_unwrap(iblk)) and IsInt(_unwrap(n)/_unwrap(oblk)) and
                ForAll(ipat, i->IsInt(i) and i < n/iblk) and ForAll(opat, i->IsInt(i) and i < n/oblk)),
            [_unwrap(n), h, _unwrap(oblk), opat, _unwrap(iblk), ipat, false]),
        (n,h,oblk,opat,iblk,ipat,isFreqData) -> Checked(IsPosIntSym(n), IsPosIntSym(oblk), IsList(opat), IsPosIntSym(iblk), IsList(ipat),
            AnySyms(n,iblk,oblk) or (IsInt(_unwrap(n)/_unwrap(iblk)) and IsInt(_unwrap(n)/_unwrap(oblk)) and
                ForAll(ipat, i->IsInt(i) and i < n/iblk) and ForAll(opat, i->IsInt(i) and i < n/oblk)),
            [_unwrap(n), h, _unwrap(oblk), opat, _unwrap(iblk), ipat, isFreqData])
        ],
    dims      := self >> [self.params[3]*Length(self.params[4]), self.params[5]*Length(self.params[6])],
    terminate := self >> let(size := self.params[1],
        taps := When(self.params[7], MatSPL(IPRDFT(self.params[1], -1)) * When(IsList(self.params[2]), self.params[2], self.params[2].list), self.params[2]),
        oblk := self.params[3], opat := self.params[4], iblk := self.params[5], ipat := self.params[6],
        Tensor(Mat(List(opat, i->BasisVec(size/oblk, i))), I(oblk)).terminate() *
        spiral.transforms.filtering.Circulant(size, taps, -size).terminate() *
        Tensor(Mat(List(ipat, i->BasisVec(size/iblk, i))).transpose(), I(iblk)).terminate()),
    isReal    := self >> true,
    normalizedArithCost := self >> let(n := self.params[1], IntDouble(5 * n * d_log(n) / d_log(2))),
    TType := TReal
));

NewRulesFor(IOPrunedRConv, rec(
    IOPrunedRConv_PrunedIPRDFT_PrunedPRDFT := rec(
       forTransposition := false,
       applicable :=  (self, nt) >> not nt.hasTags(),
       children := nt -> [[ PrunedIPRDFT(nt.params[1], 1, nt.params[3], nt.params[4]),
                            PrunedPRDFT(nt.params[1], -1, nt.params[5], nt.params[6]) ]],
       apply := (nt, C, cnt) -> Cond(ObjId(nt.params[2]) = var, # a variable is frequency domain data -> FData
                                Checked(nt.params[7], C[1] * RCDiag(FData(nt.params[2])) * C[2]), # if data is provided it needs to be frequency domain
                                nt.params[7],   # if params[7] is true then data is frequency domain data
                                C[1] *  RCDiag(When(IsList(nt.params[2]), FList(TReal, nt.params[2]), nt.params[2])) * C[2],
                                let(cxfftdiag := 1/nt.params[1]*ComplexFFT(List(nt.params[2].tolist(), i->ComplexAny(_unwrap(i)))),
                                    rcdiag := cxfftdiag{[1..Rows(PRDFT(nt.params[1], -1))/2]},
                                    rlist := Flat(List(rcdiag, i->[ReComplex(_unwrap(i)), ImComplex(_unwrap(i))])),
                                    C[1] * RCDiag(FList(TReal, rlist)) * C[2]))
    ),
));

_complexify := l -> List([1..Length(l)/2], i->l[2*i-1] + E(4) * l[2*i]);

#=========================================================================================================================================================
Declare(IOPrunedMDRConv);
Class(IOPrunedMDRConv, TaggedNonTerminal, rec(
    abbrevs := [
        (n,h,oblk,opat,iblk,ipat) -> Checked(ForAll(n, IsPosIntSym), IsPosIntSym(oblk), ForAll(opat, IsList), IsPosIntSym(iblk), ForAll(ipat, IsList),
            [_unwrap(n), h, _unwrap(oblk), opat, _unwrap(iblk), ipat, false]),
        (n,h,oblk,opat,iblk,ipat,isFreqData) -> Checked(ForAll(n, IsPosIntSym), IsPosIntSym(oblk), ForAll(opat, IsList), IsPosIntSym(iblk), ForAll(ipat, IsList),
            [_unwrap(n), h, _unwrap(oblk), opat, _unwrap(iblk), ipat, isFreqData])
        ],
    dims      := self >> [self.params[3]*Product(self.params[4], i->Length(i)), self.params[5]*Product(self.params[6], i->Length(i))],
    isReal    := self >> true,
    normalizedArithCost := self >> let(n := self.params[1], IntDouble(5 * n * d_log(n) / d_log(2))),
    TType := TReal,
    terminate := self >> When(self.params[7],
                            let(nlist := self.params[1],
                                n := nlist[1],
                                nfreq := n/2+1,
                                nrem := Drop(nlist, 1),
                                idft := Tensor(IPRDFT(n, -1), I(Product(nrem))) *
                                    Tensor(I(nfreq), L(2*Product(nrem), 2)) * Tensor(I(nfreq), RC(MDDFT(nrem, -1))),
                                tlist := MatSPL(idft) *
                                    Cond(IsList(self.params[2]), self.params[2], IsLambda(self.params[2]), self.params[2].tolist(), IsBound(self.params[2].list), self.params[2].list),
                                IOPrunedMDRConv(self.params[1], FList(TReal, tlist), self.params[3], self.params[4],
                                    self.params[5], self.params[6], false).terminate()
                            ),
                            let(nlist := self.params[1],
                                oblk := self.params[3], opat_md := self.params[4],
                                iblk := self.params[5], ipat_md := self.params[6],
                                scat3d := Tensor(List(Zip2(ipat_md, nlist), j->let(ipat := j[1], size := j[2],
                                    Tensor(Mat(List(ipat, i->BasisVec(size/iblk, i))).transpose(), I(iblk)).terminate()))::[Mat([[1], [0]])]),
                                gath3d := Tensor(List(Zip2(opat_md, nlist), j->let(opat := j[1], size := j[2],
                                    Tensor(Mat(List(opat, i->BasisVec(size/oblk, i))), I(oblk)).terminate()))::[Mat([[1, 0]])]),

                                dft3dr := RC(MDDFT(nlist, -1)),
                                idft3dr := RC(MDDFT(nlist, 1)),
                                gfd := List(1/Product(nlist) * MatSPL(MDDFT(nlist, 1)) * self.params[2].list, i->ComplexAny(_unwrap(i))),
                                gdiagr := RC(Diag(gfd)),
                                t := gath3d * idft3dr * gdiagr * dft3dr * scat3d,
                                t.terminate()
                            )
                        ),
    hashAs := self >> ApplyFunc(ObjId(self), [self.params[1], fUnk(self.params[2].range(), self.params[2].domain())]::Drop(self.params, 2))
));

NewRulesFor(IOPrunedMDRConv, rec(
    IOPrunedMDRConv_time2freq := rec(
       forTransposition := false,
       applicable :=  (self, nt) >> not nt.hasTags() and not nt.params[7],
       children := nt -> let(htime := nt.params[2].list,
                            nlist := nt.params[1],
                            dft := 1/Product(nlist) * MatSPL(MDDFT(nlist, 1)),
                            nfreq := nlist[1]/2+1,
                            hfreq := dft*htime,
                            hfreqc := List(hfreq, ComplexAny),
                            hfreqce := hfreqc{[1..nfreq*Product(Drop(nlist,1))]},
                            hfreqr := Flat(List(hfreqce, i->[ReComplex(i), ImComplex(i)])),
                            [[IOPrunedMDRConv(nt.params[1], FList(TReal, hfreqr), nt.params[3], nt.params[4], nt.params[5], nt.params[6], true)]]
                        ),
       apply := (nt, C, cnt) -> C[1]
    ),

    ## GPU/TITAN V Hockney algotithm variant
    ## 2-trip, 5-step, ZYX ====================================================
    IOPrunedMDRConv_3D_2trip_zyx_freqdata := rec(
       forTransposition := false,
       applicable :=  (self, nt) >> not nt.hasTags() and Length(nt.params[1]) = 3 and nt.params[7],
       children := nt -> let(nlist := nt.params[1],
                            diag := nt.params[2],
                            oblk := nt.params[3],
                            opats := nt.params[4],
                            iblk := nt.params[5],
                            ipats := nt.params[6],
                            nfreq := nlist[1]/2+1,
                            i := Ind(nfreq*nlist[2]),
                            hfunc := Cond(ObjId(diag) = Lambda,
                                let(j := Ind(nlist[3]),
                                    # Lambda(j, cxpack(diag.at(2*(j + i*nlist[3])), diag.at(2*(j + i*nlist[3])+1)))
                                    pos := i +j*nfreq*nlist[2],
                                    Lambda(j, cxpack(diag.at(2*pos), diag.at(2*pos+1)))
                                ),
                                ObjId(diag) = fUnk,
                                fUnk(TComplex, nlist[3]),
                                let(list := nt.params[2].list,  # here we assume FList(TReal, [...])
                                    clist := List([1..Length(list)/2], i->Cplx(list[2*i-1], list[2*i])),
                                    fc := FList(TComplex, clist),
                                    gf := fTensor(fBase(i), fId(nlist[3])),
                                    fCompose(fc, gf)
                                )
                            ),
                            [[ PrunedPRDFT(nlist[1], -1, iblk, ipats[1]),  # stage 1: PRDFT z
                                PrunedDFT(nlist[2], -1, iblk, ipats[2]),    # stage 2: DFT y
                                IOPrunedConv(nlist[3], hfunc, oblk, opats[3], iblk, ipats[3], true), # stage 3+4+5: complex conv in x
                                PrunedIDFT(nlist[2], 1, oblk, opats[2]), # stage 6: iDFT in y
                                PrunedIPRDFT(nlist[1], 1, oblk, opats[1]),   # stage 7: iPRDFT in z
                                InfoNt(i)
                            ]]),

       apply := (nt, C, cnt) -> let(prdft1d := C[1],
                                    pdft1d := C[2],
                                    iopconv := C[3],
                                    ipdft1d := C[4],
                                    iprdft1d := C[5],
                                    i := cnt[6].params[1],
                                    nlist := nt.params[1],
                                    n1 := nlist[1],
                                    nfreq := nlist[1]/2+1,
                                    n2 := nlist[2],
                                    n3 := nlist[3],
                                    oblk := nt.params[3],
                                    opats := nt.params[4],
                                    iblk := nt.params[5],
                                    ipats := nt.params[6],
                                    ns1 := iblk * Length(ipats[1]),
                                    ns2 := iblk * Length(ipats[2]),
                                    ns3 := iblk * Length(ipats[3]),
                                    nd1 := oblk * Length(opats[1]),
                                    nd2 := oblk * Length(opats[2]),
                                    nd3 := oblk * Length(opats[3]),
                                    stage1 := L(2*nfreq*ns3*ns2, ns3) * Tensor(I(ns2), Tensor(L(2*nfreq, 2) * prdft1d, I(ns3))) * Tensor(L(ns2*ns1, ns2), I(ns3)),
                                    stage2 := Tensor(I(ns3), Tensor(RC(pdft1d), I(nfreq))),
                                    pp := Tensor(L(ns3*n2*nfreq, n2*nfreq), I(2)) * Tensor(I(ns3), L(2*nfreq*n2, nfreq)),
                                    ppi := Tensor(I(nd3), L(2*nfreq*n2, 2*n2)) * Tensor(L(nd3*n2*nfreq, nd3), I(2)),
                                    stage543 := ppi * IDirSum(i, RC(iopconv)) * pp,
                                    stage76 := Tensor(L(nd2*nd1, nd1), I(nd3)) * Grp(Tensor((Tensor(I(nd2), iprdft1d * L(2*nfreq, nfreq)) *
                                        Tensor(RC(ipdft1d), I(nfreq))), I(nd3)) * L(2*nfreq*nd3*n2, 2*nfreq*n2)),
                                    conv3dr := stage76 * stage543 * stage2 * stage1,
                                    conv3dr
                            ),
    ),

    ## X dimension last (convolution)
    ## 5-step, ZYX ====================================================
    IOPrunedMDRConv_3D_5step_zyx_freqdata := rec(
       forTransposition := false,
       applicable :=  (self, nt) >> not nt.hasTags() and Length(nt.params[1]) = 3 and nt.params[7],
       children := nt -> let(nlist := nt.params[1],
                            diag := nt.params[2],
                            oblk := nt.params[3],
                            opats := nt.params[4],
                            iblk := nt.params[5],
                            ipats := nt.params[6],
                            nfreq := nlist[1]/2+1,
                            i := Ind(nfreq*nlist[2]),
                            hfunc := Cond(ObjId(diag) = Lambda,
                                let(j := Ind(nlist[3]),
                                    Lambda(j, cxpack(diag.at(2*(j + i*nlist[3])), diag.at(2*(j + i*nlist[3])+1)))
                                ),
                                ObjId(diag) = fUnk,
                                fUnk(TComplex, nlist[3]),
                                let(list := nt.params[2].list,  # here we assume FList(TReal, [...])
                                    clist := List([1..Length(list)/2], i->Cplx(list[2*i-1], list[2*i])),
                                    fc := FList(TComplex, clist),
                                    gf := fTensor(fBase(i), fId(nlist[3])),
                                    fCompose(fc, gf)
                                )
                            ),
                            [[ PrunedPRDFT(nlist[1], -1, iblk, ipats[1]),  # stage 1: PRDFT z
                                PrunedDFT(nlist[2], -1, iblk, ipats[2]),    # stage 2: DFT y
                                IOPrunedConv(nlist[3], hfunc, oblk, opats[3], iblk, ipats[3], true), # stage 3+4+5: complex conv in x
                                PrunedIDFT(nlist[2], 1, oblk, opats[2]), # stage 6: iDFT in y
                                PrunedIPRDFT(nlist[1], 1, oblk, opats[1]),   # stage 7: iPRDFT in z
                                InfoNt(i)
                            ]]),

       apply := (nt, C, cnt) -> let(prdft1d := Grp(C[1]),
                                    pdft1d := Grp(C[2]),
                                    iopconv := Grp(C[3]),
                                    ipdft1d := Grp(C[4]),
                                    iprdft1d := Grp(C[5]),
                                    i := cnt[6].params[1],
                                    nlist := nt.params[1],
                                    n1 := nlist[1],
                                    nfreq := nlist[1]/2+1,
                                    n2 := nlist[2],
                                    n3 := nlist[3],
                                    oblk := nt.params[3],
                                    opats := nt.params[4],
                                    iblk := nt.params[5],
                                    ipats := nt.params[6],
                                    ns1 := iblk * Length(ipats[1]),
                                    ns2 := iblk * Length(ipats[2]),
                                    ns3 := iblk * Length(ipats[3]),
                                    nd1 := oblk * Length(opats[1]),
                                    nd2 := oblk * Length(opats[2]),
                                    nd3 := oblk * Length(opats[3]),
                                    stage1 := Tensor(I(nfreq), L(2*ns2*ns3, ns2*ns3)) * Tensor(prdft1d, I(ns2*ns3)),
                                    stage2 := RC(Tensor(I(nfreq), pdft1d, I(ns3))),
                                    stage543c := IDirSum(i, iopconv),
                                    stage543 := RC(stage543c),
                                    stage6 := RC(Tensor(I(nfreq), ipdft1d, I(nd3))),
                                    stage7 := Tensor(iprdft1d, I(nd2*nd3)) * Tensor(I(nfreq), L(2*nd2*nd3, 2)),
                                    conv3dr := stage7 * stage6 * stage543 * stage2 * stage1,
                                    conv3dr
                            ),
    ),
    ## 5-step, YZX ====================================================
    IOPrunedMDRConv_3D_5step_yzx_freqdata := rec(
       forTransposition := false,
       applicable :=  (self, nt) >> false, #not nt.hasTags() and Length(nt.params[1]) = 3 and nt.params[7], Breaks!
       children := nt -> let(nlist := nt.params[1],
                            diag := nt.params[2],
                            oblk := nt.params[3],
                            opats := nt.params[4],
                            iblk := nt.params[5],
                            ipats := nt.params[6],
                            nfreq := nlist[1]/2+1,
                            i := Ind(nfreq*nlist[2]),
                            hfunc := Cond(ObjId(diag) = Lambda,
                                let(j := Ind(nlist[3]),
                                    Lambda(j, cxpack(diag.at(2*(j + i*nlist[3])), diag.at(2*(j + i*nlist[3])+1)))
                                ),
                                ObjId(diag) = fUnk,
                                fUnk(TComplex, nlist[3]),
                                let(list := nt.params[2].list,  # here we assume FList(TReal, [...])
                                    clist := List([1..Length(list)/2], i->Cplx(list[2*i-1], list[2*i])),
                                    fc := FList(TComplex, clist),
                                    gf := fTensor(fBase(i), fId(nlist[3])),
                                    fCompose(fc, gf)
                                )
                            ),
                            [[ PrunedPRDFT(nlist[2], -1, iblk, ipats[2]),  # stage 1: PRDFT y
                                PrunedDFT(nlist[1], -1, iblk, ipats[1]),    # stage 2: DFT z
                                IOPrunedConv(nlist[3], hfunc, oblk, opats[3], iblk, ipats[3], true), # stage 3+4+5: complex conv in x
                                PrunedIDFT(nlist[1], 1, oblk, opats[1]), # stage 6: iDFT in z
                                PrunedIPRDFT(nlist[2], 1, oblk, opats[2]),   # stage 7: iPRDFT in y
                                InfoNt(i)
                            ]]),

       apply := (nt, C, cnt) -> let(prdft1d := C[1],
                                    pdft1d := C[2],
                                    iopconv := C[3],
                                    ipdft1d := C[4],
                                    iprdft1d := C[5],
                                    i := cnt[6].params[1],
                                    nlist := nt.params[1],
                                    n1 := nlist[1],
                                    nfreq := nlist[2]/2+1,
                                    n2 := nlist[2],
                                    n3 := nlist[3],
                                    oblk := nt.params[3],
                                    opats := nt.params[4],
                                    iblk := nt.params[5],
                                    ipats := nt.params[6],
                                    ns1 := iblk * Length(ipats[1]),
                                    ns2 := iblk * Length(ipats[2]),
                                    ns3 := iblk * Length(ipats[3]),
                                    nd1 := oblk * Length(opats[1]),
                                    nd2 := oblk * Length(opats[2]),
                                    nd3 := oblk * Length(opats[3]),
                                    stage1 := Tensor(I(nfreq*ns2), L(2*ns3, ns3)) * Tensor(I(ns2), prdft1d, I(ns3)),
                                    stage2 := RC(Tensor(pdft1d, I(nfreq*ns3))),
                                    stage543c := IDirSum(i, iopconv),
                                    stage543 := RC(stage543c),
                                    stage6 := RC(Tensor(ipdft1d, I(nfreq*nd3))),
                                    stage7 := Tensor(I(nd2), iprdft1d, I(nd3)) * Tensor(I(nfreq*nd2), L(2*nd3, 2)),
                                    conv3dr := stage7 * stage6 * stage543 * stage2 * stage1,
                                    conv3dr
                            ),
    ),
    ## 3-step, (ZY)X ====================================================
    IOPrunedMDRConv_3D_3step_zyx_freqdata := rec(
       forTransposition := false,
       applicable :=  (self, nt) >> not nt.hasTags() and Length(nt.params[1]) = 3 and nt.params[7],
       children := nt -> let(nlist := nt.params[1],
                            diag := nt.params[2],
                            oblk := nt.params[3],
                            opats := nt.params[4],
                            iblk := nt.params[5],
                            ipats := nt.params[6],
                            nfreq := nlist[1]/2+1,
                            i := Ind(nfreq*nlist[2]),
                            hfunc := Cond(ObjId(diag) = Lambda,
                                let(j := Ind(nlist[3]),
                                    Lambda(j, cxpack(diag.at(2*(j + i*nlist[3])), diag.at(2*(j + i*nlist[3])+1)))
                                ),
                                ObjId(diag) = fUnk,
                                fUnk(TComplex, nlist[3]),
                                let(list := nt.params[2].list,  # here we assume FList(TReal, [...])
                                    clist := List([1..Length(list)/2], i->Cplx(list[2*i-1], list[2*i])),
                                    fc := FList(TComplex, clist),
                                    gf := fTensor(fBase(i), fId(nlist[3])),
                                    fCompose(fc, gf)
                                )
                            ),
                            [[ PrunedPRDFT(nlist[1], -1, iblk, ipats[1]),  # stage 1a: PRDFT z
                                PrunedDFT(nlist[2], -1, iblk, ipats[2]),    # stage 1b: DFT y
                                IOPrunedConv(nlist[3], hfunc, oblk, opats[3], iblk, ipats[3], true), # stage 3+4+5: complex conv in x
                                PrunedIDFT(nlist[2], 1, oblk, opats[2]), # stage 6b: iDFT in y
                                PrunedIPRDFT(nlist[1], 1, oblk, opats[1]),   # stage 6a: iPRDFT in z
                                InfoNt(i)
                            ]]),

       apply := (nt, C, cnt) -> let(prdft1d := Grp(C[1]),
                                    pdft1d := Grp(C[2]),
                                    iopconv := Grp(C[3]),
                                    ipdft1d := Grp(C[4]),
                                    iprdft1d := Grp(C[5]),
                                    i := cnt[6].params[1],
                                    nlist := nt.params[1],
                                    n1 := nlist[1],
                                    nfreq := nlist[1]/2+1,
                                    n2 := nlist[2],
                                    n3 := nlist[3],
                                    oblk := nt.params[3],
                                    opats := nt.params[4],
                                    iblk := nt.params[5],
                                    ipats := nt.params[6],
                                    ns1 := iblk * Length(ipats[1]),
                                    ns2 := iblk * Length(ipats[2]),
                                    ns3 := iblk * Length(ipats[3]),
                                    nd1 := oblk * Length(opats[1]),
                                    nd2 := oblk * Length(opats[2]),
                                    nd3 := oblk * Length(opats[3]),
                                    stage1 := Tensor(I(nfreq), L(2*ns2, ns2)) * Tensor(prdft1d, I(ns2)),
                                    stage2 := RC(Tensor(I(nfreq), pdft1d)),
                                    stage12 := Tensor(I(nfreq * n2), L(2*ns3, ns3)) * Tensor(stage2 * stage1, I(ns3)),
                                    stage543c := IDirSum(i, iopconv),
                                    stage543 := RC(stage543c),
                                    stage6 := RC(Tensor(I(nfreq), ipdft1d)),
                                    stage7 := Tensor(iprdft1d, I(nd2)) * Tensor(I(nfreq), L(2*nd2, 2)),
                                    stage67 := Tensor(stage7 * stage6, I(nd3)) * Tensor(I(nfreq * n2), L(2*nd3, 2)),
                                    conv3dr := stage67 * stage543 * stage12,
                                    conv3dr
                            ),
    ),
    ## 3-step, (YZ)X ====================================================
    IOPrunedMDRConv_3D_3step_yzx_freqdata := rec(
       forTransposition := false,
       applicable :=  (self, nt) >> false, #not nt.hasTags() and Length(nt.params[1]) = 3 and nt.params[7], Breaks!
       children := nt -> let(nlist := nt.params[1],
                            diag := nt.params[2],
                            oblk := nt.params[3],
                            opats := nt.params[4],
                            iblk := nt.params[5],
                            ipats := nt.params[6],
                            nfreq := nlist[1]/2+1,
                            i := Ind(nfreq*nlist[2]),
                            hfunc := Cond(ObjId(diag) = Lambda,
                                let(j := Ind(nlist[3]),
                                    Lambda(j, cxpack(diag.at(2*(j + i*nlist[3])), diag.at(2*(j + i*nlist[3])+1)))
                                ),
                                ObjId(diag) = fUnk,
                                fUnk(TComplex, nlist[3]),
                                let(list := nt.params[2].list,  # here we assume FList(TReal, [...])
                                    clist := List([1..Length(list)/2], i->Cplx(list[2*i-1], list[2*i])),
                                    fc := FList(TComplex, clist),
                                    gf := fTensor(fBase(i), fId(nlist[3])),
                                    fCompose(fc, gf)
                                )
                            ),
                            [[ PrunedPRDFT(nlist[2], -1, iblk, ipats[2]),  # stage 1a: PRDFT y
                                PrunedDFT(nlist[1], -1, iblk, ipats[1]),    # stage 1b: DFT z
                                IOPrunedConv(nlist[3], hfunc, oblk, opats[3], iblk, ipats[3], true), # stage 3+4+5: complex conv in x
                                PrunedIDFT(nlist[1], 1, oblk, opats[1]), # stage 6b: iDFT in z
                                PrunedIPRDFT(nlist[2], 1, oblk, opats[2]),   # stage 6a: iPRDFT in y
                                InfoNt(i)
                            ]]),

       apply := (nt, C, cnt) -> let(prdft1d := C[1],
                                    pdft1d := C[2],
                                    iopconv := C[3],
                                    ipdft1d := C[4],
                                    iprdft1d := C[5],
                                    i := cnt[6].params[1],
                                    nlist := nt.params[1],
                                    n1 := nlist[1],
                                    nfreq := nlist[2]/2+1,
                                    n2 := nlist[2],
                                    n3 := nlist[3],
                                    oblk := nt.params[3],
                                    opats := nt.params[4],
                                    iblk := nt.params[5],
                                    ipats := nt.params[6],
                                    ns1 := iblk * Length(ipats[1]),
                                    ns2 := iblk * Length(ipats[2]),
                                    ns3 := iblk * Length(ipats[3]),
                                    nd1 := oblk * Length(opats[1]),
                                    nd2 := oblk * Length(opats[2]),
                                    nd3 := oblk * Length(opats[3]),
                                    stage1 := Tensor(I(ns1), prdft1d),
                                    stage2 := RC(Tensor(pdft1d, I(nfreq))),
                                    stage12 := Tensor(I(nfreq * n2), L(2*ns3, ns3)) * Tensor(stage2 * stage1, I(ns3)),
                                    stage543c := IDirSum(i, iopconv),
                                    stage543 := RC(stage543c),
                                    stage6 := RC(Tensor(ipdft1d, I(nfreq))),
                                    stage7 := Tensor(I(nd1), iprdft1d),
                                    stage67 := Tensor(stage7 * stage6, I(nd3)) * Tensor(I(nfreq * n2), L(2*nd3, 2)),
                                    conv3dr := stage67 * stage543 * stage12,
                                    conv3dr
                            ),
    ),
#======================

    IOPrunedMDRConv_2D_2trip_xy_freqdata := rec(
       forTransposition := true,
       applicable :=  (self, nt) >> not nt.hasTags() and Length(nt.params[1]) = 2 and nt.params[7],
       children := nt -> let(nlist := nt.params[1],
                            diag := nt.params[2],
                            oblk := nt.params[3],
                            opats := nt.params[4],
                            iblk := nt.params[5],
                            ipats := nt.params[6],
                            nfreq := nlist[2]/2+1,
                            i := Ind(nfreq),
                            hfunc := Cond(ObjId(diag) = Lambda,
                                let(j := Ind(nlist[1]),
                                    pos := i + j*nfreq,
                                    Lambda(j, cxpack(diag.at(2*pos), diag.at(2*pos+1)))
                                ),
                                ObjId(diag) = fUnk,
                                fUnk(TComplex, nlist[1]),
                                let(list := nt.params[1].list,  # here we assume FList(TReal, [...])
                                    clist := List([1..Length(list)/2], i->Cplx(list[2*i-1], list[2*i])),
                                    fc := FList(TComplex, clist),
                                    gf := fTensor(fBase(i), fId(nlist[1])),
                                    fCompose(fc, gf)
                                )
                            ),
                            [[ PrunedPRDFT(nlist[2], -1, iblk, ipats[2]),  # stage 1: PRDFT x
                                IOPrunedConv(nlist[1], hfunc, oblk, opats[1], iblk, ipats[1], true), # stage 2+3+4: complex conv in y
                                PrunedIPRDFT(nlist[2], 1, oblk, opats[1]),   # stage 5: iPRDFT in x
                                InfoNt(i)
                            ]]),

       apply := (nt, C, cnt) -> let(prdft1d := C[1],
                                    iopconv := C[2],
                                    iprdft1d := C[3],
                                    i := cnt[4].params[1],
                                    nlist := nt.params[1],
                                    n1 := nlist[1],
                                    nfreq := nlist[1]/2+1,
                                    n2 := nlist[2],
                                    oblk := nt.params[3],
                                    opats := nt.params[4],
                                    iblk := nt.params[5],
                                    ipats := nt.params[6],
                                    ns1 := iblk * Length(ipats[1]),
                                    ns2 := iblk * Length(ipats[2]),
                                    nd1 := oblk * Length(opats[1]),
                                    nd2 := oblk * Length(opats[2]),

#                                    stage1 := L(2*nfreq*ns3*ns2, ns3) * Tensor(I(ns2), Tensor(L(2*nfreq, 2) * prdft1d, I(ns3))) * Tensor(L(ns2*ns1, ns2), I(ns3)),
#                                    stage2 := Tensor(I(ns3), Tensor(RC(pdft1d), I(nfreq))),
#                                    pp := Tensor(L(ns3*n2*nfreq, n2*nfreq), I(2)) * Tensor(I(ns3), L(2*nfreq*n2, nfreq)),
#                                    ppi := Tensor(I(nd3), L(2*nfreq*n2, 2*n2)) * Tensor(L(nd3*n2*nfreq, nd3), I(2)),
#                                    stage543 := ppi * IDirSum(i, RC(iopconv)) * pp,
#                                    stage76 := Tensor(L(nd2*nd1, nd1), I(nd3)) * Grp(Tensor((Tensor(I(nd2), iprdft1d * L(2*nfreq, nfreq)) *
#                                        Tensor(RC(ipdft1d), I(nfreq))), I(nd3)) * L(2*nfreq*nd3*n2, 2*nfreq*n2)),

                                    stage1 := Tensor(I(ns2), prdft1d),
                                    pp := Tensor(L(ns2*nfreq, nfreq), I(2)),
                                    ppi := Tensor(L(nd2*nfreq, nd2), I(2)),
                                    stage432 := IDirSum(i, RC(iopconv)),
                                    stage5 :=  Tensor(I(nd2), iprdft1d),

                                    conv2dr := Grp(stage5 * ppi) * stage432 * Grp(pp * stage1),
                                    conv2dr
                            ),
    ),
#     IOPrunedMDRConv_2D_2trip_yx_freqdata := rec(
#       forTransposition := true,
#       applicable :=  (self, nt) >> not nt.hasTags() and Length(nt.params[1]) = 2 and nt.params[7],
#       children := nt -> let(nlist := nt.params[1],
#                            diag := nt.params[2],
#                            oblk := nt.params[3],
#                            opats := nt.params[4],
#                            iblk := nt.params[5],
#                            ipats := nt.params[6],
#                            nfreq := nlist[1]/2+1,
#                            i := Ind(nfreq),
#                            hfunc := Cond(ObjId(diag) = Lambda,
#                                let(j := Ind(nlist[2]),
#                                    pos := i + j*nfreq,
#                                    Lambda(j, cxpack(diag.at(2*pos), diag.at(2*pos+1)))
#                                ),
#                                ObjId(diag) = fUnk,
#                                fUnk(TComplex, nlist[2]),
#                                let(list := nt.params[2].list,  # here we assume FList(TReal, [...])
#                                    clist := List([1..Length(list)/2], i->Cplx(list[2*i-1], list[2*i])),
#                                    fc := FList(TComplex, clist),
#                                    gf := fTensor(fBase(i), fId(nlist[2])),
#                                    fCompose(fc, gf)
#                                )
#                            ),
#                            [[ PrunedPRDFT(nlist[1], -1, iblk, ipats[1]),  # stage 1: PRDFT y
#                                IOPrunedConv(nlist[2], hfunc, oblk, opats[2], iblk, ipats[2], true), # stage 2+3+4: complex conv in x
#                                PrunedIPRDFT(nlist[1], 1, oblk, opats[1]),   # stage 5: iPRDFT in y
#                                InfoNt(i)
#                            ]]),
#
#       apply := (nt, C, cnt) -> let(prdft1d := C[1],
#                                    iopconv := C[2],
#                                    iprdft1d := C[3],
#                                    i := cnt[4].params[1],
#                                    nlist := nt.params[1],
#                                    n1 := nlist[1],
#                                    nfreq := nlist[1]/2+1,
#                                    n2 := nlist[2],
#                                    oblk := nt.params[3],
#                                    opats := nt.params[4],
#                                    iblk := nt.params[5],
#                                    ipats := nt.params[6],
#                                    ns1 := iblk * Length(ipats[1]),
#                                    ns2 := iblk * Length(ipats[2]),
#                                    nd1 := oblk * Length(opats[1]),
#                                    nd2 := oblk * Length(opats[2]),
#
##                                    stage1 := L(2*nfreq*ns3*ns2, ns3) * Tensor(I(ns2), Tensor(L(2*nfreq, 2) * prdft1d, I(ns3))) * Tensor(L(ns2*ns1, ns2), I(ns3)),
##                                    stage2 := Tensor(I(ns3), Tensor(RC(pdft1d), I(nfreq))),
##                                    pp := Tensor(L(ns3*n2*nfreq, n2*nfreq), I(2)) * Tensor(I(ns3), L(2*nfreq*n2, nfreq)),
##                                    ppi := Tensor(I(nd3), L(2*nfreq*n2, 2*n2)) * Tensor(L(nd3*n2*nfreq, nd3), I(2)),
##                                    stage543 := ppi * IDirSum(i, RC(iopconv)) * pp,
##                                    stage76 := Tensor(L(nd2*nd1, nd1), I(nd3)) * Grp(Tensor((Tensor(I(nd2), iprdft1d * L(2*nfreq, nfreq)) *
##                                        Tensor(RC(ipdft1d), I(nfreq))), I(nd3)) * L(2*nfreq*nd3*n2, 2*nfreq*n2)),
#
#                                    stage1 := L(2*nfreq*ns3*ns2, ns3) * Tensor(I(ns2), Tensor(L(2*nfreq, 2) * prdft1d, I(ns3))) * Tensor(L(ns2*ns1, ns2), I(ns3)),
#                                    pp := Tensor(L(ns3*n2*nfreq, n2*nfreq), I(2)) * Tensor(I(ns3), L(2*nfreq*n2, nfreq)),
#                                    ppi := Tensor(I(nd3), L(2*nfreq*n2, 2*n2)) * Tensor(L(nd3*n2*nfreq, nd3), I(2)),
#                                    stage432 := ppi * IDirSum(i, RC(iopconv)) * pp,
#                                    stage5 := Tensor(L(nd2*nd1, nd1), I(nd3)) * Grp(Tensor((Tensor(I(nd2), iprdft1d * L(2*nfreq, nfreq)) *
#                                        Tensor(RC(ipdft1d), I(nfreq))), I(nd3)) * L(2*nfreq*nd3*n2, 2*nfreq*n2)),
#
#                                    conv2dr := stage5 * stage432 * stage1,
#                                    conv2dr
#                            ),
#    ),

));


text:
# Built-in
LcmInt(12, 18);
# 36

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

NewRulesFor(DFT, rec(
    DFT_PRDFT := rec(
    switch := false,
    applicable := t -> Rows(t) >= 3 and t.getTags()=[],
    children := t -> [[ PRDFT1(Rows(t),t.params[2]) ]],
    apply := (t, ch, nonterms) -> let(n := Rows(t),
        bc := Mat([[1, E(4)], [1, -E(4)]]),
        DirectSum(
        Mat([[1,0]]),
        LIJ(n-1).transpose() * Cond(IsEvenInt(n),
            DirectSum(Tensor(I((n-2)/2), bc), Mat([[1,0]])),
            Tensor(I((n-1)/2), bc))) *
        ch[1])

    )
));

NewRulesFor(DFT2, rec(
    DFT2_PRDFT2 := rec(
    switch := false,
    applicable := t -> Rows(t) >= 3,
    children := t -> [[ PRDFT2(Rows(t), t.params[2]) ]],
    apply := (t, ch, nonterms) -> let(n := Rows(t), sgn := (-1)^t.params[2],
        bc := Mat([[1, E(4)], sgn*[1, -E(4)]]),
        DirectSum(
        Mat([[1,0]]),
        LIJ(n-1).transpose() * Cond(IsEvenInt(n),
            DirectSum(Tensor(I((n-2)/2), bc), Mat([[0,E(4)]])),
            Tensor(I((n-1)/2), bc))) *
        ch[1])

    )
));

NewRulesFor(DFT3, rec(
    DFT3_PRDFT3 := rec(
    switch := false,
    applicable := t -> Rows(t) >= 3,
    children := t -> [[ PRDFT3(Rows(t), t.params[2]) ]],
    apply := (t, ch, nonterms) -> let(n := Rows(t),
        bc := Mat([[1, E(4)], [1, -E(4)]]),
        LIJ(n).transpose() *
        Cond(IsEvenInt(n),
         Tensor(I(n/2), bc),
         DirectSum(Tensor(I((n-1)/2), bc), Mat([[1,0]]))) *
        ch[1])

    )
));

NewRulesFor(DFT4, rec(
    DFT4_PRDFT4 := rec(
    switch := false,
    applicable := t -> Rows(t) >= 3,
    children := t -> [[ PRDFT4(Rows(t),t.params[2]) ]],
    apply := (t, ch, nonterms) -> let(n := Rows(t), sgn := (-1)^t.params[2], jj := E(4)^t.params[2],
        bc := Mat([[1, E(4)], sgn*[1, -E(4)]]),
        LIJ(n).transpose() *
        Cond(IsEvenInt(n),
         Tensor(I(n/2), bc),
                 # last element could be real or imaginary based on
                 # rotation (t.params[2]) (but not both), so we scale both slots
         DirectSum(Tensor(I((n-1)/2), bc), Mat([[jj,jj]]))) *
        ch[1])

    )
));

ext: gi


text:
#
# usage: (!!!Achtung: This will overwrite all manual modifications)
#
# Read("~/Workspace/Chevalley.gap/init.gi"); Read(Filename(home_dir,"load.gi"));
# Read(Filename(e8_dir_char3,"init_char3.gi")); Read(Filename(e8_dir_char3,"gen_files.gi"));
#

for l in labels do
    f:=Filename(e8_dir_char3,Concatenation(l,".gi"));
    #PrintTo(f,l);
    PrintTo(f,"#\n");
    AppendTo(f,"# usage:\n");
    AppendTo(f,"# Read(\"~/Workspace/Chevalley.gap/init.gi\"); Read(Filename(home_dir,\"load.gi\")); ");
    AppendTo(f,"Read(Filename(e8_dir_char3,\"init_char3.gi\")); Read(Filename(e8_dir_char3,\"");
    AppendTo(f,Concatenation(l,".gi"));
    AppendTo(f,"\"));\n");
    AppendTo(f,"#\n");
    AppendTo(f,"# Read(Filename(e8_dir_char3,\"init_char3.gi\")); Read(Filename(e8_dir_char3,\"");
    AppendTo(f,Concatenation(l,".gi"));
    AppendTo(f,"\"));\n\n");

    AppendTo(f,"label:=\"");
    AppendTo(f,l);
    AppendTo(f,"\";\n");

    AppendTo(f,"orb_nr:=Position(labels,label);\n\n");
    AppendTo(f,"orb:=AllClasses(orbs)[orb_nr];\n");
    AppendTo(f,"info:=handleClassShort(orb);\n\n");

    AppendTo(f,"Print(\"Consider the class \",Label(orb),\" in characteristic \",Characteristic(ring(sys)),\":\\n\");\n");
    AppendTo(f,"Print(\"\\tBorel representative \\n\\t\",coefficients(BorelRep(orb)),\"\\n\");\n");
    AppendTo(f,"Print(\"\\tconnected C_U(u) \\n\\t\",coefficients(info[1]),\"\\n\");\n");
    AppendTo(f,"Print(\"\\tconnected C_U(u) in Levi \\n\\t\",coefficients(info[2]),\"\\n\");\n");
    AppendTo(f,"Print(\"\\tconnected double C_U(u) \\n\\t\",coefficients(info[3]),\"\\n\");\n");
od;


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

NewRulesFor(DFT, rec(

    DFT_PD := rec(
        forTransposition := false,
        maxSize          := 13,

        core := meth(self, N, k, root, is_real)
            local coeffs, l1, l2, n, circ, toep;
            coeffs := List([0..N-2], x->E(N)^(k*root^(-x) mod N));
            n := N-1;
            l1 := coeffs{[1..n/2]};
            l2 := coeffs{[n/2+1..n]};
            circ := 1/2 * (l1 + l2);
            toep := Reversed(1/2 * Concat(Drop(l2 - l1, 1), l1 - l2));
            if is_real then toep := toep / E(4); fi;
            return
            DirectSum(I(1),
                  transforms.filtering.Circulant(circ).terminate(),
                  transforms.filtering.Toeplitz(toep).terminate());
        end,

        A := N -> let(nn := (N-1)/2, ones := m -> fConst(m, 1),
             VStack(HStack(RowVec(ones(nn+1)), O(1, nn)),
               HStack(VStack(-2*ColVec(ones(nn)), O(nn, 1)), I(2*nn)))),

        applicable     := (self, nt) >> nt.params[1] > 2 and nt.params[1] <= self.maxSize and IsPrime(nt.params[1]) and not nt.hasTags(),

        apply := (self,nt,C,cnt) >> let(N := nt.params[1], k := nt.params[2], root := PrimitiveRootMod(N),
            Scat(RR(N, 1, root)) *
            DirectSum(I(1), Tensor(F(2), I((N-1)/2))) *
            Mat(MatSPL(
                self.core(N, k, root, false) *
                self.A(N))) *
            DirectSum(I(1), Tensor(F(2), I((N-1)/2))) *
    #        DirectSum(I(1), OS(N-1, -1)) *
            Gath(RR(N, 1, 1/root mod N))
        )
#D    isApplicable     := (self, P) >> P[1] > 2 and P[1] <= self.maxSize and IsPrime(P[1]) and P[3] = [],
#D
#D    rule := (self,P,C) >> let(N := P[1], k := P[2], root := PrimitiveRootMod(N),
#D        Scat(RR(N, 1, root)) *
#D        DirectSum(I(1), Tensor(F(2), I((N-1)/2))) *
#D        Mat(MatSPL(
#D            self.core(N, k, root, false) *
#D            self.A(N))) *
#D        DirectSum(I(1), Tensor(F(2), I((N-1)/2))) *
#D#        DirectSum(I(1), OS(N-1, -1)) *
#D        Gath(RR(N, 1, 1/root mod N))
#D    )
    )
));


text:
# Constants
BindGlobal("Graph6ByteInteger", 62);
BindGlobal("Graph6WordInteger", 258047);
BindGlobal("Graph6DwordInteger", 68719476735);
BindGlobal("Graph6MaxInteger", Graph6DwordInteger);
BindGlobal("Graph6BitVector", List([0..5], i -> 2^(5-i)));
BindGlobal("Graph6IntVector", List([0..5], i -> 64^(5-i)));

# graph6 integer encoding
BindGlobal("Graph6EncodeInteger", function(x)
    if not IsInt(x) or x < 0 or x >= Graph6MaxInteger then
        Error("not an integer in allowed range");
        return fail;
    fi;
    if x <= Graph6ByteInteger then
        return [x+63];
    elif x <= Graph6WordInteger then
        return [63, Int(x/4096), Int(x/64), x] mod 64 + 63;
    else
        return Concatenation([126, 126],
                List(Graph6IntVector, y -> Int(x/y) mod 64 + 63));
    fi;
end);

# graph6 integer decoding
BindGlobal("Graph6DecodeInteger", function(i, s)
    if s[i] <> 126 then
        return [i+1, s[i] - 63];
    elif s[i+1] <> 126 then
        return [i+4, (s{[i+1..i+3]} - 63)*Graph6IntVector{[4..6]}];
    else
        return [i+8, (s{[i+2..i+7]} - 63)*Graph6IntVector];
    fi;
end);

# Convert an integer to a bitstring of length w
BindGlobal("IntToBits", function(x, w)
    return List([1..w], i -> Int(x / 2^(w-i)) mod 2);
end);

# Convert a bitstring to a list of integers given the width w
BindGlobal("BitsToInt", function(b, w)
    local out, i;
    out := [];
    i := 0;
    while i+w <= Length(b) do
        Add(out, b{[i+1..i+w]} * List([1..w], j -> 2^(w-j)));
        i := i+w;
    od;
    return out;
end);

# Convert a bitstring to a graph6 string with padding p
BindGlobal("BitsToString", function(b, p)
    local c;
    c := Concatenation(b, ListWithIdenticalEntries(-Length(b) mod 6, p));
    return List([0,6..Length(c)-6], i -> c{[i+1..i+6]}*Graph6BitVector + 63);
end);

# Convert a graph6 string to a bitstring
BindGlobal("StringToBits", function(s)
    return Concatenation(List(s, x -> List(Graph6BitVector, y -> Int((x-63)/y) mod 2)));
end);

# graph6 string
BindGlobal("Graph6String", function(G)
    local A, s;
    A := CollapsedAdjacencyMat(Group(()), G);
    s := Concatenation(Graph6EncodeInteger(G.order),
        BitsToString(Concatenation(List([1..G.order], i -> A[i]{[1..i-1]})), 0));
    return List(s, CharInt);
end);

# Read graph from graph6 string
BindGlobal("GraphFromGraph6String", function(s)
    local A, b, n, i, j, t;
    s := List(s, IntChar);
    t := Graph6DecodeInteger(1, s);
    i := t[1];
    n := t[2];
    b := StringToBits(s{[i..Length(s)]});
    A := NullMat(n, n);
    i := 1;
    for j in [2..n] do
        A[j]{[1..j-1]} := b{[i..i+j-2]};
        A{[1..j-1]}[j] := b{[i..i+j-2]};
        i := i+j-1;
    od;
    return AdjFunGraph([1..n], MatrixAdjacency(A));
end);

# sparse6 string
BindGlobal("Sparse6String", function(G)
    local i, j, b, c, v, w;
    w := Log2Int(G.order-1)+1;
    b := [];
    v := 1;
    for i in [1..G.order] do
        for j in [1..i] do
            if IsVertexPairEdge(G, i, j) then
                if i = v+1 then
                    c := 1;
                else
                    if i > v then
                        Add(b, 1);
                        Append(b, IntToBits(i-1, w));
                    fi;
                    c := 0;
                fi;
                Add(b, c);
                Append(b, IntToBits(j-1, w));
                v := i;
            fi;
        od;
    od;
    if G.order in [2,4,8,16] and (-Length(b)) mod 6 > w
            and Length(Adjacency(G, G.order-1)) > 0
            and Length(Adjacency(G, G.order)) = 0 then
        Add(b, 0);
    fi;
    return List(Concatenation([58], Graph6EncodeInteger(G.order),
                                BitsToString(b, 1)), CharInt);
end);

# Read graph from sparse6 string
BindGlobal("GraphFromSparse6String", function(s)
    local A, b, t, i, m, n, v, w, c, x, y, z;
    if s[1] = ':' then
        s := s{[2..Length(s)]};
    fi;
    s := List(s, IntChar);
    t := Graph6DecodeInteger(1, s);
    i := t[1];
    n := t[2];
    w := Log2Int(n-1)+1;
    b := StringToBits(s{[i..Length(s)]});
    A := List([1..n], j -> []);
    i := 1;
    v := 1;
    z := BitsToInt(b, w+1);
    m := 2^w;
    for y in z do
        c := Int(y/m);
        x := y mod m + 1;
        if c = 1 then
            v := v+1;
        fi;
        if v > n or x > n then
            break;
        fi;
        if x > v then
            v := x;
        else
            Add(A[x], v);
            Add(A[v], x);
        fi;
    od;
    return AdjFunGraph([1..n], ListAdjacency(A));
end);

# sparse6 string
BindGlobal("IncrementalSparse6String", function(G, H)
    local i, j, b, c, v, w;
    if G.order <> H.order then
        Error("graphs have different orders");
        return fail;
    fi;
    w := Log2Int(G.order-1)+1;
    b := [];
    v := 1;
    for i in [1..G.order] do
        for j in [1..i] do
            if IsVertexPairEdge(G, i, j) <> IsVertexPairEdge(H, i, j) then
                if i = v+1 then
                    c := 1;
                else
                    if i > v then
                        Add(b, 1);
                        Append(b, IntToBits(i-1, w));
                    fi;
                    c := 0;
                fi;
                Add(b, c);
                Append(b, IntToBits(j-1, w));
                v := i;
            fi;
        od;
    od;
    if G.order in [2,4,8,16] and (-Length(b)) mod 6 > w
            and Adjacency(G, G.order-1) = Adjacency(H, G.order-1)
            and Adjacency(G, G.order) = Adjacency(H, G.order) then
        Add(b, 0);
    fi;
    return List(Concatenation([59], BitsToString(b, 1)), CharInt);
end);

# Read graph from incremental sparse6 string
BindGlobal("GraphFromIncrementalSparse6String", function(s, H)
    local A, B, b, i, m, v, w, c, x, y, z;
    if s[1] = ';' then
        s := s{[2..Length(s)]};
    fi;
    B := List([1..H.order], j -> Adjacency(H, j));
    A := List(B, ShallowCopy);
    s := List(s, IntChar);
    w := Log2Int(H.order-1)+1;
    b := StringToBits(s);
    i := 1;
    v := 1;
    z := BitsToInt(b, w+1);
    m := 2^w;
    for y in z do
        c := Int(y/m);
        x := y mod m + 1;
        if c = 1 then
            v := v+1;
        fi;
        if v > H.order or x > H.order then
            break;
        fi;
        if x > v then
            v := x;
        elif x in A[v] then
            Remove(A[x], Position(A[x], v));
            if x <> v then
                Remove(A[v], Position(A[v], x));
            fi;
        else
            Add(A[x], v);
            Add(A[v], x);
        fi;
    od;
    return AdjFunGraph([1..H.order], ListAdjacency(A));
end);

# auto6 string
BindGlobal("Auto6String", function(G)
    local b, g, h, i, l, p, r, s, w, x, sch;
    g := GeneratorsOfGroup(G.group);
    l := [1..Length(g)];
    i := 0;
    for p in [1..Length(g)] do
        if g[p] = () then
            i := i+1;
        else
            l[p] := l[p] - i;
        fi;
    od;
    sch := List(G.schreierVector, function(y)
                                    if y < 0 then
                                        return y;
                                    else
                                        return l[y];
                                    fi;
                                  end);
    g := Filtered(g, y -> y <> ());
    w := Log2Int(G.order-1)+1;
    r := Length(G.representatives);
    b := IntToBits(r, w);
    for i in [1..r] do
        Append(b, IntToBits(G.representatives[i]-1, w));
        Append(b, IntToBits(Length(G.adjacencies[i]), w));
        for x in G.adjacencies[i] do
            Append(b, IntToBits(x-1, w));
        od;
    od;
    for h in g do
        for i in [1..G.order] do
            Append(b, IntToBits(i^h - 1, w));
        od;
    od;
    if Length(g) > 1 then
        s := Log2Int(Length(g)) + 1;
        for i in [1..G.order] do
            if sch[i] < 0 then
                Append(b, IntToBits(0, s));
            else
                Append(b, IntToBits(sch[i], s));
            fi;
        od;
    fi;
    return List(Concatenation([33], Graph6EncodeInteger(G.order),
                Graph6EncodeInteger(Length(g)), BitsToString(b, 0)), CharInt);
end);

# Read graph from auto6 string
BindGlobal("GraphFromAuto6String", function(s)
    local G, A, R, S, a, b, d, r, t, g, i, j, k, n, v, w;
    if s[1] = '!' then
        s := s{[2..Length(s)]};
    fi;
    s := List(s, IntChar);
    t := Graph6DecodeInteger(1, s);
    i := t[1];
    n := t[2];
    w := Log2Int(n-1)+1;
    t := Graph6DecodeInteger(i, s);
    i := t[1];
    g := t[2];
    b := StringToBits(s{[i..Length(s)]});
    d := BitsToInt(b, w);
    r := d[1];
    if r = 0 then
        r := n;
    fi;
    i := 2;
    R := [];
    A := [];
    for j in [1..r] do
        Add(R, d[i]+1);
        Add(A, d{[i+2..i+d[i+1]+1]}+1);
        i := i + d[i+1] + 2;
    od;
    if g = 0 then
        G := ();
        v := List([1..n], x -> -Position(R, x));
    else
        G := [];
        for j in [1..g] do
            Add(G, PermList(d{[i..i+n-1]}+1));
            i := i+n;
        od;
        if g = 1 then
            v := List([1..n], function(x)
                local y;
                y := Position(R, x);
                if y = fail then
                    return 1;
                else
                    return -y;
                fi;
            end);
        else
            k := w*(i-1);
            t := Log2Int(g)+1;
            v := BitsToInt(b{[k+1..k+n*t]}, t);
            for j in [1..r] do
                v[R[j]] := -j;
            od;
        fi;
    fi;
    return rec(
        isGraph := true,
        order := n,
        names := [1..n],
        representatives := Immutable(R),
        adjacencies := A,
        group := Group(G),
        schreierVector := Immutable(v)
    );
end);

max_stars_repo_name: jaanos/gap-graphs

max_issues_repo_name: jaanos/gap-graphs

max_forks_repo_name: jaanos/gap-graphs


text:
Enqueue := function(v, x)
   Add(v[1], x);
end;

Dequeue := function(v)
   local n, x;
   n := Size(v[2]);
   if n = 0 then
      v[2] := Reversed(v[1]);
      v[1] := [ ];
      n := Size(v[2]);
      if n = 0 then
         return fail;
      fi;
   fi;
   return Remove(v[2], n);
end;

# a new queue
v := [[], []];

Enqueue(v, 3);
Enqueue(v, 4);
Enqueue(v, 5);
Dequeue(v);
# 3
Enqueue(v, 6);
Dequeue(v);
# 4
Dequeue(v);
# 5
Dequeue(v);
# 6
Dequeue(v);
# fail

max_stars_repo_name: djgoku/RosettaCodeData

max_issues_repo_name: djgoku/RosettaCodeData

max_forks_repo_name: djgoku/RosettaCodeData


text:
doors := function(n)
  local a,j,s;
  a := [ ];
  for j in [1 .. n] do
    a[j] := 0;
  od;
  for s in [1 .. n] do
    j := s;
    while j <= n do
      a[j] := 1 - a[j];
      j := j + s;
    od;
  od;
  return Filtered([1 .. n], j -> a[j] = 1);
end;

doors(100);
# [ 1, 4, 9, 16, 25, 36, 49, 64, 81, 100 ]


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

Class(arctan, AutoFoldExp, rec(
  ev := self >> self._ev(self.args).ev(),
  computeType := self >> TReal,
));

Class(TArcTan, Tagged_tSPL_Container, rec(
    abbrevs :=  [ () -> []],
    dims := self >> [1, 2],
    transpose := self >> CopyFields(self, rec(transposed := not self.transposed)),
    isReal := True
));

NewRulesFor(TArcTan, rec(
    TArcTan_Base := rec(
        applicable := True,
        forTransposition := false,
        apply := (t, C, Nonterms) -> let( x := var.fresh_t("r", TDouble), y := var.fresh_t("r", TDouble),
            OLCompose(BinOp(1, Lambda([x,y], arctan(x,y)))))
    )
));

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:
InputTextFile("input.txt");
s := ReadAll(f);;  # two semicolons to hide the result, which may be long
CloseStream(f);


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

sys:=AlgebraicU(ChevalleyAdj("F",4,GF(2)));
w:=Witt(1,sys);
Print("Created Witt ",w,"\n");
Print("\t[",IsWitt(w)=true,"] IsWitt(w)=true\n");


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

#F ExpFunc - ExpDiag element function.
#F
Class(ExpFunc, Lambda, rec(
    _isExpFunc := true,
    # this is not a ranked index function, not related to GT,
    # GT should not mess with it in the first place
    downRank := (self, loopid, ind) >> self,
    downRankFull := (self, inds) >> self,
    upRank := self >> self,
    upRankBy := (self, n) >> self,
    rotate := (self, n) >> self,
    split := (self, loopid, inner_its, outer_its) >>  self,
));

IsExpFunc := (obj) -> IsRec(obj) and IsBound(obj._isExpFunc) and obj._isExpFunc;

#F ExpDiag(<func>, <d>)
#F
#F ExpDiag(ExpFunc( [a,b], sub(a,b) ), 3)
#F
#F corresponds to matrix:
#F   [ [ 1 -1  0  0  0  0 ],
#F     [ 0  0  1 -1  0  0 ],
#F     [ 0  0  0  0  1 -1 ] ]
#F

Class(ExpDiag, BaseMat, rec(
    _short_print := true,
    abbrevs := [ (L, d) -> Checked( IsExpFunc(L) and IsPosIntSym(d), [L, d] )],

    new := (self, L, d) >> SPL(WithBases(self, rec(element := L, TType := L.range(), d := d))).setDims(),

    #-----------------------------------------------------------------------
    dims := self >> [self.d, self.d*Length(self.element.vars)],
    #-----------------------------------------------------------------------
    isPermutation := self >> false,
    #-----------------------------------------------------------------------
    isReal := self >> let(
	t := self.element.range(),
        tt := Cond(IsVecT(t), t.t, t),
	tt <> TComplex and ObjId(tt)<>T_Complex),
    #-----------------------------------------------------------------------
    toAMat         := abstract(),
    #-----------------------------------------------------------------------
    transpose      := abstract(),
    conjTranspose  := abstract(),
    inverse        := abstract(),
    #-----------------------------------------------------------------------
    arithmeticCost := abstract(),

    sums           := self >> self,
    rChildren      := (self) >> [self.element, self.d],
    rSetChild      := rSetChildFields("element", "d"),
    area           := self >> self.d
));


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

# Auxiliary Matrices
# ==================
#
# This file contains useful matrices to be used in rules
#

IsRealMat := mat -> ForAll(mat, r -> ForAll(r, IsRealNumber));

ReMat := mat -> MapMat(mat, e -> re(e));

# Permutations
# ------------

#F
#F GenPerm0(generator_func, size)
#F   returns a permutation generated by applying generator_func to
#F   list [0..size-1]
#F
GenPerm0 := function(size, generator_func)
     Constraint(IsInt(size) and size > 0);
     return PermList(List([0..size-1], generator_func) + 1);
end;

#F
#F GenPerm1(generator_func, size)
#F   returns a permutation generated by applying generator_func to
#F   list [1..size]
#F
GenPerm1 := function(size, generator_func)
     Constraint(IsInt(size) and size > 0);
     return PermList(List([1..size], generator_func));
end;

#F
#F LexMap(r,s)
#F   returns list of lexicographical ordering pairs
#F
LexMap := function(r,s)
   local i,j, result;
   result:=[];
   for i in [0..r-1] do
      for j in [0..s-1] do
         Append(result, [[i,j]]);
      od;
   od;
   return result;
end;

#F
#F ToeplitzMat([..elements..])
#F    Returns a Toeplitz matrix given by element list
#F
ToeplitzMat := function(elements)
    local i, j, n, Ls, mat;
    Constraint(IsList(elements));
    Constraint(Length(elements) > 0);
    Constraint(IsOddInt(Length(elements)));
    n := (Length(elements)+1)/2;
    mat := [ ];
    for i in [1..n] do
        Ls := List([1..n], k->0);
	for j in [1..n] do
	    Ls[j] := elements[n+i-j];
	od;
	Add(mat,Ls);
    od;
    return mat;
end;

# -----------------------------------------------------------------------------
# Matrices used in some breakdown rules
# NOTE: these should be moved to transforms/*
# -----------------------------------------------------------------------------

# perm2( <n> )
#   (x_0, .., x_n-1) --> (x_0, x_n-1, x_1, x_n-2, ..)
# For even <n> it is L(n,n/2) * IJ(n,n/2);
# for odd <n>    OddStride(n,(n+1)/2) * DirectSum(I((n+1)/2), J((n-1)/2))
#
perm2 := n -> When(IsEvenInt(n), M(n,n/2), LIJ(n));

# perm3( <n> )
#   (x_0, .., x_n-1) --> (x_0, x_n/2, x_1, x_n/2+1, ..),
# can be expressed by the previous two:
#   perm3(n) = perm2(n) * (I_n/2 dirsum J(n/2))
# for 2|n it is the stride permutation L^n_n/2
perm3 := n -> When(IsEvenInt(n), L(n, n/2), OS(n, (n+1)/2));

# perm4( <n> )
#   permutes a direct sum of 2x2 blocks
#   into star form, first block outermost etc.
#   for n even it is equal to perm2
perm4 := n -> When(IsEvenInt(n), LIJ(n), Z(n).transpose()*LIJ(n));

# perm5( <n> )
#   reverse the odd indices:
#   (x_0, x_1, .., x_2k) --> (x_0, x_2k-1, x_2, x_2k-3, .., x_1, x_2k)
perm5 := n -> LIJ(n) * L(n,2);

# perm6( <n> )
#   a permutation for rule DCT4_Iterative.
perm6 := SymFunc("perm6", function ( n )
  local k, p, i;

  k := LogInt(n, 2);
  p := PermSPL(perm4(n));
  for i in [1..k - 2] do
    p :=
      p *
      TensorProductPerm(
        [ 2^(i - 1), n/2^(i - 1) ],
        [ (),
          DirectSumPerm(
            [ n/2^i, n/2^i ],
            [ PermSPL(perm4(n/2^i)),
              PermSPL(perm4(n/2^i)) ^ PermSPL(J(n/2^i)) ]
          ) ]
      );
  od;

  return Perm(p, n);
end);

# perm7( <n> )
#   a permutation for converting an RDFT into a DCT2.
perm7 := SymFunc("perm7", function ( n )
  local k, L, i;

  k := (n-1)/2;
  L := [ ];
  for i in [0..k] do
    Add(L, (-1)^(i + k + 1)*i + k);
  od;
  for i in [k + 1..n - 1] do
    Add(L, (-1)^(i + k + 1)*(n - i) + k);
  od;

  return Perm(PermListList(L, [0..n-1]), n);
end);

# Block Matrices
# --------------

# blocks1( <n> )
#   returns an spl that represents a direct sum of
#   DFT_2's in star form. If n is odd then it contains
#   a 1x1 block in the middle.
blocks1 := n -> When(IsEvenInt(n),
    Tensor(I(Int(n/2)), F(2)) ^ perm4(n),
    DirectSum(I(1), Tensor(I(Int(n/2)), F(2))) ^ perm4(n));

blocks4 := (n, k) ->
    DirectSum(I(1), Tensor(I(Int(n/2)), F(2) * Diag([1, k*E(4)]))) ^ perm4(n);

# Sum Matrices
# ------------

# n x n matrix
# diagonal and upper diagonal 1's
sums1 := SymFunc("sums1", function ( n )
  local M, i;

  M := [ ];
  for i in [1..n] do
    Add(M, [i, i, 1]);
  od;
  for i in [1..n-1] do
    Add(M, [i, i+1, 1]);
  od;
  return Sparse(M);
end);

# n x 2n matrix
# [ identity, opposite identity ]
sums2 := SymFunc("sums2", n -> HStack(I(n), J(n)));

# n+1 x 3*n/2 matrix, e.g.,
#[ [   1,   0,   0,   0,   0,   0,   0,   0,   1,   0,   0,   0 ],
#  [   0,   1,   0,   0,   0,   0,   0,   0,   0,   1,   0,   0 ],
#  [   0,   0,   1,   0,   0,   0,   0,   0,   0,   0,   1,   0 ],
#  [   0,   0,   0,   1,   0,   0,   0,   0,   0,   0,   0,   1 ],
#  [   0,   0,   0,   0,   1,   0,   0,   0,   0,   0,   0,   0 ],
#  [   0,   0,   0,   0,   0,   1,   0,   0,   0,   0,   0,  -1 ],
#  [   0,   0,   0,   0,   0,   0,   1,   0,   0,   0,  -1,   0 ],
#  [   0,   0,   0,   0,   0,   0,   0,   1,   0,  -1,   0,   0 ],
#  [   1,   0,   0,   0,   0,   0,   0,   0,  -1,   0,   0,   0 ] ]
sums3 := SymFunc("sums3", function ( n )
  local M, i;

  M := [ ];
  for i in [1..n] do
    Add(M, [i, i, 1]);
  od;
  for i in [1..n/2] do
    Add(M, [i, n + i, 1]);
    Add(M, [n/2 + 1 + i, 3*n/2 - i + 1, -1]);
  od;
  Add(M, [n+1, 1, 1]);

  return Sparse(M);
end);

# n x 2n matrix
# [ identity, -opposite identity ]
sums4 := n -> HStack(I(n), -J(n));

# 2n x n matrix
# [ opposite lower diagonal 1;
#   identity
# ]

# n x 3n/2 matrix, e.g.,
#[ [  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0 ],
#  [  0,  1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  1 ],
#  [  0,  0,  1,  0,  0,  0,  0,  0,  0,  0,  1,  0 ],
#  [  0,  0,  0,  1,  0,  0,  0,  0,  0,  1,  0,  0 ],
#  [  0,  0,  0,  0,  0,  0,  0,  0,  1,  0,  0,  0 ],
#  [  0,  0,  0,  0,  0,  1,  0,  0,  0,  1,  0,  0 ],
#  [  0,  0,  0,  0,  0,  0,  1,  0,  0,  0,  1,  0 ],
#  [  0,  0,  0,  0,  0,  0,  0,  1,  0,  0,  0,  1 ] ]
sums5 := SymFunc("sums5", function ( n )
  local M, i;

  M := [ ];
  for i in [2..n/2] do
    Add(M, [i, i, 1]);
    Add(M, [i, 3*n/2 - i + 2, 1]);
  od;
  Add(M, [n/2 + 1, n + 1, 1]);
  for i in [1..n/2 - 1] do
    Add(M, [n/2 + 1 + i, n/2 + 1 + i, 1]);
    Add(M, [n/2 + 1 + i, n + i + 1, 1]);
  od;

  return Sparse(M);
end);

# a factorization of a matrix like
#   [ [   1/2,     0,     0,     0,     0,     0,     0,     0 ],
#     [  -1/2,     1,     0,     0,     0,     0,     0,     0 ],
#     [   1/2,    -1,     1,     0,     0,     0,     0,     0 ],
#     [  -1/2,     1,    -1,     1,     0,     0,     0,     0 ],
#     [   1/2,    -1,     1,    -1,     1,     0,     0,     0 ],
#     [  -1/2,     1,    -1,     1,    -1,     1,     0,     0 ],
#     [   1/2,    -1,     1,    -1,     1,    -1,     1,     0 ],
#     [  -1/2,     1,    -1,     1,    -1,     1,    -1,     1 ] ]
# into 1 mult and n - 1 subs. If the first column is divided by 1/2,
# then we obtain the transpose - inverse of sums1.
sums6 := SymFunc("sums6", function ( n )
  local T, M, M1, i;

  # small cases
  if n = 1 then
    return Diag([1/2]);
  fi;
  if n = 2 then
    return Mat([[1, 0], [-1, 1]]) * Diag([1/2, 1]);
  fi;

  T    := List([1..n], i -> 1);
  T[1] := 1/2;
  T    := Diag(T);
  M1   := Mat([[1, 0], [-1, 1]]);
  T    := DirectSum(M1, I(n - 2)) * T;
  for i in [3..n - 1] do
    T := DirectSum(I(i - 2), M1, I(n - i)) * T;
  od;
  return DirectSum(I(n - 2), M1) * T;
end);

# sparse1( <n>, <r> )
#   returns an spl for the matrix
#     I_m  -Z_m
#           D_m
#   where m = n/2, Z_m has ones on the lower opposite diagonal,
#   and D_m = cos(r*pi) * diag(1, 2, ..., 2)
#   This matrix occurs in Steidl's rule for DCT2/3.
sparse1 := (n,r) -> CheckedD(IsInt(n), CanBeScalar(r) or IsScalar(r), n mod 2 = 0,
    let(m := n/2,
	d := DirectSum(Diag(cospi(r)), Scale(2*cospi(r), I(n/2-1))),

	BlockMat(
	    [ [I(m), DirectSum(O(1), -J(m-1))],
	      [O(m),          d            ] ])));

# Copy Matrices
# -------------

# 2n x n+1 matrix, e.g.,
#[ [  1,  0,  0,  0,  0 ],
#  [  0,  1,  0,  0,  0 ],
#  [  0,  0,  1,  0,  0 ],
#  [  0,  0,  0,  1,  0 ],
#  [  0,  0,  0,  0,  1 ],
#  [  0,  0,  0,  1,  0 ],
#  [  0,  0,  1,  0,  0 ],
#  [  0,  1,  0,  0,  0 ] ]

copy1 := SymFunc("copy1", function ( n )
  local M, i;

  M := [ ];
  for i in [1..n + 1] do
    Add(M, [i, i, 1]);
  od;
  for i in [1..n - 1] do
    Add(M, [n + 1 + i, n - i + 1, 1]);
  od;

  return Sparse(M);
end);

# 2n x n matrix, e.g.,
#[ [   1,   0,   0,   0 ],
#  [   0,   1,   0,   0 ],
#  [   0,   0,   1,   0 ],
#  [   0,   0,   0,   1 ],
#  [   0,   0,   0,   0 ],
#  [   0,   0,   0,  -1 ],
#  [   0,   0,  -1,   0 ],
#  [   0,  -1,   0,   0 ] ]

copy2 := SymFunc("copy2", function ( n )
  local M, i;

  M := [ ];
  for i in [1..n] do
    Add(M, [i, i, 1]);
  od;
  for i in [1..n - 1] do
    Add(M, [n + 1 + i, n - i + 1, -1]);
  od;

  return Sparse(M);
end);

# Monomial Matrices
# -----------------

# monomial1 ( <n> )
#   Let m = n/4. Then monomial1(<n>) returns an spl of type "sparse"
#   that represents
#   I_m+1  0   0    0
#    0     0   J_m  0
#    0     0   0   -J'_m
#    0   I_m-1 0    0
# where the prime indicates that the entry in the last row
# has opposite sign

monomial1 := SymFunc("monomial1", function ( n )
  local m, L, i;

  m := n/4;
  L := [ ];
  for i in [1..m+1] do
    Add(L, [i, i, 1]);
  od;
  for i in [1..m] do
    Add(L, [m+1+i, 3*m+1-i, 1]);
  od;
  for i in [1..m-1] do
    Add(L, [2*m+1+i, 4*m+1-i, -1]);
  od;

  # the exception
  Add(L, [3*m+1, 3*m+1, 1]);
  for i in [1..m-1] do
    Add(L, [3*m+1+i, m+1+i, 1]);
  od;

  return Sparse(L);
end);


text:
# This file has been generated by the GAP build system,
# do not edit manually!
GAParch=FAKE-GAP-ARCH
GAP_ABI=64
GAP_HPCGAP=no

GAP_VERSION="4.12dev"
GAP_BUILD_VERSION="4.12dev"
GAP_LIBTOOL_CURRENT=8
GAP_LIBTOOL_AGE=0
GAP_KERNEL_MAJOR_VERSION=8
GAP_KERNEL_MINOR_VERSION=0

GAP_BIN_DIR="/workspace/destdir/share/gap"
GAP_LIB_DIR="/workspace/destdir/share/gap"

GAP="/workspace/destdir/bin/gap"
GAC="/workspace/destdir/share/gap/gac"

GAP_CC="cc "
GAP_CXX="c++ -std=gnu++11 "
GAP_CFLAGS="-g -O2"
GAP_CXXFLAGS="-g -O2"
GAP_CPPFLAGS="-I/workspace/destdir/include/gap -DUSE_JULIA_GC=1 -fPIC "
GAP_LDFLAGS="-L/workspace/destdir/lib "
GAP_LIBS="-lgmp"

max_stars_repo_name: MichelJuillard/Yggdrasil

max_issues_repo_name: MichelJuillard/Yggdrasil

max_forks_repo_name: MichelJuillard/Yggdrasil


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

#F VO1dsJ(<n>, <v>) -- vector construct for I(1) \dirsum J(n-1)
#F
Class(VO1dsJ, BaseMat, SumsBase, rec(
    new := (self, n, v) >> SPL(WithBases(self, rec(n := n, v := v))).setDims(),
    dims := self >> [self.n, self.n],
    #-----------------------------------------------------------------------
    rChildren := self >> [self.n, self.v],
    rSetChild := rSetChildFields("n", "v"),
    #-----------------------------------------------------------------------
    transpose := self >> self,
    conjTranspose := self >> self,
    area := self >> self.n * self.n,
    #-----------------------------------------------------------------------
    toAMat := self >> AMatMat(MatSPL(DirectSum(O(1), J(self.n-1))))
));

#F VIxJ2(<v>) -- vector construct for I(v/2) \tensor J(2)
#F
Class(VIxJ2, BaseMat, SumsBase, rec(
    new := (self, v) >> SPL(WithBases(self, rec(v := v))).setDims(),
    dims := self >> [self.v, self.v],
    #-----------------------------------------------------------------------
    isReal := True,
    rChildren := self >> [self.v],
    rSetChild := rSetChildFields("v"),
    #-----------------------------------------------------------------------
    transpose := self >> self,
    conjTranspose := self >> self,
    area := self >> self.v * self.v,
    #-----------------------------------------------------------------------
    toAMat := self >> AMatMat(MatSPL(Tensor(I(self.v/2), J(2))))
));

#F For v > 2
#F    RCVIxJ2(v) == RC(VIxJ2(v/2))
#F               == I(v/4) \tensor J(2) \tensor I(2)
#F For v=2
#F    RCVIxJ2(2) == J(2)
#F
Class(RCVIxJ2, BaseMat, SumsBase, rec(
    new := (self, v) >> SPL(WithBases(self, rec(v := v))).setDims(),
    dims := self >> [self.v, self.v],
    #-----------------------------------------------------------------------
    isReal := True,
    rChildren := self >> [self.v],
    rSetChild := rSetChildFields("v"),
    #-----------------------------------------------------------------------
    transpose := self >> self,
    conjTranspose := self >> self,
    area := self >> self.v * self.v,
    #-----------------------------------------------------------------------
    toAMat := self >> Cond(self.v=2, J(2).toAMat(), # YSV: NOTE: check if this is correct
	                             Tensor(I(self.v/4), J(2), I(2)).toAMat())
));

#F VJxI(<m>, <v>) -- vector construct for J(m) \tensor I(v/m)
#F
Class(VJxI, BaseMat, SumsBase, rec(
    new := (self, m, v) >> SPL(WithBases(self, rec(m := m, v := v))).setDims(),
    dims := self >> [self.v, self.v],
    #-----------------------------------------------------------------------
    isReal := True,
    rChildren := self >> [self.m, self.v],
    rSetChild := rSetChildFields("m", "v"),
    #-----------------------------------------------------------------------
    transpose := self >> self,
    conjTranspose := self >> self,
    area := self >> self.v * self.v,
    #-----------------------------------------------------------------------
    toAMat := self >> AMatMat(MatSPL(Tensor(J(self.m), I(self.v/self.m))))
));