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

# helps translate VTensorInd(Scat(f), i) -> VScat_sv(f \tensor_i fId(i.range))
_fTensorInd := function(fnc, idx)
    local ds, irng, fdom, nidx, lambda;
    irng := idx.range;

    if not (idx in fnc.free()) then
        return fTensor(fnc, fId(irng));
    fi;

    fdom := fnc.domain();
    nidx := Ind(irng*fdom);

    lambda := Lambda(nidx, imod(nidx, irng)+V(irng)*Lambda(idx, fnc.at(idiv(nidx, irng))).at(imod(nidx, irng))).setRange(V(fnc.range())*V(irng));
    return lambda;
end;

#prevent VRC Rules from introducing VRCL/VRCR in the wrong spot
DontNeedSpecialVRC := e->not ForAny(e.children(), f->IsBound(f.needSpecialVRC) and f.needSpecialVRC());

# Test for VRCs that, for any construct, has a needSeparateVRC set
anyNeedSeparateVRC := e->ForAny(e.children(), f->IsBound(f.needSeparateVRC) and f.needSeparateVRC());

# Test for VRCs that for all constructs, cannot change data formats
TotallyNonChangeable := function(e)
    if ObjId(e)=Compose or ObjId(e)=ComposeDists then
      return(ForAll(e.children(), f->IsBound(f.totallyCannotChangeDataFormat) and f.totallyCannotChangeDataFormat()));
    else
      return(IsBound(e.totallyCannotChangeDataFormat) and e.totallyCannotChangeDataFormat());
    fi;
end;

# Test for VRCs that have at least one construct that cannot change data formats
HasNonchangeable := function(e)
    if ObjId(e)=Compose or ObjId(e)=ComposeDists then
      return(ForAny(e.children(), f->IsBound(f.cannotChangeDataFormat) and f.cannotChangeDataFormat()));
    else
      return(IsBound(e.cannotChangeDataFormat) and e.cannotChangeDataFormat());
    fi;
end;

NoSpecialOrNonchangeableVRC := e->(DontNeedSpecialVRC(e) and not HasNonchangeable(e));
NoSpecialAndChangeableVRC   := e->(DontNeedSpecialVRC(e) and     HasNonchangeable(e));

HandleCannotChangeDataFormatVRC := function(vrc, ch, v)
    local remch, ele, left, right, i, vrc1, vrc2;

    remch := Compose(Drop(ch, 1));
    #Error("BP: cannotChangeDataFormatVRC");

    if vrc=VRC then
        # NOTE: How do we handle this?
        # If ALL children cannot change data format, simply distribute the VRC amongst them

        if TotallyNonChangeable(Compose(ch)) then
            return( Compose(VRC(ch[1], v), VRC(remch,v)) );
        fi;

        # if any of the children NOTE: this is inefficient: what we really
        # should to is to bunch together the children that needSeparateVRC, and
        # the children that are regular, and handle them separately. For now,
        # this hack works because we should have no cases where such children
        # will be mixed together.

        if anyNeedSeparateVRC(Compose(ch)) then
            return( Compose(VRC(ch[1], v), VRC(remch,v)) );
        fi;

        # If not, we must do a VRCR/VRCL split
        # Find the rightmost child that can change dataformats. That and
        # everything right of that becomes a VRCL. Everything left of that is a VRCR.

        # Bad (but good:)) HACK: for things of the form S*A*G, where S and G
        # cannot change data formats, stuff the VRC into A, and hope that A is
        # big enough to handle it.

        if Length(ch)=3 and TotallyNonChangeable(ch[1]) and TotallyNonChangeable(ch[3]) then
          return(Compose( VRC(ch[1], v), VRC(ch[2], v), VRC(ch[3], v) ));
        fi;

        for i in Reversed([1..Length(ch)]) do
          if not TotallyNonChangeable(ch[i]) then
            # We found the rightmost child
            #Error("VRC: BP");
            left  := Compose(List([1..(i-1)], e->ch[e]));
            right := Compose(List([i..Length(ch)], e->ch[e]));
            #Error("VRC: BP");
            return( Compose( VRCR(left, v), VRCL(right, v) ) );
          fi;
        od;

        #Error("BP: cannotChangeDataformatVRC: VRC needs to be split into VRCL/VRCR");
    fi;

    if vrc=VRCLR then
        # This should really be handled as if the cannotChangeDataformat doesn't exist.
        #Error("VRCLR: BP");
        vrc1 := [[VRCL, VRCR], [VRCLR, VRCLR]];
        vrc2 := When(NeedInterleavedRight(ch[1]) or NeedInterleavedLeft(ch[2]), vrc1[1], vrc1[2]);
        #Error("VRCLR: BP");
        return( Compose(vrc2[1](ch[1], v), vrc2[2](Compose(Drop(ch, 1)), v)) );
    fi;

    if vrc=VRCL then
        # Find the first child from the right that is capable of making a
        # format change. We then have 2 cases.

        # Note: we really have 3 cases to be general: VRCL -> (VRCLR, VRCL, VRC)
        # Sometimes, the LR won't exist, and sometimes, the VRC won't exist
        # But it's easier to code if we only always break down to 2 cases

        for i in Reversed([1..Length(ch)]) do
          if not TotallyNonChangeable(ch[i]) then
            # We found the child that will do the dataformat change (VRCL)
            # We now have 2 cases.
            #   Case 1: VRCL -> (VRCL,  VRC)   (found child is not rightmost)
            #   Case 2: VRCL -> (VRCLR, VRCL)  (found child is rightmost)

            if i <> Length(ch) then # Case 1
              left  := Compose(List([1..i],              e->ch[e]));
              right := Compose(List([(i+1)..Length(ch)], e->ch[e]));
              #Error("BP:VRCL1\n, ---------VRCL(left)---------", left, "\n---------VRC(right)--------", right);
              return( Compose( VRCL(left, v), VRC(right, v) ) );
            else                    # Case 2
              left  := Compose(List([1..(i-1)], e->ch[e]));
              right := ch[i];
              #Error("BP:VRCL2\n, ---------VRCLR(left)---------", left, "\n---------VRCL(right)--------", right);
              return( Compose( VRCLR(left, v), VRCL(right, v) ) );
            fi;
          fi;
        od;
        Error("cannotChangeDataFormatVRC/L: Didn't find an appropriate child");
    fi;

    if vrc=VRCR then
        # Find the first child from the left that is capable of making a
        # format change. We then have 2 cases.

        for i in [1..Length(ch)] do
          if not TotallyNonChangeable(ch[i]) then
            # We found the child that will do the dataformat change (VRCR)
            # We now have 2 cases.
            #   Case 1: VRCR -> (VRC,  VRCR)   (found child is not leftmost)
            #   Case 2: VRCR -> (VRCR, VRCLR)  (found child is leftmost)

            if i <> 1 then # Case 1
              left  := Compose(List([1..(i-1)],      e->ch[e]));
              right := Compose(List([i..Length(ch)], e->ch[e]));
              #Error("BP:VRCR1\n, ---------VRC(left)---------", left, "\n---------VRCR(right)--------", right);
              return( Compose( VRC(left, v), VRCR(right, v) ) );
            else                    # Case 2
              left  := ch[i];
              right := Compose(List([(i+1)..Length(ch)], e->ch[e]));
              #Error("BP:VRCR2\n, ---------VRCR(left)---------", left, "\n---------VRCLR(right)--------", right);
              return( Compose( VRCR(left, v), VRCLR(right, v) ) );
            fi;
          fi;
        od;
        Error("cannotChangeDataFormatVRC/R: Didn't find an appropriate child");

    fi;

    Error("cannotChangeDataFormat: I don't know what to do with this one. None of the cases matched.");

#            vrc=VRCL and     HasNonchangeable(ch1) and not HasNonchangeable(ch2), Compose(VRCLR(ch1, v), VRCL(remch, v) ),
#            vrc=VRCL and not HasNonchangeable(ch1) and     HasNonchangeable(ch2), Compose(VRCL(ch1, v),  VRC(remch,v)   ),
#            vrc=VRCL and     HasNonchangeable(ch1) and     HasNonchangeable(ch2),
#                When(NeedInterleavedRight(ch1) or NeedInterleavedLeft(ch2), Compose(VRCL(I(ch1.dims()[1]),v),  VRC(ch,v)),
#                                                                            ComposeVRCLR(ch,v), VRCL(I(ch1.dims()[1]),v)),
#
#            vrc=VRCR and     HasNonchangeable(ch1) and not HasNonchangeable(ch2), Compose(VRC(ch1, v),   VRCR(remch, v) ),
#            vrc=VRCR and not HasNonchangeable(ch1) and     HasNonchangeable(ch2), Compose(VRCR(ch1, v),  VRCLR(remch,v) ),
#            vrc=VRCR and     HasNonchangeable(ch1) and     HasNonchangeable(ch2),
#                When(NeedInterleavedRight(ch1) or NeedInterleavedLeft(ch2), Compose(VRCL(I(ch1.dims()[1]),v),  VRC(ch,v)),
#                                                                            ComposeVRCLR(ch,v), VRCL(I(ch1.dims()[1]),v))
#          )
#        )
#    ),

end;

Scat.needInterleavedLeft:=False;
Scat.needInterleavedRight:=False;

# NOTE: Why were these set to False???
Gath.needInterleavedLeft:=True;
Gath.needInterleavedRight:=True;

VBlk.needInterleavedLeft:=False;
VBlk.needInterleavedRight:=False;

# NOTE: remove Diag here
Diag.needInterleavedLeft:=False;
Diag.needInterleavedRight:=False;

Inplace.needInterleavedRight := self >> self.child(1).needInterleavedRight();
Inplace.needInterleavedLeft := self >> self.child(1).needInterleavedLeft();
Inplace.cannotChangeDataFormat := self >> self.child(1).cannotChangeDataFormat();
Inplace.totallyCannotChangeDataFormat := self >> self.child(1).totallyCannotChangeDataFormat();

BB.needInterleavedRight := self >> self.child(1).needInterleavedRight();
BB.needInterleavedLeft := self >> self.child(1).needInterleavedLeft();
BB.cannotChangeDataFormat := self >> self.child(1).cannotChangeDataFormat();
BB.totallyCannotChangeDataFormat := self >> self.child(1).totallyCannotChangeDataFormat();

VContainer.needInterleavedRight := self >> self.child(1).needInterleavedRight();
VContainer.needInterleavedLeft := self >> self.child(1).needInterleavedLeft();
VContainer.cannotChangeDataFormat := self >> self.child(1).cannotChangeDataFormat();
VContainer.totallyCannotChangeDataFormat := self >> self.child(1).totallyCannotChangeDataFormat();

NoDiagPullin.needInterleavedRight := self >> self.child(1).needInterleavedRight();
NoDiagPullin.needInterleavedLeft := self >> self.child(1).needInterleavedLeft();
NoDiagPullin.cannotChangeDataFormat := self >> self.child(1).cannotChangeDataFormat();
NoDiagPullin.totallyCannotChangeDataFormat := self >> self.child(1).totallyCannotChangeDataFormat();

NoDiagPullinLeft.needInterleavedRight := self >> self.child(1).needInterleavedRight();
NoDiagPullinLeft.needInterleavedLeft := self >> self.child(1).needInterleavedLeft();
NoDiagPullinLeft.cannotChangeDataFormat := self >> self.child(1).cannotChangeDataFormat();
NoDiagPullinLeft.totallyCannotChangeDataFormat := self >> self.child(1).totallyCannotChangeDataFormat();

NoDiagPullinRight.needInterleavedRight := self >> self.child(1).needInterleavedRight();
NoDiagPullinRight.needInterleavedLeft := self >> self.child(1).needInterleavedLeft();
NoDiagPullinRight.cannotChangeDataFormat := self >> self.child(1).cannotChangeDataFormat();
NoDiagPullinRight.totallyCannotChangeDataFormat := self >> self.child(1).totallyCannotChangeDataFormat();

SymSPL.needInterleavedRight := self >> self.child(1).needInterleavedRight();
SymSPL.needInterleavedLeft := self >> self.child(1).needInterleavedLeft();
SymSPL.cannotChangeDataFormat := self >> self.child(1).cannotChangeDataFormat();
SymSPL.totallyCannotChangeDataFormat := self >> self.child(1).totallyCannotChangeDataFormat();

VGath_sv.rcVariant := RCVGath_sv;
VScat_sv.rcVariant := RCVScat_sv;

RCVGath_sv.vecVariant := (f,v)->VGath(f,v);
RCVScat_sv.vecVariant := (f,v)->VScat(f,v);

ScatGath.needInterleavedLeft:=True;
ScatGath.needInterleavedRight:=True;

SMPBarrier.needInterleavedLeft := (self) >> self.child(1).needInterleavedLeft();
SMPBarrier.needInterleavedRight := (self) >> self.child(1).needInterleavedRight();

_VRCFamily := [VRC, VRCL, VRCR, VRCLR];

