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

Import(approx);

SubMat := (mat, ri, ci, rsize, csize) -> let(
    d := DimensionsMat(mat),
    nrows := d[1],
    ncols := d[2],
    List(mat{[ri .. Min2(ri+rsize-1, nrows)]},
	row -> row{[ci .. Min2(ci+csize-1, ncols)]})
);

GetBlocksMat := (mat, rblk, cblk) -> Checked(IsPosInt(rblk), IsPosInt(cblk), let(
    d := DimensionsMat(mat),
    nrows := d[1],
    ncols := d[2],
    List([0 .. CeilingRat(nrows/rblk)-1], rb ->
	List([0 .. CeilingRat(ncols/cblk)-1], cb ->
	    SubMat(mat, 1+rb*rblk, 1+cb*cblk, rblk, cblk)))
));

# Cyclic diagonal vectorization:
VectorizedZMatSPL := ( isa, spl ) -> Checked( Rows(spl)=Cols(spl) and isa.v = Rows(spl), let(
    mat := Cond(IsMat(spl), spl, MatSPL(spl)),
    VBlk( [List( [0 .. isa.v-1], j -> isa.t.value(
                List( [0 .. isa.v-1], i -> mat[1 + i][1 + (i-j) mod isa.v])
         ))],
         isa.v)
));

# Straight diagonal vectorization:
VectorizedDMatSPL := ( isa, spl ) -> Checked(Rows(spl) <= isa.v, let(
    mat   := Cond(IsMat(spl), spl, MatSPL(spl)),
    r     := Rows(mat),
    c     := Cols(mat),
    padd  := (offs, l) -> Replicate(offs, 0) :: l :: Replicate(isa.v - offs - Length(l), 0),
    VBlk( [ List( [0 .. r-2], i -> isa.t.value(padd(r-i-1, List( [0 .. _unwrap(min(i,c-1))], j -> mat[r-i+j][j+1] )))) ::
            List( [0 .. c-1], i -> isa.t.value(padd(    0, List( [0 .. When(i+r<=c, r-1, c-i-1)], j -> mat[j+1][i+j+1] )))) ],
         isa.v)
));

# Zero pad both dimensions:
ZeroPadMat := (mat, v) -> let(
    r   := Rows(mat),
    c   := Cols(mat),
    dr  := CeilingRat(r/v)*v - r,
    dc  := CeilingRat(c/v)*v - c,
    List( mat, row -> row :: Replicate(dc, 0) ) :: Replicate(dr, Replicate(c+dc, 0))
);

# Zero pad matrix by adding zero rows:
ZeroPadMatV := (mat, v) -> let(
    r   := Rows(mat),
    c   := Cols(mat),
    dr  := CeilingRat(r/v)*v - r,
    mat :: Replicate(dr, Replicate(c, 0))
);