Class(RulesVRC, RuleSet);
RewriteRules(RulesVRC, rec(
    VRC_VContainer := Rule([@(1, _VRCFamily), @(2, VContainer)],
	e -> let(cont := @(2).val,
	    VContainer(ObjId(@(1).val)(cont.child(1), @(1).val.v), cont.isa))),

    VRC_ISum := Rule([@(1, _VRCFamily), @(2, [ISum, SMPSum, SMPBarrier, SUM])], e->let(s := @(2).val,
    CopyFields(s, rec(_children := List(s.children(), c->ObjId(@(1).val)(c, @(1).val.v)),
                  dimensions := @(1).val.dimensions)))),

    VRC_Container := Rule([@(1, _VRCFamily),  @(2, [BB,Buf,Inplace,RecursStep,NoDiagPullin, NoDiagPullinLeft, NoDiagPullinRight])],
        e -> ObjId(@(2).val)(ObjId(@(1).val)(@(2).val.child(1), @(1).val.v))),

    VRC_Data := Rule([@(1, _VRCFamily), @(2, Data)],
	e -> Data(@(2).val.var, @(2).val.value, ObjId(e)(@(2).val.child(1)))),

#   VRC_SymSPL := Rule([@(1, _VRCFamily), @(2, SymSPL)],
#        e -> ObjId(@(1).val)(@(2).val.child(1), @(1).val.v)),

   VRC_Compose := Rule([@(1, _VRCFamily), @(2, [Compose,ComposeStreams]).cond(NoSpecialOrNonchangeableVRC)], # Use with VRC_ComposeNonchangeable below
        e->let(v := @(1).val.v, ch := @(2).val.children(), vrc := ObjId(@(1).val),
        vrc1 := Cond(
                 vrc=VRC,   [[VRC,  VRC ], [VRCR,  VRCL ]],
                 vrc=VRCLR, [[VRCL, VRCR], [VRCLR, VRCLR]],
                 vrc=VRCL,  [[VRCL, VRC ], [VRCLR, VRCL ]],
                 vrc=VRCR,  [[VRC,  VRCR], [VRCR,  VRCLR]]),
        vrc2 := When(NeedInterleavedRight(ch[1]) or NeedInterleavedLeft(ch[2]), vrc1[1], vrc1[2]),
            Compose(vrc2[1](ch[1], v), vrc2[2](Compose(Drop(ch, 1)), v)))),

    # NOTE: When building the format-change I() below, ch1.dims()[1] is
    #        always used. This is probably incorrect in some cases.
    VRC_ComposeNonchangeable := Rule([@(1, _VRCFamily), @(2, Compose).cond(NoSpecialAndChangeableVRC)],
        e->let(v     := @(1).val.v,
               ch    := @(2).val.children(),
               vrc   := ObjId(@(1).val),
               HandleCannotChangeDataFormatVRC(vrc, ch, v)
        )
    ),

    #VRC_ComposeStreams := Rule([@(1, VRC), @(2, ComposeStreams)],
    #    e->let(v := @(1).val.v, ch := @(2).val.children(),
    #        ComposeStreams( VRC(ch[1], v), VRC(ComposeStreams(Drop(ch,1)), v) )
    #        )
    #),

    VRCLR_VRCDiag:= Rule([@(1, _VRCFamily), @(2, VRCDiag)], e -> ObjId(@(1).val)(@(2).val.toloop().sums(), @(1).val.v)),

    VRCLR_VBlk:= Rule([@(1, VRCLR),@(2, VBlk)], e -> VBlk(RealVMatComplexVMat(@(2).val.element), @(2).val.v)),

    VRCLR_VGath_zero:= Rule([@(1, VRCLR),@(2, VGath_zero)], e -> let(g:=@(2).val,
        VRCLR(VGath(fAdd(g.N, g.n, 0), getV(g)), getV(g)))),

    VRCLR_VGathScat := Rule([@(1, VRCLR), @(2, [VGath, VScat, VScatAcc])],
        e -> ObjId(@(2).val)(fTensor(@(2).val.func, fId(2)), @(2).val.v)),

    VRCLR_GathScat := Rule([@(1, VRCLR), @(2, [Gath, Scat, ScatAcc])],
        e -> ObjId(@(2).val)(fTensor(@(2).val.func, fId(2)))),

    VRC_VGathVScat := Rule([@(1, VRC), @(2, [VGath, VScat, VScatAcc])],
        e -> ObjId(@(2).val)(fTensor(@(2).val.func, fId(2)), @(2).val.v)),

    VRCLR_VGathScat_sv := Rule([@(1, VRCLR), @(2, [VGath_sv, VScat_sv])], e->  # NOTE: acc variant missing
        @(2).val.rcVariant(@(2).val.func, @(2).val.v, @(2).val.sv, @(2).val.rem)),

    VRCLR_VTensor := Rule([@(1, VRCLR), @(2, VTensor)],
        e -> VTensor(RC(@(2).val.child(1)), @(2).val.vlen)),

    VRCLR_Diag := Rule([@(1, VRCLR), @(2, [Diag, VDiag])],
        e -> let(d := @(2).val.element, n := d.domain(), v := @(1).val.v,
            VRCDiag(VData(fCompose(RCData(d), fTensor(fId(n / v), L(2*v, 2))), v), v))),

    VRC_BlockVPerm := Rule([@(1, _VRCFamily), @(2, [BlockVPerm, BlockVPerm2])],
        e->let(bd := @(2).val.child(1), vrc := ObjId(@(1).val),
            BlockVPerm(@2.val.n, @2.val.vlen, vrc(bd, @(1).val.v),
            MatSPL(vrc(@2.val.perm, @1.val.v))))),

    fPrecompute_VData := Rule([VData, [fPrecompute, @(1)], @(2)], e -> fPrecompute(VData(@(1).val, @(2).val))),
    fPrecompute_VDup := Rule([VDup, [fPrecompute, @(1)], @(2)], e -> fPrecompute(VDup(@(1).val, @(2).val))),

    RCVScat_sv__fId := Rule(@(1,RCVScat_sv,e->IsInt(2*e.func.domain()/getV(e)) and ObjId(e.func) = fId),
        e->let(v:=getV(@(1).val), When(Cols(@(1).val)<>2*@(1).val.func.domain(),
            VGath_zero(Cols(@(1).val)/v, 2*@(1).val.func.domain()/v, v), VGath(fId(2*@(1).val.func.domain()/v), v))))
));

Class(RulesVRCTermDiag, RuleSet);
RewriteRules(RulesVRCTermDiag, rec(
    VRCLR_VDiag_x_I := Rule([@(1, VRCLR), @(2, VDiag_x_I)], e->VTensor(RC(Diag(@(2).val.element)), @(2).val.v))
));

Class(RulesVRCTerm, RulesVRCTermDiag);
RewriteRules(RulesVRCTerm, rec(
    VIxL := Rule(VIxL, (e, cx) -> e.implement(cx.opts.vector.isa)),
    VL := Rule(VL, (e, cx) -> e.implement(cx.opts.vector.isa)),
#--
### Rules for SAR -- to be moved and fixed...

    Term_VTensorInd_ScatGath := Rule([@(1, VTensorInd), @(2, ScatGath), @(3)],
        e->ScatGath(_fTensorInd(@(2).val.sfunc, @(3).val), _fTensorInd(@(2).val.gfunc, @(3).val))),

    RC_ScatGath := Rule([@@(1, RC), @(2, ScatGath)], (e,cx)->@(2).val.toloopRCVec(@(2).val.maxBkSize(), Last(cx.VContainer).isa.getV())),
##    VRCLR_ScatGath := Rule([@(1, VRCLR), @(2, ScatGath)], e->@(2).val.toloopRCVec(@(2).val.maxBkSize(), @(1).val.v)),

    RCVGathRCSVcat_sv_fTensor := Rule([@(1, [RCVGath_sv, RCVScat_sv]),
                            [@(2,fTensor), ..., [fId, @(3).cond(e -> Gcd(@(1).val.v/@(1).val.sv, EvalScalar(e)) = @(1).val.v)]]],
     e -> let(v := @(1).val.v, sv := @(1).val.sv,
          n := EvalScalar(@(3).val),   gcd := Gcd(v / sv, n),
          @(1).val.vecVariant(fTensor(DropLast(@(2).val.children(), 1), fId(2*n/gcd)), v))),

#    Id_RCV1 := Rule([@(1,[RCVScat_sv,RCVGath_sv]), @(2,fId,e->(IsValue(e.n) or IsInt(e.n)) and IsInt(EvalScalar(e.n/@(1).val.v)))], e->I(EvalScalar(@(2).val.n))),
#
#    Term__VScat_ScatGath := ARule(Compose, [@(1, VScat), @@(2, ScatGath, (e,cx)->not ForAny(_VRCFamily, i->IsBound(cx.(i.name)) and cx.(i.name) <> []))],
#        e->[@(1).val * @@(2).val.toloopVec(@@(2).val.maxBkSize(), @(1).val.v)]),
#
#    # in the precond of this rule there was a guard on VRCXX, but that was never matching due to an error. so I dont know if it was needed at all...
#    Term__ScatGath := Rule(@@(1, ScatGath, (e,cx)->not ForAny([VTensor, VTensorInd], i->IsBound(cx.(i.name)) and cx.(i.name) <> []) and
#        (IsBound(cx.VContainer) and cx.VContainer <> [])),
#        (e,cx)->@@(1).val.toloopVec(@@(1).val.maxBkSize(), Last(cx.VContainer).isa.getV())),
#
    Term_NeedInterleavedComplex_VRC := Rule([@(1, VRC), @(2, NeedInterleavedComplex)], e->RC(@(2).val.child(1))),
#    Term_NeedInterleavedComplex_RC := Rule([@(1, RC), @(2, NeedInterleavedComplex)], e->DRC(@(2).val.child(1))),
#
#    RC_GathScat_sv_fTensor := Rule([@(1, [VGath_sv, VScat_sv]), [@(2,fTensor), ..., [fId, 2]]],
#        e -> let(v := @(1).val.v, sv := @(1).val.sv,
#            @(1).val.rcVariant(fTensor(DropLast(@(2).val.children(), 1)), v, sv))),
#
    RC_XXX := Rule([@(1, RC), @(2, [VGath, VScat, VScatAcc])], e->ObjId(@(2).val)(fTensor(@(2).val.func, fId(2)), @(2).val.v)),
#    RC_XXX_sv := Rule([@(1, RC), @(2, [VGath_sv, VScat_sv])], e->@(2).val.rcVariant(@(2).val.func, @(2).val.v, @(2).val.sv)),
### end SAR
#--
#   these rules are buggy, as Scat_sv does zero padding.
#    RCVScat_sv_toVScat := Rule(@(1, RCVScat_sv, e->e.sv=1 and e.v=2), e -> VScat(e.func, e.v)),
#    RCVGath_sv_toVGath := Rule(@(1, RCVGath_sv, e->e.sv=1 and e.v=2), e -> VGath(e.func, e.v)),

#-- VPrm_x_I --

    #NOTE: Check if this rule is correct!
    #VRC_VPrm_x_I := Rule([@(1, VRC), @(2, VPrm_x_I)], e->let(p := @(2).val, v := p.v, b := p.dims()[1]/v,
    #    Compose(
    #        Tensor(I(2), VTensor(Prm(p.func), v)).sums().unroll()
    #    ))),

    # YSV: we use NoDiagPullin to prevent diagonals from going into the loops resulting from Tensor(I(2),...)
    #      because in some cases they can't be sucked in completely and one ends up with multiple Scat * Diag * Scat
    #      sequences inside a SUM, which overlap, implying that it is a SUMAcc. However, Spiral does not know
    #      that, and generates invalid code. I don't know of a better way to handle this at the moment.
    VRCLR_VPrm_x_I := Rule([@(1, [VRCLR, VRC]), @(2, VPrm_x_I)], e->let(p := @(2).val, v := p.v, b := p.dims()[1]/v,
        Compose(
            VTensor(Prm(L(2*b, b)), v),
            NoDiagPullin(Tensor(I(2), VTensor(Prm(p.func), v)).sums().unroll()),
            VTensor(Prm(L(2*b, 2)), v)
        ))),

    VRCL_VPrm_x_I := Rule([@(1, VRCL), @(2, VPrm_x_I)], e->let(p := @(2).val, v := p.v, b := p.dims()[1]/v,
        Compose(
            VIxL(b, 2, v),
            VTensor(Prm(L(2*b, b)), v),
            NoDiagPullin(Tensor(I(2), VTensor(Prm(p.func), v)).sums().unroll()),
            VTensor(Prm(L(2*b, 2)), v)
        ))),

    VRCR_VPrm_x_I := Rule([@(1, VRCR), @(2, VPrm_x_I)], e->let(p := @(2).val, v := p.v, b := p.dims()[1]/v,
        Compose(
            VTensor(Prm(L(2*b, b)), v),
            NoDiagPullin(Tensor(I(2), VTensor(Prm(p.func), v)).sums().unroll()),
            VTensor(Prm(L(2*b, 2)), v),
            VIxL(b, v, v)
        ))),

#-- VGath ------------------
    VRCL_VGath := Rule([@(1, VRCL), @(2, VGath)], e->let(v := @(2).val.v,
        Compose(
            VIxL(@(2).val.func.domain(), 2, v),
            VGath(fTensor(@(2).val.func, fId(2)), v)
        ))),

    VRCL_VScat := Rule([@(1, VRCL), @(2, [VScat, VScatAcc])], e->let(v := @(2).val.v,
        Compose(
            ObjId(@(2).val)(fTensor(@(2).val.func, fId(2)), v),
            VIxL(@(2).val.func.domain(), 2, v)
        ))),

    VRCL_VGath_sv := Rule([@(1, VRCL), @(2, VGath_sv)], e->let(v := @(2).val.v, sv := @(2).val.sv,
        Compose(
            VIxL(Rows(@(2).val)/v, 2, v),
            RCVGath_sv(@(2).val.func, v, sv, @(2).val.rem)
        ))),

    VRCL_IxVGath_pc := Rule([@(1, VRCL), @(2, IxVGath_pc)], e -> let(
	g := @(2).val, v := g.v,
        Compose(
            VIxL(g.k * _roundup(g.n, g.v) / v, 2, v),
            IxRCVGath_pc(g.k, g.N, g.n, g.ofs, v)
        ))),

    VRCLR_IxVGath_pc := Rule([@(1, VRCLR), @(2, IxVGath_pc, x->x.N*x.k mod x.v=0)], e -> let(
	g := @(2).val, v := g.v,
        Compose(
            VIxL(g.k * _roundup(g.n, v) / v, 2, v),
            IxRCVGath_pc(g.k, g.N, g.n, g.ofs, v),
            VIxL(g.k * g.N / v, v, v)
        ))),

    VRCL_VStretchGath := Rule([@(1, VRCL), @(2, VStretchGath)], e -> let(
	g := @(2).val,
	v := g.v,
	#XXX  Cond(???,
            Compose(
		VIxL(_roundup(Rows(g), v) / v, 2, v),
		RCVStretchGath(g.func, g.part, v))

	#XXX    VStretchGath(
	#XXX	fCompose(fTensor(fId(g.func.domain()/v), L(2*v, v)), fTensor(g.func, fId(2))),
	#XXX	...,
	#XXX	...)

        )),

    VRCR_VGath_zero:= Rule([@(1, VRCR),@(2, VGath_zero)], e -> let(g:=@(2).val,
        VGath_zero(2*g.N, 2*g.n, g.v) * VIxL(g.N, g.v, g.v))),

#-- VScat ------------------
    VRCR_VScat := Rule([@(1, VRCR), @(2, [VScat, VScatAcc])], e->let(v := @(2).val.v,
        Compose(
            ObjId(@(2).val)(fTensor(@(2).val.func, fId(2)), v),
            VIxL(@(2).val.func.domain(), v, v)
        ))),

    VRCR_VScat_sv := Rule([@(1, VRCR), @(2, VScat_sv)], e->let(v := @(2).val.v, sv := @(2).val.sv, # NOTE: acc variant missing
        Compose(
            RCVScat_sv(@(2).val.func, v, sv, @(2).val.rem),
            VIxL(Cols(@(2).val)/v, v, v)
        ))),

    VRCR_IxVScat_pc := Rule([@(1, VRCR), @(2, IxVScat_pc)], e->let(s:=@(2).val, v := s.v, # NOTE: acc variant missing
        Compose(
            IxRCVScat_pc(s.k, s.N, s.n, s.ofs, v),
            VIxL(s.k*_roundup(s.n, v)/v, v, v)
        ))),

    VRCLR_IxVScat_pc := Rule([@(1, VRCLR), @(2, IxVScat_pc, x->x.N*x.k mod x.v=0)], e -> let(
        s := @(2).val, v := s.v,
        Compose(
            VIxL(s.k * s.N / v, 2, v),
            IxRCVScat_pc(s.k, s.N, s.n, s.ofs, v),
            VIxL(s.k*_roundup(s.n, v)/v, v, v)
        ))),

    VRCR_VStretchScat := Rule([@(1, VRCR), @(2, VStretchScat)], e->let(s:=@(2).val, v := s.v,
        Compose(
            RCVStretchScat(s.func, s.part, v),
            VIxL(_roundup(Cols(s), v)/v, v, v)
        ))),

#----------------------------
    VRCLR_VPerm := Rule([@(1, VRCLR), @(2, VPerm)], e->let(p := @(2).val, v := p.vlen, b := p.dims()[1]/v,
        Compose(
            VTensor(Prm(L(2*b, b)), v),
            Tensor(I(2), p).sums().unroll(),
            VTensor(Prm(L(2*b, 2)), v)
        ))),

    VRCL_VPerm := Rule([@(1, VRCL), @(2, VPerm)], e->let(p := @(2).val, v := p.vlen, b := p.dims()[1]/v,
        Compose(
            VTensor(Prm(L(2*b, b)), v),
            Tensor(I(2), p).sums().unroll(),
            VTensor(Prm(L(2*b, 2)), v),
            VIxL(b, 2, v)
        ))),

    VRCR_VPerm := Rule([@(1, VRCR), @(2, VPerm)], e->let(p := @(2).val, v := p.vlen, b := p.dims()[1]/v,
        Compose(
            VIxL(b, v, v),
            VTensor(Prm(L(2*b, b)), v),
            Tensor(I(2), p).sums().unroll(),
            VTensor(Prm(L(2*b, 2)), v)
        ))),

#----------------------------
    VRCLR_Perm := Rule([@(1, VRCLR), @(2, Prm)], e->let(p := @(2).val, v := getV(@1.val), b := p.dims()[1]/v, n := @(2).val.dims()[2]/v,
        Compose(
            VIxL(n, 2, v),
            RCVScat_sv(fId(p.func.domain()), v, 1),
            RCVGath_sv(p.func, v, 1),
            VIxL(n, v, v)
        ))),

    VRC_Perm := Rule([@(1, VRC), @(2, Prm)], e->let(p := @(2).val, v := getV(@1.val), b := p.dims()[1]/v, n := @(2).val.dims()[2]/v,
            Compose (RCVScat_sv(fId(p.func.domain()), v, 1),
            RCVGath_sv(p.func, v, 1))
        )),

    VRCL_Split := Rule([@(1, VRCL), @(2, [VDiag, VTensor])], e->let(v := getV(@(2).val),
        Compose(
            VRCLR(@(2).val, v),
            VIxL(@(2).val.dims()[2]/v, 2, v)
        ))),

    VRCR_Split := Rule([@(1, VRCR), @(2, [VDiag, VTensor])], e->let(v := getV(@(2).val),
        Compose(
            VIxL(@(2).val.dims()[1]/v, v, v),
            VRCLR(@(2).val, v)
        ))),

    VRC_VTensor := Rule([@(1, VRC), @(2, [VTensor])], e->let(v := getV(@(2).val),
        Compose(
            VIxL(@(2).val.dims()[1]/v, v, v),
            VRCLR(@(2).val, v),
            VIxL(@(2).val.dims()[2]/v, 2, v)
        ))),

    VRC_VDiag := Rule([@(1, VRC), @(2, [VDiag])], e->let(
       v := getV(@(2).val), d := @(2).val.dims()[1],
       Cond(@(2).val.element.isReal(),
	   VDiag(diagTensor(@(2).val.element, fConst(TReal, 2, 1)), @(2).val.v),
	   VIxL(d/v, v, v) * VRCLR(@(2).val, v) * VIxL(d/v, 2, v)))),

    VRC_Blk := Rule([@(1, VRC), @(2, Blk)], e->RC(@(2).val)),

    VRCLR_VScat_zero := Rule([@(1, VRCLR), @(2, VScat_zero)], e->let(s:=@(2).val, VScat_zero(2*s.N, 2*s.n, s.v))),

    VRC_Gath := Rule([@(1, VRC), @(2, Gath)], e->let(g := @(2).val, v := getV(@1.val), b := g.dims()[1]/v, n := @(2).val.dims()[2]/v,
            Compose (VRCLR(VScat_sv(fId(g.func.domain()), v, 1), v),
            VRCLR(VGath_sv(g.func, v, 1), v))
        )),

    FormatPrm_Term := Rule([FormatPrm, @(1)], e -> Prm(@(1).val)),
));

Class(RulesSplitComplex, RuleSet);
RewriteRules(RulesSplitComplex, rec(
    VRCR_IxVScat_pc := ARule(Compose, [[@(1,[Prm, FormatPrm]), @(4,L,e->e.params[2]=2)], [@(2,VRCR), @(3,IxVScat_pc)]],
                e -> let(v:=@(2).val.v, c:=Cols(@(3).val), s:=@(3).val,
                    [IxVScat_pc(2*s.k, s.N, s.n, s.ofs, v), VPrm_x_I(L(2*c/v, 2), v)])),

    VRCL_IxVGath_pc := ARule(Compose, [[@(2,VRCL), @(3,IxVGath_pc)], [@(1,[Prm, FormatPrm]), @(4,L,e->e.params[2]=e.params[1]/2)]],
                e -> let(v:=@(2).val.v, r:=Rows(@(3).val), g:=@(3).val,
                    [VPrm_x_I(L(2*r/v, r/v), v), IxVGath_pc(2*g.k, g.N, g.n, g.ofs, v)])),

    VRCR_VStretchScat := ARule(Compose, [[@(1,[FormatPrm,Prm]), @(4,L,e->e.params[2]=2)], [@(2,VRCR), @(3,VStretchScat)]],
                e -> let(v:=@(2).val.v, c:=Cols(@(3).val), s:=@(3).val,
                    [VStretchScat(fTensor(fId(2), s.func), 2*s.part, v), VPrm_x_I(L(2*c/v, 2), v)])),

    VRCL_VStretchGath := ARule(Compose, [[@(2,VRCL), @(3,VStretchGath)],[@(1,[FormatPrm,Prm]), @(4,L,e->e.params[2]=e.params[1]/2)]],
                e -> let(v:=@(2).val.v, r:=Rows(@(3).val), g:=@(3).val,
                    [VPrm_x_I(L(2*r/v, r/v), v), VStretchGath(fTensor(fId(2), g.func), 2*g.part, v)])),

    VRCR_VScat_sv := ARule(Compose, [[@(1,[FormatPrm,Prm]), @(4,L,e->e.params[2]=2)], [@(2,VRCR), @(3,VScat_sv)]],
                e -> let(v:=@(2).val.v, c:=Cols(@(3).val), s:=@(3).val,
                    [VStretchScat(fTensor(fId(2), s.func), 2, v), VPrm_x_I(L(2*c/v, 2), v)])),

    VRCL_VGath_sv := ARule(Compose, [[@(2,VRCL), @(3,VGath_sv)],[@(1,[FormatPrm,Prm]), @(4,L,e->e.params[2]=e.params[1]/2)]],
                e -> let(v:=@(2).val.v, r:=Rows(@(3).val), g:=@(3).val,
                    [VPrm_x_I(L(2*r/v, r/v), v), VStretchGath(fTensor(fId(2), g.func), 2, v)])),

    VRCR_VScat := ARule(Compose, [[@(1,[Prm, FormatPrm]), @(4,L,e->e.params[2]=2)], [@(2,VRCR), @(3,[VScat, VScatAcc])]],
                e -> let(v:=@(2).val.v, c:=Cols(@(3).val), s:=@(3).val,
                    [ObjId(@(3).val)(fTensor(fId(2), s.func), v), VPrm_x_I(L(2*c/v, 2), v)])),

    VRCL_VGath := ARule(Compose, [[@(2,VRCL), @(3,VGath)],[@(1,[Prm, FormatPrm]), @(4,L,e->e.params[2]=e.params[1]/2)]],
                e -> let(v:=@(2).val.v, r:=Rows(@(3).val), g:=@(3).val,
                    [VPrm_x_I(L(2*r/v, r/v), v), VGath(fTensor(fId(2), g.func), v)])),

    VRCL_VScat := ARule(Compose, [[@(2,VRCL), @(3,[VScat, VScatAcc])],[@(1,Prm), @(4,L,e->e.params[2]=e.params[1]/2)]],
                e -> let(v:=@(2).val.v, r:=Rows(@(3).val), g:=@(3).val,
                    [VPrm_x_I(L(2*r/v, r/v), v), ObjId(@(3).val)(fTensor(fId(2), g.func), v)])),

    VRCR_VTensor1 := ARule(Compose, [[@(1, [FormatPrm, Prm]), @(4,L,e->e.params[2]=2)], [@(2,VRCR), @(3,VTensor)]],
                e -> let(v:=@(2).val.v, c:=Cols(@(3).val), s:=@(3).val,
                    [VPrm_x_I(L(2*c/v, 2), v), VRCLR(s, v)])),

    VRCR_VTensor1a := ARule(Compose, [[@(1, [FormatPrm, Prm]), @(4,L,e->e.params[2]=2)], [@(2,VRC), @(3,VTensor)]],
                e -> let(v:=@(2).val.v, c:=Cols(@(3).val), s:=@(3).val,
                    [VPrm_x_I(L(2*c/v, 2), v), VRCL(s, v)])),

    VRCL_VTensor2 := ARule(Compose, [[@(2,VRCL), @(3,VTensor)],[@(1, [FormatPrm, Prm]), @(4,L,e->e.params[2]=e.params[1]/2)]],
                e -> let(v:=@(2).val.v, c:=Cols(@(3).val), g:=@(3).val,
                    [VRCLR(g,v), VPrm_x_I(L(2*c/v, c/v), v)])),

    VRCL_VTensor2a := ARule(Compose, [[@(2,VRC), @(3,VTensor)],[@(1, [FormatPrm, Prm]), @(4,L,e->e.params[2]=e.params[1]/2)]],
                e -> let(v:=@(2).val.v, r:=Rows(@(3).val), g:=@(3).val,
                    [VRCR(g,v), VPrm_x_I(L(2*r/v, r/v), v)])),

#    VRCLR_VGath_split := Rule([@(1, VRCLR), @(2, VGath_sv)],
#                e -> let(g := @(2).val, v:= g.v, func := g.func, n := func.domain(), N:= func.range(),
#                        VGath_sv(fCompose(fTensor(L(2*N/v,2), fId(v)), fTensor(fId(2), func), fTensor(L(2*n/v,n/v), fId(v))), v, g.sv))),
#
#    VRCLR_VScat_split := Rule([@(1, VRCLR), @(2, VScat_sv)],
#                e -> let(g := @(2).val, v:= g.v, func := g.func, n := func.domain(), N:= func.range(),
#                        VScat_sv(fCompose(fTensor(L(2*N/v,2), fId(v)), fTensor(fId(2), func), fTensor(L(2*n/v,n/v), fId(v))), v, g.sv)))
#
    VGath_FormatPrm := ARule(Compose, [[@(1, VGath), [fTensor, fBase, @(2, fId)]],
            [@(3, FormatPrm), [fTensor, fId, @(4, L, e->IsInt(@(2).val.domain()/e.domain()))]]],
        e->[FormatPrm(fTensor(fId(Rows(@(1).val)/@(4).val.domain()), @(4).val)), @(1).val]),

    FormatPrm_VScat := ARule(Compose, [[@(1, FormatPrm), [fTensor, fId, @(2, L)]],
            [@(3, [VScat, VScatAcc]), [fTensor, fBase, @(4, fId, e->IsInt(e.domain()/@(2).val.domain()))]]],
        e->[@(3).val, FormatPrm(fTensor(fId(Cols(@(3).val)/@(2).val.domain()), @(2).val))]),

    VRCL_VGath_fIdxL := ARule(Compose, [[@(2,VRCL), @(3,VGath)], [@(1,FormatPrm), [@(4,fTensor), fId, @(5, L, e->e.params[1]=2*e.params[2])]]],
        e -> let(v:=@(2).val.v, r:=Rows(@(3).val), g:=@(3).val, [VGath(fTensor(g.func, fId(2)), v)])),

    VRCR_VScat_fIdxL := ARule(Compose, [[@(1, FormatPrm), [@(4,fTensor), fId, @(5, L, e->e.params[2]=2)]], [@(2,VRCR), @(3,[VScat, VScatAcc])]],
        e -> let(v:=@(2).val.v, c:=Cols(@(3).val), s:=@(3).val, [ObjId(@(3).val)(fTensor(s.func, fId(2)), v)]))

));