NewRulesFor(TMat, rec(
    TMat_Base := rec(
	applicable := (self, t) >> not t.hasTags(),
        apply := (self, t, C, Nonterms) >> Mat(t.params[1])
    ),

    TMat_Vec := rec(
	forTransposition := false,
	maxSize := 128,
	applicable := (self, t) >> t.hasTags() and t.isTag(1, AVecReg) and
	    (self.maxSize<0 or (Rows(t) <= self.maxSize and Cols(t) <= self.maxSize)),

        apply := (self, t, C, Nonterms) >> let(
	    isa := t.getTag(AVecReg, 1).isa,
	    VectorizedMatSPL(isa, t)
        )
    ),

    # vectorization of square matrix with dimensions equal to vector width using cyclic shifts
    TMat_VecZ := rec(
        forTransposition := false,
	applicable := (self, t) >> t.hasTags() and t.isTag(1, AVecReg) and
	    Rows(t) = Cols(t) and Rows(t)=t.getTag(AVecReg).v,

        apply := (self, t, C, Nonterms) >> let( isa := t.getTag(AVecReg).isa, i := Ind(isa.v),
	    VectorizedZMatSPL(isa, t)
	     * ISum(i, VScat(fBase(isa.v, i), isa.v) * VPerm(Z(isa.v,-i), isa.rotate_left(i), isa.v, 1))
	     * VGath(fId(1), isa.v)
        )
    ),

    # Vectorization of non square matrix by diagonals
    TMat_VecD := rec(
        forTransposition := false,
        maxSize := 256,
	applicable := (self, t) >> t.hasTags() and t.isTag(1, AVecReg) and
	    Rows(t) <= t.getTag(AVecReg).v and (self.maxSize<0 or Cols(t) <= self.maxSize),

        apply := (self, t, C, Nonterms) >> let(
            isa := t.getTag(AVecReg).isa,
            v   := isa.v,
            r   := Rows(t),
            c   := Cols(t),
            ru  := _roundup(c + v-1, v),
            i   := Ind(r-1),
            j   := Ind(c),
            VScat_pc( r, r, 0, v) *
            VectorizedDMatSPL(isa, t) *
            SUM(
                When(r>1,
                [ ISum(i,
                    VScat(fBase(r+c-1,     i), v) * VPerm(Z(v,v-r+i+1), isa.rotate_left(r-i-1), v, 1)
                     * VGath_pc(c, min(c, v), 0, v))], []) ::
                [ISum(j,
                    VScat(fBase(r+c-1, j+r-1), v)
                     * VGath_pc(c, min(v, c-j), j, v)).unrolledChildren()]
            )
        )
    ),

    TMat_UnrolledBlkH := rec(
        forTransposition := false,
	applicable := (self, t) >> t.hasTags() and t.isTag(1, AVecReg) and
	    Rows(t) > t.getTag(AVecReg).v,

        freedoms := (self, t) >> [[t.getTag(AVecReg).v]],

	child := (self, t, fr) >> Flat(MapMat(GetBlocksMat(ZeroPadMatV(t.params[1], t.getTag(AVecReg).v), fr[1], Cols(t)), e -> TMat(e).withTags(t.getTags()))),

        apply := (self, t, C, Nonterms) >> let(
            isa := t.getTag(AVecReg).isa,
            v   := isa.v,
            c   := Cols(t),
            r   := Rows(t),
            VScat_pc(r, r, 0, v) * SUM(
                List( [1..Length(C)], i -> VScat(fAdd(Length(C), 1, i-1), v) * C[i] ))
        )
    ),

    TMat_UnrolledBlkV := rec(
        forTransposition := false,
	applicable := (self, t) >> t.hasTags() and t.isTag(1, AVecReg) and
	    Cols(t) > t.getTag(AVecReg).v,

        freedoms := (self, t) >> [[1]],

	child := (self, t, fr) >> let(
	    v := t.getTag(AVecReg).v,
	    m := ZeroPadMat(t.params[1], v),
	    Flat(MapMat(GetBlocksMat(m, Rows(m), v), e -> TMat(e).withTags(t.getTags())))
	),

        apply := (self, t, C, Nonterms) >> let(
            isa := t.getTag(AVecReg).isa,
            v   := isa.v,
            c   := Cols(t),
            r   := Rows(t),
            h   := CeilingRat(r/v),
            VScat_pc(r, r, 0, v) *
            VBlk( MapMat(TensorProductMat( [Replicate( CeilingRat(c/v), 1)], IdentityMat(h) ), e -> When(e=0, isa.t.zero(), isa.t.one())), v)*
            SUM(
                List( [1..Length(C)], i -> VScat(fAdd(Length(C)*h, h, (i-1)*h), v) * C[i] * VGath(fAdd(Length(C), 1, i-1), v) ))*
            #HStack(C) *
            VGath_pc(c, c, 0, v)
        )
    ),

    TMat_Blk := rec(
    	forTransposition := false,

    	applicable := (self, t) >> t.hasTags() and t.isTag(1, AVecReg)
	    and (Rows(t) > self.rblk or Cols(t) > self.cblk),

	# Matrix partition
	# N - normal block (rblk x cblk), T - tall blocks, S - short blocks, C - corner bock
	#
	# N ... N T
	# . \   . T
	# .   \ . T
	# N ... N T
        # S S S S S
	#
	_NTSC_Split := meth(self, mat, rblk, cblk)
	    local r, c, nr, nc, N, T, S, C;
	    [ r,  c] := Dimensions(mat);
	    [nr, nc] := [r - r mod rblk,  c - c mod cblk];
	    N := mat{[1..nr]}{[1..nc]};
	    T := mat{[1..nr]}{[nc+1..c]};
	    S := mat{[nr+1..r]}{[1..nc]};
	    C := mat{[nr+1..r]}{[nc+1..c]};
	    return [N, T, S, C];
	end,

        apply := meth(self, t, C, Nonterms)
	    local isa, N, T, S, C, blks, dat, dvar, i, j;

    	    isa := t.getTag(AVecReg, 1).isa;

	    [N, T, S, C] := self._NTSC_Split(t.params[1], self.rblk, self.cblk);
	    blks := rec(
		N := MapMat(GetBlocksMat(N, self.rblk, self.cblk), x -> VectorizedMatSPL(isa, x)),
		T := MapMat(GetBlocksMat(T, self.rblk, self.cblk), x -> VectorizedMatSPL(isa, x)),
		S := MapMat(GetBlocksMat(S, self.rblk, self.cblk), x -> VectorizedMatSPL(isa, x))
	    );
	    dat := rec(
		N := V(List(Collect(blks.N, VBlk), x->x.element)),
		T := V(List(Collect(blks.T, VBlk), x->x.element)),
		S := V(List(Collect(blks.S, VBlk), x->x.element))
	    );
	    dvar := rec(
		N := var.fresh_t("DN", dat.N.t),
		T := var.fresh_t("DT", dat.T.t),
		S := var.fresh_t("DS", dat.S.t)
	    );

	    [i, j] := [Ind(Rows(blks.N)), Ind(Cols(blks.N))];

	    return
	    Data(dvar.N, dat.N,
	    Data(dvar.T, dat.T,
	    Data(dvar.S, dat.S,
		VStack(
		    IterVStack(i,
			HStack(
			    IterHStack(j,
				BB(SubstTopDownNR(blks.N[1][1], VBlk, x->VBlk(nth(dvar.N, i*Cols(blks.N)+j), isa.v)))
			    ),
			    BB(SubstTopDownNR(blks.T[1][1], VBlk, x->VBlk(nth(dvar.T, i), isa.v)))
			)
		    ),
		    HStack(
			IterHStack(j,
			    BB(SubstTopDownNR(blks.S[1][1], VBlk, x->VBlk(nth(dvar.S, j), isa.v)))
			),
			VectorizedMatSPL(isa, C)
		    )
		)
	    )));
	end
    )
));