Class(RulesVBlkInt, RuleSet);
RewriteRules(RulesVBlkInt, rec(
    Merge_RulesVBlkInt_VRC_rt := ARule(Compose, [@(1, VBlkInt), @(2, VRC)],
                e -> let(v1 := @(1).val, v2 := @(2).val, v := v1.v, [ VIxL(Rows(v1)/(2*v), v, v), v1.child(1), VRCL(v2.child(1), v) ])),

    Merge_RulesVBlkInt_VRC_lft := ARule(Compose, [@(2, VRC), @(1, VBlkInt)],
                e -> let(v1 := @(1).val, v2 := @(2).val, v := v1.v, [ VRCR(v2.child(1), v), v1.child(1), VIxL(Rows(v1)/(2*v), 2, v) ])),

    Merge_RulesVBlkInt_vRC_rt := ARule(Compose, [@(1, VBlkInt), @(2, vRC)],
                e -> let(v1 := @(1).val, v := v1.v, [ VIxL(Rows(v1)/(2*v), v, v), v1.child(1), VIxL(Cols(v1)/(2*v), 2, v), @(2).val ])),

    Merge_RulesVBlkInt_vRC_lft := ARule(Compose, [@(2, vRC), @(1, VBlkInt)],
                e -> let(v1 := @(1).val, v := v1.v, [ @(2).val, VIxL(Rows(v1)/(2*v), v, v), v1.child(1), VIxL(Rows(v1)/(2*v), 2, v) ])),
));


text:
InstallMethod(AlgebraicU,
              "Simple algebraic of adjoint type",
              [IsString,IsPosInt,IsRing],
              function(Type,Rank,Ring)
              local object,
              avarnames;
              object:=Objectify(NewType(NewFamily("AlgebraicUFamily"),
                                        IsAttributeStoringRep and
                                        IsChevalleyAdj and
                                        IsAlgebraicU),
                                rec());

              Settype(object,Type);
              Setrank(object,Rank);

              SetchevalleyAdj(object,ChevalleyAdj(Type,Rank,Ring));

              SetBaseRing(object,Ring);
              avarnames:=List([1..2*Length(positiveRoots(chevalleyAdj(object)))],i->Concatenation("a_{",String(i),"}"));
              Setring(object,PolynomialRing(Ring,avarnames));

              SetCharacteristic(object,Characteristic(Ring));

              SetlieAlgebra(object,SimpleLieAlgebraTypeA_G(Type,Rank,ring(object)));
              SetrootSystem(object,RootSystem(SimpleLieAlgebra(Type,Rank,Integers)));
              SetpositiveRoots(object,positiveRoots(rootSystem(object)));
              SetallRoots(object,Concatenation(
                  positiveRoots(rootSystem(object)),
                  -positiveRoots(rootSystem(object))));
              SetweylGroup(object,WeylGroup(rootSystem(object)));

              SetA(object,A_rs(object));
              SetH(object,Eta(object));
              SetN(object,N_rs(object));
              SetM(object,M_rsi(object));
              SetC(object,C_ijrs(object));

              SetName(object,Concatenation("<simple adjoint ",Type,String(Rank),
                                           " in characteristic ",String(Characteristic(Ring)),">"));

              return object;
end);

InstallMethod(AlgebraicU,
              "Simple algebraic of adjoint type",
              [IsChevalleyAdj],
              function(sys)
              local object,
              avarnames;
              object:=Objectify(NewType(NewFamily("AlgebraicUFamily"),
                                        IsAttributeStoringRep and
                                        IsChevalleyAdj and
                                        IsAlgebraicU),
                                rec());

              Settype(object,type(sys));
              Setrank(object,rank(sys));

              avarnames:=List([1..2*Length(positiveRoots(sys))],i->Concatenation("a_{",String(i),"}"));
              Setring(object,PolynomialRing(ring(sys),avarnames));

              SetCharacteristic(object,Characteristic(sys));

              SetlieAlgebra(object,SimpleLieAlgebraTypeA_G(type(object),rank(object),ring(object)));
              SetrootSystem(object,rootSystem(sys));
              SetpositiveRoots(object,positiveRoots(sys));
              SetallRoots(object,allRoots(sys));
              SetweylGroup(object,weylGroup(sys));

              SetA(object,A(sys));
              SetH(object,H(sys));
              SetN(object,N(sys));
              SetM(object,M(sys));
              SetC(object,C(sys));

              SetBaseRing(object,ring(sys));
              SetchevalleyAdj(object,sys);

              SetName(object,Concatenation("<simple adjoint ",type(object),String(rank(object)),
                                           " in characteristic ",String(Characteristic(object)),">"));

              return object;
end);

ext: gi


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

FindUnexpandableNonterminal := function(t, opts)
    local trees, res;
    Constraint(IsSPL(t));
    trees := ExpandSPL(t, opts);
    if trees=[] then return [t];
    else
        return ConcatList(trees, tr ->
	    ConcatList(tr.children, c -> FindUnexpandableNonterminal(c, opts)));
    fi;
end;

RuleTreeClass_rChildren := self >> [self.node, self.children];
RuleTreeClass_rSetChild := rSetChildFields("node", "children");

EnableRuleTreeRewriting := function(switch)
    if switch then
        RuleTreeClass.rChildren := RuleTreeClass_rChildren;
        RuleTreeClass.rSetChild := RuleTreeClass_rSetChild;
    else
        Unbind(RuleTreeClass.rChildren);
        Unbind(RuleTreeClass.rSetChild);
    fi;
end;

VerifyRulesInRuleTree := function(r, opts)
    local nt, n, bad;
    EnableRuleTreeRewriting(true);
    nt := CollectNR(r, @.cond(IsNonTerminal));
    bad := [];
    for n in nt do
        #if n.free() <> [] then
            n:=HashAsSPL(n);
        #fi;
        if ObjId(n) <> InfoNt then
            Print(Red("-------", n, "-------"), "\n");
            if VerifyRulesForSPL(n, opts) = false then
                Add(bad, n);
            fi;
        fi;
    od;
    EnableRuleTreeRewriting(true);
    Print(Red("Failed:\n    "), bad, "\n");
end;

_VerifySubRuleTrees := function(r, opts, rt_tomat_func, ind)
    local subtrees, rt, bugs, hashrt, t, innerbugs;
    EnableRuleTreeRewriting(true);
    bugs := [];
    for rt in r.children do
        t := HashAsSPL(rt.node);
	hashrt := ApplyRuleTreeSPL(rt, t, opts);
        if InfinityNormMat(MatSPL(rt_tomat_func(hashrt)) - MatSPL(hashrt.node)) > 1e-13 then
	    PrintLine(Blanks(ind), hashrt.node, RedStr("  FAIL"));
	    bugs := [hashrt];
	    innerbugs := _VerifySubRuleTrees(hashrt, opts, rt_tomat_func, ind+4);
	    return Cond(innerbugs=[], bugs, innerbugs);
	else
	    PrintLine(Blanks(ind), hashrt.node, GreenStr("  OK"));
	fi;
    od;
    EnableRuleTreeRewriting(false);
    return bugs;
end;

VerifySubRuleTreesSPL := (r, opts) -> _VerifySubRuleTrees(r, opts, SPLRuleTree, 0);

#VerifySubRuleTreesSums := (r, opts) -> _VerifySubRuleTrees(r, opts, x->SumsRuleTree(x, opts), 0);


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

#F Perm_CCS(<n>) - SPL objects for CCS -> Perm IPP data format conversion
#F
#F Returns an SPL object that converts the result of an RDFT in CCS format into
#F Perm format. Both of these formats are defined in the IPP manual.
#F IPP manual does not define Perm format for odd RDFT size, so in this case
#F we defined it to coincide with CCS.
#F
#F The output of SPIRAL's PRDFT transform is equivalent to IPP CCS format.
#F Suggested uses are:
#F     Perm_CCS(n) * PRDFT(n)
#F     IPRDFT(n) * Perm_CCS(n).transpose()
#F
Perm_CCS := n -> When(IsEvenInt(n),
    Gath(H(n+2,n,0,1))*Perm((2,n+1),n+2),
    DirectSum(Mat([[1,0]]), I(n-1))
);

#F Pack_CCS(<n>) - SPL objects for CCS -> Pack IPP data format conversion of RDFT of type 1
#F
#F Returns an SPL object that converts the result of an RDFT in CCS format into
#F Pack format. Both of these formats are defined in the IPP manual.
#F
#F The output of SPIRAL's PRDFT transform is equivalent to IPP CCS format.
#F Suggested uses are:
#F     Pack_CCS(n) * PRDFT(n)
#F     IPRDFT(n) * Pack_CCS(n).transpose()
#F
Pack_CCS := n -> When(IsEvenInt(n),
    DirectSum(Mat([[1,0]]), When(n=2, [], I(n-2)), Mat([[1,0]])),
    DirectSum(Mat([[1,0]]), I(n-1))
);

#F Pack_CCS(<n>) - SPL objects for CCS -> Pack IPP data format conversion of RDFT of type 3
#F
#F [ note that IPP does not have RDFT of type 3, we just use their format terminology ]
#F Returns an SPL object that converts the result of an RDFT in CCS format into
#F Pack format. Both of these formats are defined in the IPP manual.
#F Perm format does not make sense for PRDFT3.
#F
#F The output of SPIRAL's PRDFT3 transform is equivalent to IPP CCS format.
#F Suggested uses are:
#F     Pack_CCS3(n) * PRDFT3(n)
#F     PRDFT3(n).inverse() * Pack_CCS3(n).transpose()
#F
Pack_CCS3 := n -> When(IsEvenInt(n),
    I(n),
    DirectSum(I(n-1), Mat([[1,0]]))
);


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

# the rule set to promote nonterminals
Class(RulesFFTXPromoteNT, RuleSet);

Class(RulesFFTXPromoteNT_Cleanup, RuleSet);

#RewriteRules(RulesFFTXPromoteNT, rec(
#    IPRDFT_RCDiag_PRDFT__Circulant_Rule := Rule([Compose, @(1,IPRDFT), [@(2,RCDiag), @(4, FDataOfs, e->e.ofs = 0), @(5,I)], @(3,PRDFT)],
#        e->let(n := @(1).val.params[1], fdata := @(2).val.element,
#            Circulant(n, FDataNT(fdata.var, @(1).val), -n))
#    )
#));

_toSymList := l -> [Minimum(l)..Maximum(l)];

RewriteRules(RulesFFTXPromoteNT, rec(
#   batch1_MDXFFT := Rule(@@(1, [MDDFT, MDPRDFT, IMDPRDFT, TColMajor],
#       (e,cx) -> ForAll([TTensorI, TColMajor, TDAG, TDAGNode, HStack, VStack, Compose], i -> (not IsBound(cx.(i.name)) or cx.(i.name) = []))),
#       e-> TTensorI(@@(1).val, 1, AVec, AVec)),

    Scat_Circulant_Gath__IOPrunedRConv := ARule(Compose, [[@(1, Gath), fAdd],  @(2,Circulant), [@(3, Scat), fAdd]],
        e-> [IOPrunedRConv(@(2).val.params[1],
            FDataOfs(@(2).val.params[2].var, 2*(@(2).val.params[1]/2+1), 0),
            1, _toSymList(List(@(1).val.func.tolist(), _unwrap)),
            1, _toSymList(List(@(3).val.func.tolist(), _unwrap)), true)]),

# This is one mega promotion rule for Hockney that needs to be broken apart after MDRConv and/or PrunedMDRDFT is introduced
    Hockney_hack := ARule(Compose, [[@(6, Gath), @(8,fTensor, e->ForAll(e.children(), i->ObjId(i)=fAdd))],
                                     @(1,IMDPRDFT, e -> e.params[2] = 1), [@(2,RCDiag), @(4, FDataOfs, e->e.ofs = 0), @(5,I)],
                                     @(3,MDPRDFT, e -> e.params[2] = Product(e.params[1])-1),
                                     [@(7, Scat), @(9,fTensor, e->ForAll(e.children(), i->ObjId(i)=fAdd))]],
        e-> let(ii := Ind(Rows(@(2).val)),
            sym := @(2).val.element.var,
            symf := Lambda(ii, nth(sym,ii)),
            opat := List(@(8).val.children(), i-> _toSymList(List(i.tolist(), _unwrap))),
            ipat := List(@(9).val.children(), i-> _toSymList(List(i.tolist(), _unwrap))),
        [ IOPrunedMDRConv(@(1).val.params[1], symf, 1, opat, 1, ipat, true) ])),

# This is one mega promotion rule for Hockney that needs to be broken apart after MDRConv and/or PrunedMDRDFT is introduced
    Hockney_hack2 := ARule(Compose, [[@(6, Gath), @(8,fTensor, e->ForAll(e.children(), i->ObjId(i)=fAdd))],
                                     @(1,IMDPRDFT, e -> e.params[2] = 1), [@(2,RCDiag), @(4, FData), @(5,I)],
                                     @(3,MDPRDFT, e -> e.params[2] = Product(e.params[1])-1),
                                     [@(7, Scat), @(9,fTensor, e->ForAll(e.children(), i->ObjId(i)=fAdd))]],
        e-> let(ii := Ind(Rows(@(2).val)),
            sym := @(2).val.element.var,
            symf := Lambda(ii, nth(sym,ii)),
            opat := List(@(8).val.children(), i-> _toSymList(List(i.tolist(), _unwrap))),
            ipat := List(@(9).val.children(), i-> _toSymList(List(i.tolist(), _unwrap))),
        [ IOPrunedMDRConv(@(1).val.params[1], symf, 1, opat, 1, ipat, true) ])),

# DAG to Compose rules
    DAG_collapse1 := Rule([@(1, TDAG), ..., @(2, TDAGNode, e->ForAny(Drop(@(1).val.params[1], 1), k->k.params[3] = e.params[2])), ...],
        e->  let(nbr := Filtered(@(1).val.params[1], k->k.params[3] = @(2).val.params[2])[1],
                 ch := Filtered(Drop(@(1).val.params[1],1) , i -> i<> nbr),
                 dn := TDAGNode(nbr.params[1] * @(2).val.params[1], nbr.params[2], @(2).val.params[3]),
                 TDAG([dn]::ch))),
    DAG_remove := Rule([@(1,TDAG), [@(2, TDAGNode), @(3),
            @(4).cond(e->IsList(e) and Length(e) = 1 and let(_e := e[1].free(), Length(_e) = 1 and ObjId(_e[1]) = var and _e[1].id in ["Y", "Yptr"] and Collect(e[1], TColMaj) = [])),
            @(5).cond(e->IsList(e) and Length(e) = 1 and let(_e := e[1].free(), Length(_e) = 1 and ObjId(_e[1]) = var and _e[1].id in ["X", "Xptr"] and Collect(e[1], TColMaj) = [])),
            ...],...],
        e-> @(2).val.params[1]),

    propagate_ColMaj := Rule([@(1,TDAGNode), [TRC, @(2,MDDFT),...],
            @(3).cond(e->ObjId(e[1])=tcast and ObjId(e[1].args[1].t) = TColMaj),
            @(4).cond(e->ObjId(e[1])=tcast and ObjId(e[1].args[1].t) = TColMaj),...],
        e -> TDAGNode(TRC(TColMajor(@(2).val)), @(3).val[1].args[2], @(4).val[1].args[2])),

    Drop_TDecl := Rule(@(1, TDecl, e->ForAll(e.params[2], i->not i in e.params[1].free())),
        e -> @(1).val.params[1]),

    Drop_TTensorI := Rule(@(1, TTensorI, e -> e.params[2] = 1), e -> e.params[1])
));