text:
IsLeapYear := function(n)
  return (n mod 4 = 0) and ((n mod 100 <> 0) or (n mod 400 = 0));
end;

# alternative using built-in function
IsLeapYear := function(n)
  return DaysInYear(n) = 366;
end;

ext: gap

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


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

#F DiagFunc - base class for diagonal generating functions
#F Required fields in subclasses:
#F    .abbrevs = [ (p1)     -> [parlist1],
#F                 (p1,...) -> [parlist2] ]
#F       Mapping from various lengths of parameter lists to a
#F       fixed signature .def
#F
#F    .lambda(self)  - must return Lambda() object representing the actual
#F                     function. One should use .params field here
#F
#F    .domain(self) -- must return the domain of the function (integer or type)
#F    .range(self)  -- must return the range of the function (integer or type)
#F
#F Optional:
#F    .def(self, p1, ...)  - this method is called from the constructor
#F                           before the instance object is created.
#F                           <self> = class in this case
#F                           .def() must return a record, all fields of the record
#F                           will be copied to the resulting instance object,
#F                           which will also have the field .params := [p1, ...].
#F                           additional error checking can be performed in .def
#F                           NB: parameter lists are normalized here, using .abbrevs
Class(DiagFunc, Function, Sym, rec(
    isSPL := false,

    range  := self >> Error(ObjId(self), " does not implement .range()"),
    domain := self >> Error(ObjId(self), " does not implement .domain()"),

    tensor_op := (v,rows,cols,fl,gl) ->
        Lambda(v, fl.at(idiv(v, cols)) * gl.at(imod(v, cols))),

    sums := self >> self,
    rSetChild := meth(self, n, newChild)
        self.params[n] := newChild;
    end,

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

    def := arg -> rec(),

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

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

        result := SPL(WithBases(self, CopyFields(
                    rec(params := params),
                    ApplyFunc(self.def,params))));

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

    free := self >> Union(List(self.params, FreeVars)),
    __call__ := ~.fromDef,

    isReal := self >> not IsComplexT(self.range())
));

#F II - interval diagonal function
#F
#F II(<N>) - diagonal of N 1's
#F II(<N>, <to>) - diagonal with 1's from 0..<to>, and 0's after
#F II(<N>, <from>, <to>) - diagonal with 1's in <from>..<to>, 0's elsewhere
#F
Class(II, DiagFunc, rec(
    abbrevs := [ (N)         -> [N,0,N],
                 (N,to)      -> [N,0,to],
                 (N,from,to) -> [N,from,to] ],

    def := (N, from, to) -> Checked(
	IsPosInt0Sym(N), IsPosInt0Sym(from), IsPosInt0Sym(to), rec()),

    lambda := self >> let(
        i := Ind(self.params[1]), from := self.params[2], to := self.params[3],
        Lambda(i, cond(leq(from, i, to-1), 1, 0))),

    domain := self >> self.params[1],
    range := self >> TInt,
));

#F fConst(<t>, <N>, <c>) - constant diagonal function, <N> values of <c> of type <t>
#F
Class(fConst, DiagFunc, rec(
    abbrevs := [   (N,c) -> Checked(           IsPosInt0Sym(N), [TReal, N,c]),
                 (t,N,c) -> Checked(IsType(t), IsPosInt0Sym(N), [t, N, c]) ],
    range := self >> self.params[1],
    domain := self >> self.params[2],
    lambda := self >> let(i := Ind(self.params[2]), Lambda(i, self.params[3]))
));

#F fUnk(<t>, <N>) - dummy diagonal that denotes an "unknown" function, <N> values of type <t>
#F   To be used for .hashAs methods of transforms
#F
Class(fUnk, DiagFunc, rec(
    abbrevs := [ N -> Checked(IsIntSym(N), [TReal, N]),
	        (t, N) -> Checked(IsType(t), IsIntSym(N), [t, N]) ],
    range := self >> self.params[1],
    domain := self >> self.params[2],
    lambda := self >> let(N:=self.params[2], i:=Ind(N),
	Cond(IsIntT(self.params[1]), Lambda(i, i+1),
	                             Lambda(i, self.params[1].one()/N*(i+1))))
));

# FUnk is for backwards compatibility, use fUnk
FUnk := n -> fUnk(TReal, n);

Declare(LD);

# UD(<N>, <shift>) - upper diagonal
#
Class(UD, Sym, rec(
    def := (n,k) -> Checked(k < n, Z(n,k) * Diag(II(n, k, n))),
    transpose := self >> ApplyFunc(LD, self.params)
));

# LD(<N>, <shift>) - lower diagonal
#
Class(LD, Sym, rec(
    def := (n,k) -> Checked(k < n, Z(n,n-k) * Diag(II(n, 0, n-k))),
    transpose := self >> ApplyFunc(UD, self.params)
));

# [ 1 1
#     1 1
#       ... ]
Class(S, Sym, rec(
    def := n -> SUM(I(n), UD(n,1))));

# Class(fComputeOnline, PlaceholderExp);
Class(fComputeOnline, FuncClassOper, rec(
    range  := self >> self._children[1].range(),
    domain := self >> self._children[1].domain(),
    lambda := self >> self._children[1].lambda(),
    tolist := self >> self._children[1].tolist(),
    free   := self >> self._children[1].free(),
    ev     := self >> self._children[1].ev(),
    eval   := self >> self._children[1].eval(),
));

Declare(fPrecompute);

# fPrecompute(<func>)
#
# Marks functions that should be precomputed by either generating a
# table (opts.generateInitCode := true) or runtime initialization code
# (opts.generateInitCode := true).
#
# Actual precomputation is done by Process_fPrecompute. Which uses
# options record 'opts' to determine the precomputation strategy.
#
Class(fPrecompute, FuncClassOper, rec(
    def    := func -> rec(t := TPtr(func.range())),
    range  := self >> self._children[1].range(),
    domain := self >> self._children[1].domain(),
    lambda := self >> self._children[1].lambda(),
    tolist := self >> self._children[1].tolist(),
    free   := self >> self._children[1].free(),
    isReal := self >> self._children[1].isReal()
));

# diagAdd(<f1>, <f2>,...) : i -> f1(i) + f2(i) + ...
# Pointwise addition of several diagonal functions
Class(diagAdd, FuncClassOper, rec(
    domain := self >> self.child(1).domain(),
    range  := self >> UnifyTypes(List(self.children(), x->x.range())),
    lambda := self >> let(v := Ind(self.domain()),
	Lambda(v, ApplyFunc(add, List(self.children(), i->i.at(v))))),
));

# diagMul(<f1>, <f2>,...) : i -> f1(i) * f2(i) * ...
# Pointwise multiplication of several diagonal functions
Class(diagMul, FuncClassOper, rec(
    domain := self >> self.child(1).domain(),
    range  := self >> UnifyTypes(List(self.children(), (x) -> x.range())),
    lambda := self >> let(v := Ind(self.domain()),
	Lambda(v, ApplyFunc(mul, List(self.children(), i->i.at(v))))),
));

Class(diagCond, FuncClassOper, rec(
    domain := self >> self.child(2).domain(),
    range  := self >> UnifyTypes(List(Drop(self.children(), 1), x->x.range())),
    lambda := self >> let(v := Ind(self.domain()),
	Lambda(v, cond(self.child(1), self.child(2).at(v), self.child(3).at(v))))
));

#F fCast(<toT>, <fromT>)  -- symbolic type cast function
#F
#F Example: fCast(TReal, TInt)
#F
Class(fCast, DiagFunc, rec(
    abbrevs := [ (toT, fromT)  -> Checked(IsType(toT), IsType(fromT), [toT, fromT]) ],
    range := self >> self.params[1],
    domain := self >> self.params[2],
    lambda := self >> let(i:=TempVar(self.domain()), Lambda(i, tcast(self.range(), i))),
));

ext: gi


text:
v := [1 .. 8];

Sum(v);
# 36

Product(v);
# 40320

# You can sum or multiply the result of a function

Sum(v, n -> n^2);
# 204

Product(v, n -> 1/n);
# 1/40320


text:
fib := function(n)
  local a;
  a := [[0, 1], [1, 1]]^n;
  return a[1][2];
end;


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

########################################################################
#   rules for A x I, I x A, (A x I)L, (I x A)L
Declare(AxI_vec);