RewriteRules(RulesFFTXPromoteNT_Cleanup, rec(
# FIXME: the promotion rule does not (yet) have the guard to ensure the value of k is correct
    IPRDFT_RCDiag_PRDFT__Circulant_ARule := ARule(Compose, [@(1,IPRDFT, e -> e.params[2] = 1), [@(2,RCDiag),
            @(4, FDataOfs, e->e.ofs = 0), @(5,I)], @(3,PRDFT, e -> e.params[2] = e.params[1]-1)],
        e->let(n := @(1).val.params[1], fdata := @(2).val.element,
            [Circulant(n, FDataNT(fdata.var, @(1).val), -n) ])),

# FIXME: the promotion rule does not (yet) have the guard to ensure the value of k is correct
    IMDPRDFT_RCDiag_MDPRDFT__RConv_ARule := ARule(Compose, [@(1,IMDPRDFT, e -> e.params[2] = 1), [@(2,RCDiag), @(4, FDataOfs, e->e.ofs = 0), @(5,I)],
            @(3,MDPRDFT, e -> e.params[2] = Product(e.params[1])-1)],
        e->let(n := @(1).val.params[1], fdata := @(2).val.element,
            [MDRConv(n, fdata.var, true) ])),

# FIXME: the promotion rule does not (yet) have the guard to ensure the value of k is correct
    IMDPRDFT_Diag_MDPRDFT__RConv_ARule := ARule(Compose, [@(1,IMDPRDFT, e -> e.params[2] = 1), [@(2,Diag), [diagTensor, @(4, FDataOfs, e->e.ofs = 0), @(5,fConst, e->e.params = [ TReal, 2, 1 ])]],
            @(3,MDPRDFT, e -> e.params[2] = Product(e.params[1])-1)],
        e->let(n := @(1).val.params[1], fdata := @(4).val,
            [MDRConvR(n, fdata.var, true) ])),

# MDDFT * Scat -> PrunedMDDFT
    MDDFT_Scat__PrunedMDDFT := ARule(Compose, [@(1,MDDFT), [@(2, Scat), @(3,fTensor, e->ForAll(e.children(), i->ObjId(i)=fAdd))]],
        e -> [ PrunedMDDFT(@(1).val.params[1], @(1).val.params[2], 1,  List(@(3).val.children(), i-> _toSymList(List(i.tolist(), _unwrap))))] ),

# Gath * MDDFT -> PrunedIMDDFT
    Gath_MDDFT__PrunedIMDDFT := ARule(Compose, [[@(1, Gath), @(2,fTensor, e->ForAll(e.children(), i->ObjId(i)=fAdd))], @(3, MDDFT)],
        e -> [PrunedIMDDFT(@(3).val.params[1], @(3).val.params[2], 1,  List(@(2).val.children(), i-> _toSymList(List(i.tolist(), _unwrap))))]),

# PRDFT * Scat -> PrunedPRDFT
    PRDFT_Scat__PrunedPRDFT := ARule(Compose, [@(1,PRDFT), [@(2, Scat), @(3,fAdd)]],
        e -> [ PrunedPRDFT(@(1).val.params[1], @(1).val.params[2], 1,  _toSymList(List(@(3).val.tolist(), _unwrap)))] ),

# MDPRDFT * Scat -> PrunedMDPRDFT
    MDPRDFT_Scat__PrunedMDPRDFT := ARule(Compose, [@(1,MDPRDFT), [@(2, Scat), @(3,fTensor, e->ForAll(e.children(), i->ObjId(i)=fAdd))]],
        e -> [ PrunedMDPRDFT(@(1).val.params[1], List(@(3).val.children(), i-> _toSymList(List(i.tolist(), _unwrap))), @(1).val.params[2])] ),

# Gath * IMDPRDFT -> PrunedIMDPRDFT
    Gath_IMDPRDFT__PrunedIMDPRDFT := ARule(Compose, [[@(1, Gath), @(2,fTensor, e->ForAll(e.children(), i->ObjId(i)=fAdd))], @(3, IMDPRDFT)],
        e -> [ PrunedIMDPRDFT(@(3).val.params[1], List(@(2).val.children(), i-> _toSymList(List(i.tolist(), _unwrap))), @(3).val.params[2])] ),

# Gath * PRDFT -> PrunedIPRDFT
    Gath_PRDFT__PrunedIPRDFT := ARule(Compose, [[@(1, Gath), @(2,fAdd)], @(3, IPRDFT)],
        e -> [PrunedIPRDFT(@(3).val.params[1], @(3).val.params[2], 1,  _toSymList(List(@(2).val.tolist(), _unwrap)))]),

    Compose_TCompose := Rule(@@(1, Compose, (e,cx) -> IsBound(cx.TFCall) and Length(cx.TFCall) = 1 and cx.TFCall[1].hasTags()),
        e -> TCompose(e.children()))
));

# WarpX stuff
Class(RulesFFTXPromoteWarpX1, RuleSet);
RewriteRules(RulesFFTXPromoteWarpX1, rec(
    MultiX := Rule([@(0, TDAGNode), @(1), @(2), @(3).cond(e->IsList(e) and Length(e) = 1 and ObjId(e[1]) = nth and e[1].loc = X), ...],
        e -> TDAGNode(@(1).val * GathPtr(@(3).val[1], fId(Cols(@(1).val))), @(2).val, [ @(3).val[1].loc ]).withTags(@(0).val.getTags())),

    MultiY := Rule([@(0, TDAGNode), @(1), @(2).cond(e->IsList(e) and Length(e) = 1 and ObjId(e[1]) = nth and e[1].loc = Y), @(3), ...],
        e -> TDAGNode(ScatPtr(@(2).val[1], fId(Rows(@(1).val))) * @(1).val, [ @(2).val[1].loc ], @(3).val).withTags(@(0).val.getTags())),
));

Class(RulesFFTXPromoteWarpX2, RuleSet);
RewriteRules(RulesFFTXPromoteWarpX2, rec(
    TDAGNode_VStack1 := Rule([@(1, TDAG), ..., @(2, TDAGNode,
            e-> IsList(e.params[2]) and Length(e.params[2]) = 1 and ObjId(e.params[2][1]) = nth and e.params[2][1].loc <> Y and _unwrap(e.params[2][1].idx) = 0
                and not ObjId(e.params[3][1]) = nth),
            ...],
        e -> let(
            vr := @(2).val.params[2][1].loc,
            ch := Filtered(@(1).val.params[1], c->IsList(c.params[2]) and Length(c.params[2]) = 1 and ObjId(c.params[2][1]) = nth and c.params[2][1].loc = vr
                     and c.params[3] = @(2).val.params[3]),
            sl := Flat(SortRecordList(ch, c->_unwrap(c.params[2][1].idx))),
            spls := List(sl, i-> i.params[1]),
            vstk := VStack(spls),
            dn := TDAGNode(vstk, [vr], @(2).val.params[3]),
            nch := Filtered(@(1).val.params[1], c->(not IsList(c.params[2])) or (not Length(c.params[2]) = 1) or (not ObjId(c.params[2][1]) = nth) or (not c.params[2][1].loc = vr)
                     or (not c.params[3] = @(2).val.params[3])),
            TDAG([dn]::nch).withTags(@(1).val.tags)
        )),

    TDAGNode_VStack2 := Rule([@(1, TDAG), ..., @(2, TDAGNode,
            e-> IsList(e.params[3]) and Length(e.params[3]) = 1 and ObjId(e.params[3][1]) = nth and e.params[3][1].loc <> X and _unwrap(e.params[3][1].idx) = 0
                and not ObjId(e.params[2][1]) = nth),
            ...],
        e -> let(
            vr := @(2).val.params[3][1].loc,
            ch := Filtered(@(1).val.params[1], c->IsList(c.params[3]) and Length(c.params[3]) = 1 and ObjId(c.params[3][1]) = nth and c.params[3][1].loc = vr
                     and c.params[2] = @(2).val.params[2]),
            sl := Flat(SortRecordList(ch, c->_unwrap(c.params[3][1].idx))),
            spls := List(sl, i-> i.params[1]),
            hstk := HStack(spls),
            dn := TDAGNode(hstk,  @(2).val.params[2], [vr]),
            nch := Filtered(@(1).val.params[1], c->(not IsList(c.params[3])) or (not Length(c.params[3]) = 1) or (not ObjId(c.params[3][1]) = nth) or (not c.params[3][1].loc = vr)
                     or (not c.params[2] = @(2).val.params[2])),
            TDAG(nch::[dn]).withTags(@(1).val.tags)
        )),

   TResample_TGath := Rule(@(1, TResample, e->TResample_TGath.applicable(e)),
       e->let(l := @(1).val,
              c := TResample_TGath.children(l)[1],
              TResample_TGath.apply(l, c, c)
           )),

   TResample_TScat := Rule(@(1, TResample, e->TResample_TScat.applicable(e)),
       e->let(l := @(1).val,
              c := TResample_TScat.children(l)[1],
              TResample_TScat.apply(l, c, c)
           )),

   TResample_MD_nofrac := Rule(@(1, TResample, e->TResample_MD_nofrac.applicable(e)),
       e->let(l := @(1).val,
              c := TResample_MD_nofrac.children(l)[1],
              TResample_MD_nofrac.apply(l, c, c)
           )),

   TGath_Gath := Rule(@(1, TGath, e-> TGath_base.applicable(e)), e->ApplyRuleSPL(TGath_base, @(1).val)),

   TScat_Scat := Rule(@(1, TScat, e-> TScat_base.applicable(e)), e->ApplyRuleSPL(TScat_base, @(1).val)),

   Gath_GathPtr := ARule(Compose, [@(1, Gath, e-> ObjId(e.func) = fTensor and ForAll(e.func.children(), c->ObjId(c) = fId or (ObjId(c) = fAdd and c.params[3] = 0))), @(2, GathPtr)],
       e->[ GathPtr(@(2).val.ptr, fId(Rows(@(1).val))) ]),

   ScatPtr_Scat := ARule(Compose, [@(1, ScatPtr), @(2, Scat, e-> ObjId(e.func) = fTensor and ForAll(e.func.children(), c->ObjId(c) = fId or (ObjId(c) = fAdd and c.params[3] = 0)))],
       e->[ ScatPtr(@(1).val.ptr, fId(Cols(@(2).val))) ]),

   Expand_TResample := Rule(@(1, TResample, e->TResample_MD_frac.applicable(e)),
       e->let(trs := @(1).val,
               c := TResample_MD_frac.children(trs)[1],
               TResample_MD_frac.apply(trs, c, c)
           )),
));