NewRulesFor(TTensorI, rec(
#   vectorization cases

#   I_n x A_rxs -> I_n/v x L^rv_v(A_rxs x I_v)L^sv_s
    IxA_vec_push := rec(
        info := "IxA vec",
        forTransposition := false,
        applicable := nt -> IsParPar(nt.params) and (nt.isTag(1,AVecReg) or nt.isTag(1,AVecRegCx)),
        children := nt -> let(
            pv := nt.firstTag(),
            v := pv.v,
            d := nt.params[1].dims(),
            krnl := nt.params[1].withTags(nt.getTags()),
            [[ krnl ]]
	    # YSV: below is the original piece of code and comment. This is no longer works
	    #      and does not seem to be necessary. Above works quite well. Delete the below
	    #      as soon as everything is tested.
            #D  When(krnl.isReal(), krnl.setWrap(VWrapId), krnl.setWrap(VWrapTRC(pv.isa))) ]]
        ),  #D  thats bad but required to keep wrapping working :(
        apply := (nt, c, cnt) -> Tensor(I(nt.params[2]), c[1])

    ),

#   A x I_n -> (A x I_n/v) x I_v
    AxI_vec := rec(
        info := "A x I_n -> (A x I_n/v) x I_v",
        forTransposition := false,
        applicable := nt -> nt.hasTags() and IsVecVec(nt.params) and (nt.isTag(1,AVecReg) or nt.isTag(1,AVecRegCx)) and IsInt(nt.params[2]/nt.firstTag().v),
        children := nt -> let(
            r := nt.params[2] / nt.firstTag().v,
            isa := nt.firstTag().isa,
            [[ When(r = 1,
                When(nt.numTags() = 1,
                    nt.params[1].setWrap(VWrap(isa)),
                    nt.params[1].setWrap(Drop(nt.getTags(), 1)).setWrap(VWrap(isa))
                ),
                TTensorI(nt.params[1].setWrap(VWrap(isa)), r, AVec, AVec).withTags(Drop(nt.getTags(), 1))
            )]]
        ),
        apply := (nt, c, cnt) -> VTensor(c[1], nt.firstTag().v)
    ),
#   I_n x A_rxs -> I_n/v x L^rv_v(A_rxs x I_v)L^sv_s
    IxA_vec := rec(
        info := "IxA vec",
        forTransposition := false,
        applicable := nt -> IsParPar(nt.params) and nt.hasTags() and (nt.isTag(1,AVecReg) or nt.isTag(1,AVecRegCx)) and IsInt(nt.params[2]/nt.firstTag().v),
        children := nt -> let(
            pv := nt.getTags(),
            v := pv[1].v,
            isa := pv[1].isa,
            d := nt.params[1].dims(),
            [[
                TL(d[1]*v, v, 1, 1).withTags(pv).setWrap(VWrapId),
                When(Length(pv)=1, nt.params[1].setWrap(VWrap(isa)), nt.params[1].setpv(Drop(pv, 1)).setWrap(VWrap(isa))),
                TL(d[2]*v, d[2], 1, 1).withTags(pv).setWrap(VWrapId)
            ]]
        ),
        apply := (nt, c, cnt) -> let(
            l := nt.params[2] / nt.firstTag().v,
            A := c[1] * VTensor(c[2], nt.firstTag().v) * c[3],
            NoDiagPullin(When(l=1, A , Tensor(I(l), A)))
        )
    ),
#   (I_n x A_rxs)L^nr_n -> (I_n/v x L^rv_v(A_rxs x I_v))(L^ns/v_n/v x I_v)
    IxA_L_vec := rec(
        info := "(IxA)L vec",
        forTransposition := false,
        maxSize := -1,
        applicable := (self, nt) >>
            ObjId(nt.params[1]) <> TTensorI
            and IsParVec(nt.params)
            and nt.hasTags()
            and (nt.isTag(1,AVecReg) or nt.isTag(1,AVecRegCx))
            and IsInt(nt.params[2] / nt.firstTag().v)
            and (
                self.maxSize = -1
                or (
                    Rows(nt.params[1]) <= self.maxSize
                    and Cols(nt.params[1]) <= self.maxSize
                )
            ),
        children := nt -> let(
            pv := nt.getTags(),
            v := pv[1].v,
            isa := pv[1].isa,
            d := nt.params[1].dims(),
            [[
                TL(d[1]*v, v).withTags(pv).setWrap(VWrapId),
                When(Length(pv)=1,
                    nt.params[1].setWrap(VWrap(isa)),
                    nt.params[1].withTags(Drop(pv, 1)).setWrap(VWrap(isa))
                )
            ]]
        ),
        apply := (nt, c, cnt) -> let(
            v := nt.firstTag().v,
            m := nt.params[2] / v,
            d := nt.params[1].dims(),
            NoDiagPullinLeft(
                When(m = 1,
                    c[1] * VTensor(c[2], v),
                    Tensor(I(m), c[1] * VTensor(c[2], v)) * VTensor(L(m*d[2], m), v)
                )
            )
        ),
        switch := false

    ),
#   L^nr_r(I_n x A_rxs) -> (L^nr/v_r/v x I_v)(I_n/v x (A_rxs x I_v)L^sv_s)
    L_IxA_vec := rec(
        info := "L(IxA) vec",
        forTransposition := false,
        maxSize := -1,
        applicable := (self, nt) >>
            ObjId(nt.params[1]) <> TTensorI
            and IsVecPar(nt.params)
            and nt.hasTags()
            and (nt.isTag(1,AVecReg) or nt.isTag(1,AVecRegCx))
            and IsInt(nt.params[2] / nt.firstTag().v)
            and (
                self.maxSize = -1
                or (
                    Rows(nt.params[1]) <= self.maxSize
                    and Cols(nt.params[1]) <= self.maxSize
                )
            ),
        children := nt -> let(
            pv := nt.getTags(),
            v := pv[1].v,
            isa := pv[1].isa,
            d := nt.params[1].dims(),
            [[
                When(nt.numTags() = 1,
                    nt.params[1].setWrap(VWrap(isa)),
                    nt.params[1].withTags(Drop(pv, 1)).setWrap(VWrap(isa))
                ),
                TL(d[2]*v, d[2]).withTags(pv).setWrap(VWrapId)
            ]]
        ),
        apply := (nt, c, cnt) -> let(
            v := nt.firstTag().v,
            m := nt.params[2] / v,
            d := nt.params[1].dims(),
            NoDiagPullinRight(
                When(m = 1,
                    VTensor(c[1], v) * c[2],
                    VTensor(L(m*d[1], d[1]), v) * Tensor(I(m), VTensor(c[1], v) * c[2])
                )
            )
        ),

        switch := false
    ),
############################################################################
#   L splitting variant
#   (I_n x A_rxs) L^ns_n
    IxA_L_split_vec := rec(
        info := "split (I_n x A_rxs) L^ns_n",
        forTransposition := false,
        applicable := nt -> IsParVec(nt.params) and nt.hasTags() and (nt.isTag(1,AVecReg) or nt.isTag(1,AVecRegCx)),
        children := nt -> [[ TTensorI(nt.params[1], nt.params[2], AVec, AVec).withTags(nt.getTags()) ]],
        apply := (nt, c, cnt) -> let(d := nt.params[1].dims()[1],
            L(d * nt.params[2], nt.params[2]) * c[1]
        ),
        switch := false
    ),
#   L^nr_r (I_n x A_rxs)
    L_IxA_split_vec := rec(
        info := "split L^nr_r (I_n x A_rxs)",
        forTransposition := false,
        applicable := nt -> IsVecPar(nt.params) and (nt.isTag(1,AVecReg) or nt.isTag(1,AVecRegCx)),
        children := nt -> [[ TTensorI(nt.params[1], nt.params[2], AVec, AVec).withTags(nt.getTags()) ]],
        apply := (nt, c, cnt) -> let(d := nt.params[1].dims()[2],
            c[1] * L(d * nt.params[2], d)
        ),
        switch := false
    ),
    IxA_split_vec := rec(
        info := "split (I_n x A_rxs)",
        forTransposition := false,
        applicable := nt -> IsParPar(nt.params) and (nt.isTag(1,AVecReg) or nt.isTag(1,AVecRegCx)),
        children := nt -> [[ TTensorI(nt.params[1], nt.params[2], AVec, AVec).withTags(nt.getTags()) ]],
        apply := (nt, c, cnt) -> let(d := nt.params[1].dims(),
            L(d[1] * nt.params[2], nt.params[2]) * c[1] * L(d[2] * nt.params[2], d[2])
        ),
        switch := false
    ),
    IxA_conj_vec := rec(
        info := "split (I_n x A_rxs)",
        forTransposition := false,
        applicable := nt -> IsParPar(nt.params) and (nt.isTag(1,AVecReg) or nt.isTag(1,AVecRegCx)),
        children := nt -> [[ TTensorI(nt.params[1], nt.params[2], AVec, AVec).withTags(nt.getTags()) ]],
        apply := (nt, c, cnt) -> let(d := nt.params[1].dims(),
            ConjLR(c[1], L(d[1] * nt.params[2], nt.params[2]), L(d[2] * nt.params[2], nt.params[2]))
        ),
        switch := false
    ),

############################################################################
#   for now, these rules are for unrolled cod only
#   non-vref TTensorI
#   A x I_n, v\nmid n
    AxI_svec := rec(
        info := "A x I_n for v\nmid n",
        forTransposition := false,
        applicable := nt ->
            IsVecVec(nt.params)
            and (nt.isTag(1,AVecReg) or nt.isTag(1,AVecRegCx))
            #and let( d := nt.params[1].dims(), d[1] = d[2] )
            and not IsInt(nt.params[2] / nt.firstTag().v),

        children := nt -> let(
            v := nt.firstTag().v,
            r := v * QuoInt(nt.params[2] + v - 1, v),
            [[ TTensorI(nt.params[1], r, AVec, AVec).withTags(nt.getTags()) ]]
        ),
        apply := (nt, c, cnt) -> let(
            d := nt.params[1].dims(),
            v := nt.firstTag().v,
            IxVScat_pc(d[1], nt.params[2], nt.params[2], 0, v)
            * c[1]
            * IxVGath_pc(d[2], nt.params[2], nt.params[2], 0, v)
        ),
    ),
#   (I_m x A_n )L^mn_m, v\nmid m
    IxA_L_svec := rec(
        info := "(I_m x A_n )L^mn_m, v\nmid m",
        forTransposition := false,
        applicable := nt -> IsParVec(nt.params) and (nt.isTag(1,AVecReg) or nt.isTag(1,AVecRegCx))
            and let(
                v := nt.firstTag().v,
                d := nt.params[1].dims(),
                d[1] = d[2]
                and (
                    not IsInt(nt.params[2] / v)
                    or (not IsInt(d[1] / v))
                )
            ),
        children := nt -> let(
            v := nt.firstTag().v,
            d := nt.params[1].dims(),
            r := v * QuoInt(nt.params[2]+v-1, v),
            s := v * QuoInt(d[1]+v-1, v),
            [[
                TL(r*s, r, 1, 1).withTags(nt.getTags()),
                TTensorI(nt.params[1], r, AVec, AVec).withTags(nt.getTags())
            ]]
        ),
        apply := (nt, C, cnt) -> let(
            d:=nt.params[1].dims()[1],
            v := nt.firstTag().v,
            c1c := Cols(C[1])/v, c1r := Rows(C[1])/v,
            c2c := Cols(C[2])/v, c2r := Rows(C[2])/v,
            C0 := IxVScat_pc(nt.params[2],d,d,0,v),
            c0c := Cols(C0)/v,
            C0
#           * VGath(fAdd(c1r,c0c,0), v)
            * VGath_zero(c1r,c0c, v)
#           * VTensor(HStack(I(c0c), O(c0c,c1r-c0c)),v)
            * C[1]
            * VScat_zero(c1c,c2r,v)
#           * VTensor(VStack(I(c2r), O(c1c-c2r,c2r)),v)
            * C[2]
            * IxVGath_pc(d,nt.params[2],nt.params[2],0,v)
        )
    ),
    TTensorI_oddvloop := rec(
        switch  := false,
        maxSize := false,
        filter  := e->true,
        applicable := nt ->
            IsVecVec(nt.params)
            and (nt.isTag(1,AVecReg) or nt.isTag(1,AVecRegCx))
            and not IsInt(nt.params[2] / nt.firstTag().v),

        children := nt -> let(
            v := nt.firstTag().v,
           [[ TTensorI(nt.params[1], v, AVec, AVec).withTags(nt.getTags()) ]]
        ),

        apply := (t, C, Nonterms) -> let(
  	        v := t.firstTag().v,
            kernel := C[1],
            m := Rows(t.params[1]),
            n := Cols(t.params[1]),
            k := t.params[2],
            kd := _rounddown(k,v).v,
            kr := k-kd,
            i := Ind(kd/v),
            #---------------------
            SUM(
                ISum(i,
                    IxVScat_pc(n, k, v, v*i, v) *
                    kernel *
                    IxVGath_pc(m, k, v, v*i, v)
                ),
                IxVScat_pc(n, k, k-kd, kd, v) *
                kernel *
                IxVGath_pc(m, k, k-kd, kd, v)
            )
        )
    )
));


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