Class(RulesFFTXPromoteWarpX3, RuleSet);
RewriteRules(RulesFFTXPromoteWarpX3, rec(
    Merge_TTensorI_VStack := ARule(Compose, [ @(1, TTensorI, e->e.params[3] = APar and e.params[4] = APar),
        @(2, VStack, e->@(1).val.params[2] = Length(e.children()) and ForAll(e.children(), e->Rows(e) = Cols(@(1).val.params[1]))) ],
        e->[VStack(List(@(2).val.children(), c->@(1).val.params[1] * c))]),

    Merge_HStack_TTensorI := ARule(Compose, [ @(1, HStack),
        @(2, TTensorI, e->e.params[3] = APar and e.params[4] = APar and
            e.params[2] = Length(@(1).val.children()) and ForAll(@(1).val.children(), j->Cols(j) = Rows(e.params[1]))) ],
        e->[HStack(List(@(1).val.children(), c->c * @(2).val.params[1]))]),

# this rule is WRONG: the sign needs to be inverse of each other!
    MDPRDFT_IMDPRDFT := ARule(Compose, [@(1, MDPRDFT), @(2, IMDPRDFT, e->Length(@(1).val.params[1]) = Length(e.params[1]) and Mod(@(1).val.params[2] + e.params[2], Product(e.params[1])) = 0 and
        ForAll([1..Length(e.params[1])], j->e.params[1][j] = @(1).val.params[1][j]))],
        e->[Diag(fConst(TReal, Rows(@(1).val), Product(@(1).val.params[1])))]),

    Diag_RCDiag := ARule(Compose, [@(1, Diag, e->ObjId(e.element) = fConst and e.element.range() = TReal), @(2, RCDiag)],
        e -> [RCDiag(diagMul(@(1).val.element, @(2).val.element))]),

    RCDiag_Diag := ARule(Compose, [@(1, RCDiag), @(2, Diag, e->ObjId(e.element) = fConst and e.element.range() = TReal)],
        e -> [RCDiag(diagMul(@(2).val.element, @(1).val.element))]),

 # these rules need a bit more guards to not misfire
    Peel_HStack_RCDiag := Rule(@(1, HStack, e->ForAll(e.children(), c->let(cc := c.children(), Length(cc) = 3 and
            ObjId(cc[1]) = ScatPtr and ObjId(cc[2]) = IMDPRDFT and ObjId(cc[3]) = RCDiag  and ObjId(cc[3].element) = RCData))),
        e -> let(cc := @(1).val.children(), i := Ind(Length(cc)), c1 := cc[1].children()[1], c2 := cc[1].children()[2],
#                Error(),
                fts := Flat(List(cc, c->Collect(c, fTensor))),
#                rng := List(fts, f->List([1..Length(f.children())], ii->f.child(ii).range())),
                _rng := Flat(List(fts, f->List([1..Length(f.children())], ii->f.child(ii).range()))),
#                rngf := FData(rng),
                _rngf := FData(_rng),
                dom := fts[1].child(1).domain(),
#                fts2 := fTensor(List([0..Length(fts[1].children())-1],  ii->fAdd(nth(rngf.at(i), ii), dom, 0))),
                _fts2 := fTensor(List([0..Length(fts[1].children())-1],  ii->fAdd(_rngf.at(Length(fts[1].children()) * i + ii), dom, 0))),
            TNoDiagPullinRight(TIterHStack(TCompose([ScatPtr(nth(c1.ptr.loc, i), _fts2), c2]), i)) *
                TNoPullLeft(RCDiag(RCData(diagDirsum(List(cc, c->Last(c.children()).element.func))))))),

    Peel_RCDiag_VStack := Rule(@(1, VStack, e->ForAll(e.children(), c->let(cc := Reversed(c.children()), Length(cc) in [2, 3] and
            ObjId(cc[1]) = GathPtr and ObjId(cc[2]) = MDPRDFT and When(Length(cc) = 3, ObjId(cc[3]) = RCDiag and ObjId(cc[3].element) = RCData, true)))),
        e-> let(cc := @(1).val.children(), i := Ind(Length(cc)), rcc := Reversed(cc[1].children()),
                fts := Flat(List(cc, c->Collect(c, fTensor))),
#                rng := List(fts, f->List([1..Length(f.children())], ii->f.child(ii).range())),
                _rng := Flat(List(fts, f->List([1..Length(f.children())], ii->f.child(ii).range()))),
#                rngf := FData(rng),
                _rngf := FData(_rng),
                dom := fts[1].child(1).domain(),
#                fts2 := fTensor(List([0..Length(fts[1].children())-1],  ii->fAdd(nth(rngf.at(i), ii), dom, 0))),
                _fts2 := fTensor(List([0..Length(fts[1].children())-1],  ii->fAdd(_rngf.at(Length(fts[1].children()) * i + ii), dom, 0))),

            TNoPullRight(RCDiag(RCData(diagDirsum(List(cc, c->When(ObjId(c.children()[1]) = RCDiag, c.children()[1].element.func,
                fConst(TComplex, Rows(rcc[2])/2, V(1)))))))) *
            TNoDiagPullinLeft(TIterVStack(TCompose([rcc[2], GathPtr(nth(rcc[1].ptr.loc, i), _fts2)]), i))
        )),

    Fuse_Diag := ARule(Compose, [ @(1, TNoPullLeft, e->ObjId(e.params[1]) = RCDiag), @(2, TRC), @(3, TNoPullRight, e->ObjId(e.params[1]) = RCDiag)],
        e->[ TRC(TGrp(TCompose([Diag(@(1).val.params[1].element.func), @(2).val.params[1], Diag(@(3).val.params[1].element.func)])))]),

    Promote_Compose := Rule(@(1, Compose, e->ForAll(e.children(), c->ObjId(c) in [TRC, TNoDiagPullinLeft, TNoDiagPullinRight])),
        e->TCompose(@(1).val.children())),

    flatten_I := Rule(@(1, Tensor, e->ForAll(e.children(), c->ObjId(c)=I)), e->I(Rows(@(1).val))),

));

RewriteRules(RulesSums, rec(
 PullInRightScatPtr := ARule( Compose,
       [ @(1, [ScatPtr, TNoPullLeft]),
         @(2, [RecursStep, Grp, BB, SUM, Buf, ISum, Data, COND, TNoDiagPullin, TNoDiagPullinLeft, TNoDiagPullinRight, NeedInterleavedComplex, SIMTISum]) ],
    e -> [ CopyFields(@(2).val, rec(
             _children :=  List(@(2).val._children, c -> @(1).val * c),
             dimensions := [Rows(@(1).val), Cols(@(2).val)] )) ]),

 PullInLeftGathPtr := ARule( Compose,
       [ @(1, [RecursStep, Grp, BB, SUM, SUMAcc, Buf, ISum, ISumAcc, Data, COND, TNoDiagPullin, TNoDiagPullinLeft, TNoDiagPullinRight, NeedInterleavedComplex, SIMTISum]),
         @(2, [GathPtr, TNoPullRight]) ],
    e -> [ CopyFields(@(1).val, rec(
                _children := List(@(1).val._children, c -> c * @(2).val),
                dimensions := [Rows(@(1).val), Cols(@(2).val)] )) ]),
 ComposeGathGathPtr := ARule(Compose, [ @(1, Gath), @(2, GathPtr) ], # o 1-> 2->
     e -> [ GathPtr(@(2).val.ptr, fCompose(@(2).val.func, @(1).val.func)) ]),

 ComposeScatPtrScat := ARule(Compose, [ @(1, ScatPtr), @(2, Scat) ], # <-1 <-2 o
     e -> [ ScatPtr(@(1).val.ptr, fCompose(@(1).val.func, @(2).val.func)) ]),

));

# IOPrunedRConv(@(2).val.params[1], @(2).val.params[2], 1, _toSymList(List(@(1).val.func.tolist(), _unwrap)), 1, _toSymList(List(@(3).val.func.tolist(), _unwrap)));
# IOPrunedRConv(@(2).val.params[1], FDataOfs(@(2).val.params[2].var, 2*(@(2).val.params[1]/2+1), 0), 1, _toSymList(List(@(1).val.func.tolist(), _unwrap)), 1, _toSymList(List(@(3).val.func.tolist(), _unwrap)), true);

max_stars_repo_name: mikefranusich/spiral-package-fftx

max_issues_repo_name: mikefranusich/spiral-package-fftx

max_forks_repo_name: mikefranusich/spiral-package-fftx


text:
if <condition> then
    <statements>
elif <condition> then
    <statements>
else
    <statements>
fi;

ext: gap

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text:
Reversed("abcdef");
# "fedcba"


text:
local reps,tori,labels,I,II,ords,blocks;

#......Ross
reps:=[[[1,1],[2,1],[3,1],[4,1],[5,1]],    # D_5
       [[5,1],[4,1],[9,1],[7,1],[1,1]]     # D_5(a_1)
       ];

I:=[[],[]];

#
# The parabolic where the representative is distinguished
#
II:=[[],
     []];

labels:=["D_5",
         "D_5(a_1)"
         ];

tori:=[
       [],
       []
    ];

blocks:=[];

ords:=[];

return [reps,tori,labels,I,II,ords,blocks];


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

##################################
# I currently have a streaming algorithm (from Nikara et al.) for two variations of DCT-2.
#
# First, we have Spiral's standard DCT-2.  This is called TDCT2, defined:
#    [ cos(k*(l+1/2)*pi/n) | k,l = 0...n-1 ]
#
# Then, I have added a scaled version, defined:
#    [ Diag(1/sqrt(2), 1, 1, ...  ] * [ cos(k*(l+1/2)*pi/n) | k,l = 0...n-1 ]
# which matches the definition in the Nikara paper.

# We use essentially the same algorithm for both; the only change is in one diagonal.

# Using opts := InitStreamUnrollHw(); will enable these rules.

# Other DCT-2 definitions also have a sqrt(2/n) scaling factor in front of the whole matrix.

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

Declare(Sc_DCT2_func);

#F Sc_DCT2(<n>) - Scaled Discrete Cosine Transform, Type II, non-terminal
#F Definition: (n x n)-matrix [ Diag(1/sqrt(2), 1, 1, ...  ] * [ cos(k*(l+1/2)*pi/n) | k,l = 0...n-1 ]
#F Example:    DCT2(8)
Class(Sc_DCT2, TaggedNonTerminal, rec(

    abbrevs := [
    (n)       -> Checked(IsPosIntSym(n),
        [_unwrap(n)]),
    ],

    hashAs := self >> ObjId(self)(self.params[1], 1).withTags(self.getTags()),

    dims := self >> [ self.params[1], self.params[1] ],

    isReal := self >> true,

    terminate := self >> Mat(Sc_DCT2_func(self.params[1])),
    transpose := self >> Mat(TransposedMat(Sc_DCT2_func(self.params[1]))),
));

# hadamard_func(n,i): n-point Hadamard permutation, position i.
# Defined as in:
#    Z. Wang, Pruning the fast discrete Cosine transform, IEEE
#    Tr. Communications 39(5), May 1991, 640-643.
hadamard_func := function(n,i)
   if ((n=1) and (i=0)) then return 0; fi;
   if (imod(n,2) <> 0) then return -1;  fi; # ERROR
   if (imod(i,2) = 0) then return hadamard_func(n/2, i/2); fi;

   return n-1-hadamard_func(n/2, (i-1)/2);

end;

# An Exp wrapper for hadamard_func.  Used so Spiral does not try to evaluate the expression
# until Process_fPrecompute is called.

Class(hadamard_func_exp, Exp, rec(
    ev := self >> hadamard_func(self.args[1].ev(), self.args[2].ev())
));

# Had(n): Hadamard permutation on n points.
# See:
#    Z. Wang, Pruning the fast discrete Cosine transform, IEEE
#    Tr. Communications 39(5), May 1991, 640-643.
Class(Had, PermClass, rec(
#    exportSymbol := self>>self.name,
#    exportParams := self>>self.params,
#    export := Sym.export,

    def := (n) -> Checked(
        IsPosIntSym(n),
        rec(size := n)),

    lambda := self >> let(
        n := self.params[1], i := Ind(n),
        lt := List([0..n-1], it->hadamard_func(n,it)),
        FList(TInt, lt).lambda()
    ),

    transpose := self >> Error("Transpoed Hadamard permutation not currently supported."),

    # Only symmetric for size 2.
    isSymmetric := self >> (self.params[1] = 2)
));

# Don't think I need this anymore.
# tSPL Hadamard permutation on n points
# See:
#    Z. Wang, Pruning the fast discrete Cosine transform, IEEE
#    Tr. Communications 39(5), May 1991, 640-643.
Class(THad, Tagged_tSPL_Container, rec(
    abbrevs :=  [ size -> [size] ],

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

    terminate := self >> Had(self.params[1]),
    transpose := self >> Error("Transposed Hadamard permutation not currently supported."),
    isReal := self >> true,
    isSymmetric := self >> self.params[1] = 2,

    # We can verify this with:
    # PermMatrixToBits(MatSPL(Had(n)))
    permBits := self >> let(n:=self.params[1], logn:=Log2Int(n),
        (MatSPL(DirectSum(J(logn-1), O(1,1))) + MatSPL(J(logn))) *
        GF(2).one
    )

));

########################################
## The following are helper functions used in the streaming DCT2
## diagonal given in
##    Nikara et al., Discrete cosine and sine transforms--regular
##    algorithms and pipeline architectures, Signal Processing 86,
##    2006.

mu := (s,i) -> Cond(imod(i,2^s)=0, 0, 1);
tau := (i,s) -> Cond(s=i, 0, 1);

Class(mu_exp, Exp, rec(
   ev := self >> mu(self.args[1].ev(), self.args[2].ev())
));

Class(tau_exp, Exp, rec(
   ev := self >> tau(self.args[1].ev(), self.args[2].ev())
));

dct_diag_f := (k,i,s) -> (imod(i,2) + (1-tau_exp(0,i)) * (1-tau_exp(k-1,s)));

#dct_diag_d := i -> let(K := 2^(Log2Int(i)), t := i-K,
#    h := hadamard_func(K, t),
#    cospi((h+1/2)/(2*K)));

dct_diag_d := i -> let(K := 2^(floor(fdiv(log(i), log(2)))), t := i-K,
    h := hadamard_func_exp(K, t),
    cospi(fdiv((h+1/2), (2*K))));

dct_sc_diag_g := (k,i,s) -> ((2^(mu_exp(s, floor(fdiv(i,2))))) *
    (dct_diag_d(2^(k-s-1) + floor(fdiv(i, (2^(s+1)))))))^dct_diag_f(k,i,s);

dct_unsc_diag_g := (k,i,s) -> ((2^(mu_exp(s, floor(fdiv(i,2))))) *
    (dct_diag_d(2^(k-s-1) + floor(fdiv(i, (2^(s+1)))))))^imod(i,2);

## This represents the diagonal matrix used in the streaming
## DCT2 algorithm. The problem is of size 2^k, and s represents
## the iteration
#Str_DCT2_Diag := (k, s) -> Diag(List([0..((2^k)-1)],
#    i-> dct_diag_g(k, i, s)));

Class(Str_Sc_DCT2_Diag, DiagFunc, rec(
    abbrevs := [(k, s) -> [k, s]],
    def := (k, s) -> rec(size := 2^k),
    lambda := self >> let(k := self.params[1], s := self.params[2], i := Ind(2^k), Lambda(i, dct_sc_diag_g(k, i, s))),
    range := self >> TReal,
));