Class(VectorizedBaseMat, BaseMat, SumsBase, rec(
    transpose := self >> Inherit(self, rec(transposed := not self.transposed)),
    #-----------------------------------------------------------------------
    print := (self,i,is) >> Print(self.name, "(", self.n, ", ", self.v,")",
	self.printA(), When(self.transposed, ".transpose()", "")),
));

#F VS(<n>, <v>)  --  vectorized version of doubly-diagonal matrix
#F
#F Represents the matrix that has 1s on the main diagonal and upper shifted by 1 diagonal
#F
#F   <n> - matrix size
#F   <v> - vector length, does not affect matrix shape, only needed for Codegen
#F
#F Example:
#F
#F spiral> PrintMat(MatSPL(VS(8,4)));
#F [ [ 1, 1,  ,  ,  ,  ,  ,   ],
#F   [  , 1, 1,  ,  ,  ,  ,   ],
#F   [  ,  , 1, 1,  ,  ,  ,   ],
#F   [  ,  ,  , 1, 1,  ,  ,   ],
#F   [  ,  ,  ,  , 1, 1,  ,   ],
#F   [  ,  ,  ,  ,  , 1, 1,   ],
#F   [  ,  ,  ,  ,  ,  , 1, 1 ],
#F   [  ,  ,  ,  ,  ,  ,  , 1 ] ]
#F
Class(VS, VectorizedBaseMat, rec(
    new := (self, n, v) >> SPL(WithBases(self,
        rec(n := n, v := v, transposed := false))).setDims(),
    #-----------------------------------------------------------------------
    dims := self >> [self.n, self.n],
    #-----------------------------------------------------------------------
    rChildren := self >> [self.n, self.v],
    rSetChild := rSetChildFields("n", "v"),
    #-----------------------------------------------------------------------
    area := self >> self.n * self.n,
    #-----------------------------------------------------------------------
    toAMat := self >> let(s:=S(self.n), AMatMat(MatSPL(When(self.transposed, s.transpose(), s))))
));

Declare(VLD);

#F VUD(<n>, <v>)  --  vectorized upper diagonal matrix
#F
#F VUD(n,v) is equivalent to UD(n,1)
#F
#F   <n> -- matrix size
#F   <v> -- vector length, does not affect matrix shape, only needed for Codegen
#F
#F spiral> PrintMat(MatSPL(VUD(8,4)));
#F [ [  , 1,  ,  ,  ,  ,  ,   ],
#F   [  ,  , 1,  ,  ,  ,  ,   ],
#F   [  ,  ,  , 1,  ,  ,  ,   ],
#F   [  ,  ,  ,  , 1,  ,  ,   ],
#F   [  ,  ,  ,  ,  , 1,  ,   ],
#F   [  ,  ,  ,  ,  ,  , 1,   ],
#F   [  ,  ,  ,  ,  ,  ,  , 1 ],
#F   [  ,  ,  ,  ,  ,  ,  ,   ] ]
#F
Class(VUD, VectorizedBaseMat, rec(
    new := (self, n, v) >> SPL(WithBases(self, rec(n := n, v := v))).setDims(),
    #-----------------------------------------------------------------------
    transpose := (self) >> VLD(self.n, self.v),
    #-----------------------------------------------------------------------
    dims := self >> [self.n, self.n],
    #-----------------------------------------------------------------------
    rChildren := self >> [self.n, self.v],
    rSetChild := rSetChildFields("n", "v"),
    #-----------------------------------------------------------------------
    area := self >> self.n,
    #-----------------------------------------------------------------------
    toAMat := self >> UD(self.n).toAMat()
));