Class(Str_DCT2_Diag, DiagFunc, rec(
    abbrevs := [(k, s) -> [k, s]],
    def := (k, s) -> rec(size := 2^k),
    lambda := self >> let(k := self.params[1], s := self.params[2], i := Ind(2^k), Lambda(i, dct_unsc_diag_g(k, i, s))),
    range := self >> TReal,
));

Class(Str_DCT2_Perm, PermClass, rec(
    def := (k) -> Checked(
        IsPosIntSym(k),
        rec(size := 2^k)),

    # Sigh, this is broken it seems.
    lambda := self >> let(
        k := self.params[1],
	fCompose(Reversed(List([0..(k-2)], i-> let(rsize := 2^k - 2^(k-i),

            fCompose(
		fTensor(fId(2^i), L(2^(k-i), 2^(k-i-1))),
		fDirsum(fId(2^(k-i)), fTensor(fId(rsize/4), fDirsum(fId(2), J(2))))
	    )

	)))).lambda()),

    transpose := self >> self,
    isSymmetric := self >> true,

));

Class(Str_DCT2_Perm_tspl, TaggedNonTerminal, rec(

    abbrevs := [(k)      -> [k]],

    dims := self >> [ 2^self.params[1], 2^self.params[1] ],

    terminate := self >> let(k := self.params[1],
	    Compose(List([0..(k-2)], i-> let(rsize := 2^k - 2^(k-i),
            DirectSum(I(2^(k-i)),
                Tensor(I(rsize/4), DirectSum(I(2), J(2)))
            ) *
            Tensor(I(2^i), L(2^(k-i), 2^(k-i-1)))
	    )
        ))
    ),

    isReal := self >> true,

    print := meth(self, indent, indentStep)
        local lparams, mparams;
        if not IsBound(self.params) then Print(self.name); return; fi;
        Print(self.name, "(");
        if IsList(self.params) then
            lparams := Filtered(self.params, i->not (IsList(i) and i=[]));
            mparams := Filtered(lparams, i->not (IsBool(i) and not i));
            DoForAllButLast(mparams, x -> Print(x, ", "));
            Print(Last(mparams));
        else
            Print(self.params);
        fi;
    Print(")", When(self.transposed, ".transpose()", ""));
    end,

));

Str_DCT2_M := (n,s) -> let(l := Ind(n/2), TTensorInd(COND(eq(imod(l, 2^s), 0), I(2), Mat([[1,0],[-1,1]])), l, APar, APar));

Str_DCT2_H := (n,s) -> let(k := Log2Int(n), l := Ind(n/4), Cond(s=0, I(n), TTensorInd(COND(eq(imod(l, 2^(s-1)),0), I(4), Mat([[1,0,0,0],[0,0,0,1],[0,0,1,0],[0,1,0,0]])), l, APar, APar)));

#Sc_DCT2_func := (N) -> List([0..N-1], m-> List([0..N-1], n-> ((Sqrt(2/N) * Cond(m=0, 1/Sqrt(2), 1) * CosPi(m*(n+1/2)/N)))));
Sc_DCT2_func := (N) -> List([0..N-1], m-> List([0..N-1], n-> (Cond(m=0, 1/Sqrt(2), 1) * CosPi(m*(n+1/2)/N))));

NewRulesFor(Sc_DCT2, rec(
    Sc_DCT2_Stream := rec(
        forTransposition := false,
        applicable := (self, nt) >> nt.params[1] > 4 and nt.isTag(1, AStream) and nt.firstTag().bs >= 4,
        children := (self, t) >> let(
            n := t.params[1],
            k := Log2Int(n),
            [[ TCompose(Concatenation(
                  [TPrm(Str_DCT2_Perm(k))],
                  Reversed(List([1..k-1], s-> TCompose([Str_DCT2_M(n, s), TDiag(fPrecompute(Str_Sc_DCT2_Diag(k, s))), Str_DCT2_H(n, s), TTensorI(F(2), n/2, APar, APar), TTensorI(TPrm(L(2^(s+1), 2^s)), 2^(k-s-1), APar, APar),
#TL(2^(s+1), 2^s, 2^(k-s-1), 1)
]))),
                  [TDiag(fPrecompute(Str_Sc_DCT2_Diag(k, 0))), TTensorI(F(2), n/2, APar, APar), TPrm(Had(n))])).withTags(t.getTags())
            ]]
         ),

        apply := (t,c,nt) -> c[1],
    )
));

NewRulesFor(TDCT2, rec(
    DCT2_Stream := rec(
        forTransposition := false,
        applicable := (self, nt) >> nt.params[1] > 4 and nt.isTag(1, AStream) and nt.firstTag().bs >= 4,
        children := (self, t) >> let(
            n := t.params[1],
            k := Log2Int(n),
            [[ TCompose(Concatenation(
                  [Cond(t.firstTag().bs = t.params[1],
			  TPrm(Str_DCT2_Perm(k)),
			  TPrm(Str_DCT2_Perm_tspl(k))
		      )
		      ],
                  Reversed(List([1..k-1], s-> TCompose([Str_DCT2_M(n, s), TDiag(fPrecompute(Str_DCT2_Diag(k, s))), Str_DCT2_H(n, s), TTensorI(F(2), n/2, APar, APar), TTensorI(TPrm(L(2^(s+1), 2^s)), 2^(k-s-1), APar, APar),
]))),
                  [TDiag(fPrecompute(Str_DCT2_Diag(k, 0))), TTensorI(F(2), n/2, APar, APar), TPrm(Had(n))])).withTags(t.getTags())
            ]]
         ),

        apply := (t,c,nt) -> c[1],
    )
));


text:
>1apf
GVPCLCDSDGPRPRGNTLSGILWFYPSGCPS2GWH1NCKAHGPNIGWCCKK2
>1ahl
GVSCLCDSDGPSVRGNTLSGTLWLYPSGCPS2GWH1NCKAHGPTIGWCCKQ2
>1atx
GAACLCKSDGPNTRGNSMSGTIWVF2GCPS2GWN1NCEGRA1IIGYCCKQ2
>1sh1
1AACKCDDEGPDIRTAPLTGTVDLG2SCNA2GWE1KCASYYTIIADCCRKKK
>1bds
AAPCFCSGKP7GRGDLWILRGTCPGGYGYTSNCYK2WPNICCYPH2

max_stars_repo_name: JoseEspinosa/nf-benchmark

max_issues_repo_name: JoseEspinosa/nf-benchmark

max_forks_repo_name: JoseEspinosa/nf-benchmark


text:
%Terminals
    DollarSign ::= '$'
    Percent ::= '%'
    _

    a    b    c    d    e    f    g    h    i    j    k    l    m
    n    o    p    q    r    s    t    u    v    w    x    y    z
%End