#F VLD(<n>, <v>)  --  vectorized lower diagonal matrix
#F
#F VLD(n,v) is equivalent to LD(n,1)
#F
#F   <n> -- matrix size
#F   <v> -- vector length, does not affect matrix shape, only needed for Codegen
#F
#F spiral> PrintMat(MatSPL(VLD(8,4)));
#F [ [  ,  ,  ,  ,  ,  ,  ,   ],
#F   [ 1,  ,  ,  ,  ,  ,  ,   ],
#F   [  , 1,  ,  ,  ,  ,  ,   ],
#F   [  ,  , 1,  ,  ,  ,  ,   ],
#F   [  ,  ,  , 1,  ,  ,  ,   ],
#F   [  ,  ,  ,  , 1,  ,  ,   ],
#F   [  ,  ,  ,  ,  , 1,  ,   ],
#F   [  ,  ,  ,  ,  ,  , 1,   ] ]
#F
Class(VLD, VectorizedBaseMat, rec(
    new := (self, n, v) >> SPL(WithBases(self, rec(n := n, v := v))).setDims(),
    #-----------------------------------------------------------------------
    transpose := (self) >> VUD(self.n, self.v),
    #-----------------------------------------------------------------------
    dims := self >> [self.n, self.n],
    #-----------------------------------------------------------------------
    rChildren := self >> [self.n, self.v],
    rSetChild := rSetChildFields("n", "v"),
    #-----------------------------------------------------------------------
    area := self >> self.n,
    #-----------------------------------------------------------------------
    toAMat := self >> LD(self.n, 1).toAMat()
));


text:
IsLuhn := function(n)
	local c, d, i, j, r;
	d := "0123456789";
	j := 1;
	r := 0;
	for c in Reversed(String(n)) do
		i := Position(d, c);
		if i = fail then
			continue;
		fi;
		i := j*(i - 1);
		r := r + QuoInt(i, 10) + RemInt(i, 10);
		j := 3 - j;
	od;
	return RemInt(r, 10) = 0;
end;

List([49927398716, 49927398717, 1234567812345678, 1234567812345670], IsLuhn);
# [ true, false, false, true ]

# Will also work on strings, and will skip non-digits
IsLuhn("4-992-739-871-6");
# true