%Headers
    /.
        final static int tokenKind[] = new int[128];
        static
        {
            tokenKind['$'] = $sym_type.$prefix$DollarSign$suffix$;
            tokenKind['%'] = $sym_type.$prefix$Percent$suffix$;
            tokenKind['_'] = $sym_type.$prefix$_$suffix$;

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

        final int getKind(int c)
        {
            return ((c & 0xFFFFFF80) == 0 /* 0 <= c < 128? */ ? tokenKind[c] : 0);
        }
    ./
%End

max_stars_repo_name: kuafuwang/LPG-VScode

max_issues_repo_name: kuafuwang/LPG-VScode

max_forks_repo_name: kuafuwang/LPG-VScode


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

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


text:
a := [1 .. 4];
b := ShallowCopy(a);

# Apply and replace values
Apply(a, n -> n*n);
a;
# [ 1, 4, 9, 16 ]

# Apply and don't change values
List(b, n -> n*n);
# [ 1, 4, 9, 16 ]

# Apply and don't return anything (only side effects)
Perform(b, Display);
1
2
3
4

b;
# [ 1 .. 4 ]


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

Declare(ScatPtr);

Class(GathPtr, Gath, rec(
    rChildren := self >> [self.ptr, self.func],
    rSetChild := rSetChildFields("ptr", "func"),
    new := (self, ptr, func) >> SPL(WithBases(self, rec(
        ptr := ptr,
      	func := Checked(IsFunction(func) or IsFuncExp(func), func)))).setDims()
));

Class(ScatPtr, Scat, rec(
    rChildren := self >> [self.ptr, self.func],
    rSetChild := rSetChildFields("ptr", "func"),
    new := (self, ptr, func) >> SPL(WithBases(self, rec(
        ptr := ptr,
	    func := Checked(IsFunction(func) or IsFuncExp(func), func)))).setDims()
));

Class(OO, O);

max_stars_repo_name: franzfranchetti/spiral-package-fftx

max_issues_repo_name: franzfranchetti/spiral-package-fftx

max_forks_repo_name: franzfranchetti/spiral-package-fftx


text:
%Headers
    --
    -- Additional methods for the action class not provided in the template
    --

	-- [cwd] Template provided by LPG defines a constructor that uses an Option
	--       class that does not exist in LPG Runtime.  Deleted this constructor

    /.
        final void makeToken(int kind)
        {
            int startOffset = getLeftSpan(),
                endOffset = getRightSpan();
            makeToken(startOffset, endOffset, kind);
            if (printTokens) printValue(startOffset, endOffset);
        }

        final void makeComment(int kind)
        {
            int startOffset = getLeftSpan(),
                endOffset = getRightSpan();
            super.getIPrsStream().makeAdjunct(startOffset, endOffset, kind);
        }

        final void skipToken()
        {
            if (printTokens) printValue(getLeftSpan(), getRightSpan());
        }

        final void checkForKeyWord()
        {
            int startOffset = getLeftSpan(),
                endOffset = getRightSpan(),
            kwKind = kwLexer.lexer(startOffset, endOffset);
            makeToken(startOffset, endOffset, kwKind);
            if (printTokens) printValue(startOffset, endOffset);
        }

        final void printValue(int startOffset, int endOffset)
        {
            String s = new String(getInputChars(), startOffset, endOffset - startOffset + 1);
            System.out.print(s);
        }

        //
        //
        //
        public final static int tokenKind[] =
        {
            Char_CtlCharNotWS,    // 000    0x00
            Char_CtlCharNotWS,    // 001    0x01
            Char_CtlCharNotWS,    // 002    0x02
            Char_CtlCharNotWS,    // 003    0x03
            Char_CtlCharNotWS,    // 004    0x04
            Char_CtlCharNotWS,    // 005    0x05
            Char_CtlCharNotWS,    // 006    0x06
            Char_CtlCharNotWS,    // 007    0x07
            Char_CtlCharNotWS,    // 008    0x08
            Char_HT,              // 009    0x09
            Char_LF,              // 010    0x0A
            Char_CtlCharNotWS,    // 011    0x0B
            Char_FF,              // 012    0x0C
            Char_CR,              // 013    0x0D
            Char_CtlCharNotWS,    // 014    0x0E
            Char_CtlCharNotWS,    // 015    0x0F
            Char_CtlCharNotWS,    // 016    0x10
            Char_CtlCharNotWS,    // 017    0x11
            Char_CtlCharNotWS,    // 018    0x12
            Char_CtlCharNotWS,    // 019    0x13
            Char_CtlCharNotWS,    // 020    0x14
            Char_CtlCharNotWS,    // 021    0x15
            Char_CtlCharNotWS,    // 022    0x16
            Char_CtlCharNotWS,    // 023    0x17
            Char_CtlCharNotWS,    // 024    0x18
            Char_CtlCharNotWS,    // 025    0x19
            Char_CtlCharNotWS,    // 026    0x1A
            Char_CtlCharNotWS,    // 027    0x1B
            Char_CtlCharNotWS,    // 028    0x1C
            Char_CtlCharNotWS,    // 029    0x1D
            Char_CtlCharNotWS,    // 030    0x1E
            Char_CtlCharNotWS,    // 031    0x1F
            Char_Space,           // 032    0x20
            Char_Exclamation,     // 033    0x21
            Char_DoubleQuote,     // 034    0x22
            Char_Sharp,           // 035    0x23
            Char_DollarSign,      // 036    0x24
            Char_Percent,         // 037    0x25
            Char_Ampersand,       // 038    0x26
            Char_SingleQuote,     // 039    0x27
            Char_LeftParen,       // 040    0x28
            Char_RightParen,      // 041    0x29
            Char_Star,            // 042    0x2A
            Char_Plus,            // 043    0x2B
            Char_Comma,           // 044    0x2C
            Char_Minus,           // 045    0x2D
            Char_Dot,             // 046    0x2E
            Char_Slash,           // 047    0x2F
            Char_0,               // 048    0x30
            Char_1,               // 049    0x31
            Char_2,               // 050    0x32
            Char_3,               // 051    0x33
            Char_4,               // 052    0x34
            Char_5,               // 053    0x35
            Char_6,               // 054    0x36
            Char_7,               // 055    0x37
            Char_8,               // 056    0x38
            Char_9,               // 057    0x39
            Char_Colon,           // 058    0x3A
            Char_SemiColon,       // 059    0x3B
            Char_LessThan,        // 060    0x3C
            Char_Equal,           // 061    0x3D
            Char_GreaterThan,     // 062    0x3E
            Char_QuestionMark,    // 063    0x3F
            Char_AtSign,          // 064    0x40
            Char_A,               // 065    0x41
            Char_B,               // 066    0x42
            Char_C,               // 067    0x43
            Char_D,               // 068    0x44
            Char_E,               // 069    0x45
            Char_F,               // 070    0x46
            Char_G,               // 071    0x47
            Char_H,               // 072    0x48
            Char_I,               // 073    0x49
            Char_J,               // 074    0x4A
            Char_K,               // 075    0x4B
            Char_L,               // 076    0x4C
            Char_M,               // 077    0x4D
            Char_N,               // 078    0x4E
            Char_O,               // 079    0x4F
            Char_P,               // 080    0x50
            Char_Q,               // 081    0x51
            Char_R,               // 082    0x52
            Char_S,               // 083    0x53
            Char_T,               // 084    0x54
            Char_U,               // 085    0x55
            Char_V,               // 086    0x56
            Char_W,               // 087    0x57
            Char_X,               // 088    0x58
            Char_Y,               // 089    0x59
            Char_Z,               // 090    0x5A
            Char_LeftBracket,     // 091    0x5B
            Char_BackSlash,       // 092    0x5C
            Char_RightBracket,    // 093    0x5D
            Char_Caret,           // 094    0x5E
            Char__,               // 095    0x5F
            Char_BackQuote,       // 096    0x60
            Char_a,               // 097    0x61
            Char_b,               // 098    0x62
            Char_c,               // 099    0x63
            Char_d,               // 100    0x64
            Char_e,               // 101    0x65
            Char_f,               // 102    0x66
            Char_g,               // 103    0x67
            Char_h,               // 104    0x68
            Char_i,               // 105    0x69
            Char_j,               // 106    0x6A
            Char_k,               // 107    0x6B
            Char_l,               // 108    0x6C
            Char_m,               // 109    0x6D
            Char_n,               // 110    0x6E
            Char_o,               // 111    0x6F
            Char_p,               // 112    0x70
            Char_q,               // 113    0x71
            Char_r,               // 114    0x72
            Char_s,               // 115    0x73
            Char_t,               // 116    0x74
            Char_u,               // 117    0x75
            Char_v,               // 118    0x76
            Char_w,               // 119    0x77
            Char_x,               // 120    0x78
            Char_y,               // 121    0x79
            Char_z,               // 122    0x7A
            Char_LeftBrace,       // 123    0x7B
            Char_VerticalBar,     // 124    0x7C
            Char_RightBrace,      // 125    0x7D
            Char_Tilde,           // 126    0x7E
            Char_CtlCharNotWS,    // 127    0x7F

            Char_Acute,           // for the acute accent 0xb4
            Char_AfterASCIINotAcute,  // for all chars in range 0x80..0xfffe excluding the acute accent
            Char_EOF              // for '\uffff' or 65535
        };

        @Override public final int getKind(int i)  // Classify character at ith location
        {
            char c = (i >= getStreamLength() ? '\uffff' : getCharValue(i));
            return (c < 128)? // ASCII Character
                      tokenKind[c] :
                      (c == '\uffff')?
                           Char_EOF :
                           (c == '\u00b4')?
                           Char_Acute :
                               Char_AfterASCIINotAcute;
        }
    ./
%End


text:
Declare(CUFFTCall);

Class(CUFFTCall, BaseMat, SumsBase, rec(
    dims := self >> self.L.dims(),
    isReal := self >> self.L.isReal(),
    #-----------------------------------------------------------------------
    rChildren := self >> [self.L, self.codegen],
    rSetChild := rSetChildFields("L", "codegen"),
    #-----------------------------------------------------------------------
    new := (self, L, codegen) >> SPL(WithBases(self,
        rec(L   := L,
            codegen   := CopyFields(codegen, rec(
                plan_var := var.fresh_t("plan", TSym("cufftHandle")),
                size_var := var.fresh_t("size", TInt),
                data := (self, c) >> data(self.size_var, V(self.N), c),
                init := self >> call(rec(id:="cufftPlanMany"), addrof(self.plan_var), V(1), addrof(self.size_var), addrof(self.size_var), V(1), V(self.M),
    	            addrof(self.size_var), V((self.N * self.K) / (self.p1 * self.p2)), V(1), "CUFFT_Z2Z", V((self.N * self.K) / (self.p1 * self.p2))),
            )),
            dimensions     := L.dims())
    )),

    #-----------------------------------------------------------------------
    transpose := self >> CUFFTCall(self.L.transpose(), self.codegen),
    #-----------------------------------------------------------------------

    print := (self,i,is) >> Print(self.name, "(",
        self.L.print(i+is,is), ", <codegen>)"),
    #-----------------------------------------------------------------------
    toAMat := self >> self.L.toAMat(),
    #-----------------------------------------------------------------------
    isPermutation := False
));

RewriteRules(RulesRC, rec(
    RC_CUFFTCall := Rule([RC, @(1, CUFFTCall)], e -> CUFFTCall(RC(@(1).val.L), @(1).val.codegen)),
));

CudaCodegen.CUFFTCall := (self, o, y, x, opts) >> simt_block(
    call(rec(id := "cufftExecZ2Z", codegen := o.codegen), o.codegen.plan_var,
        tcast(TPtr(TSym("cufftDoubleComplex")), x), tcast(TPtr(TSym("cufftDoubleComplex")), y), When(o.codegen.k = 1, "CUFFT_INVERSE", "CUFFT_FORWARD")));


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

################################################################################
##
#F pseudounitary_sort(<list>,<list>,<list>) . . . . sorts S and T to give pseudo
##. . . . . . . . . . . . . . .  . . . . . . .unitary data. FPdims are computed.
##
InstallGlobalFunction(pseudounitary_sort, function(S,T,Simples)
  local FPdim,FPdimC,x,r,i,FPRow,FPRow_idx,j,is_real,is_pos;
  #
  # look along the diagonal for 1's and check that the corresponding row is real
  # and positive
  FPRow:=[];
  r := Size(Simples);
  for i in [1..r] do
    if S[i][i] = 1 then
      is_real := true;
      is_pos := true;
      for j in [1..r] do
        if ImaginaryPart(S[i][j]) <> 0 then
          is_real := false;
          break;
        fi;
        if S[i][j] <= 0 then
          is_pos := false;
          break;
        fi;
      od;
      if is_real and is_pos then
        Add(FPRow,ShallowCopy(S[i]));
      fi;;
    fi;
  od;

  # find out which row is largest
  Sort(FPRow);
  FPRow:=Reversed(FPRow);
  FPRow_idx:=Position(S,FPRow[1]);

  # switch 1 with FPRow_idx
  if FPRow_idx <> 1 then
    S{[1,FPRow_idx]}:=S{[FPRow_idx,1]};
    S:=TransposedMatMutable(S);
    S{[1,FPRow_idx]}:=S{[FPRow_idx,1]};
    S:=TransposedMatMutable(S);
    Simples{[1,FPRow_idx]}:=Simples{[FPRow_idx,1]};
    T{[1,FPRow_idx]}:=T{[FPRow_idx,1]};
  fi;

  # get the FPdimension
  FPdim:=ShallowCopy(S[1]);
  FPdimC:=Sum(List(FPdim,x->x^2));

  return [Simples,S,T,FPdim,FPdimC];
end);

################################################################################
##
#F dimension_sort(<list>,<list>,<list>,<list>) . . . . sorts data by increasing
#                                               dimension.
##
InstallGlobalFunction(dimension_sort, function(S,T,Simples,FPdim)
  local made_switch, x, r;
  made_switch := true;
  r:=Size(FPdim);
  while made_switch do
    made_switch := false;
    for x in [3..r] do
      if FPdim[x] < FPdim[x-1] then
        S{[x,x-1]} := S{[x-1,x]};
        S:=TransposedMatMutable(S);
        S{[x,x-1]} := S{[x-1,x]};
        S:=TransposedMatMutable(S);
        T{[x,x-1]} := T{[x-1,x]};
        FPdim{[x,x-1]} := FPdim{[x-1,x]};
        Simples{[x,x-1]} := Simples{[x-1,x]};
        made_switch := true;
      fi;
    od;
  od;
  return [Simples,S,T,FPdim];
end);
#E pseudounitary_sort.gi . . . . . . . . . . . . . . . . . . . . . . . ends here

max_stars_repo_name: pnnl/GroupTheoretical

max_issues_repo_name: pnnl/GroupTheoretical

max_forks_repo_name: pnnl/GroupTheoretical


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

benchMACRO := function()
    return rec(
        2xf := rec(
            wht := _defaultSizes(s->doSimdWht(s, MACRO_2xf, rec(verify := true, oddSizes := false, svct := true)), List([3..8], i->2^i)),
            1d := rec(
                dft_ic := rec(
                    small := rec(
                        real := _defaultSizes(s->doSimdDft(s, MACRO_2xf, rec(verify := true, globalUnrolling := 10000, tsplBluestein:=true,
                          PRDFT:=true, URDFT:= true, CT := true, minCost := true, tsplPFA := true, tsplCT := true, oddSizes := true, PD := true, PFA := true,
                          RealRader := true, Rader := true, realVect := true, cplxVect := false, stdTTensor := false, interleavedComplex := true)), [2..64]),
                        cmplx := _defaultSizes(s->doSimdDft(s, MACRO_2xf, rec(verify := true, globalUnrolling := 10000, tsplBluestein:=true,
                          PRDFT:=true, URDFT:= true, CT := true, minCost := true, tsplPFA := true, tsplCT := true, oddSizes := true, PD := true, PFA := true,
                          RealRader := true, Rader := true, realVect := false, cplxVect := true, stdTTensor := false, interleavedComplex := true)), [2..64])
                    ),
                    medium := _defaultSizes(s->doSimdDft(s, MACRO_2xf, rec(verify := true, globalUnrolling := 128, tsplRader:=false, tsplBluestein:=false, tsplPFA:=false, oddSizes:=false,
                            interleavedComplex := true, cplxVect := false, realVect := true)),
                          List([2..16], i->2^i)),
                    medium_cx := _defaultSizes(s->doSimdDft(s, MACRO_2xf, rec(globalUnrolling := 64, tsplRader:=false, tsplBluestein:=false, tsplPFA:=false, oddSizes:=false,
                            interleavedComplex := true, cplxVect := true, realVect := false)),
                          List([1..16], i->2^i))
                ),
                dft_sc := _defaultSizes(s->doSimdDft(s, MACRO_2xf, rec(verify := true, globalUnrolling := 10000, tsplRader:=false, tsplBluestein:=false, tsplPFA:=false, oddSizes:=false, interleavedComplex := false)), 8*[1..16])
            )
        ),

        4xf := rec(
            wht := _defaultSizes(s->doSimdWht(s, MACRO_4xf, rec(verify := true, oddSizes := false, svct := true)), List([4..8], i->2^i)),
            1d := rec(
                dft_ic := rec(
                    small := rec(
                        real := _defaultSizes(s->doSimdDft(s, MACRO_4xf, rec(verify := true, globalUnrolling := 10000, tsplBluestein:=true,
                          PRDFT:=true, URDFT:= true, CT := true, minCost := true, tsplPFA := true, tsplCT := true, oddSizes := true, PD := true, PFA := true,
                          RealRader := true, Rader := true, realVect := true, cplxVect := false, stdTTensor := false, interleavedComplex := true)), [2..64]),
                        cmplx := _defaultSizes(s->doSimdDft(s, MACRO_4xf, rec(verify := true, globalUnrolling := 10000, tsplBluestein:=true,
                          PRDFT:=true, URDFT:= true, CT := true, minCost := true, tsplPFA := true, tsplCT := true, oddSizes := true, PD := true, PFA := true,
                          RealRader := true, Rader := true, realVect := false, cplxVect := true, stdTTensor := false, interleavedComplex := true)), [2..64])
                    ),
                    medium := _defaultSizes(s->doSimdDft(s, MACRO_4xf, rec(verify := true, globalUnrolling := 128, tsplRader:=false, tsplBluestein:=false, tsplPFA:=false, oddSizes:=false,
                            interleavedComplex := true, cplxVect := false, realVect := true)),
                          List([4..16], i->2^i)),
                    medium_cx := _defaultSizes(s->doSimdDft(s, MACRO_4xf, rec(globalUnrolling := 64, tsplRader:=false, tsplBluestein:=false, tsplPFA:=false, oddSizes:=false,
                            interleavedComplex := true, cplxVect := true, realVect := false)),
                          List([2..16], i->2^i))
                ),
                dft_sc := _defaultSizes(s->doSimdDft(s, MACRO_4xf, rec(verify := true, globalUnrolling := 10000, tsplRader:=false, tsplBluestein:=false, tsplPFA:=false, oddSizes:=false, interleavedComplex := false)), 16*[1..8])
            )
        ),
        8xf := rec(
            wht := _defaultSizes(s->doSimdWht(s, MACRO_8xf, rec(verify := true, oddSizes := false, svct := true)), List([6..8], i->2^i)),
            1d := rec(
                dft_ic := rec(
                    small := rec(
                        real := _defaultSizes(s->doSimdDft(s, MACRO_8xf, rec(verify := true, globalUnrolling := 10000, tsplBluestein:=true,
                          PRDFT:=true, URDFT:= true, CT := true, minCost := true, tsplPFA := true, tsplCT := true, oddSizes := true, PD := true, PFA := true,
                          RealRader := true, Rader := true, realVect := true, cplxVect := false, stdTTensor := false, interleavedComplex := true)), [2..64]),
                        cmplx := _defaultSizes(s->doSimdDft(s, MACRO_8xf, rec(verify := true, globalUnrolling := 10000, tsplBluestein:=true,
                          PRDFT:=true, URDFT:= true, CT := true, minCost := true, tsplPFA := true, tsplCT := true, oddSizes := true, PD := true, PFA := true,
                          RealRader := true, Rader := true, realVect := false, cplxVect := true, stdTTensor := false, interleavedComplex := true)), [2..64])
                    ),
                    medium := _defaultSizes(s->doSimdDft(s, MACRO_8xf, rec(verify := true, globalUnrolling := 128, tsplRader:=false, tsplBluestein:=false, tsplPFA:=false, oddSizes:=false,
                            interleavedComplex := true, cplxVect := false, realVect := true)),
                          List([6..16], i->2^i)),
                    medium_cx := _defaultSizes(s->doSimdDft(s, MACRO_8xf, rec(globalUnrolling := 128, tsplRader:=false, tsplBluestein:=false, tsplPFA:=false, oddSizes:=false,
                            interleavedComplex := true, cplxVect := true, realVect := false)),
                          List([4..16], i->2^i))
                ),
                dft_sc := _defaultSizes(s->doSimdDft(s, MACRO_8xf, rec(verify := true, globalUnrolling := 10000, tsplRader:=false, tsplBluestein:=false, tsplPFA:=false, oddSizes:=false, interleavedComplex := false)), 64*[1..4])
            )
        )
    );
end;